[gmx-users] Production Run mdp file for Umbrella Sampling

[Nidhin Thomas]

Hello GROMACS Users,

I am trying to run an umbrella sampling simulation of a lipid bilayer-protein system. The objective of the simulation is to obtain the free energy difference in inserting a protein into the bilayer. I created the systems with different protein insertion depths using CHARMM-GUI. I obtained the mdp files from CHARMM-GUI as well. However, I modified the mdp files based on Dr. Lemkul’s umbrella sampling simulation tutorial. If I understand correctly, the only constraint for each of my simulations should be just the relative distance of PROT and MEMB groups (reaction coordinate).

But I have few doubts regarding the mdp file and would like to clarify them before the simulation. 

1. Do I have to maintain the restraints on the reference group (MEMB) during production run as well ? (I assume that Dr. Lemkul maintained the restraints on chain  B during production run as well)

2. Does this command "pull_coord1_start   = yes" make sure that PROT com and MEMB com are constrained production runs?

3. Would this "pull_coord1_k  = 1000" be the constraint constant during the production run?

I have attached the mdp file below.

It would be a great help if anyone could help me out here.

integrator             = md
;define                  = -DPOSRES_MEMB
dt                          = 0.002
nsteps                  = 10000000
nstlog                  = 1000
nstxout                 = 5000
nstvout                 = 5000
nstfout                 = 5000
nstcalcenergy           = 100
nstenergy               = 1000
;
cutoff-scheme           = Verlet
nstlist                 = 20
rlist                   = 1.2
coulombtype             = pme
rcoulomb                = 1.2
vdwtype                 = Cut-off
vdw-modifier            = Force-switch
rvdw_switch             = 1.0
rvdw                    = 1.2
;
tcoupl                  = Nose-Hoover
tc_grps                 = PROT   MEMB   SOL_ION
tau_t                   = 1.0    1.0    1.0
ref_t                   = 303.15 303.15 303.15
;
pcoupl                  = Parrinello-Rahman
pcoupltype              = semiisotropic
tau_p                   = 5.0
compressibility         = 4.5e-5  4.5e-5
ref_p                   = 1.0     1.0
;
constraints             = h-bonds
constraint_algorithm    = LINCS
continuation            = yes
;
nstcomm                 = 100
comm_mode               = linear
comm_grps               = PROT   MEMB   SOL_ION
;
refcoord_scaling        = com

; Pull code             
pull                    = yes       
pull_ngroups            = 2
pull_ncoords            = 1
pull_group1_name        = MEMB
pull_group2_name        = PROT
pull_coord1_type        = umbrella      ; harmonic biasing force
pull_coord1_geometry    = distance      ; simple distance increase
pull_coord1_groups      = 1 2
pull_coord1_dim         = N N Y
pull_coord1_rate        = 0.0
pull_coord1_k           = 1000          ; kJ mol^-1 nm^-2
pull_coord1_start       = yes           ; define initial COM distance > 0

Thanks,

Nidhin Thomas