GROMACS version: 2025.3
GROMACS modification: No, just the MPI library.
I’m currently researching the dimerization of a fragment of zinc-finger protein (calculations are made for two coupled identical chains of 68 amino acids), and I want to calculate PMF for coupling with AWH.
As a reaction coordinate was chosen the dictance between centers of masses of two aminoacids that are connected with a hydrogen bond together.
My main reference is DNA base pair opening tutorial, in which plot for AWH_pullx.xvg looks like 1) (on the picture).
Covering is reached in less than 2000 ps, and the simulation goes through reaction coordinate as it should.
However, in my case, it looks like 2) (awh1-dim1-force-constant = 400000, awh1-dim1-diffusion = 5e-5):
Or like 3) (awh1-dim1-force-constant = 400000, awh1-dim1-diffusion = 4e-3):
Bigger force constant caused simulation to crash, and i feel like it isn’t the problem, because covering interval is reached, but for some reason not in a way i want it to. Not to mention that the constant is already huge, comparing to similar projects by other people.
Why is covering in my system reached so difficult? Should I keep trying different force constants? Or should I try multidimensional reaction coordinate?
Here’s my mdp file, not much changed from the one from tutorial:
title = AWH simulation
integrator = md
dt = 0.002
nsteps = 50000000
nstlog = 1000
nstenergy = 1000
nstxout-compressed = 2000
; Settings that make sure we run with parameters in harmony with the selected force-field
coulombtype = pme
rcoulomb = 1.0
fourierspacing = 0.121 ; grid spacing for FFT
;
vdwtype = cut-off
vdw-modifier = potential-shift-verlet
rvdw = 1.0
dispCorr = enerpres
constraints = h-bonds
; Pressure coupling is on
pcoupl = parrinello-rahman
pcoupl-type = semiisotropic ; isotropic in the x and y direction, but different in the z direction
tau_p = 5.0
ref_p = 1.0 1.0
compressibility = 4.5e-5 4.5e-5
nstpcouple = 10
; Keep system temperature fluctuating physically correct
tcoupl = v-rescale
tc-grps = system
tau_t = 0.5
ref_t = 300
gen-vel = no
gen-temp = 300
gen-seed = -1
periodic-molecules = yes ; for systems with molecules that couple to themselves through the periodic boundary conditions, this requires a slower PBC algorithm and molecules are not made whole in the output.
pull = yes ; The reaction coordinate (RC) is defined using pull coordinates.
pull-ngroups = 2 ; The number of atom groups needed to define the pull coordinate.
pull-ncoords = 1 ; Number of pull coordinates.
pull-nstxout = 5000 ; Step interval to output the coordinate values to the pullx.xvg.
pull-nstfout = 5000 ; Step interval to output the applied force (skip here).
pull-group1-name = G20 ; Name of pull group 1 corresponding to an entry in an index file.
pull-group2-name = V22 ; Same, but for group 2.
pull-coord1-groups = 1 2 ; Which groups define coordinate 1? Here, groups 1 and 2.
pull-coord1-geometry = distance ; How is the coordinate defined? Here by the COM distance.
pull-coord1-type = external-potential ; Apply the bias using an external module.
pull-coord1-potential-provider = AWH ; The external module is called AWH!
awh = yes ; AWH on.
awh-nstout = 50000 ; Step interval for writing awh*.xvg files.
awh-nbias = 1 ; One bias, could have multiple.
awh1-ndim = 1 ; Dimensionality of the RC, each dimension per pull coordinate.
awh1-dim1-coord-index = 1 ; Map RC dimension to pull coordinate index (here 1–>1)
awh1-dim1-start = 0.3 ; Sampling interval min value (nm)
awh1-dim1-end = 5.1 ; Sampling interval max value (nm)
awh1-dim1-force-constant = 400000 ; Force constant of the harmonic potential (kJ/(mol*nm^2))