AMBER Archive (2008)

Subject: Re: AMBER: What's making my system lose total energy?

From: Thomas Cheatham (tec3_at_utah.edu)
Date: Mon Mar 17 2008 - 22:36:30 CDT

> NSTEP = 100 TIME(PS) = 25182.850 TEMP(K) = 307.65 PRESS = 0.0
> Etot = 4401.8340 EKtot = 2659.7186 EPtot = 1742.1154
> BOND = 417.5512 ANGLE = 1388.9901 DIHED = 601.5975
> 1-4 NB = 363.9760 1-4 EEL = 16825.8254 VDWAALS = -2246.5884
> EELEC = -15609.2364 EHBOND = 0.0000 RESTRAINT = 0.0000
> Ewald error estimate: 0.7533E-05
>
> ...
>
> NSTEP = 106400 TIME(PS) = 25289.150 TEMP(K) = 296.16 PRESS = 0.0
> Etot = 4398.9124 EKtot = 2560.4280 EPtot = 1838.4844
> BOND = 407.1341 ANGLE = 1415.2345 DIHED = 607.2740
> 1-4 NB = 348.7719 1-4 EEL = 16857.0617 VDWAALS = -2249.7184
> EELEC = -15547.2735 EHBOND = 0.0000 RESTRAINT = 0.0000
> Ewald error estimate: 0.2851E-04

In addition to SHAKE tolerance and worries about NSCM as mentioned by Dr.
Walker, it is not readily possible from snapshots of the energy to
understand what is happening, at least for me, especially over a 100 ps
window... I see a large rise in potential energy with increases in some
components and decreases in others. What I would recommend is looking at
ns or longer averages or the time course over the whole simulation.
While you are at it, I also suggest running with constant pressure, but a
coupling time that is longer than the simulation; this would help see what
is happening to pressure in the system. This is not a particularly large
system so large fluctuations are expected as the fluctuations are larger
as the system size gets smaller. Also, for a small system as this, you
can bump up the Ewald; set grid size < 1 A (power of two) and reduce
DSUM_TOL by an order of magnitude or so.

What I would also check is if you have any large force constants for
bonds/angles that may require a smaller time step.

> stable levels of potential energy. That's when I found that it was this issue
> of losing TOTAL energy all along. Can anyone help?

Finally, are you convinced that the system is in equilibrium and has
settled down? Although TIP3P water relaxes in 50-100 ps, this does not
mean that all systems do. Proteins, for example, take considerably
longer... What I would do is try the same protocol with a box of water
and see if this shows the same problems; if not, this suggests things are
OK. If so, perhaps more digging is necessary (or force constants need to
be lowered). This is one reason why I advocate NPT with (very) weak
coupling; I believe we can conserve energy so like the ability to modulate
the temp/pressure for things out of equilibrium. However, I do not want
to start a religious war on ensembles, conservation or coupling.

-- tec3

-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu