Subject: chloroform solvent simulation difficulties

From: Chas Simpson <[email-masked]>
Date: Mon, 9 Oct 2000 11:43:12 +0200 (SAST)

Hi All,

I am trying to simulate (with MD) a chloroform solution using periodic
boundary conditions in an NVE ensemble. The simulation of chloroform is
well documented in literature but I am still unable to get it to stablise
correctly.  I heat the system up to 288K from 200K and equilibrate in an
NVT ensemlbe at 288K with no problems.  When I try and equilibrate in an
NVE ensemble, the temperature jumps to about 350K no matter what I try. (I
am running these simulations using CHARMM and using SHAKE to constrain the
motion of the hydrogen atoms).

My first question is how can I reduce the temperature in an NVE
simulation, ie not just transfer the kinetic energy to potential energy
but to actually reduce the kinetic energy term..... if that's
possible/legal(!!)? Or alternatively is there something obvious that I am
doing wrong?

My second question is I am using particle mesh Ewald summations to treat
the long-range nonbonded interactions. Charmm complains and says that I
must specify the frequency for stopping the translation and rotation of
the molecules during the simulation. Why is this necessary and what
frequency should be used? And finally charmm automatically stops the
translation and rotation of the molecules during heating. Why is this
necessary?

Thanking you in advance for any help (which will be GREATLY
appreciated!!). 

Kind regards

Chas

[email-masked]

---

From: Michael Brunsteiner <[email-masked]>
Date: Mon, 09 Oct 2000 13:11:27 +0100

Chas Simpson wrote:
> NVE ensemble, the temperature jumps to about 350K no matter what I try.
> (I am running these simulations using CHARMM and using SHAKE to constrain the
> simulation, 

in an NVE simluation the temperature should depend (in the long run) on
the 
density only (and on the interaction potential parameters of course, but
i assume
you use the standard non-bonded parameters provided by charmm anyway
...)

> ie not just transfer the kinetic energy to potential energy
> possible/legal(!!)?

you are not supposed to/cannot do that ... the ratio between kinetic and
potential
energy is a constant for a given system in equilibrium (get yourself a 
textbook on statistical thermodynamics (MacQuarrie is good) and look for 
equipartition, energy ...)

> Or alternatively is there something obvious that I am
> doing wrong?

no ... send me you input script ...

> translation and rotation of the molecules during heating. Why is this

it is ALWAYS good to do this in the beginning of a simulation ...
otherwise your thermostat may end up measuring the wrong temperature
(apart from other more sublime problems) ... imagine your system
uniformly 
translating through (periodic) space at a constant velocity in addition
to 
the thermal motions .... your calculated kinetic energy would be rubbish 
... (maybe thats also the reason for the temperature problem you
describe 
above ...) concerning the frequency ... usually it should do to do it
only once
(i.e. in the beginning) ... so for a simulation with 100000 timesteps
the frquency would be 1/100000 ... but in any case you should check 
whether charmm really does it ... it should say so in the output ...

mic

==============================================================
Man invented language to satisfy his deep need to complain.
                -- Lily Tomlin

--------------------------------------------------------------
Michael Brunsteiner, Centre for Theoretical and Computational 
Chemistry, University College London
mailto:[email-masked], http://www.ucl.ac.uk/~uccambr/
--------------------------------------------------------------

---

From: Rick Venable <[email-masked]>
Date: Mon, 9 Oct 2000 11:56:08 -0400

> possible/legal(!!)? Or alternatively is there something obvious that I am

I recommend running some NPT simulations first, to ensure that your
chosen density is consistent with the molecular potential.  Also, be
careful about system size; a rule of thumb is that you need 500
particles or more to get good stability.  Your unit cell edge should
probably be on the order of 40 A as opposed to 20 A, for instance.

You should extract and analyze and/or plot the system Hamiltonian for
both your NVT and NVE simulations; for NVE this is ETOT (KE+PE) which
MUST be constant (flat) for a valid NVE simulation.  The Hamiltonian for
NVT depends on the method chosen; it's the term labeled XTLT for NVT
done with the extended system (CPT LEAPFROG), and the Total Energy term
printed in the optional data file (IUNO) for the NOSE method.  Again, it
should be flat, i.e. no apparent upward or downward drift; if not, it
indicates a problem with energy conservation.

As noted, there are a couple different ways of doing NVT simulations in
CHARMM, and they are not particularly straightforward.  An alternative
is to use the more classic equilibration within CHARMM (non-zero
IEQFRQ), using Gaussian assignment of velocities.  I typically use
something like

	ieqfrq 500 iasors 1 ichecw 1 twindh 5.0 twindl -5.0 

This would assign velocities every 500 steps if T is more than 5 deg
away from the target (FINALT), while the other 499 steps are propagated
using standard NVE Verlet.  This won't help if your volume (i.e.
density) is not correct, however.

> translation and rotation of the molecules during heating. Why is this

The net translation and rotation of the entire system are removed, and
not that of individual molecules; the contributions from these 6 global
degrees of freedom are redistributed over the remaining degrees of
freedom within the system.  Because of consistent errors probably
related to shifting charges to a uniform grid for the mesh, it is
necessary to do this regularly when using PM Ewald within CHARMM.  
Otherwise, a published phenomenon termed "the flying ice cube" can
occur, where the entire system is being translated through space.  This
is done automatically when velocities are assigned every IHTFRQ (or
IEQFRQ) steps, but NTRFRQ must be set explicitly for use with PM Ewald.  
As noted in the .doc file, the step frequency should be less than IPRFRQ
to avoid an error exit.  I've used

	ntrfrq 2500 iprfrq 5000

for simulations using PM Ewald electrosatics.

--
Rick Venable                  =====\     |=|    "Eschew Obfuscation"
FDA/CBER Biophysics Lab       |____/     |=|
Bethesda, MD  U.S.A.          |   \    / |=|  ( Not an official statement or
[email-masked]          |    \  /  |=|    position of the FDA; for that,
[email-masked]              \/   |=|    see   http://www.fda.gov  )