[gmx-users] PMF from pull code, unexpected results

[Michael Brunsteiner]

Dear all,

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
unexpected results:

The PMFs look very differnt depending on whether I use "pulldim  Y Y Y" 
or "pulldim  "N N Y":
http://www.brunsteiner.net/tbut-pmf.gif
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 
constraints.

With "pulldim  Y Y Y" the RDF (option nolog in g_wahm) looks, in fact, 
like a NON-normalized rdf:
http://www.brunsteiner.net/tbut-rdf.gif
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 
files)
http://www.brunsteiner.net/tbut-f.gif
http://www.brunsteiner.net/tbut-z.gif
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
results:
http://www.brunsteiner.net/2at-pmf.gif
http://www.brunsteiner.net/2at-rdf.gif
http://www.brunsteiner.net/2at-z.gif
http://www.brunsteiner.net/2at-f.gif

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

my specs:
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!

cheers,
Michael



=============
mdp file:

integrator      = sd        ; this is better than md/NHT for systems with very 
few DOFs
tau_t           = 2.0 2.0
ref_t           = 290 290   ; a bit lower than 298 since sd with default 
parameters
                            ; 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 
distance;
                                    ; 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