Subject: Re: vmd-l: Re: CCL:NAMD failure on Berendsen pressure scheme

From: Ioana Cozmuta <ioana-.at.-nas.nasa.gov>
Date: Thu, 31 Jul 2003 11:03:24 -0700 (PDT)

Hi Vlad,

In the input file you've sent:

#***temperature regulation for equilibration step (NPT) using
#***temperature coupling
temperature      100.0
tcouple          on
tcoupletemp      300.0
tcouplefile      rnaWI_ceq2.pdb  ### friction coef.=1.0
tcouplecol       B

you assign a value of 100K for the initial temperature (that will be used
to set the initial velocities of your atoms). Right after that you couple
your system to 300K. If you did a previous heating/equilibration to bring
your system up to 300K I think it is "cleaner" to continue with the
velocities in the .vel file produced in the heating/equilibration run. The
reason is that those velocities correspond to an equilibrated system while
by resetting the temperature in your NPT run to 300K does not necesarily
give you an equilibrated distribution of velocities for your system.
Setting it to 100 I think it makes even less sense. I am basing my
comments on the input file you've sent to the list.
Basically not continuing with the velocities written in the .vel file in
the previous run and using the temperature setting to initialize
velocities your first ps in the NPT run will be used to reequilibrate the
system (depending on the size and complexity of the system the
equilibration time varies).

I did not say that 1fs is not fine, I just suggested that you could try to
increase the integration time step of the EOM to 2fs specially in
combination with SHAKE. But then again, this was only a suggestion and the
final choice is of course yours.

There is a tutorial at
http://www.ks.uiuc.edu/Research/namd/tutorial/NCSA2002/hands-on/
where it explains how to generate the files where you define the values
for the restraining force constant. I remember I've tried the 0.5 value
and I encountered no problems. But I've always generated the files via
vmd (so maybe there is a format problem?)

1. If you look into the NAMD src directory in Controller.C
there are two definitions, one for PRESSURE and one for GPRESSURE

./Controller.C:  iout << "PRESSURE: " << step << " "
./Controller.C:  <<PRESSUREFACTOR * pressure << "\n"
./Controller.C:  << "GPRESSURE: " <<step << " "
./Controller.C:  << PRESSUREFACTOR * groupPressure << "\n" <<
endi;

PRESSURE is defined as pressure*PRESSUREFACTOR
GPRESSURE is defined as groupPressure*PRESSUREFACTOR

PRESSUREFACTOR is defined in the SimParameters.C
depending on the type of pressure coupling used (Berendsen coupling or
Langevin)

./SimParameters.C:     berendsenPressureTarget /= PRESSUREFACTOR;
./SimParameters.C:     langevinPistonTarget /= PRESSUREFACTOR;

and is read from your input file as berendsenPresureTarget or
langevinPistonTarget.

Pressure and groupPressure are calculated acording to the virial theorem
just that the sum is split into a fast term and a slow one (and I think
this is implemented this way because NAMD offers the possibility to use
MTS methods). Thus,to calculate "pressure" you use:

./Controller.C:      pressure = pressure_normal + pressure_nbond + pressure_slow;
./Controller.C:      pressure_normal = virial_normal / volume;
./Controller.C:        pressure_nbond = virial_nbond / volume;
./Controller.C:        pressure_slow = virial_slow / volume;

and groupPressure:

./Controller.C:      groupPressure_normal = ( virial_normal - intVirial_normal ) / volume;
./Controller.C:        groupPressure_nbond = ( virial_nbond - intVirial_nbond ) / volume;
./Controller.C:        groupPressure_slow = ( virial_slow - intVirial_slow ) / volume;
./Controller.C:  groupPressure = groupPressure_normal +groupPressure_nbond +
./Controller.C:                                        groupPressure_slow;

To my understanding pressure is the instantaneous value calculated at each
step for the pressure while GPRESSURE is the pressure difference between
the instantaneous pressure and the pressure achieved after the reduction
is applied (it is assumed that you are trying to equilibrate your system
to the pressure target which is read as PRESSUREFACTOR). If the pressure
in your system is too low then the factor is used to increase the
instantaneous values, if it is too high compared to the target the factor
is used to decrease the instantaneous values.
If anyone has more comments on this (more complete I mean) please enter
this conversation!

2. I did write my own perl script to get out of the NAMD output some
averages, values for density etc. It's not MSc-CS art but I can attach it
to an e-mail if you and other people think it is useful.
I do not have though the RMSD values but I could put that in as well. It's
actually a good suggestion.

Ioana

On Thu, 31 Jul 2003 Vlad.Cojocaru-.at.-mpi-bpc.mpg.de wrote:

> Dear Ioana,
>    Thank you very much for the reply. this is a second equilibration step
> (first NVT and this one is NPT so the timestep 1.0 is still ok...the
> production runs I will run with 2.0) The "temperature=300.0" and not 100 and I
> believe you dont need to set "velocities=file" if I set the temp to 300. I
> traced the problem..it was something to do with the constraints applied to
> hold the solute fixed.
> In Amber I use:
> Group input for holding the solute fixed
>   500.0
> RES 1 142
> END
> END
> So I tried to apply the same in NAMD by using 500.0 for the force constant in
> the B column of the file for restraints...somehow if I changed this 500.000 to
> 1.000 in that file thigs work ok.  And the values of the properties are
> comparable to the amber output.
> If you have a bit of time could you tell me:
> 1. what is GPRESSURE in the NAMD output?
> 2. How can find the averages and the RMS deviations (T, P, Energy) in the
> output of NAMD (something simi lar to amber "AVERAGES AFTER N STEPS" and "RMSD
> AFTER N STEPS")..I am intersted especiall in Temp, Pressure, Total energy,
> Volume and Density.
> Thank you very much again,
> Ceau
> Vlad
>
>
>
>
> -------------------
> > Hi Vlad,
> >
> > 1. You say that you did an NVT equilibration previous to the NPT dynamics
> > you want to run. But there is nowhere in your script a command where you
> > load the velocities from the file generated during equilibration
> > (velocities   file.vel). Why do you heat up your system if you do not use
> > in the end the velocities of the solvent? I do not understand why you are
> > going back to Temperature 100.
> > 2. I usually also read the unit cell information from the xsc file from
> > the previous equilibration
> > extendedSystem   file.xsc
> > 3. You should be able to use a
> > timestep of 2fs. I am not sure what water model you are using but with
> > Shake/Rattle on this should work (it works for me).
> > 4. Here is what I use for the Berendsen pressure coupling:
> >
> > berendsenPressure on
> > berendsenPressureTarget 1.01325
> > berendsenPressureCompressibility 4.9e-5
> > berendsenPressureRelaxationTime 500
> > berendsenPressureFreq 1
> >
> > 5. I also define the other parameters (even if they are =1)
> > fullElectFrequency 1
> > nonbondedFreq      1
> > stepspercycle    1
> >
> > 6. As a general rule, if something is not working, maybe you should look
> > back at your heating procedure and minimization, plot out some parameters
> > (T, Epotential, Etotal, Evdw..etc). Read Berendsen's paper about the way
> > the temperature and pressure coupling works.
> >
> > 7. Read careful the information listed in the header of the NAMD output
> > file (INFO), NAMD is sensitive to the name of the variables you define in
> > your input file. And may read them wrong or not at all.If they are
> > critical your simulation will not start at all but it may be that they are
> > not critical for the simulation to run and still their value will be the
> > default one
> >
> >  Hope this helps,
> > Ioana
> >
> >
> >
> >
> >
> > On Wed, 30 Jul 2003, Vlad Cojocaru wrote:
> >
> > > Dear fellows,
> > >  If you have some experience with namd could somebody explain why this
> > > is happening (see error message at the bottom of the message)? I am
> > > using the following script:. It is a constant pressure equilibration
> > > using Berendsen algorithm for temp and pressure. This equilibration step
> > > follows after a minim and a const Vol. eq in which the temperature is
> > > adjusted to 300 (also with temp coupling scheme). The protocol below
> > > works fine in AMBER but maybe there some differences between the values
> > > used by the 2 programs that I didnt figure out. The first eq. at const
> > > Vol works fine both in Amber and Namd so I guess there is something
> > > about pressure here. If you want more details please let me know!.
> > >   I would appreciate any advice. Thank you very much ,
> > > Vlad
> > >
> > > # *** AMBER force field *** ---------------------------------
> > > amber            on
> > > parmfile         rnaWI.top
> > > coordinates      rnaWI_mdeq2.pdb.coor
> > > readexclusions   yes
> > > exclude          scaled1-4
> > > 1-4scaling       0.83333333   #=1/1.2
> > > scnb             2
> > >
> > > #*** approximations for nonbonded interactions ------------
> > > switching        off
> > > cutoff           9
> > > stepspercycle    1
> > >
> > > #***equilibration parameters
> > > minimization     off
> > > numsteps         10000
> > > timestep         1.0
> > >
> > > #***SHAKE use
> > > rigidbonds all
> > >
> > > #***temperature regulation for equilibration step (NPT) using
> > > #***temperature coupling
> > > temperature      100.0
> > > tcouple          on
> > > tcoupletemp      300.0
> > > tcouplefile      rnaWI_ceq2.pdb  ### friction coef.=1.0
> > > tcouplecol       B
> > >
> > >
> > > #***constant pressure equilibration using Berendsen bath coupling
> > > usegrouppressure                         on
> > > useflexiblecell                          on
> > > berendsenpressure                        on
> > > berendsenpressuretarget                  1.0
> > > berendsenpressurecompressibility         0.0000446
> > > berendsenpressurerelaxationtime          500
> > > berendsenpressurefreq                    100
> > >
> > > #***PME
> > > PME         on
> > > cellBasisVector1    56.21    0    0
> > > cellBasisVector2    0    58.53    0
> > > cellBasisVector3    0    0    57.89
> > > cellOrigin        0   0   0
> > > PMEGridSizeX        50
> > > PMEGridSizeY        55
> > > PMEGridSizeZ        50
> > > XSTfile         rnaWI_mdeq3.XST
> > > XSTfreq         200
> > >
> > > #***constraints (holding the solute fixed)
> > > constraints      on
> > > consexp          2
> > > consref          rnaWI_mdeq2.pdb.coor
> > > conskfile        rnaWI_res500.pdb
> > > conskcol         B
> > >
> > > #***outputnames
> > > outputname       rnaWI_mdeq3.pdb
> > > restartname      rnaWI_mdeq3.rst
> > > dcdfile          rnaWI_mdeq3.dcd
> > > dcdfreq          200
> > > veldcdfile       rnaWI_mdeq3.vel.dcd
> > > veldcdfreq       200
> > > restartfreq      200
> > > outputenergies   100
> > > outputpressure   100
> > > binaryoutput     off
> > > binaryrestart    off
> > >
> > >
> > >
> > > ###---error-----------------
> > >
> > > RESSURE: 300 76798.7 -40631.2 58693.1 -13378.6 29200.3 -5152.76 69223.5
> > > -42523.7 81493.1
> > > GPRESSURE: 300 76687.2 -40836.2 58854.3 -13914.7 29538.5 -4262.26
> > > 68552.2 -42832.5 81484.6
> > > ENERGY:    300 2002.7469   2400.4670   550.4454    0.0000
> > > -61779.2755 9089.4231   60866.5902  0.0000      31247.2457  44377.6430
> > > 1056.1747   62497.3627  62570.0965  185386.7344 -1757.5130  -1600.7691
> > > ERROR: Constraint failure in RATTLE algorithm for atom 121!
> > > ERROR: Constraint failure; simulation has become unstable.
> > > ERROR: Exiting prematurely.
> > > ==========================================
> > >
> > > --
> > > Vlad Cojocaru
> > > Max Planck Institute for Biophysical Chemistry
> > > Department: 060
> > > Am Fassberg 11, 37077 Goettingen, Germany
> > > tel: ++49-551-201.1327
> > > e-mail: Vlad.Cojocaru-.at.-mpi-bpc.mpg.de
> > > home tel: ++49-551-9963204
> > >
> > >
> > >
> > >
> > >
> > > -= This is automatically added to each message by mailing script =-
> > > To send e-mail to subscribers of CCL put the string CCL: on your Subject:
> line
> > > and send your message to:  CHEMISTRY-.at.-ccl.net
> > >
> > > Send your subscription/unsubscription requests to:
> CHEMISTRY-REQUEST-.at.-ccl.net
> > > HOME Page: http://www.ccl.net   | Jobs Page: http://www.ccl.net/jobs
> > >
> > > If your mail is bouncing from CCL.NET domain send it to the maintainer:
> > > Jan Labanowski,  jkl-.at.-ccl.net (read about it on CCL Home Page)
> > > -+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
> > >
> > >
> > >
> > >
> > >
> > >
> >
>