[gmx-users] refcoord_scaling

[Mark Abraham]

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.
>
>
>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20111107/67f15677/attachment.html>