Subject: P convergence in NPT MD runs - summary

From: "Gustavo A. Mercier Jr" <[email-masked]>
Date: Mon, 11 Nov 1996 16:28:05 -0500 (EST)

Hi, CClers!

A few weeks ago I posted a message regarding convergence
on the pressure in NPT simulations. The
answers I received are below.

In a nutshell, the pressure will not converge as tightly as
the temperature. As explained by D. Case (and is detailed in the
amber web paged mentioned below), the fluctuations will be large,
ca. +/-100-150bars. Mean values close to +/-(1-10) bars can be
obtained over long periods, ca. 100 ps.

Bye!

***** Answers *******

*******************
If I understand correctly, what you are seeing is what is expected.  At
equilibrium,  the "instantaneous" pressure will have what may seem like
large fluctuations.  These will never go away, and are related to the
compressibility of the solvent.

This is in the AMBER FAQ list: see

   http://www.amber.ucsf.edu/amber/amber.html

go to "answers to commonly asked questions", and look for
  "constant pressure".  This contains a disucssion of the behavior that
one should expect, which sounds not too far from what you are seeing.

 ....hope this helps...

 ...dave case
**********
Hi Gustavo, 

        I have also had some problems like yours. I solved it by the following 
manner:

        1.- I slightly increased the temperature for a few picoseconds, with 
the followinfg coupling constants: temperature 0.05  and pressure 0.5 ps.

        2.- After this, I reduced the temperature to the original one, and
I used a coupling constants of 0.1 and 0.5 for temperature and pressure
respectivelly. 

        If this does not work, I advise you to run a few picoseconds at 
constant volumen, followed of the steps 1 and 2 such as was described above.

Good luck,

Javier

*******************************************************************************
Dr. Jose Javier Lopez Cascales     | Office phone:Int-34-68-307100,ext:2190   *
Departamento de Quimica-Fisica     | Fax:    Int-34-68-364148                 *
Fac. de Quimica, Uni. de Murcia,   | e-mail: [email-masked]                 *
Campus de Espinardo, 30071 Murcia, | www -> http://leonardo.fcu.um.es/~javier/*
Spain                              |                                          *
*******************************************************************************
*************
Gustavo,

As for the accuracy of the pressure calculation, it is always
difficult to get very accurate value, because pressure is
calculated by the virial theorem and it is calculated by
subtracting a big value from a big value. The difference
between the attractive and repulsive potential curve determines
the value and makes the value not accurate. Moreover, you are
using NPT ensemble and the particle system interacts with one
dynamical system which represents external pressure. Therefore,
simple oscillation may appear in the monitoring values. I think
that if you add more non-linearity into the interaction between
particle system and external pressure system, the equilibriation
will be obtained more quickly.

Best Regards,

-----------------------------------------------------------------------
Masa Katagiri
Queen's University of Belfast                                          
Atomistic Simulation Group               http://titus.phy.qub.ac.uk
School of Mathematics and Physics        Tel:+44 1232 273528
Belfast BT7 1NN                          Fax:+44 1232 241958
Northern Ireland                    
                                  E-mail:[email-masked]
                                         [email-masked]
-----------------------------------------------------------------------
********
Hi Gustavo,

    I have observed the same kind of behaviour. Pressure does converge more slowly
and fluctuates more widely than Temperature. You could try tighter coupling; i.e.
a smaller TAUP value. Note, though, tighter coupling for the Temperature (TAUTP)
can lead to anomalous MD; I don't know whether the same thing can happen when the
TAUP coupling constant is made very small.

Cheers,

Vidana C. Epa    

Biomolecular Research Institute,
343 Royal Parade,
Parkville, Vic. 3052,
AUSTRALIA.
Tel:   (61) - 3 - 9342 - 4300
Fax:   (61) - 3 - 9342 - 4301
email: [email-masked] 
********
Hi Gustav,

as for the pressure fluctuations, if I were you I should not worry to much
about that. The pressure is known to fluctuate enormously, and only the
pressures averaged over sufficiently long periods (10-100ps) will have
moderate fluctuations. So nothing is wrong in the system, except that if you
have to tight coupling (and I think 0.1ps is already below the limit; I'd
recommend 0.4ps for the pressure as the lower limit), you will not generate
boltzmann-weigthed trajectories. The same holds for the T-coupling; do not
couple to tightly. 

Hope this is of some help,

Andrepeter Heiner.
=================================================================
| Andreas P. Heiner          | voice : +358-9-456 5105          |
| VTT/Biotechnology and      | fax   : +358-9-455 2103          |
|     Food Research          | e-mail: [email-masked] |
| Biologinkuja 1,            |                                  |
| PO Box 1500, Espoo         |                                  |
| FIN-02044 VTT Finland      |                                  | 
=================================================================

-------------------------------------------------------------------------

Use a longer coupling constant for the pressure (0.5 or 1.0 ps)
The fluctuations in the pressure are also dependent on the way the virial
is computed: based on atoms, or based on molecules. The methods are in
principle the same and do give the same averages. The atom based virial
gives larger fluctuations. In our GROMACS code we have atom based virial
and use at least 1.0 ps for the pressure coupling.

David van der Spoel
---------------------------------------------------------
EMAIL:	[email-masked]
WWW:	http://rugmd0.chem.rug.nl/~spoel
PHONE:	31-50-3634327	FAX: 31-50-3634800
MAIL:	Nijenborgh 4, 9747 AG Groningen, The Netherlands.
---------------------------------------------------------

---------------------------------------------------------

The pressure will never converge as you expect - please
see amber41/Questions/pressure, or for the latest version,
the version on our web site

	http://www.amber.ucsf.edu/amber/

Bill Ross

---------------------------------------------------------

Gustavo,
	It is normal for the instananeous pressure to fluctuate wildly 
when one calculates the pressure from the "atomic" virial (there is a 
discussion in Allen and Tildesley about different ways to calculate 
pressure).  If you average over 100 ps or so, then the AVERAGE pressure 
should be within a few atmospheres of your target.
	Before you invest a lot of CPU in this project, you might want to 
read a paper that Bernie Brooks, Rich Pastor and I published in JCP last 
year.  JCP, 103, p 4613 (September 15, 1995 issue).  There are some 
serious problems using the Berendsen method with an inhomogenous system, 
especially is SHAKE is used to maintain rigid water.

Scott

---------------------------------------------------------

   Gustavo,

< personal material deleted :-) >

      As regards your scientific question, I think that big fluctuations
 are a fact of life with constant pressure simulations - we use the 
 same algorithm, and see the same behavior. Have you tried varying the
 coupling factor? I doubt that would accomplish a great deal, but I have
 never done (or seen) a systematic study. I cannot tell you off the top of
 my head if these big fluctuations are theoretically reasonable for a small
 ensemble of particles, or if the algorithm is to blame.

     Regards,

     Randy

All opinions expressed here are mine, not my employer's

///////////////////////////////////////////////////////////////////////// 
\\ Randy J. Zauhar, PhD             | E-mail: [email-masked]        //
\\ Tripos, Inc.                     |       : [email-masked] //
\\ 1699 S. Hanley Rd., Suite 303    |  Phone: (314) 647-1099 Ext. 3382 //
\\ St. Louis, MO 63144              |                                  //
/////////////////////////////////////////////////////////////////////////
**                                                                     **
**  "If you have conceptions of things that you can have no conception **
**   of, then the conception and the thing appear to co-incide."       **
**   --- C.G. Jung                                                     **
*************************************************************************
*************
Hola Gus,
Perdon por la demora por si alguien te contesto primero.
La respuesta facil es:
no hay problema con las fluctuaciones de la presion, y son solo producto de
la compresibilidad del agua (TIP3).

Aca va lo que dicen los expertos sobre el tema:

        I assume that the pressure should stablize at approximately 1.

Except 1) the physics of pressure does not lend itself to instantaneous
evaluation/correction and 2) the algorithm used by Amber (scaling the
box & translating everything, and then applying SHAKE or letting bond,
angle, etc. potentials rectify things) is not the best available. Having
most values in the range -100..100 is the best I ever hope for w/ my 8-10000
atom dna/water systems.

Note by Dave Case:

    -- pressure fluctuations are large because the pressure-volume
isotherm of a liquid is extremely steep: it takes a large pressure
change to yield a small volume change, and conversely: a volume that is
instantaneously slightly out of equilibrium yields a large pressure.
The fluctuations in pressure are then related to this: in a thermodynamic
sense (e.g. Landau & Lifshitz, "Statistical Physics" sec. 114):

             mean value of (delta p)**2 = -kT(dp/dV)

where the derivative is evaluated at constant S (!).    Rearranging the
above equation, approximating (dp/dV) at constant S by 1/(dV/dp) at
constant T, yields

           (delta p)**2 = kT/(V*beta)  or delta p ~= sqrt(1/beta)

where beta is isothermal compressibility and I have set kT/V to 1 bar
in the final equation.  If we use the macroscopic value for beta for
water (roughly 4 x 10**-5 bar**-1) we obtain an rms pressure fluctation
microscopic ideas: pressure fluctuations in very small volumes would
probably be larger than this.

Having said this, it must be noted that the procedure used in Amber
(based on ideas of Berendsen) does not generate a simple thermodynamic
ensemble, and (in addition) should not be expected to yield fluctuations
that match the thermodynamic formula cited above.   It will be interesting
to see if anyone reports a systematic study of these numbers as a function
of the type of simulation being done (especially on the volume of the
simulated system.)

---
************

 ... I have done a large amount of dynamics with periodic boundaries
[in charmm]. I believe you can reasonably use about 1 to 100 atmospheres
without any real perturbation. My experimental work in my job indicates
little pressure effects on most systems of common interest. Besides, the
RMS deviation of the pressure in the system is related to the
compressability of the system and thus there is a finite and large
variation in the internal pressure.  The internal pressure of the system
will be a strong function of the temperature at which the dynamics are
conducted.  One approach is to use a boundary
which gives a water density very close to standard tables for the given
temperature.  I have noticed that sometimes this can be incorrect as the
second phase (or solute) can have too low of a density.  Just recently I have
taken a new approach (at least it is new for me).  Since standard CHARMm
output does not give any idea of the internal pressure or pressure
fluctuations for constant volume dynamics, I have started looking at the
potential energy of the system as the boundary is decreased.  That is,
about 5 ps of dynamics are conducted and the average potential energy at the
end is recorded.  The periodic boundary is then reduced by 0.1 Angstroms and
the dynamics restarted for another 5 ps.  The procedure is repeated for
for several steps of boundary reduction.  One will observe a sharp increase
in potential energy as the system reaches a pressure of about 1 - 10
atmospheres.  This procedure assumes that the initial periodic boundary
was too large for the system being studied.  In addition, I should mention
that the size of my systems range from about 40 Ang on a cubic edge to
60 Ang. on a cubic edge.  You may have to take smaller step sizes than I do.

---
 User dialogue with Dave Case:

> in Sander output files (standard name: mdout and mden) values
> of the pressure are reported. I'm wandering in which unit
> this values are given. I would gess "mbar", but maybe
> somebody knows it exactly. A futher problem with the pressure is
> that I have got negative values. What does negative values mean?
> Is the pressure given as relative pressure?

Pressures are reported in bars; negative pressures indicate that
the system would like to contract its volume, i.e. that an outward
force on the walls of the container would be required to keep the system
at its current volume; positive pressure is the opposite: the system
would like to increase its volume, and an inward-directed external force
would be required to keep the system at its current volume.

> As the reported pressure is the absolute pressure I would suggest to
> interpret the values in the following way:
>
> If the instantaneous pressure is greater than PRES0 the system
> wants to increase is volume (overpressure); if the instantaneous pressure
> is lower than PRES0 but greater than zero the system wants to contract its
> volume (subpressure); an instantaneous pressure  below zero
> is an artifact due to the way the virial term is calculated.

I must say that I don't agree with your interpretation....

The pressure is computed from the virial in the standard way.  The only
influence of PRES0 is on how the volume of the system will respond: if
the instantaneous pressure is greater than PRES0, the MD algorithm in
Amber will slowly increase the volume of the system; vice-versa if the
instantaneous pressure is below PRES0.  The *interpretation* of the
computed pressure is independent of the value of PRES0.

There is nothing "wrong" with negative pressures, and they need not arise
from errors or truncations in the virial calculation. The equation of
state of a rare gas fluid (for example) has substantial regions of density
and temperature where the equilibrium pressure is negative [see, for
example, chapter 14 of McQuarrie's "Statistical Mechanics", or other
textbooks.]  One of the tests I made (many years ago) to verify the correctness
of the way Amber calculates pressure terms was to reproduce these equations
of state for a Lennard-Jones 6-12 fluid.  [One can hope that bugs have
not crept in since then!]
-- 
                            Adrian E. Roitberg
========================================================================
    National Institute for   | Phone: (301) 975-4469. Fax (301) 330-3447
   Standards and Technology  | E-mail : [email-masked]
Bioprocess Engineering Group |
    Building 222, A-353      | URL : http://www.chem.nwu.edu/~adrian
  Gaithersburg, MD 20899     |
========================================================================

-- 
                                      ("`-/")_.-'"``-._
Gustavo A. Mercier,Jr.,MD,PhD         (. . `) -._    )-;-,_() 
Division of Nuclear Medicine          (v_,)'  _  )`-.\  ``-
Dept. of Radiology                    _;- _,-_/ / ((,'
University of Pennsylvania           ((,.-'  ((,/
3400 Spruce St.                  [email-masked]
Philadelphia, PA 19104         215-662-3069/3091 fax: 215-349-5843