Huge energies after energy minimization, only in larger simulations

Dear LAMMPS users,

Hello, I am having an issue with minimizing my energy. I don't think the issue is specifically with my minimization command, though. I can successfully minimize energy only if I'm using a small number of molecules: When I create a box of 27 (3X3X3) molecules (Tetrolic acid: http://www.chemspider.com/Chemical-Structure.61810.html), the minimization step completes just fine, minimizing the energy from -178.857129425 to -262.013061474 Kcal/mol (I am using real units).
However, this simulation is unstable. Particularly the hydrogens on the OH groups begin flying around at exceptional speeds, quickly raising the temperature to almost 2000K before one of them leaves and cannot be found. It may be worth noting that, even though I run 100 timesteps with random velocities before minimization to get rid of symmetry, the minimized configuration (as seen in VMD) is very symmetrical, with all molecules ending up in the same orientation with equal spacing. Before failing, several atoms begin overlapping (but they approach each other very slowly until they are nearly in the same position). I'm not sure if this has anything to do with my problems minimizing larger simulations.

In simulations with 4X4X4 molecules and larger, the energy begins at some sort of reasonable value. Afterwards, however, the total energy is on the order of 10^16 or 10^17 Kcal/mol, at which point LAMMPS obviously refuses to run the simulation. I had this problem once before, but fixed it and the cause must be different this time. The first time I ran into this, the issue was that my bond coefficients, pair coefficients, etc. were in a different order in my input script than they were in my data file. That is no longer the case and I am fairly positive that all those values are correct and in the right order at this time. If anybody else has had this issue, do you have any idea what my mistake might be? In case you would like to run my simulation, I've posted some text files online. Copy and paste them into a text editor and make sure to name the data files 'data.TTA' so that the input script can call it correctly.

data file with 27 TTA molecules: (This one should get through the minimization step just fine, but fail before completing the simulation (5000 timesteps to start)
http://justpaste.it/j3rv

data file with 64 TTA molecules: (This one should not minimize correctly. Energy before minimization: ~ -440 kcal/mol. After: ~ 3.567e17 kcal/mol
http://justpaste.it/j3rz

Input file:
http://justpaste.it/j3ru
*NOTE: depending on your processor, you may need to change or remove the first line. My computer has an i7 (which has 4 cores, 8 threads) so I tried to make use of all of them, but you could just remove that line. The simulation does not take very long.

If you need any more information, let me know. Thanks very much in advance for the help!
N.H.

If you minimize, set all the velocities to 0.0 or to a very low Temp, then run NVE, what happens?
Is your timestep small enough? Is eneregy conserved? Does the pressure

or temperature jump to high values?

Steve

Dear LAMMPS users,

Hello, I am having an issue with minimizing my energy. I don't think the issue is specifically with my minimization command, though. I can successfully minimize energy only if I'm using a small number of molecules: When I create a box of 27 (3X3X3) molecules (Tetrolic acid: http://www.chemspider.com/Chemical-Structure.61810.html), the minimization step completes just fine, minimizing the energy from -178.857129425 to -262.013061474 Kcal/mol (I am using real units).
However, this simulation is unstable. Particularly the hydrogens on the OH groups begin flying around at exceptional speeds, quickly raising the temperature to almost 2000K before one of them leaves and cannot be found. It may be worth noting that, even though I run 100 timesteps with random velocities before minimization to get rid of symmetry, the minimized configuration (as seen in VMD) is very symmetrical, with all molecules ending up in the same orientation with equal spacing. Before failing, several atoms begin overlapping (but they approach each other very slowly until they are nearly in the same position). I'm not sure if this has anything to do with my problems minimizing larger simulations.

my guess is that this is due to your ridiculously short LJ cutoff of
0.2 \AA. charges may get much too close and then atoms start
accelerating like crazy.

axel.

Thank you for your reply, Steve. I tried reducing the temperature before minimization and "successfully" minimized a 4X4X4 box of TTA, which I had been unable to do previously. I put "successfully" in quotes because, while it didn't fail, the total energy barely changed (from -583.3 kcal/mol to -591.7 kcal/mol), but LAMMPS used far fewer than the maximum allowed timesteps for minimization. In addition, the resulting minimized structure was again completely symmetrical. In fact, the behavior of all 64 molecules was almost identical for the first few hundred timesteps. I think that's a result of dropping the temperature and then scaling the velocities back up again after minimization.

I did try running NVE instead of NVT as well, but the result was pretty much the same. The timestep I used was 0.25 fs, and making it smaller did not affect the simulation, so I assume it's fine. Particles do not "jump" across the screen at the beginning of the simulation (when viewed in VMD), but move smoothly, until it gets out of control (like I said, though, it's only really the Hydrogen on the OH group of every molecule that's moving around significantly). Energy is not conserved. The temperature increases apparently exponentially from the beginning of the simulation, and pressure varies wildly. Using this last example, which is identical to my previous script except for scaling the temperature to 0K before minimization, running NVE instead of NVT, and still using 0.25 fs timesteps, I got the following outputs for temperature and pressure every 200 timesteps:

step 0: T=300 K, P=1398 atm TotEng=-20 kcal/mol
step 200: T=316 K, p=-616 atm TotEng=314 kcal/mol
step 400: T=631 K, p=-15,615 atm TotEng=3,813 kcal/mol
step 600: T=15,082 K, p=20,232 atm TotEng=48,194 kcal/mol
step 800: T=124,926 K, p=-139,963 atm TotEng=459,193 kcal/mol
Followed by an error (Dump dcd of non-matching # of atoms).

I hope this helps. It really looks like that one hydrogen is the problem, but I'm out of troubleshooting ideas. I can't figure out what else might be wrong. Thanks for any help/suggestions in advance!

Best,
N.H.

Thank you for your reply, Steve. I tried reducing the temperature before minimization and "successfully" minimized a 4X4X4 box of TTA, which I had been unable to do previously. I put "successfully" in quotes because, while it didn't fail, the total energy barely changed (from -583.3 kcal/mol to -591.7 kcal/mol), but LAMMPS used far fewer than the maximum allowed timesteps for minimization. In addition, the resulting minimized structure was again completely symmetrical. In fact, the behavior of all 64 molecules was almost identical for the first few hundred timesteps. I think that's a result of dropping the temperature and then scaling the velocities back up again after minimization.

I did try running NVE instead of NVT as well, but the result was pretty much the same. The timestep I used was 0.25 fs, and making it smaller did not affect the simulation, so I assume it's fine. Particles do not "jump" across the screen at the beginning of the simulation (when viewed in VMD), but move smoothly, until it gets out of control (like I said, though, it's only really the Hydrogen on the OH group of every molecule that's moving around significantly). Energy is not conserved. The temperature increases apparently exponentially from the beginning of the simulation, and pressure varies wildly. Using this last example, which is identical to my previous script except for scaling the temperature to 0K before minimization, running NVE instead of NVT, and still using 0.25 fs timesteps, I got the following outputs for temperature and pressure every 200 timesteps:

step 0: T=300 K, P=1398 atm TotEng=-20 kcal/mol
step 200: T=316 K, p=-616 atm TotEng=314 kcal/mol
step 400: T=631 K, p=-15,615 atm TotEng=3,813 kcal/mol
step 600: T=15,082 K, p=20,232 atm TotEng=48,194 kcal/mol
step 800: T=124,926 K, p=-139,963 atm TotEng=459,193 kcal/mol
Followed by an error (Dump dcd of non-matching # of atoms).

I hope this helps. It really looks like that one hydrogen is the problem, but I'm out of troubleshooting ideas. I can't figure out what else might be wrong. Thanks for any help/suggestions in advance!

i already said that the problem is consistent with your unreasonably
short LJ cutoff. this allows the point charges to get too close to
each other and at 0.2 \AA your time step is far too large and the LJ
potential far too steep so that your atoms will shoot through the
system. since the hydrogen atoms don't have an LJ sphere around them,
they are particularly vulnerable. normally the LJ repulsion will keep
the heavier atoms about 2-3 \AA apart, which leaves ample room for the
point charges on the hydrogens to not get too close to other point
charges and create a "coulomb catastrophe". however, with a cutoff of
0.2 \AA this is not given and they can get arbitrarily close. it is
just a matter of chance whether they hit a heavy atom or not.

change the pair_style command to lj/cut/coul/long 9.0
and you'll be fine. well, 9.0 angstrom LJ cutoff used to be a
reasonably choice many, many moons ago. these days 12 to 15 \AA is
considered a reasonable compromise between accuracy and speed. you
probably also want to improve kspace accuracy and use pppm instead of
ewald, especially for your larger systems.

that is what concerns the unusual (or not) acceleration. there are two
other drastic problems in your input script:

1) you set the neighbor list update to be far too infrequent for your
system with
neigh_modify every 10 delay 20 check yes
this will create loads of "dangerous builds" atoms will move far too
far between reneighboring leading to very inconsistent force
computations. it is good practice to check your output for this and
adjust neighbor list settings to have no dangerous builds.

2) you time integrate atoms twice when adding the nvt fix. LAMMPS
warns you about this, too.

axel.