Strange behavior of global scalar

Hello,
  I am trying to run a simulation of 3-site acetonitrile using the TRAPPE potential with the fix rigid/npt command. In order to check that the simulation is working, I checked the Lammps manual for the conserved quantity, which I believe should be the sum of the global scalar and the potential energy of my system. I printed both, and while the potential energy quickly levels off, the global scalar does not. It is also multiple orders of magnitude larger than the potential energy of the system (~-1400 kcal/mol for PE and ~700000 kcal/mol for f_1), which I believe should not be the case. The strange thing is that all other quantities that I print seem to level off relatively quickly, just the global scalar does not. I have attached the input, and dada files as well as plots of the potential energy, global scalar (both plotted from 0.4ns - 5ns, since the first 0.4 skew the graph and make it hard to see the increase), and density (plotted for the full 5ns). Any suggestions on what I might be doing wrong (or if everything is all right and I am just interpreting this wrong) would be appreciated.
Thanks,
  Tal Aharon

data.aceto3site (90.9 KB)

in.hirata (2.99 KB)

density.png

global_scalar.png

pe.png

Hello,
I am trying to run a simulation of 3-site acetonitrile using the TRAPPE potential with the fix rigid/npt command. In order to check that the simulation is working, I checked the Lammps manual for the conserved quantity, which I believe should be the sum of the global scalar and the potential energy of my system. I printed both, and while the potential energy quickly levels off, the global scalar does not. It is also multiple orders of magnitude larger than the potential energy of the system (~-1400 kcal/mol for PE and ~700000 kcal/mol for f_1), which I believe should not be the case. The strange thing is that all other quantities that I print seem to level off relatively quickly, just the global scalar does not. I have attached the input, and dada files as well as plots of the potential energy, global scalar (both plotted from 0.4ns - 5ns, since the first 0.4 skew the graph and make it hard to see the increase), and density (plotted for the full 5ns). Any suggestions on what I might be doing wrong (or if everything is all right and I am just interpreting this wrong) would be appreciated.

​you should check, if your simulation can conserve energy with fix rigid/nve or fix rigid or fix rigid/small without any thermostats (after initial equilibration). it may be, that your simulation is not conserving energy and then the thermostat would be compensating ​it and thus you’d have a continuous transfer of energy into the heat bath.

​axel.​

Dear Dr. Kohlmeyer,
I apologize, after further discussion of my data, I was informed that I had been looking at the wrong value to begin with. Some others who i’ve been working with suggested that the conserved quantity is the total energy via the sum of pe and ke, which should not include contributions from the thermostat and barostat. Indeed, a plot of the sum of the two shows that it does in fact stabilize very quickly. As I am new to running MD simulations, I am not sure what is the correct way to tell if my simulation has run correctly, but the lack of drift in the total energy does seem promising. Thank you for your time.
Sincerely,
Tal Aharon

Dear Dr. Kohlmeyer,
I apologize, after further discussion of my data, I was informed that I had been looking at the wrong value to begin with. Some others who i’ve been working with suggested that the conserved quantity is the total energy via the sum of pe and ke, which should not include contributions from the thermostat and barostat.

that is not correct. it only applies to using fix nve without a thermostat.

Indeed, a plot of the sum of the two shows that it does in fact stabilize very quickly. As I am new to running MD simulations, I am not sure what is the correct way to tell if my simulation has run correctly, but the lack of drift in the total energy does seem promising. Thank you for your time.

no, it does not. like it mentioned before, this may be just hiding an
incorrect time integration by removing the drift instantaneous by the
thermostat.

the conserved quantity for an NVT (not NVE) ensemble is the total
potential energy plus the total kinetic energy plus the energy
differential transferred into or from the heat bath.

axel.

PS.: let me repeat how important it is to do a test for energy
conservation without any thermostat. that is the test you should do
first (after equilibrating the system, where energy conservation is
irrelevant). that will allow you to find good and suitable settings
from length of timestep, cufoffs, neighborlist updates and more. only
after this is sorted out, you should look at how total energy is
conserved when including a thermostat/barostat. otherwise the
thermostat/barostat can hide bad simulations settings and you may end
of with bogus trajectories without noticing.

axel.

Dear Dr. Kohlmeyer,
  I apologize for my the delay since my last response.
  I have recently acquired the results of an NVE simulation for the previously mentioned system. I ran a 1 ns nve equilibration calculation, and then took the coordinates from the end of that to run a 5ns nve calculation. I’ve attached the data file, input file, an image of the plotted total energy (from the etotal call in thermo_style), as well as a plot of the total energy without the first 50 steps, so I could check for drift. The plotted fit shows the slope of the data without the first 50 steps. It seems that with no thermostats, the system energy is conserved, but that does not explain why when I run with a thermostat I get such strange behavior of the thermostat energy (as you may recall from before, it continuously increases). Thank you for your time and assistance.
Sincerely,
  Tal Aharon

data.aceto3site (90.9 KB)

Total_energy.png

Drift_Slope.png

in.hirata (1.07 KB)

ATT00001.txt (4.76 KB)