I have an atomistic polymer system using OPLS obtained from tinker via the moltemplate distribution that I’m using. After reviewing the pair coefficients and not finding anything unusual, I’m not sure why the dynamics are too fast or so stiff that I’m getting an “out of range atoms” error as soon as the simulation starts. The error doesn’t come up if I use a timestep of 0.1 fs, but it does take considerably longer to run.
My input files are attached.
Any thoughts on what it is in my system that is causing this? Any chance it could be related to the hydrogen atoms having both pair coefficients (epsilon and sigma) equal to 0? I tried changing epsilon to 0.1, still crashed.
There are a number of contributing causes, not just a single reason.
That could be one contributing factor.
Changing only epsilon doesn’t really make a difference. When sigma is also 0. That would still turn off the interaction. But also changing sigma would be changing the force field, so that is not really a good option. Instead you should check with the relevant publication for OPLS/AA that you want to apply what is the correct atom type to be used for those hydrogen atoms and what parameters should be used. If they turn out to be 0.0 / 0.0, then this is correct. You may need to look for an alternate force field for your system then. If not, you have to update your input deck accordingly.
There are multiple other factors to consider:
your 1fs time step is very large for a system with hydrogen atoms (very light and thus very fast). In combination with the rather stiff force constants for bonds (when compared to “soft” Lennard-Jones interactions), your time integration for the bonds is likely diverging. For sustaining 1fs you usually would have to use fix shake to make the bonds involving hydrogens rigid.
your KSpace convergence with 1.0e-3 is on the low side. That leads to significant errors on forces and thus can also contribute to out of range errors.
you are running your simulations at unusually high temperatures for molecular systems. That would also lead to instabilities and larger per-time step displacements due to the higher kinetic energy. So that would often require a further reduction of the time step.
So lowering the time step is a must. For the non-bonded interactions, however, a choice of 0.5fs should suffice. The trouble is with the bonds. There you have two choices:
apply fix shake: e.g. with fix rigidh all shake 1.0e-4 20 1000 m 1.0
use r-RESPA: e.g. with run style respa 2 5 bond 1 pair 2
Keep monitoring your energies with frequent thermo output until you have a better handle on the system. From what it looks like, you need to have your system undergo a significant reorganization. So it also would be better to not put the entire simulation into one giant input file, but rather run the various stages separately and write out data files. That would allow you to tweak settings for individual steps, if needed for stability, without having to re-run everything and once your system is closer to equilibrium, you will likely be able to raise the timestep somewhat.