[gmx-users] Problem in PMF simulation

[Shi Li]

Dear Gromacs users,

I posted this question a few days ago, but still couldn’t solve it. This problem has bothered my for a long time, so I am reposting this for help. 

I am currently doing some PMF simulations for two small molecular types in vacuum. I applied an energy minimization, a NVT equilibrium and a production run on my system. Following are the mdp files for each. I used a commend in my run script to change the Window in pull_coord1_init in the pull code for each sample window. 

Energy minimization
---
integrator               = steep
nsteps                   = 20000

nstenergy                = 500
nstlog                   = 500
nstxout-compressed       = 1000

cutoff-scheme            = group

coulombtype              = Cut-off
rcoulomb                 = 0

vdwtype                  = Cut-off
rvdw                     = 0
rlist                    = 0
pbc                      = no

; Pull code
pull             = umbrella
pull_geometry    = distance
pull_start       = no
pull-ncoords     = 1
pull_group1_name = ADT1
pull_group2_name = ADT2
pull-coord1-groups = 1 2
pull_coord1_init = WINDOW
pull_coord1_rate = 0.0
pull_coord1_k = 5000

—
NVT
----
integrator               = md
dt                       = 0.002     ; 2 fs
nsteps                   = 500000    ; 1.0 ns

nstenergy                = 200
nstlog                   = 2000
nstxout-compressed       = 10000

continuation             = yes
constraint-algorithm     = lincs
constraints              = h-bonds

cutoff-scheme            = group
rlist                    = 0

rcoulomb                 = 0

rvdw                     = 0
pbc                      = no
tcoupl                   = V-rescale 
tc-grps                  = System
tau-t                    = 2.0
ref-t                    = 298.15
nhchainlength            = 1

; Pull code
pull             = umbrella
pull_geometry    = distance
;pull_dim         = N N Y
pull_start       = no
pull-ncoords     = 1
pull_ngroups     = 2
pull_group1_name = ADT1
pull_group2_name = ADT2
pull-coord1-groups = 1 2
pull_coord1_init = WINDOW
pull_coord1_rate = 0.0
pull_coord1_k = 5000

——
Production run:
----
integrator               = md
dt                       = 0.002   ; 2 fs
nsteps                   = 10000000   ; 20.0 ns

nstenergy                = 5000
nstlog                   = 5000
nstxout-compressed       = 5000

continuation             = yes
constraint-algorithm     = lincs
constraints              = h-bonds

cutoff-scheme            = Group

rcoulomb                 = 0
rlist                    = 0

rvdw                     = 0
pbc                      = no

tcoupl                   = V-rescale; Nose-Hoover
tc-grps                  = System
tau-t                    = 2.0
ref-t                    = 298.15
nhchainlength            = 1

; Pull code
pull             = umbrella
pull_geometry    = distance
pull_start       = no
pull-ncoords     = 1
pull_ngroups     = 2
pull_group1_name = ADT1
pull_group2_name = ADT2
pull-coord1-groups = 1 2
pull_coord1_init = WINDOW
pull_coord1_rate = 0.0
pull_coord1_k = 5000

———
I have a very strange profile and histogram. I uploaded to these links.

PMF: https://www.dropbox.com/s/i10hn1p30w2j71r/F-anti-adt.eps?dl=0
HISTO: https://www.dropbox.com/s/1xi9i4cj60g97n2/F-anti-histo.eps?dl=0

The molecule is about 1.5 nm long, so the smooth profile before 1.5 nm looks fine, but I don’t understand why the profile had this fluctuation after 1.5 nm and the histogram seems to have some split peaks after this length. Can anyone tell me why this happened? Is there something wrong in my mdp file or something wrong in my topology? I have used this similar mdp files for my other types of molecular systems and the results were fine. So it is really confused me now. 

I appreciate for any help! Thanks!

Shi