[gmx-users] Combining position restraints with Parrinello-Rahman pressure coupling leads to instabilities
Naba
nabajyoti.goswami at gmail.com
Tue Aug 21 08:33:24 CEST 2018
Dear Dr. Justin,
Thanks a lot for suggestions and hints.
I issued:
gmx dump -s npt.tpr
to check the velocities.
Following are last few lines after execution of the above command:
v[120488]={-8.24012e-01, -3.10189e-01, 3.53090e-01}
v[120489]={-5.70083e-02, 6.79958e-01, 1.43155e+00}
v[120490]={ 4.60961e-01, 8.73407e-02, -1.26087e+00}
v[120491]={ 3.43934e-01, -3.55264e-02, 5.41598e-02}
v[120492]={-2.71835e-01, 2.20839e-01, -2.95827e-01}
v[120493]={ 5.59665e-02, 6.66364e-01, -2.97680e-01}
v[120494]={ 1.22493e-01, 5.81446e-01, 1.32195e-01}
Group statistics
T-Coupling : 33740 86755 (total 120495 atoms)
Energy Mon. : 120495 (total 120495 atoms)
Acceleration: 120495 (total 120495 atoms)
Freeze : 120495 (total 120495 atoms)
User1 : 120495 (total 120495 atoms)
User2 : 120495 (total 120495 atoms)
VCM : 33740 86755 (total 120495 atoms)
Compressed X: 120495 (total 120495 atoms)
Or. Res. Fit: 120495 (total 120495 atoms)
QMMM : 120495 (total 120495 atoms)
-----
Using -maxwarn 1 in grompp, I produced the tpr file for Berendsen barostate
also and got the same output. So, I can proceed with Berendsen pressure
coupling, am I right?
Thanks & regards,
Naba
On 8/20/18 6:52 AM, Naba wrote:
> > Dear Gromacs users and developers,
> >
> > I am using Gromacs 2018.2.
> > Following the membrane protein simulation tutorial, I am planning to run
> > long simulations of a tetramer that needs larger lipid bilayer than 128
> > lipids as described in the tutorial. So, I have replicated the
> > pre-equilibrated POPC bilayer of 128 lipids from lipidbook (
> > https://lipidbook.bioch.ox.ac.uk/lipid/) to get a larger bilayer of 512
> > lipids using genconf. Moreover, I have extended gromos54a7_lipid.ff
> > forcefield to Berger lipid parameters for long membrane protein
> > simulations. I have set 10 ns for NVT and 20 ns for NPT equilibration.
> > NVT simulation ran successfully, but when I reached to NPT equilibration,
> > grompp showed me a note like the following:
> > "You are combining position restraints with Parrinello-Rahman pressure
> > coupling, which can lead to instabilities. If you really want to
> combine
> > position restraints with pressure coupling, we suggest to use
> Berendsen
> > pressure coupling instead."
>
> I think this note is a bit imprecise. grompp is guessing that the
> combination of Parrinello-Rahman and restraints implies equilibration,
> and it should really say so, or otherwise do a more rigorous check to
> see if velocities are being generated. In principle, there shouldn't be
> anything wrong with this approach, but Parrinello-Rahman *usually* isn't
> the best choice for equilibration.
>
> When I used Berendsen pressure coupling, grompp terminated with a warning
> > as:
> > "Using Berendsen pressure coupling invalidates the true ensemble for the
> > thermostat"
>
> This only matters for production runs. The Berendsen method should never
> be used for actual data collection. For equilibration, it's perfectly
> fine because you're going to throw this time out, anyway.
>
> -Justin
>
> > What to do with this? Can I proceed with the note using Parrinello-Rahman
> > pressure coupling? Please help.
> >
> > Following is the mdp parameters I have used:
> >
> > title = NPT Equilibration for KALP15-DPPC
> > define = -DPOSRES ; position restrain the protein
> > ; Run parameters
> > integrator = md ; leap-frog integrator
> > nsteps = 10000000 ; 2 * 10000000 = 20000 ps (20 ns)
> > dt = 0.002 ; 2 fs
> > cutoff-scheme = verlet
> > ; Output control
> > nstxout = 2500 ; save coordinates every 5 ps
> > nstvout = 2500 ; save velocities every 5 ps
> > nstenergy = 2500 ; save energies every 5 ps
> > nstlog = 2500 ; update log file every 5 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 cels
> > nstlist = 5 ; 10 fs
> > rlist = 1.2 ; short-range neighborlist cutoff (in nm)
> > rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
> > rvdw = 1.2 ; 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 = Nose-Hoover ; More accurate
> thermostat
> > tc-grps = Protein_POPC Water_and_ions ; two coupling groups -
> > more accurate
> > tau_t = 0.5 0.5 ; time constant, in ps
> > ref_t = 300 300 ; reference temperature,
> > one for each group, in K
> > ; Pressure coupling is on
> > pcoupl = Berendsen ; Pressure coupling on in NPT
> > pcoupltype = semiisotropic ; uniform scaling of x-y box vectors,
> > independent z
> > tau_p = 5.0 ; time constant, in ps
> > ref_p = 1.0 1.0 ; reference pressure, x-y, z (in bar)
> > compressibility = 4.5e-5 4.5e-5 ; isothermal compressibility, 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
> > ; COM motion removal
> > ; These options remove motion of the protein/bilayer relative to the
> > solvent/ions
> > nstcomm = 1
> > comm-mode = Linear
> > comm-grps = Protein_POPC Water_and_ions
> > ; Scale COM of reference coordinates
> > refcoord_scaling = com
> > nstcalcenergy = 1
> > nhchainlength = 1
> >
> >
> > Thank in advance.
> >
> > Regards,
> >
> > Nabajyoti Goswami
> >
> > Research Associate
> > Bioinformatics Infrastructure Facility
> > Department of Animal Biotechnology
> > College of Veterinary Science
> > Khanapara,Guwahati 781022
> > Assam, India
>
> --
> ==================================================
>
> Justin A. Lemkul, Ph.D.
> Assistant Professor
> Virginia Tech Department of Biochemistry
>
> 303 Engel Hall
> 340 West Campus Dr.
> Blacksburg, VA 24061
>
> jalemkul at vt.edu | (540) 231-3129
> http://www.thelemkullab.com
>
> ==================================================
> <https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users>
> <http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List>
> <http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List>
More information about the gromacs.org_gmx-users
mailing list