[gmx-users] Questions about parameters
Justin Lemkul
jalemkul at vt.edu
Fri Mar 18 13:08:26 CET 2016
On 3/17/16 10:44 AM, abhishek khetan wrote:
> 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.
If there are bonds defined in a molecule, it won't "disintegrate."
> 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 =
The normal leap-frog integrator does support Nose-Hoover.
> 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
Neither position restraints nor constraints are necessary for this. You should
have normal bonds defined in your topology. Maybe you should revisit some of
your assumptions, because at this point it seems you're sort of randomly
combining terminology in the hopes that something will work. Not a wise approach.
> 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.
>
With the md integrator, Nose-Hoover and Parrinello-Rahman are supported. I use
them all the time. Likely grompp told you that generating velocities in
combination with Nose-Hoover is not stable and ended with a warning. This does
not mean that you can't use N-H with md. This just means you should not use it
in your initial equilibration because it's not stable.
>> 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 ?
>
http://www.gromacs.org/Documentation/Terminology/Pressure
http://pubs.acs.org/doi/abs/10.1021/ct300688p
https://github.com/shirtsgroup/checkensemble
>> 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
Well, a restraint has to be defined as the deviation from some position. The
reference coordinates are taken from grompp -c unless you provide -r (which is
only for more interesting cases like flat-bottom restraints).
> 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.
>
For a simple system of pure liquid, you likely shouldn't be bothering with any
restraints at all.
> 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
>
--
==================================================
Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow
Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201
jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul
==================================================
More information about the gromacs.org_gmx-users
mailing list