PMF from SMD simulations using constant force option

User discussions
, ,
1

GROMACS version:
GROMACS modification: Yes/No
Here post your question
Hello Gromacs users,
I would like to perform steered MD simulations using gromacs. I have followed umbrella sampling tutorial in which first the protein is pulled to generate conformations which then ran and analysed using wham. I want to pull small molecule through a channel, for which I want to use option constant force for my pulling, which gives two output xvg files force.xvg, position.xvg. My question is, how to get PMF profile from these output files?

2

Are you doing umbrella sampling or steered MD? Steered MD is often used as a preparatory step for US, but you cannot really get any reasonable free energy estimate from a single steered MD run (even with Jarzynski you’d need at least dozens, if not thousands).

3

Thank you so much for reply. I am doing steered MD simulations, pulling a small molecule through a channel. I have the following question,

  1. My reaction coordinate is com of small molecule and com of dummy kept at just below the center of the channel. The initial distance between them is 38 Angstorms, i have used following pull.mdp file
    integrator = md
    dt = 0.002
    nsteps = 1500000
    nstxout-compressed = 500
    ;nstvout = 50000
    ;nstfout = 50000
    nstcalcenergy = 5000
    nstenergy = 5000
    nstlog = 5000
    ;
    cutoff-scheme = Verlet
    nstlist = 20
    rlist = 1.2
    vdwtype = Cut-off
    vdw-modifier = Force-switch
    rvdw_switch = 1.0
    rvdw = 1.2
    coulombtype = PME
    rcoulomb = 1.2
    ;
    tcoupl = nose-hoover
    tc_grps = Protein wat_cap
    tau_t = 1.0 1.0
    ref_t = 310 310
    constraints = h-bonds
    constraint_algorithm = LINCS
    continuation = yes
    ;; Pull code
    pull = yes
    pull-coord1-dim = N Y N
    pull_ncoords = 1 ; only one reaction coordinate
    pull_ngroups = 2 ; one c6 and three dummies
    pull_coord1_groups = 1 2
    pull_group1_name = HEP
    pull_group2_name = DUM

pull_coord1_type = umbrella ; harmonic potential
pull_coord1_geometry = distance ; simple angle increase
pull-nstfout = 10
pull-nstxout = 10
pull-print-components = no
pull-print-ref-value = yes
pull_coord1_start = yes ; define initial COM distance > 0
pull_coord1_rate = -0.001 ; 0.01 nm per ps = 10 nm per ns
pull_coord1_k = 1000 ; kJ mol^-1 nm^-2
pull-coord1-vec = 0 -1 0
pull-pbc-ref-prev-step-com = yes
according to the above mdp file, the reaction coordinate should cover 3 nm but when I use wham to get pmf. I get pmf from 3.8 nm to 4.2 nm, rest of the reaction coordinate is not covered? Could you please let me know why is it not covering entire reaction coordinate?

4

Your vector is towards negative values of Y, and your pulling rate is negative, so you end up pulling the molecule towards higher values of Y. Remove one of the minus signs and the direction should be fine.

But again, steered MD by itself does not give you the PMF, make sure you run a free energy method such as US after your SMD.

5

With steered MD the only option to get a PMF is using Jarzynski’s equality. You need (very) many simulations though.

6

Thank you so much. I am running simulation and will let you know

7

Thank you so much for a reply. If I Run 20 simulations using the above pull code for 1 ns each. after obtaining force.xvg, pos.xvg files as output. How can we convert using an equation? could you please provide stepwise details or direct me to any of the papers which contain step-by-step details?