Subject: Re: CCL:Velocity re-assignment

From: Rick Venable <[email-masked]>
Date: Sat, 17 Nov 2001 11:55:26 -0500

There are several methods and reasons for modifying MD velocities:

+ periodic assignment during equilibration for NVE dynamics
+ periodic scaling during equilibration for NVE dynamics
+ changed every timestep, coupled to a piston for NVT dynamics

NVE and NVT refer to constant volume ensembles, with explicit conserved
system Hamiltonians; velocities are not changed during NVE dynamics, but
100 ps or more of equilibration with velocity adjustment is usually
needed to get a system w/o large temperature drift.  One should not
include the equilibration time period in the calculation of molecular
properties.

The use of a piston to "extend" the system (adding degrees of freedom)
is often used to provide true NVT dynamics; the potential and kinetic
energy of the piston are included in the system Hamiltonian.  The same
method can be applied to NPT (const T,P) and NPH (const P) ensembles in
a rigorous way, from the viewpoint of statistical mechanics.

I should note that a method termed "weak coupling", popularized by HC
Berendsen's implementation, has also been used for constant T
simulations; a theoretical issue with that approach is that the ensemble
isn't known, and hence neither is the Hamiltonian.  The extended system
approach is preferred, due to the rigorous stat mech ensemble.  There's
a more complete discussion of these issues in:

"Computer Simulation of Liquids", MP Allen and DJ Tildesley, Clarendon
Press, Oxford UK (1987; pbk 1989,1990).

Feller, Zhang, Pastor & Brooks, JCP, 103, 4613 (1995)

The JCP paper discusses the extended system implementation for constant
T and/or constant P in CHARMM, and compares it to the weak coupling
algorithm.

On Fri, 16 Nov 2001, Ibrahim Moustafa wrote:
>     I'm doing MD on my protein structure; in some recent work on
> simulation using MD at const. temp. I've seen people re-assign velocities
> every certain Intervals.
> 
>     I want to know excactly how to re-assign velocities periodically?
>     Is it by changing the seed value randomly periodically? or they use
>     another way to do that?
> 
>      What are the benefits for the simulation from velocity re-assignment?
>      Doesn't this drift the trajectory everytime velocities being
>      reassigned?? Is this good for simulation?
> 
>     I'd be appreciated if anyone calrify this point.  
> 
>      I use Tripos software to do the simulation.

=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=
Rick Venable
FDA/CBER/OVRR Biophysics Lab
1401 Rockville Pike    HFM-419
Rockville, MD  20852-1448  U.S.A.
(301) 496-1905   [email-masked]
ALT email:  [email-masked]
=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=

---

From: David van der Spoel <[email-masked]>
Date: Sat, 17 Nov 2001 20:17:26 +0100 (CET)

On Sat, 17 Nov 2001, Rick Venable wrote:

Without wanting to start a holy war...
The Berendsen coupling does not give a well known stat mech ensemble for
intermediate values of the coupling parameter tau_T, however in the
limiting cases (tau_T -> 0, or tau_T -> oo) it is identical to well known
ensembles. The clear advatange of the algorithm is that it works with
minimal effect on the dynamics of the system (velocity modification < 1%).

>+ changed every timestep, coupled to a piston for NVT dynamics
These algorithme do influence the dynamics of the system considerably
because the trajectories are not continuous. Besides, they dont yield a
known stat mech ensemble of your system, but of your system plus a
piston... So you can ask yourself what does that mean...

The real thing to do is avoid having to use a temperature scaling
algorithm. Equilibrate at the temperature you need using any algorithm,
then take care to conserve energy, by avoiding cut-offs. In periodic
systems you will have to use Ewald summation or PME for the Coulomb
interactions, and a smooth (long) cut-off for the non bondeds.
In non periodic systems you will have to include all interactions.

Groeten, David.
________________________________________________________________________
Dr. David van der Spoel, 	Biomedical center, Dept. of Biochemistry
Husargatan 3, Box 576,  	75123 Uppsala, Sweden
phone:	46 18 471 4205		fax: 46 18 511 755
[email-masked]	[email-masked]   http://zorn.bmc.uu.se/~spoel
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

---

From: Rick Venable <[email-masked]>
Date: Sat, 17 Nov 2001 14:50:36 -0500

On Sat, 17 Nov 2001, David van der Spoel wrote:
> The Berendsen coupling does not give a well known stat mech ensemble
> for intermediate values of the coupling parameter tau_T, however in
> the limiting cases (tau_T -> 0, or tau_T -> oo) it is identical to
> well known ensembles. The clear advatange of the algorithm is that
> it works with minimal effect on the dynamics of the system (velocity
> modification < 1%).

However, one does not typically use values close to zero or infinity, so
in general, the ensemble isn't known.  By taking care with the piston
parameters, the extended systems also have minimal effect on the
dynamics of the system.  The early work by people such as HC Andersen
and WG Hoover, and later studies by ML Klein and others have shown this.

> piston... So you can ask yourself what does that mean...

As I stated clearly and you have omitted, it is not appropriate to
include the equlibration period when computing properties from MD
simulations.  The problem is, as you indicate, that the periodic 
velocity adjustment makes the dynamics discontinuous.

The use of a piston for const T does correspond to the NVT ensemble; or
are you claiming that Allen, Tildesely, Andersen, Hoover, Klein, etc.
are wrong?

> In non periodic systems you will have to include all interactions.

I agree about equilibration (and said so in my orignal message), and
also about the use of Ewald methods, especially for true periodic
systems.  However, a protein in a cube (etc.) of water with periodic
images isn't really a true periodic system, because the effective
concentration is much higher than typical solution experimantal work.  
In most cases, one doesn't want the field effect of the protein in the
adjoining image cell, and a cutoff may be needed to exclude it.

---

From: Ibrahim Moustafa <[email-masked]>
Date: Sat, 17 Nov 2001 21:17:14 +0000 (GMT)

The debate on which algorithm is better in MD (NTV) is still going;
 and the use of the extended method is more rigorous stast. mech.

    Actually, I have to use weak coupling method coz this algorithm is
already implemented in SYBYL while the extended algorithm is not.

    I still need to know: given that i use weak coupling algorithm for
MD-NVT with coupling factor 0.1 ps, how can I re-assign the velocity in a
proper way such that not biasing the simulation results; is that by
changing the seed value periodically by entering different seed values in
the input for the simulation? if so, does this cause any harm to the
simulation??!!

      Thanks for all contributed to calarify this point.

Name         :Ibrahim M.Moustafa
Mail address :Center for Biomolecular Science,
              BIOMOLECULAR SCIENCE BUILDING,
              North Haugh,St-Andrews,
              Fife,KY16 9ST
              Scotland,U.K.
         Tel :+44(0)1334-467257
      E-mail :[email-masked]

On Sat, 17 Nov 2001, Rick Venable wrote:

> -= This is automatically added to each message by mailing script =-
> [email-masked] -- To Everybody  | [email-masked] -- To Admins
> [email-masked] -- HELP CHEMISTRY or HELP SEARCH
> [email-masked] -- archive search    |    Gopher: gopher.ccl.net 70
> Ftp: ftp.ccl.net  |  WWW: http://www.ccl.net/chemistry/   | Jan: [email-masked]
> 
> 
> 
> 
> 
>

---

From: David van der Spoel <[email-masked]>
Date: Sat, 17 Nov 2001 22:54:57 +0100 (CET)

I was basically recalling what is in my understanding in this paper:
>@Article{Morishita2000a,
>  author =      {T. Morishita},
>  title =       {Fluctuation formulas in molecular-dynamics simulations
>with the weak coupling heat bath},
>  journal =     {J. Chem. Phys.},
>  year =        2000,
>  volume =      113,
>  pages =       {2976-2982}
>}

I am not claiming that the older papers by Klein, Hoover etc. are wrong,
but if you use a piston to keep the temperature this becomes part of the
system, and hence you are not simulating the system you think you do. As
far as I know most coupling methods influence the temperature 
fluctuations (necessary for computing heat capacities etc.), which can
be quite easily seen by taking the fourier transform of the temperature
fluctuations. These usually have a peak corresponding to the coupling
constant of your piston.

>are you claiming that Allen, Tildesely, Andersen, Hoover, Klein, etc.
See above.

>adjoining image cell, and a cutoff may be needed to exclude it.
One should always ask: "how do I minimise the artefacts on my system of
interest." PME I think is less bad than cut-offs... One day when I am
bored I will do the simulations with PME with ever increasing box size,
and check the interaction between the protein and its image. See how far I
have to go to simulate a protein at infinite dilution.

On a side note, under physiological conditions (i.e. in a cell) the
protein density is quite high. So why not simulate three or four proteins
in a single simulation box with PME? The different proteins will average
out any serious "crystal" artefacts due to PME and you can write four
papers based on one extremely realistic simulation...

Groeten, David.