AWH assertion failed error

GROMACS version: 2023.1
GROMACS modification: No

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:

  1. Lowered dt to 0.001 from 0.002
  2. Ran my minimization and equilibration steps with a double precision build of gromacs
  3. Lowered my emtol parameter during minimization to 100 from 1000
  4. 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.

Any insights are appreciated, thanks!

Could you post your AWH mdp parameters?

1 Like

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)
;

My guess is that a minimum distance of 0.1 nm leads to extremely high energies and forces which could cause instabilities in the simulation.

1 Like