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

[Justin Lemkul]

On 7/6/16 11:13 AM, Shi Li wrote:
> 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?

Don't use position restraints.  This is probably completely undermining your 
entire process.

> 2. Is the profile.xvg file generated from WHAM directly plotable? Or should I do some modification before ploting it?

It's the result you need.

> 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?
>

WHAM assigns a zero value to the left-most window and determines all other free 
energies relative to it (this is in the WHAM paper, which is where you should 
start before doing anything else!).  You can adjust the zero location with the 
-zprof0 option if it makes more sense to assign the zero position elsewhere 
(e.g. longest distance).

-Justin

-- 
==================================================

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==================================================