Hello! I am trying to run an AWH simulation of a DNA/RNA hybrid molecule. My pull groups are two atoms in the backbone of the RNA in adjacent nucleotides, the reaction co-ord is the distance between these atoms. My goal is to create a pmf curve along this co-ordinate. The issue is that after I run this simulation for some time, I get the error:
Program: gmx mdrun, version 2023.1
Source file: src/gromacs/applied_forces/awh/pointstate.h (line 355)
Function: gmx::PointState::updateFreeEnergy(const gmx::BiasParams&, double)::<lambda()>
MPI rank: 0 (out of 2)
Assertion failed:
Condition: std::abs(freeEnergy_) < detail::c_largePositiveExponent
Very large free energy differences or badly normalized free energy in AWH
update.
There are several things I’ve tried to fix this issue:
Lowered dt to 0.001 from 0.002
Ran my minimization and equilibration steps with a double precision build of gromacs
Lowered my emtol parameter during minimization to 100 from 1000
added vdw-modifier = Potential-shift to my simulation mdp
Otherwise I am using the mdp files produced by charmm-gui and the awh tutorial on the gromacs site. I am using amber force fields but that hasn’t been an issue before.
No matter the parameters I run it with, the awh simulation will always run for about 70 ns before encountering this error. When i run the same system without awh I can reach 200 ns no problem. I know 70 ns is plenty, im just concerned ive set something up wrong in the awh simulation which is causing the error.
Here is the mdp as I would like it to work and how I originally tried the simulation, thank you for any comments!
integrator = md
dt = 0.002 ;
nsteps = 100000000
nstxtcout = 50000
nstvout = 50000
nstfout = 50000
nstcalcenergy = 100
nstenergy = 1000
nstlog = 1000
;
cutoff-scheme = Verlet
nstlist = 20
vdwtype = Cut-off
vdw-modifier = None ;
DispCorr = EnerPres
rvdw = 0.9
rlist = 0.9
rcoulomb = 0.9
coulombtype = PME
;
tcoupl = Nose-Hoover
tc_grps = SOLU SOLV
tau_t = 1.0 1.0
ref_t = 303.15 303.15
;
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 = SOLU SOLV
;
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 = 0 ; Step interval to output the applied force (skip here).
pull-group1-name = Oxygen43 ; Name of pull group 1 corresponding to an entry in an index file.
pull-group2-name = Phosphate44 ; 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.10 ; Sampling interval min value (nm)
awh1-dim1-end = 0.60 ; Sampling interval max value (nm)
awh1-dim1-force-constant = 128000 ; Force constant of the harmonic potential (kJ/(mol*nm^2))
awh1-dim1-diffusion = 5e-5 ; Estimate of the diffusion (nm^2/ps)
;