PBC and torchani/openMM

wiederm commented 5 months ago

Hi Peter,

As @wardhaddadin1 mentioned yesterday in the openmm/EXS call, the neighbor list implemented in Torchani seems to assume that input coordinates are wrapped in the simulation box.

I have tested this with torchani (5 Angstrom box length, pbc set in the torchani input options) and a 2 particle system in which I moved one of the particles outside the defined simulation cell. I have attached the plots showing the trajectory of the particles (the box is indicated in black, note the different x,y,z scales) and the calculated energy. The displacement outside of the initial box results in incorrect energies if distances are > box_length + cutoff.

I have also tested this with openmm-ml, with a 2 particle system running an MD simulation. The box length was 10 Angstrom, and the cutoff of ANI was 5 Angstrom. I have attached the input file and the script. The trajectory and the energy calculation seem to indicate that openmm-ml is also affected by this and the results are not correct. I have attached a plot of the potential energy at each time step. The trajectory can be reproduced using the script.

test.zip

wiederm commented 5 months ago

I have repeated the openmm-ml simulation to make sure that I can reproduce this with a reduced box size and increased simulation time: for the plot below the cell length was set to 6 Angstrom and simulation time was increased to 5 ps, otherwise I have used the script attached above.

chrisiacovella commented 5 months ago

I was able to replicate effectively the same behavior with my own script. Digging through the code, I think I found the issue.

In openmm-ml, the condition to define if the force is periodic looks at "getPeriodicBoxVectors" from the topology. Both system and topology provide functions to setPeriodicBoxVectors. In my script I had initially set these in the system (which of course should lead to this condition being false based on the code below), but setting these box vectors in the topology, it did not change the overall behavior.

https://github.com/openmm/openmm-ml/blob/688b1fc9528905633f13f3f0168b0c3a2f17e11f/openmmml/models/anipotential.py#L132

is_periodic = (topology.getPeriodicBoxVectors() is not None) or system.usesPeriodicBoundaryConditions()

In both cases, if I to call system.usesPeriodicBoundaryConditions() it is false. My understanding is this looks at the registered forces to determine if there is periodicity. So if I look at system.getForce(0).usesPeriodicBoundaryConditions() it is also false.

I just grabbed the source of openmm-ml, and put in a few print statements, and I'm pretty sure in my case, it was an order of operations.

I initialized the molecule topology, then the potential, then set the box vectors in the topology. So basically the factory setting up ani2x only saw the topology without the box vectors set; it never looks at the topology again to see if this has changed.

Granted, the way the script is set up that Marcus posted, that shouldn't be an issue (since the box vectors got set when reading the pdb file), so I'm a bit confused at what is going on there. I will upload my script and info about which versions I'm running in a little bit.