NPT in schnetpack.md

[sujona]

Hello, I have trained a PaiNN model of water, with inputs in Å, eV, eV/Å, and eV/Å^2 for the stress. In ASE both NVT and NPT simulations seem fine, but it's slow, I believe because of the ASENeighborList(?). In schnetpack.md, NVT works and is much faster, but in NPT the volume just expands to some unreasonable size. Is there something obvious to consider here? I'm using the current pip version. I have attached the run script.

npt-schnetmd.py

[jnsLs]

I can't spot an obvious problem in your script.
Since you are using a neighborlist with cutoff buffer, the neighborlist calculation should not be that expensive.
However, I would recommend using MatScipyNeighborList or VesinNeighborList as base_nbl.
How large is your system and how long do the respective calculations take?
Have you tried using the CLI for your simulations, and if so, have you been able to reproduce the problem?

[sujona]

In ASE's MD there was no skin thickness, didn't know that was possible. But I have now tried using the SkinNeighborList wrapper explicilty, as well as the neighbor lists that you suggest, and neither gave considerable speed improvement in ASE. It seems like what I thought was a constant time delay was simply quite poor scaling with CPU cores. A smaller NN runns faster, but still with poor CPU scaling, for either neighbor list (the ASE list is indeed a bit slower, which is noticeable with the smaller NN). The system is a cubic 18 Å box with 216 water molecules.

In SchnetPack's NPT (also on CPU), I tried replacing the initial momenta (300K Maxwell-Boltzmann) with equillibrated velocities from an ASE run, but the same expansion occurs. If I adjust the barostat time constant to 1000 (from 100) the expansion is slower, but still reaches 26 Å after 10 000 steps (5 ps).

Didn't try the CLI yet, I'll try to get to that soon!

[sujona]

I think I have some clue about the issue. If I set the pressure to 1000 bar, the cell fluctuates around the correct value. The NN is trained on eV and Angstrom, so returns the stress in eV/Å**3. An eV is about 100 kJ/mol, and a nm is 10 Å, where kJ/mol and nm are the internal units of schnetpack.md. So there seems to be a unit conversion error for the stress in this case. The SchNetPackCalculator is defined with the eV and Angstrom units keys.

Maybe related to virial / stress / pressure conversion when dividing by the volume in nm vs Angstrom?

On the other hand, 1 bar and 0 bar should give basically the same density so scaling target_pressure makes little sense, 1000*0 is still 0. Basically calculated pressure should be around 0. Maybe the issue could be in the kinetic contribution...

[sujona]

OK, the cell expansion in NPT in schnetpack.md happens when using TorchNeighborList. The other lists seem to work fine but a time step takes much longer.

[jnsLs]

Previously you mentioned that there must be a unit conversion error. The units are converted in the SpkCalculator and therefore should be the same for all nbh lists. Since most of the nbh lists seem to work, as far as I understand, it's unlikely that there is something wrong with the unit conversion.

Another remark: In the script you uploaded initially, you are using the schnetpack md package with the calculator and nbh list from the same package. You are using the NeighborListMD which has a skin thickness (cutoff_shell). If you however want to run the md in ASE I would recommend the calculator from our ase_interface.py

[sujona]

About the remark: the script i posted is only for schnetpack.md. I use a completely different script for running in ASE where I use the spk.interfaces.SpkCalculator, and that works fine, I think it's all in order. It is a bit slower than schnetpack.md with a TorchNeighborList, but also schnetpack.md is slower with the other neighbor lists. So it's OK!

About the NPT expansion with schnetpack.md and TorchNeighborList:
I agree, it doesn't seem like a unit error, as I wrote in the end of that message, a pressure of 0 is the same in all units, and should work fine for my system, but it didn't. The neighbor list is the deciding factor as far as I can tell.