[gmx-users] PMF from pull code, unexpected results
mbx0009 at yahoo.com
Mon Feb 21 19:47:44 CET 2011
I would like to calculate the potential of mean force between two molecules
in aqueous solution using the pull code. For a start i performed a number
of calulations with a couple of very simple model systems with settings
loosely based on the example in the gmx tutorial section (see mdp file included
at the bottom). My box is 8x8x8 nm**3, pmf is calculated from some r_min up to
approx 3.8 nm. Performing simulations with two t-buntanol molecules in implicit
solvent, anlyzing the output (x and f*.xvg files) with g_wham, i get rather
The PMFs look very differnt depending on whether I use "pulldim Y Y Y"
or "pulldim "N N Y":
Since we look at two neutral molecules with a permanent dipole
moment one would expect the PMF to be negative up to a distance
at which the two molecules come into contact, but with "pulldim Y Y Y"
the PMF is positive, i.e., repulsive, throughout.
With "pulldim N N Y" I constrain the two molecules in two
dimensions (by freezing the central atom in the x and y dimensions only,
BTW any ideas how else i could achieve that?)
One could argue that a combination of frozen atoms and pull code might be
problematic, but i freeze and pull in different dimensions, so this should
be OK, and, more importantly: with "pulldim N N Y" i get results that
look much more reasonable than with "pulldim "Y Y Y" and NO additional
With "pulldim Y Y Y" the RDF (option nolog in g_wahm) looks, in fact,
like a NON-normalized rdf:
and, interestingly if i normalize it myself (by dividing through z**2 before
taking the logarithm) the resulting PMF looks much more reasonable.
Although what I see suggests that with "pulldim Y Y Y" the RDF just doesn't
normalized properly, this issue seems to be more involved:
Looking at the average forces and positions (the content of the f*xvg and x*xvg
suggests that there's already something wrong in the mdrun output
meaning that the problem is further upstream and not connected to anything
done by g_wham.
I repeated the calculations with an even simpler system (2 single atoms
that only interact via a LJ potential) to get qualitatively the same
A few more remarks:
1) Convergence does not seem to be an issue here as i extended
some of the calculations to include 35+ windows, 2 ns each, and
the PMF remains the same nearly quantitatively.
2) The length of my reaction coordinates is always shorter
than half the box length.
3) I've compared calculations cut-off vs PME, to get very similar
4) If I use "pulldim Y Y Y" with no additional constraints but with
"comm_mode = Angular" i get results somewhere inbetween the two cases
Gromacs version: 4.5.3
Linux 2.6.35-23-generic Ubuntu x86_64
gcc version 4.4.5
I am not sure whether i overlooked something in my input,
or whether there's a problem with code.
I'd be grateful for any ideas/suggestions!
integrator = sd ; this is better than md/NHT for systems with very
tau_t = 2.0 2.0
ref_t = 290 290 ; a bit lower than 298 since sd with default
; has problems controlling the temperature.
tc_grps = p1 p2 ; also tried System here, no difference
dt = 0.002
tinit = 0
nsteps = 100000 ; played with that, results seem to have converged
nstxout = 1000
nstvout = 0
nstfout = 1000
nstenergy = 100
constraint_algorithm = lincs
constraints = hbonds
continuation = no
comm_mode = Linear
nstcomm = 1
pbc = xyz ; also tried "pbc = no" same result
nstlist = 10
ns_type = grid
rlist = 4.0
rcoulomb = 4.0
rvdw = 4.0
coulombtype = Shift
vdwtype = Shift
epsilon_r = 80
pull = umbrella
pull_geometry = distance
pull_dim = N N Y ; or "Y Y Y"
pull_start = yes
pull_ngroups = 1
pull_group0 = p1
pull_group1 = p2
pull_rate1 = 0.0
pull_k1 = 500 ; i let this number vary depending on the pull
; also played with the numbers, same result
pull_nstxout = 100
pull_nstfout = 100
pull_pbcatom0 = 1 ; my system geometry is so that the molecules
never leave the box
pull_pbcatom1 = 16 ; or cross a cell boundary.
freezegrps = c1 c2 ; with pulldim Y Y Y these lines are
freezedim = Y Y N Y Y N ; removed or commented out
energygrps = p1 p2
More information about the gromacs.org_gmx-users