GROMACS version: 2023
GROMACS modification: No
Dear all,
I am doing a WHAM analysis of an umbrella sampling of an enzyme-inhibitor complex and I have some questions.
I did SMD simulations to pull the inhibitor from the enzyme along the X-axis with pull_coord1_k = 1000 kJ mol^-1 nm^-2, at different pulling velocities 0.1, 0.05, 0.02, 0.01, 0.005 Å/ps. I got frames along the pulling axis for US sampling during 10 ns. Then, I run the wham command as:
gmx_mpi wham -it tpr-files.dat -ix pullx-files.dat -o profile_pullx -hist -unit kCal -b 1000 -min 2.375790 -max 5 -nBootstrap 100
Below is a figure of the histograms and the PMF for each pulling velocity.
At first glance, it seems that the pulling velocities of 0.1 and 0.005 Å/ps give the best histogram overlapping. The ∆G was calculated manually as the difference between the minimum and the maximum value of the PMF values (∆G = min(PMF)-max(PMF)).
1 - Does the wham command give the ∆G value directly or should it be calculated as before manually?
2 - Although the ∆G values are reasonable, these are 2x the expected value ≈ 14 kcal/mol. I am aware that this analysis can be difficult and tricky. Could you give some tips to improve this calculation?
My mdp file for the US simulations is:
title = Umbrella pulling simulation
define = -DPOSRES_X ; Restraining the enzyme
; Run parameters
integrator = md
dt = 0.002
tinit = 0
nsteps = 5000000 ; 10 ns
nstcomm = 10
; Output parameters
nstxout-compressed = 5000 ; every 10 ps
nstenergy = 5000
; Bond parameters
constraint_algorithm = lincs
constraints = all-bonds
continuation = yes
; Single-range cutoff scheme
cutoff-scheme = Verlet
nstlist = 20
ns_type = grid
rlist = 1.4
rcoulomb = 1.4
rvdw = 1.4
; PME electrostatics parameters
coulombtype = PME
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 = Nose-Hoover
tc_grps = Protein Non-Protein
tau_t = 1.0 1.0
ref_t = 310 310
; Pressure coupling is on
Pcoupl = Parrinello-Rahman
pcoupltype = isotropic
tau_p = 1.0
compressibility = 4.5e-5
ref_p = 1.0
refcoord_scaling = com
; Generate velocities is off
gen_vel = no
; Periodic boundary conditions are on in all directions
pbc = xyz
; Long-range dispersion correction
DispCorr = EnerPres
; Pull code
pull = yes
pull_ncoords = 1 ; only one reaction coordinate
pull_ngroups = 2 ; two groups defining one reaction coordinate
pull_group1_name = emzyme
pull_group2_name = inhibitor
pull_coord1_type = umbrella ; harmonic potential
pull_coord1_geometry = distance ; simple distance increase
pull_coord1_dim = Y N N
pull_coord1_groups = 1 2
pull_coord1_start = yes ; define initial COM distance > 0
pull_coord1_rate = 0.0 ; restrain in place
pull_coord1_k = 1000 ; kJ mol^-1 nm^-2
I am considering repeating the calculations with a different pulling constant to improve the spacing of the initial configurations for the US.
Thanks in advance.