[gmx-users] thermostats, canonical distribution, and stochastic dynamics

[Michael Brunsteiner]

Dear all,

I encountered problems when trying to keep
a simple system at constant temperature, and would
appreciate any help/advice!

I try to calculate a PMF between two relatively simple
particles (about spherical, each 60 atoms) in solution with
umbrella sampling. For now i model the solvent as primitive
model electrolyte (vacuum + epsilonr=80), and what i see is
this:

trying out the following thermostats:
* leapfrog + nose hoover
* velocity verlet + 10 nose hoover chains
* v-rescale
in ALL cases, i get what seems to be
a non-canonical distribution:

Since the distance between the two particles is restrained
via a harmonic potential one would expect (at least for large
distances) to get an approximately Gaussian distribution of
the distance between the two particles, centred near the
equilibrium distance of the harmonic restraint. However, I 
keep seeing a bi-modal distribution (as one would expect 
from a harmonic oscillator with no thermostat) which, I believe,
can be taken as evidence for a non-canonical phase space
coverage (and probably this also screws up the WHAM
analysis)

My system is relatively simple - pathologically so one might
argue - but then, what I recall from articles I read, especially
about Nose Hoover chains, it seems that I could expect getting
a proper canonical distribution even for much simpler systems
with no more than a few degrees of freedom (and my system
has >300 degrees of freedom)

With Langevin dynamics (instead of md) the distance
distribution DOES become Gaussian, but now I encounter
a different problem: More so than with MD convergence
becomes an issue: especially with small tau_t (<5 ps). 
E.g. to cover a distance of about 3 nm it takes at least 44 windows
with >=1 ns each to get a reasonably converged PMF. So I wonder:
i) is there any limit to the size of tau_t I am supposed to use? or 
ii) is there any other way to accelerate convergence?

FTR: I also tried using two tc_grps (one for each particle) 
but, at least qualitatively, i see the same results.

thanks in advance for any help!

regards,
Michael



===
mdp-file:

integrator      = sd               ; or md, or md-vv-avek
tau_t           = 0.2 0.2          ; varies from 0.2 to 10.0 for sd
ref_t           = 300 300
tc_grps         = p1 p2            ; or just System
dt              = 0.001
tinit           = 0
nsteps          = 400000           ; or more
nstxtcout       = 0
nstxout         = 1000
nstvout         = 0
nstfout         = 0
nstenergy       = 100
constraints     = h-bonds
comm_mode       = Linear
nstcomm         = 1
pbc             = no
nstlist         = 0
ns_type         = simple
rlist           = 5.0
rcoulomb        = 5.0
rvdw            = 5.0
coulombtype     = Cut-off
vdwtype         = User          ; repulsive only LJ pot
epsilon_r       = 80
pull            = umbrella
pull_geometry   = distance
pull_dim        = N N Y
pull_start      = no
pull_ngroups    = 1
pull_group0     = p1
pull_group1     = p2
pull_rate1      = 0.0
pull_k1         = kkk            ; varies from 500-5000
pull_nstxout    = 100
pull_nstfout    = 100
pull_pbcatom0   =  61
pull_pbcatom1   = 122
pull_init1      = ddd             ; varies from 1-4
freezegrps      = c1 c2
freezedim       = Y Y N Y Y N     ; this is to keep the particles in 1-D
energygrps      = p1 p2