[gmx-users] Problem in umbrella sampling for PMF calculation

[Shi Li]

Dear GMX users,

I am a beginner of Gromacs and I am trying to calculate the potential of mean force in my simulation system. In my simulation box I have two small molecules(70 atoms each). I am fixing one molecule and moving the other approaching to it along Z axis. I use following .mdp code to push these two molecules together. 
-------------------
title       = Umbrella pulling simulation 
define      = -DPOSRES_B
; Run parameters
integrator  = md
dt          = 0.002
tinit       = 0
nsteps      = 250000    ; 500 ps
nstcomm     = 10
; Output parameters
nstxout     = 5000     
nstvout     = 5000 
nstfout     = 500
nstxtcout   = 500      
nstenergy   = 500
; Bond parameters
constraint_algorithm    = lincs
constraints             = all-bonds
continuation            = yes     
; Single-range cutoff scheme
cutoff-scheme = group
nstlist     = 1 
ns_type     = grid
nstcalclr    =1 
rlist       = 0
rcoulomb    = 0
rvdw        = 0
;; PME electrostatics parameters
coulombtype     = cut-off 
fourierspacing  = 0.12
fourier_nx      = 0
fourier_ny      = 0
fourier_nz      = 0
pme_order       = 4
ewald_rtol      = 1e-5
optimize_fft    = yes
; Berendsen temperature coupling is on in two groups
Tcoupl      = Berendsen 
tc_grps     = System  
tau_t       = 0.5       
ref_t       = 310       
; Pull code
pull             = umbrella
pull_geometry    = distance;
pull_dim         = N N Y
pull_start       = yes
pull-ncoords     = 1
pull_ngroups     = 2
pull_group1_name = core2; moving molecule
pull_group2_name = core1; fixed molecule
pull-coord1-groups = 2 1
pull_coord1_init = 3
pull_coord1_rate = 0.005
pull_coord1_k = 100
pull_nstxout    = 1000      ; every 2 ps
pull_nstfout    = 1000      ; every 2 ps
;-----------

After the pushing process, I selected the configurations(0.2 nm each) from the result and did the equilibrium and umbrella sampling with following mdp files.

(Equilibrium)
title       = Umbrella pulling simulation 
define      = -DPOSRES
; Run parameters
integrator  = md
dt          = 0.002
tinit       = 0
nsteps      = 5000   
nstcomm     = 10
; Output parameters
nstxout     = 500      ;
nstvout     = 500 
nstfout     = 500
nstxtcout   = 500       s
nstenergy   = 500
; Bond parameters
constraint_algorithm    = lincs
constraints             = all-bonds
continuation            = no       
; Single-range cutoff scheme
nstlist     = 1
cutoff-scheme = group 
ns_type     = grid
rlist       = 0	
rcoulomb    = 0
rvdw        = 0
;; PME electrostatics parameters
coulombtype     = cut-off
fourierspacing  = 0.12
fourier_nx      = 0
fourier_ny      = 0
fourier_nz      = 0
pme_order       = 4
ewald_rtol      = 1e-5
optimize_fft    = yes
; Berendsen temperature coupling is on in two groups
Tcoupl      = Berendse
tc_grps     = System   
tau_t       = 0.5       
ref_t       = 310       
; Pull code
pull             = umbrella
pull_geometry    = distance 
pull_dim         = N N Y
pull_start       = yes
pull-ncoords     = 1
pull_ngroups     = 2 
pull_group1_name = core2
pull_group2_name = core1
pull-coord1-groups = 1 2 
pull_coord1_rate = 0.0
pull_coord1_k = 100 
pull_nstxout    = 1000     
pull_nstfout    = 1000      
;
----
(Umbrella sampling):
title       = Umbrella pulling simulation 
define      = -DPOSRES
; Run parameters
integrator  = md
dt          = 0.001
tinit       = 0
nsteps      = 5000000  
nstcomm     = 10
; Output parameters
nstxout     = 50000   
nstvout     = 50000 
nstfout     = 5000
nstxtcout   = 5000       
nstenergy   = 5000
; Bond parameters
constraint_algorithm    = lincs
constraints             = all-bonds
continuation            = yes       
; Single-range cutoff scheme
cutoff-scheme = group
nstlist     = 5 
ns_type     = grid
rlist       = 0
rcoulomb    = 0
rvdw        = 0
; PME electrostatics parameters
coulombtype     = cut-off 
fourierspacing  = 0.12
fourier_nx      = 0
fourier_ny      = 0
fourier_nz      = 0
pme_order       = 4
ewald_rtol      = 1e-5
optimize_fft    = yes
; Berendsen temperature coupling is on in two groups
Tcoupl      = Berendsen 
tc_grps     = System   
tau_t       = 0.5 
ref_t       = 310
pull             = umbrella
pull_geometry    = distance
pull_dim         = Y Y Y
pull_start       = yes
pull-ncoords     = 1
pull_ngroups     = 2 
pull_group1_name = core2
pull_group2_name = core1
pull-coord1-groups =  1 2
pull_coord1_rate = 0.0
pull_coord1_k = 100 
pull_nstxout    = 1000    
pull_nstfout    = 1000     

----------
After these, I did WHAM to generate the profile file and plot it. I got a curve went from 0 to negative as the distance increase, which didn't seem to be correct.

My questions are:
1, What is wrong with my input mdp files? What should I change in the mdp files for a simple simulation system in vacuum?
2. Is the profile.xvg file generated from WHAM directly plotable? Or should I do some modification before ploting it?
3. I thought the PMF should go close to zero as the distance increase, should I manual set a reference distance value to zero and modify the data so it is more realistic?


Thank you all very much!
Shi