<html>
<head>
<meta content="text/html; charset=ISO-8859-1"
http-equiv="Content-Type">
</head>
<body bgcolor="#FFFFFF" text="#000000">
On 7/11/2011 4:53 PM, khuchtumur bumerdene wrote:
<blockquote
cite="mid:CAGYCuvPgRbXjUEST4a22Wx-abMtnq5zVu+LJcNRh_nyhc3q2QA@mail.gmail.com"
type="cite">Hi,<br>
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. <br>
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:<br>
<br>
WARNING 1 [file npt.mdp]:<br>
You are using pressure coupling with absolute position
restraints, this<br>
will give artifacts. Use the refcoord_scaling option.<br>
</blockquote>
<br>
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 :-)<br>
<br>
<blockquote
cite="mid:CAGYCuvPgRbXjUEST4a22Wx-abMtnq5zVu+LJcNRh_nyhc3q2QA@mail.gmail.com"
type="cite"><br>
My npt.mdp file is as follows (straight from the tutorial at this
stage):<br>
<br>
define = -DPOSRES ; position restrain the protein<br>
; Run parameters<br>
integrator = md ; leap-frog integrator<br>
nsteps = 150000 ; 2 * 150000 = 300 ps<br>
dt = 0.002 ; 2 fs<br>
; Output control<br>
nstxout = 100 ; save coordinates every 0.2 ps<br>
nstvout = 100 ; save velocities every 0.2 ps<br>
nstenergy = 100 ; save energies every 0.2 ps<br>
nstlog = 100 ; update log file every 0.2 ps<br>
; Bond parameters<br>
continuation = yes ; Restarting after NVT<br>
constraint_algorithm = lincs ; holonomic constraints<br>
constraints = all-bonds ; all bonds (even heavy atom-H
bonds) constrained<br>
lincs_iter = 1 ; accuracy of LINCS<br>
lincs_order = 4 ; also related to accuracy<br>
; Neighborsearching<br>
ns_type = grid ; search neighboring grid cells<br>
nstlist = 5 ; 10 fs<br>
rlist = 1.0 ; short-range neighborlist cutoff
(in nm)<br>
rcoulomb = 1.0 ; short-range electrostatic cutoff
(in nm)<br>
rvdw = 1.0 ; short-range van der Waals cutoff
(in nm)<br>
; Electrostatics<br>
coulombtype = PME ; Particle Mesh Ewald for
long-range electrostatics<br>
pme_order = 4 ; cubic interpolation<br>
fourierspacing = 0.16 ; grid spacing for FFT<br>
; Temperature coupling is on<br>
tcoupl = V-rescale ; modified Berendsen thermostat<br>
tc-grps = Protein Non-Protein ; two coupling groups -
more accurate<br>
tau_t = 0.1 0.1 ; time constant, in ps<br>
ref_t = 300 300 ; reference temperature, one for
each group, in K<br>
; Pressure coupling is on<br>
pcoupl = Parrinello-Rahman ; Pressure coupling on in
NPT<br>
pcoupltype = isotropic ; uniform scaling of box vectors<br>
tau_p = 2.0 ; time constant, in ps<br>
ref_p = 1.0 ; reference pressure, in bar<br>
compressibility = 4.5e-5 ; isothermal compressibility of
water, bar^-1<br>
; Periodic boundary conditions<br>
pbc = xyz ; 3-D PBC<br>
; Dispersion correction<br>
DispCorr = EnerPres ; account for cut-off vdW scheme<br>
; Velocity generation<br>
gen_vel = no ; Velocity generation is off<br>
<br>
<br>
The only change I made is the nsteps as my box is bigger and thus
has more waters in it. <br>
<br>
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:<br>
<br>
<a moz-do-not-send="true" name="pc"><dt><b>refcoord_scaling:</b></dt>
<dd>
<dl compact="compact">
<dt><b>no</b></dt>
<dd>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.</dd>
<dt><b>all</b></dt>
<dd>The reference coordinates are scaled with the scaling
matrix of the pressure coupling.</dd>
<dt><b>com</b></dt>
<dd>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.
</dd>
</dl>
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. <br>
</dd>
</a></blockquote>
<br>
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.<br>
<br>
<blockquote
cite="mid:CAGYCuvPgRbXjUEST4a22Wx-abMtnq5zVu+LJcNRh_nyhc3q2QA@mail.gmail.com"
type="cite"><a moz-do-not-send="true" name="pc"><dd>
<br>
</dd>
</a>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. <br>
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.<br>
</blockquote>
<br>
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.<br>
<br>
<blockquote
cite="mid:CAGYCuvPgRbXjUEST4a22Wx-abMtnq5zVu+LJcNRh_nyhc3q2QA@mail.gmail.com"
type="cite">
<br>
So the questions are, <br>
- 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?<br>
- 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?)<br>
</blockquote>
<br>
Probably your observations are too short to be statistically
significant. The error analysis from g_energy may help you assess
this.<br>
<br>
Mark<br>
<br>
<blockquote
cite="mid:CAGYCuvPgRbXjUEST4a22Wx-abMtnq5zVu+LJcNRh_nyhc3q2QA@mail.gmail.com"
type="cite">
- Or have I made some other careless mistake?<br>
<br>
If you need more information, I'd be happy to provide them. <br>
<br>
<br>
<fieldset class="mimeAttachmentHeader"></fieldset>
<br>
</blockquote>
<br>
</body>
</html>