[gmx-users] Umbrella Sampling on Lipid Bilayers, Pulling Issue and Umbrella Sampling Window Simulation Issues.

[Alexander Phillips]

Hello Everyone,

I am currently running a CG MARTINI simulation on a bilayer of POPC Lipid
in Water. I also have a single vitamin E in the mixture aswell (inside of
the POPC bilayer). I am trying to perform umbrella sampling by pulling the
vitamin E and measuring it's binding energy to POPC. I am basing most of my
work off of the tutorial made by Justin Lemkul. The plan for my umbrella
sampling is to pull the vitamin E vertically out of the membrane  (from the
center of the membrane to some distance into the water). However I am
having issues with performing umbrella sampling. I believe the issue is one
of two things (or perhaps a mixture of both):
1) The simulation where I pull the Vitamin E out of POPC is happenning too
quickly and I cannot gather enough frames to get valid input windows for
the umbrella sampling portion.
2) The umbrella sampling portion is not keeping the Vitamin E around the
starting portion and it falls back into the membrane.

I will attach the mdp files I used for both, and hopefully someone can give
me some advice about how to do this / fix it.
I am running in Gromacs 4.6.5 so a few things are changed from the files
given in Justin's tutorial, however I believe ultimately it's about the
same.

Here is the mdp for the Pulling:


integrator         = md
dt                        = 0.03
nsteps               = 250000
nstcomm          = 10

comm-grps      =
nstxout             = 5000
nstvout             = 5000
nstfout              = 5000
nstlog                = 1000  ; Output frequency for energies to log file
nstenergy         = 100   ; Output frequency for energies to energy file
nstxtcout          = 1000  ; Output frequency for .xtc file
xtc_precision   = 5
xtc-grps             =

energygrps       = Lipids Solvent PULLING_ATOC
nstlist                  = 10
ns_type               = grid
pbc                      = xyz
rlist                       = 1.4
coulombtype    = Shift  ;Reaction_field (for use with Verlet-pairlist) ;PME
(especially with polarizable water)
rcoulomb_switch          = 0.0
rcoulomb           = 1.2
epsilon_r            = 15   ; 2.5 (with polarizable water)
vdw_type            = Shift  ;cutoff (for use with Verlet-pairlist)
rvdw_switch      = 0.9
rvdw                    = 1.2  ;1.1 (for use with Verlet-pairlist)

tcoupl                 = v-rescale
tc-grps                = Lipids Solvent PULLING_ATOC
tau_t                    = 1.0  1.0 1.0
ref_t                     = 310 310 310
Pcoupl                = parrinello-rahman
Pcoupltype        = semiisotropic
tau_p                   = 12.0 12.0 12.0 ;parrinello-rahman is more stable
with larger tau-p, DdJ, 20130422
compressibility = 3e-4 3e-4 3e-4
ref_p                    = 1.0 1.0 1.0

gen_vel               = no
gen_temp           = 310
gen_seed            = 473529

constraints         = none
constraint_algorithm  = Lincs
continuation      = no
lincs_order          = 4
lincs_warnangle  = 30

; Pull Code
pull = umbrella
pull-geometry = direction-periodic
pull-start = yes
pull-nstxout = 10
pull-nstfout = 10
pull-ngroups = 1
pull-group0 = Lipids
pull-group1 = PULLING_ATOC
pull-vec1 = 0 0 1
pull-rate1 = 0.01
pull-k1 = 20




---------------------------------------------------------------- End of MDP
1 -----------------------

The second MDP file (the one used for the Umbrella Sampling production run)
is listed here

integrator             = md
dt                            = 0.002
nsteps                   = 10000000
nstcomm              = 10
comm-grps           =

nstxout                  = 50000
nstvout                  = 50000
nstfout                   = 5000
nstlog                     = 1000  ; Output frequency for energies to log
file
nstenergy              = 5000  ; Output frequency for energies to energy
file
nstxtcout               = 5000  ; Output frequency for .xtc file
xtc_precision        = 100
xtc-grps                 =
energygrps           = Lipids Solvent PULLING_ATOC

nstlist                = 5
ns_type             = grid
pbc                    = xyz
rlist                    = 1.4

coulombtype       = Shift  ;Reaction_field (for use with Verlet-pairlist)
;PME (especially with polarizable water)
rcoulomb_switch = 0.0
rcoulomb              = 1.2
epsilon_r               = 15   ; 2.5 (with polarizable water)
vdw_type              = Shift  ;cutoff (for use with Verlet-pairlist)
rvdw_switch         = 0.9
rvdw                       = 1.2  ;1.1 (for use with Verlet-pairlist)

tcoupl                    = v-rescale
tc-grps                   = Lipids Solvent PULLING_ATOC
tau_t                       = 1.0  1.0 1.0
ref_t                        = 310 310 310
Pcoupl                   = parrinello-rahman
Pcoupltype           = semiisotropic
tau_p                      = 12.0 12.0 12.0 ;parrinello-rahman is more
stable with larger tau-p, DdJ, 20130422
compressibility    = 3e-4 3e-4 3e-4
ref_p                       = 1.0 1.0 1.0

gen_vel                      = no
gen_temp                 = 310
gen_seed                  = 473529

constraints                   = none
constraint_algorithm = Lincs
continuation                = no
lincs_order                    = 4
lincs_warnangle          = 30

; Pull Code
pull                     = umbrella
pull-geometry = direction-periodic
pull-start           = yes
pull-nstxout     = 10
pull-nstfout      = 10
pull-ngroups    = 1
pull-group0      = Lipids
pull-group1      = PULLING_ATOC
pull-vec1            = 0 0 1
pull-rate1           = 0
pull-k1                 = 1

---------------------------- END of MDP 2 -------------------


If you have any questions, please ask and I will try my best to answer
them, Thank you,
Alex