VMD-L Mailing List

From: Ioana Cozmuta (ioana_at_nas.nasa.gov)
Date: Thu Jul 31 2003 - 13:03:24 CDT

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_osc.edu (read about it on CCL Home Page)
> > > -+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
> > >
> > >
> > >
> > >
> > >
> > >
> >
>