[gmx-users] Questions about parameters

abhishek khetan askhetan at gmail.com
Thu Mar 17 15:44:35 CET 2016


Dear Justin,

Thanks for your insightful replies. I decided that I will not use position
restraints but rather all-bond constraints because I want to allow the
molecules to move freely and reorient without getting disintegrated,
although they did not look like disintegrating even without constraints.
Distance restraints might be an even better option in this regard but I
have to understand it better. Therefore, I have still some more questions
as I am totally new to MD and gromacs. They are as follows:

> 1a. Is it okay (in sense of physical accuracy) to change the tcoupl and
integrator when going from NVT to NPT? -- Without knowing what errors you
were getting, it's hard to say.  Just about all combinations should be
supported so I don't know why this is necessary.  But in  any case, any
time you change an algorithm, you need to allow sufficient time for
relaxation.
==> I am allowing for 1000 ps NVT followed by 1000 ps of NPT, which seems
okay enough for relaxation. When I use integrator = md, then I cannot use
tcoupl = nose-hoover. When I use integrator = md-vv with tcoupl =
nose-hoover, I cannot use pcoupl = MTTK because it gives the error - "MTTK
not compatible with lincs -- use shake instead." and "Constraints are not
implemented with MTTK pressure control". So I cannot use MTTK with
constraints. As, I explained, I want to allow movement of molecules without
disintegrating them, therefore I want to use constraints instead of
absolute position restraints. Then, the only option I am left with is to
use pcoupl = parrinello-rahman, for which I cannot use integrator = md-vv,
which means I cannot use tcoupl = nose-hoover and therefore I have to use
tcoupl = v-rescale, which is the closest thing to nose-hoover. Seen this
way, most (desirable and physically accurate) combinations are not
possible.

> 2.d The volume after NPT with no constraints was 2.20x2.20x2.20 nm3 and
the volume after NPT with contraints = all bonds was 2.22x2.22x2.22 nm3.
They are so close to the original, which should I trust more in terms of
physical reality (essentially same as question 2a.)? -- The volume of a
single snapshot means nothing.  The ensemble average is probably
indistinguishable, but that's what you should check.  Such a tiny box is
going  to be prone to massive pressure fluctuations, though, so anything
related to  pressure or density is probably unreliable, especially with
Parrinello-Rahman.
==> You're spot on. The pressure fluctuations are massive(-200 to 800
bars). But then, can it not find its equilibrium volume because of the
pressure coupling I used ? Is increasing the box size with several more
molecules the only way to achieve reasonable pressure and density
fluctuations? Could you please guide me to the page/tutorial/link where I
can learn more about the ensemble average and look at all the samples. Is
there a theory article on this ?

> 3a. What is way position restraints work? I understand that they
introduce a heavy energy penalty on the movement of atoms, but do they
apply this penalty on the absolute deviation of the atoms' positions or on
the deviation of the atom's positions with respect to the centre of mass of
the respective molecule? -- Depends on the what option you choose for
refcoord_scaling.
> 4a. How does the refcoord_scaling work? when i use 'com', does it mean
that the coordinate are scaled with respect to the COM of the whole system
in a way that even bond lengths within the individual molecules change ? --
This is explained in the manual.
==> The gromacs manual 5.1.2 doesn't not have the string 'refcoord' or
'refcoord_scaling'. When I look up the mdp options online, then it says
'refcoord_scaling = com'
would  "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". My question is
- whats the meaning of reference coordinates? to what 'reference' are these
coordinates referenced ? How are they different from the absolute
coordinates ? In case I understand it correctly, these are the coordinates
referenced to the COM of the entire system and this kind of
refcoord_scaling doesn't change the bond lengths but just moves the COMs of
the - which to me is also what I want. However, if I use 'refcoord_scaling
= all', does this mean the bond lengths will also change ? When I use
'refcoord_scaling=com' along with position restraints, does it mean that
the molecule's COM are scaled freely but the internal positions are
restrained with heavy energy penalty? - This would be ideal and awesome,
because right now I am using constraints instead of restraints. The manual
doesn't define these terms clearly, therefore, it would be very kind if you
could please clarify this.

Last question: The last set of lines from the *log from my NPT simulation
were:
-----------------------------------------------------------------------------------------------
   Energies (kJ/mol)
          Angle Ryckaert-Bell.          LJ-14     Coulomb-14        LJ (SR)
    1.91627e+03    3.78194e+02    5.34471e+02    2.70351e+03   -3.44494e+02
  Disper. corr.   Coulomb (SR)   Coul. recip.      LJ recip.      Potential
   -1.46846e-06   -1.93059e+03    1.47488e+02   -1.80265e+03    1.60220e+03
    Kinetic En.   Total Energy    Temperature Pres. DC (bar) Pressure (bar)
    2.70185e+03    4.30405e+03    3.08163e+02   -2.25676e-06    4.09580e+02
   Constr. rmsd
    8.13024e-07

    <======  ###############  ==>
    <====  A V E R A G E S  ====>
    <==  ###############  ======>

    Statistics over 1000001 steps using 10001 frames

   Energies (kJ/mol)
          Angle Ryckaert-Bell.          LJ-14     Coulomb-14        LJ (SR)
    1.86921e+03    3.71064e+02    5.32660e+02    2.71914e+03   -3.02512e+02
  Disper. corr.   Coulomb (SR)   Coul. recip.      LJ recip.      Potential
   -1.43373e-06   -1.93579e+03    1.42894e+02   -1.77019e+03    1.62647e+03
    Kinetic En.   Total Energy    Temperature Pres. DC (bar) Pressure (bar)
    2.63258e+03    4.25906e+03    3.00263e+02   -2.15237e-06    9.24047e+00
   Constr. rmsd
    0.00000e+00

          Box-X          Box-Y          Box-Z
    2.22960e+00    2.22960e+00    2.22960e+00

   Total Virial (kJ/mol)
    8.82221e+02   -2.19096e+00    3.19987e+00
   -2.19097e+00    8.64520e+02    7.99956e-01
    3.19995e+00    7.99981e-01    8.80886e+02

   Pressure (bar)
   -1.04437e+01    7.15563e+00   -1.23345e+01
    7.15566e+00    4.29864e+01   -1.89927e+00
   -1.23348e+01   -1.89935e+00   -4.82135e+00


-----------------------------------------------------------------------------------------------

As you can see, the potential energy is positive. Does this not mean that
these simulations are totally unphysical ? What can I do to remedy this ??
Thanks in advance for your detailed answers.

--
MfG,
abhishek


More information about the gromacs.org_gmx-users mailing list