CF(Constant force)-SMD(Steered Molecular Dynamics)

User discussions
1

Hi,
I have done cf-smd but I am not sure about the correct mdp parameters required for doing cf-smd these are my used mdp parameters.

define = -DPOSRES -DPOSRES_FC_BB=400.0 -DPOSRES_FC_SC=40.0
integrator = md
dt = 0.02
nsteps = 5000000
nstxtcout = 500
nstvout = 500
nstfout = 500
nstcalcenergy = 10
nstenergy = 100
nstlog = 100
;
cutoff-scheme = Verlet
nstlist = 20
vdwtype = Cut-off
vdw-modifier = Force-switch
DispCorr = EnerPres
rvdw = 0.9
rvdw_switch = 0.7
rlist = 0.9
rcoulomb = 0.9
coulombtype = PME
;
tcoupl = Nose-Hoover
tc_grps = Protein non-Protein
tau_t = 1.0 1.0
ref_t = 300 300
;
pcoupl = Parrinello-Rahman
pcoupltype = isotropic
tau_p = 5.0
compressibility = 4.5e-5
ref_p = 1.0
;
constraints = h-bonds
constraint_algorithm = LINCS
continuation = yes
;
nstcomm = 100
comm_mode = linear
comm_grps = Protein non-Protein
;
pbc = xyz
;
gen_vel = no
;
refcoord_scaling = com
pull = yes
pull_ncoords = 1 ; only one reaction coordinate
pull_ngroups = 2 ; two groups defining one reaction coordinate
pull_group1_name = UNL
pull_group2_name = Protein
pull_coord1_type = constant-force
pull_coord1_geometry = distance ; simple distance increase
pull_coord1_dim = N N Y ; pull along y
pull-coord1-vec = (0.0 0.0 -1.0) ;
pull_coord1_groups = 1 2 ; groups 1 (Chain A) and 2 (Chain B) define the reaction coordinate
pull_coord1_start = yes ; define initial COM distance > 0
pull_coord1_k = 10000 ; kJ mol^-1 nm^-1
pull-group2-pbcatom = 2637
pull-pbc-ref-prev-step-com = yes

Please suggest me if anything I have missed.

2

Most settings look good, but I would suggest to use geometry=direction, especially since you use a negative component in your pull vector.

But why do you want to do constant force pulling?

3

I am trying to do CF-SMD i.e constant force steered molecular dynamics (using constant force parameter) as I want to calculate dissociation constant of a protein system. I have read an article in which dissociation constant is calculated by CF-SMD(https://doi.org/10.1021/acs.jcim.2c01529).
Thank you!

4

I don’t know what you want to dissociate, but such calculations are often very problematic. Adding a biasing force often causes a lot of friction and you can end up with unphysical pathways. A good check to do is running the process in reverse from the obtained end state to see if you get close to the initial state again.

I have no experience with constant-force SMD. You would need to run multiple replicates and use Jarzynski’s equality to estimate the free-energy difference.

Another option is to use the AWH method. This samples back and forward. In many cases of a/dissociation studies, the results from AWH willl not converge, indicating that it is very difficult of impossible to sample the transition on time scales accessible by MD.