Subject: Summary: scaling velocities in MD

From: [email-masked] (Teerakiat Kerdchaoren)
Date: Tue, 7 Dec 1993 18:18:27 +0100 (NFT)

Dear computational chemistry mailing-list members,

	A couple of weeks ago, I post a question about scaling 
velocities in MD to the net. The original messages and
the responses are herewith enclosed. I would like to express
my special thanks to them whom sent me those replies. Indeed,
I found all of them valuable and helpful.
	Wish you nice X'mas.

					  	warm regards,	

						Teerakiat
						[email-masked]

---------------------- original messages -----------------------
 
Dear CCL fellows, 
	I am writing an MD program and here comes some questions.
	To my experience, the average temperature tends to raise
up gradually in a simulation resulting in the increasing of total
energy. 
        1) Is it possible to apply a "scaling temperature" algorithm
to control the temperature range ( my MD is NVE ) ? To my idea,
scaling the velocities cause the deviation of the p-q coordinate 
( generalized momentum-generalized coordinates ) in the phase space.
The trajectories are therefore deviated, is it true ? How to scale
both velocities and coordinates to keep the correct trajectories and
p-q coordinates in the phase space ?
	2) Is it possible to apply an "external bath" in which is 
generally used with NVT ensemble to an NVE ensemble ?
	
	Thanks a lot in advance.
						best wishes,
						T.K.Noy
------------------------  response -------------------------------

From: "Michael E. Coltrin" <[email-masked]>

We published a method for simulating the effects of a thermal
reservoir on lattice atoms. The method ensures the proper flow
of energy between the dyanamical atoms in your trajectory study
and a surrounding bath (if the bath is at a different temperature).
The method also will ensure that the total energy will not "drift"
as you described.  

Briefly, the method periodically resets a portion of the velocity
vector of each atom at the edge of the dynamical zone (but doesn't
reset the position). The method is very easy to implement in
a trajectory code.

The reference is M. E. Riley, M. E. Coltrin, and D. J. Diestler,
J. Chem. Phys., 88, 5934 (1988).

I hope this will be useful to you.

Mike Coltrin
-------------------------------------------------------------------------

From: [email-masked] (Craig W. Burkhart)

You may want to try the canonical ensemble dynamics of Nose' and Hoover.
Their method is inherently isothermal, whereas the standard methods
are isokinetic. That is a subtle, but important, difference.

The appropriate refs are:

S. Nose', J. Chem. Phys. 81, 511(1984).
W.H. Hoover, Phys. Rev. A 31, 1695(1985).

Hope this helps.

Craig W. Burkhart, Ph.D                
--------------------------------------------------------------------------

From: [email-masked] (Anatoli Belonoshko)

   I have written quite a few MD programs.
   To my experience, the average temperature does not tend to rise
 gradually, unless you
 1) have not "equilibrated" your system properly
 2) use inappropriate numerical scheme solving Hamilton equation
 3) use too large timestep
 4) still have some other bugs

    Usual practice for NVE ensemble is to equilibrate system scaling
 velocities to maintain the temperature, and then to stop the scaling.
 NVE is rather famous on inability to provide "exact" temperature.
 On my experience, the scaling during "averaging" stage does not cause
 sigificant errors, but it is better still to avoid such a scaling.

    "external bath" technique is applicable only to NVT ensemble by
 definition. Constant energy simulation can not be realized when
 system have an energy external source.

    Best regards,

    Anatoly Belonoshko
    Uppsala University, Sweden
    e-mail: [email-masked]

--------------------------------------------------------------------------
From: [email-masked] (K.V. Damodaran)

Yes, the average temperature tends to rise in the beginning of the 
simulation. BUT, usually the starting configuration has a very large 
potential energy and the increase in temperature is due to the strong 
forces, a consequence of which is the sharp decrease in potential 
energy. Thus in NEV simulations the total energy will not be conserved
very well in this region. But this part of the trajectory is discarded
anyway.

Velocity scaling applied for this part of the trajectory scales down
the velocities by SQRT(T-To), when the temperature goes beyond some
tolerance. NEV simulations do not require velocity scaling when you
actually collect the data. So you probably don't have to be concerned
about the deviations of the p-q coordinate in the phase space. Scaling
the coordinates will amount to using an external pressure bath.

If you apply an external temp. bath to an NVE ensemble, it will become
an NVT ensemble!! See the book by Allen & Tildesley for details.

damodaran
[email-masked]
------------------------------------------------------------------------
From: [email-masked] (David C. Doherty)
Subject: Re: Molecular Dynamics questions
Status: RO

>        2) Is it possible to apply an "external bath" in which is 
>generally used with NVT ensemble to an NVE ensemble ?

A method that many use is found in Berendsen et al., 
J. Chem. Phys., 81, 3684-3690 (1984).  This includes a 
pressure bath also.

---
David C. Doherty
Minnesota Supercomputer Center, Inc.
[email-masked]
-------------------------------------------------------------------------
If you are running an NVE simulation, then the temperature should
not rise.  This indicates problems with your simulations.  Check
if the forces are being calculated correctly.  Also, if there are
no other problems, your time step may be too large.  See if this
systematic change in the energy persists even if you reduce the
time step.  You could add a temperature controller but you should
first make sure your simulation runs correctly in NVE.  Hope this
helps.

Dilip Paithankar
Chemical Engineering, Purdue University, West Lafayette, IN 47907
Telephone: 317/495-6910 (home), 317/494-4055 (office)
           317/494-0805 (FAX).
-------------------------------------------------------------------
I am writing a MD-program myself and I encountered the same problem.
The only solution I found was to restrict the run to a sensible length
(say: energy conservation 1 : 1.0d(6-7)).

I would be interested in the other suggestions you get.
Please report to the net.

Yours sincerely,

Ferenc Molnar
EMail (SMTP):  [email-masked]
               [email-masked]
               [email-masked]
---------------------------------------------------------------------------
From: [email-masked]

	One of the really nice features of the molecular dynamics technique is
that there is usually a "conserved quantity" (such as the Hamiltonian) which
makes it relatively easy to detect if you have an error in your code.  In your
case,
if you are carrying out an NVE calculation, then the total energy MUST be
conserved.
The fact that you do not observe this suggests one of two things to me:
	(1)  You have a bug in your code.   This may be something as trivial as
an actual
		coding error, or something more subtle, such as having too
great a
		distance between points in a potential look-up table.
	(2)  The system you are trying to simulate is in an unfavorable
configuration,
		resulting in large forces on the particles.  If this is the
case, you will
		need very small time steps in order to obtain total energy
conservation.
The situation described in (2) should disappear if the system is allowed to run
on for
a while.  From your problem description, I suspect (1) is likely to be your
problem.
I suggest you use the conservation of total energy to ensure that your code is
really
correct, rather than to try cover up the problem with some sort of velocity
scaling.

Happy Bug Hunting!

	Dr. James P. Bareman
	Member of Research Staff
	Xerox Research Centre of Canada.
-------------------------------------------------------------------------
From: JEREMY WALTON <[email-masked]>

The traditional way to get the NVE simulation running at (close to) the 
temperature you want is to have an equilibration period where the velocities
are all uniformly scaled so the system is forced to have the kinetic energy
appropriate to the required temperature.  Following this period (typically
a few thousand timesteps), scaling is switched off and the averaging over
configurations can be started.  The temperature should be monitored during this
averaging period (it's one of the things you calculate as an average); if it
drifts too much, then maybe the equilibration period wasn't long enough, or
maybe your timestep is too big.

There are a number of algorithms for this.  Two which I recall are that of 
Berendsen and that of Nose.  Check out Allen & Tildesley's "Computer Simulation
of Liquids" (OUP, 1987), where all this is described in some detail.

Hope this helps,

Jeremy
----------------------------------------------------------------------------
From: <[email-masked]>

I have a comment that is a "tangental" answer. In fitting several 
variable temperature NMR simulations where you have an equilbrium
situation between several species of various populations, I noted 
that I could not obtain satisfactory fits of my spectra by simply 
increasing the rate of exchange between the various species (as 
dictated by the Eyring equation) until I began to increase the 
populations of the minor species with increasing temperature. 

Further, the best fits were obtained when I could track the 
popolations of each species at several temperatures and generate
a Van't Hoff plot that could be used to "predict" the populations
for spectra where overlap prevented population elucidation. 

If your MD model is varying in temperature, a more realistic 
model would include an increasing of minor conformation populations
with increasing temperature (as dictated by the Van't Hoff equation).

What do others think.

Jay Brown
----------------------------------------------------------------------
From: "Gerald Loeffler" <[email-masked]>

there is a 'classical' paper by by Berendsen et. al. that served as
a guideline for me, when I implemented temperature-scaling in my MD-program:
	Berendsen,Postma,van Gunsteren,DiNola,Haak:"Molecular dynamics
	with coupling to an external bath",J.Chem.Phys.,Vol.81,No.8,1984

	good luck,
	Gerald
-----------------------------------------------------------------------
From: Rob Hooft <[email-masked]>

In general, if your E rises in a NVE you should try to

 1) reduce the time step.
 2) If a cutoff radius is used, increase it.

Rising energy in NVE is due to problems with integration accuracy, so
the only real way to decrease the rising is to improve the integration
accuracy. 
------------------------------------------------------------------------
From: Tom Connor Bishop <[email-masked]>

i'm sure you've had more than one person tell you that
something must be wrong if you're total energy increases.
this is the standard test to see if the simlation is accuarate 
ie  the law of conservation of energy must hold.
often times temperatures will increase but this must be coupled with a
decrease in potential energy!!!
either something is wrong with your code, or you are not running the simulation
with parameters appropriate to the system. (note that for non-equilib systems  
the time step must be shorter in order to obtain numberical stability)

just my 
-------------------------------------------------------------------------