[gmx-users] Questions about parameters

[Justin Lemkul]

On 3/16/16 3:13 PM, abhishek khetan wrote:
> Dear gmxers,
>
> After some basic simulations of a box of a non-aqueous solvent, I want to
> know your opinion about the meaning of some parameters, the values I have
> used for these, and how they affect the accuracy of simulations and their
> speed (which at this moment is rather unimportant for me).
>
> I started with a box of 2.21x2.21x2.21 nm3 with (4x4x4=) 64 dimethoxyethane
> or DME molecules. The box was created so as to match the exact experimental
> density to begin with. First I wanted to do the NVT ensemble with Maxwell
> like distribution so I went for a total of 1000 ps md-vv integrator and
> nose-hoover thermostat with:
> dt     = 0.0005
> tcoupl = nose-hoover
> tau-t  = 5.0
> ref-t  = 300
>
> Next I did a NPT ensemble for 1000 ps with the md integrator, but with many
> parameters changed as:
> tcoupl                   = v-rescale
> tau-t                    = 1.0
> ref-t                    = 300
> pcoupl                   = parrinello-rahman
> pcoupltype               = isotropic
> tau-p                    = 2.0
> compressibility          = 1.2e-4
> ref-p                    = 1.0
> refcoord_scaling         = com
>
> I needed to change the parameters in order to make the simulations work,
> otherwise there were errors which wouldn't let it start, because certain
> tcoupl were not compatible with certain pcoupl or integrators. I *DID NOT
> USE position restraints*. My NPT simulations converged to a final box
> volume very close to the original. More details in my questions, which 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.

> 1b. How does the value of tau-t and tau-p affect my accuracy and speed? In
> my opinion, over a long time simulation, they shouldn't effect accuracy,
> but only speed. I am asking this to confirm if I can change them as I want
> in order to prevent the simulations from blowing up.
>

It has no effect on speed; it's just a value that goes into a function.  For 
accuracy, you can use ensemble checking tools, e.g. those published by the 
Shirts group.  I suspect they have very little impact overall.

>
> 2a. In one set (NVT followed by NPT) of simulations used constraints = all
> bonds and in another set I did not use anything. The simulations with
> constraints  = all-bonds is a bit slower. What effect do they have on the
> accuracy. In the most ideal case of a god given force-field, I should be
> able to get by without using any constraints, or ?

Constraints are routinely applied to bonds involving hydrogen to allow dt = 1 fs 
or larger.  People commonly constrain all bonds, though strictly speaking that 
shouldn't be done for certain force fields; the errors should be quite small in 
magnitude.

> 2b. If I do not use any constraints, then does gromacs still ensure that
> the individual solvent molecules do not disintegrate ?

GROMACS does what you tell it.  The integrity of the physical model depends on 
the topology and the algorithms you chose.  Bonds are bonds, and they can't be 
broken or formed, but the simulation can crash if you choose to do something 
physical unstable.

> 2c. When I use constraints = all-bonds, then does gromacs still allow for
> the harmonic vibration of the individual bonds within a single molecule?

No, by definition.

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

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

> 3b. If the absolute positions are restrained, then does this also not mean
> that in essence the indivdual molecules can never have Maxwell type
> velocity distribution? Would it not be better to have the latter kind of
> position restraint where the positions are penalized on their deviations
> from the centre of mass of individual molecules but the molecules can have
> non-penalized movement.
> 3c. How do these restraints affect my accuracy?

Related to both questions, the application of restraints is for specific 
purposes; what purpose would they have in your simulation?

> 3d. If I do not restrain the positions of some of the atoms in the
> molecule, does gromacs still ensure that the molecules do not disintegrate?

Again, GROMACS only does what you tell it.

> 3e. Is there any reason why in the itp file the restraints are of values
> 1000 and 0 only. Would something in the middle, like 500, make sense?
>

Try it and see.  The generic values are just that - generic and likely to work 
in most cases.

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

> 4b. How important is the value of compressibility? for many common
> solvents, this value is not available and I am choosing the value of the
> most similar specie. Would that affect the accuracy a lot for liquids,
> which have the same order of values (10^-4 to 10^-5 bar^-) ?
>

The value has little impact on anything other than the barostat response time.

-Justin

-- 
==================================================

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

==================================================