[gmx-users] Questions about parameters

[abhishek khetan]

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


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 ?
2b. If I do not use any constraints, then does gromacs still ensure that
the individual solvent molecules do not disintegrate ?
2c. When I use constraints = all-bonds, then does gromacs still allow for
the harmonic vibration of the individual bonds within a single molecule?
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.)?

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

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 ?
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^-) ?


-- 
MfG,
abhishek