On 7/11/2011 4:53 PM, khuchtumur bumerdene wrote:
Hi,
I have a question about refcoord_scaling option and its role in
pressure coupling. What I'm currently doing is just running the
lysozyme tutorial on my own protein.
I've so far had successful runs when equilibrating on my own machine,
which is currently running gmx4.5.3. NPT came to a stable average of
~1.0 bar. But when I do it on the cluster with gmx4.5.5, grompp gives
the following error:
WARNING 1 [file npt.mdp]:
You are using pressure coupling with absolute position restraints, this
will give artifacts. Use the refcoord_scaling option.
That's a warning, not an error. You need to apply judgement about
whether a warning is significant. Possibly this warning was introduced
between 4.5.3 and 4.5.5. If so, then the reason for the warning message
still applied to the 4.5.3 run and ignorance was bliss :-)
My npt.mdp file is as follows (straight from the tutorial at this stage):
define = -DPOSRES ; position restrain the protein
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 150000 ; 2 * 150000 = 300 ps
dt = 0.002 ; 2 fs
; Output control
nstxout = 100 ; save coordinates every 0.2 ps
nstvout = 100 ; save velocities every 0.2 ps
nstenergy = 100 ; save energies every 0.2 ps
nstlog = 100 ; update log file every 0.2 ps
; Bond parameters
continuation = yes ; Restarting after NVT
constraint_algorithm = lincs ; holonomic constraints
constraints = all-bonds ; all bonds (even heavy atom-H bonds)
constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
ns_type = grid ; search neighboring grid cells
nstlist = 5 ; 10 fs
rlist = 1.0 ; short-range neighborlist cutoff (in nm)
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range
electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein Non-Protein ; two coupling groups - more
accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 300 300 ; reference temperature, one for each
group, in K
; Pressure coupling is on
pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water,
bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; Velocity generation is off
The only change I made is the nsteps as my box is bigger and thus has
more waters in it.
So after reading the warning I tried to find some info on the
refcoord_scaling, which the only information found by my feeble
attempts were from the mdp options and the manual, which states:
*refcoord_scaling:*
*no*
The reference coordinates for position restraints are not
modified. Note that with this option the virial and pressure
will depend on the absolute positions of the reference
coordinates.
*all*
The reference coordinates are scaled with the scaling matrix
of the pressure coupling.
*com*
Scale the center of mass of the reference coordinates with the
scaling matrix of the pressure coupling. The vectors of each
reference coordinate to the center of mass are not scaled.
Only one COM is used, even when there are multiple molecules
with position restraints. For calculating the COM of the
reference coordinates in the starting configuration, periodic
boundary conditions are not taken into account.
I don't quite understand this completely, is the option there to
change the coordinates of the solute as the box size changes,
where the 'all' option changes the coordinates of every atom and
the com only changes the coordinates of the solute com? Please
correct me on this.
Under NPT the box size changes each step. You are using position
restraints to a pre-defined set of reference coordinates. This option
allows you to choose how those *reference coordinates* should change
when the box size changes (respectively do not scale them at all, scale
them all, or scale their COM but leave their internal geometry fixed).
Position restraints are then applied using the updated reference
coordinates.
So before posting this question I ran a few npt equilibrations with
the all and com options chosen with parrinello-rahman from 100 ps, to
300 ps, and to 1000 ps. All with 10 different initial velocities from
the nvt. The reason for increasing time is that the equilibration did
not reach the target 1 bar and ended at averages of 4-7 bars for the
10 equilibrations. The oscillation was comparable to the graph seen on
the lysozyme tutorial between 0-400.
So I tried running berendsen instead of parrinello rahman as the
manual says it should help me get to the target pressure better and I
can follow up with parrinello-rahman. But the pressure still didn't
reach the target pressure.
Pressure has large fluctuations and doesn't converge fast. How long you
need depends on the system size. Berendsen before PR is a very good idea
if your initial system is far from equilbrium.
So the questions are,
- is my naive understanding of the refcoord_scaling correct? if not
will you help explain it to me or point me to a link that could help me?
- Also, when pressure oscillations are at 100s of bars, is 3-6 bar
difference from 1 bar going to make any difference? (it should.
shouldn't it?)
Probably your observations are too short to be statistically
significant. The error analysis from g_energy may help you assess this.
Mark
- Or have I made some other careless mistake?
If you need more information, I'd be happy to provide them.
--
gmx-users mailing list gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
