Simulations Become Unstable With an Increase in Bead Count

[Harry Sullivan]

Hello!

I am currently trying to simulate bulk benzene with a tabulated potential using i-pi and lammps. I have run the simulation with a few different bead counts to see the effect. However once the bead count gets high enough the potential energy spikes and causes the benzene rings to fold in on themselves. Is this expected behavior? I ran these calculations on two different seeds to see if it was just a bad initial velocity sampling, however the behavior persisted. A plot of this is attached. Looking at the RDFs of my simulations, I feel that an increase in bead count beyond 32 is necessary for the hydrogen. The fluctuations of the bonds don't seem to converge to a steady solution. A plot of the RDFs is provided as well. Is my error due to an improper force field or am I missing something within my input scripts? Is there any reason this divergence happens at 36 beads? Any help would be greatly appreciated. I also provided a link to a google drive containing my input scripts and potential. Inside each folder there is an input.unix.xml for i-pi, and within the sub folder SOPR I have my lammps input script alongside the tabulated potential I was using. I used SLURM to run these calculations so I have some slurm input scripts as well.

https://drive.google.com/drive/folders/1RfALeHbMvPI2TJUuHcGSijUXIrX77WKn?usp=sharing

Thanks in advanced.

Winston.

[Michele Ceriotti]

Hello. Problem seems to be related to the time step. Although the normal mode integrator allows running MD roughly with the same time step as a classical MD trajectory, you can encounter resonance problems (given that this is a form of multiple time step integrator), particularly with large bead numbers. 

I tried running your examples with a time step of 0.5fs and there is no divergence or instability. You could probably also stabilize simulations at 1fs by tuning the thermostat, using a BAOAB integrator, or the Cayley integrator (which is not super-well tested, but is meant to address exactly this kind of problems). 

Hope this helps. At worst, you can prob get away with 0.7fs which is not that much slower than what you were using. 

MC.

[Winston Sullivan]

This was great advice. Thanks for giving some insight on the details as to why this error happens.

Best, Winston.