Advice needed: AWH and Pull code debugging

GROMACS version: 2023.2
GROMACS modification: No
Here post your question
Hi,

I’m working on pulling a peptide into a lipid bilayer using AWH and the pull code in GROMACS. Below is my complete configuration for the pull and AWH settings:

; Pull code settings
pull = yes
pull-ngroups = 2
pull-ncoords = 1
pull-nstxout = 5000
pull-nstfout = 0
pull-pbc-ref-prev-step-com = yes
pull-group1-name = Pept_middle
pull-group2-name = MEMB
pull-group1-pbcatom = 139
pull-group2-pbcatom = 33614
pull-coord1-k = 1000
pull-coord1-groups = 1 2
pull-coord1-geometry = cylinder
pull-cylinder-r = 1.5
pull-coord1-type = external-potential
pull-coord1-potential-provider = AWH
pull-coord1-vec = 0 0 -1
pull-coord1-start = yes
pull-coord1-rate = 0

; AWH settings
awh = yes
awh-potential = convolved
awh-nstout = 50000
awh-nbias = 1
awh1-ndim = 1
awh1-dim1-coord-index = 1
awh1-dim1-start = -2.6 ; Sampling interval min value (nm)
awh1-dim1-end = 2.6 ; Sampling interval max value (nm)
awh1-dim1-force-constant = 10000 ; Force constant of the harmonic potential (kJ/(mol*nm^2))
awh1-dim1-diffusion = 5e-5 ; Estimate of the diffusion (nm^2/ps)
awh1-dim1-coord-provider = pull
When I run mdrun, I encounter the following error:
Assertion failed:
Condition: weightSum > 0
zero probability weight when updating AWH probability weights.

I tried adjusting the awh1-dim1-start and awh1-dim1-end parameters to -1.6 and 1.6 nm, but then I received a different error:

Coordinate 1 of an AWH bias has a value 2.642145 which is more than 10 sigma
out of the AWH range of [-1.600000, 1.600000]. You seem to have an unstable
reaction coordinate setup or an unequilibrated system.

I’m new to this, and I’m having trouble to understand what might be going wrong. Any guidance or suggestions would be greatly appreciated.

Thank you!

It seems like the initial value of your reaction coordinate is 2.642. If so, the you either need to make the AWH range include this value or you need to run a short initial run with only pulling to get the value into the AWH interval.

I see now that you got an assertion failure. That should never happen. Could you provide us with the input to run grompp, so we can find and fix this issue?

Thank you for your response.

When I attempted to include the reaction coordinate in the awh1-dim1-start and awh1-dim1-end parameters, I encountered an assertion failure. I also tried running a short run, but that didn’t work either, and I ended up getting the same errors again.

Please find attached the config file and the input files required to run grompp:

step8_awh.mdp (2.3 KB)

topol.top (885 Bytes)

Do you need the .gro, .tpr, and .ndx files? These file formats are not supported in this forum and are too large for Pastebin. I can send them to you via email if you need them.

Yes, I would need all files those files. You can mail them to hess at kth.se

With a debug build I see a division by zero far before the AWH issue. This is caused by your incorrect assignment of pull groups. With geometry=cylinder, the membrane should be group 1, not group 2. We should add a check for non-zero cylinder group mass. I expect fixing this will also make the assertion failure go away.

Hi Berk,

I wanted to extend a big thank you for your help! After adjusting the pull groups, everything is now working perfectly. I’ve shared my .mdp file below, so if anyone else runs into a similar issue, they can use it as a helpful reference. I really appreciate the support!

; General settings
integrator = md
dt = 0.002
nsteps = 500000000

; Output settings
nstxout = 50000
nstvout = 50000
nstfout = 50000
nstcalcenergy = 100
nstenergy = 1000
nstlog = 1000

; Nonbonded settings
cutoff-scheme = Verlet
nstlist = 20
rlist = 1.4
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw_switch = 1.0
rvdw = 1.4
coulombtype = PME
rcoulomb = 1.4
fourierspacing = 0.12
dlb = yes

; Thermostat and barostat settings
tcoupl = Nose-Hoover
tc_grps = system
tau_t = 1.0
ref_t = 310

pcoupl = Parrinello-Rahman
pcoupltype = semiisotropic
tau_p = 5.0
compressibility = 4.5e-5 4.5e-5
ref_p = 1.0 1.0

; Constraints
constraints = h-bonds
constraint_algorithm = LINCS
continuation = yes

; Pull code settings
pull = yes
pull-ngroups = 2
pull-ncoords = 1
pull-nstxout = 5000
pull-nstfout = 0
pull-pbc-ref-prev-step-com = yes
pull-group1-name = MEMB
pull-group2-name = Pept_middle
pull-group1-pbcatom = 33614
pull-group2-pbcatom = 139
pull-coord1-k = 3000
pull-coord1-groups = 1 2
pull-coord1-geometry = cylinder
pull-cylinder-r = 1.5
pull-coord1-type = external-potential
pull-coord1-potential-provider = AWH
pull-coord1-vec = 0 0 -1
pull-coord1-dim = N N Y
pull-coord1-start = yes
pull-coord1-rate = 0

; AWH settings
awh = yes
awh-potential = convolved
awh-nstout = 50000
awh-nbias = 1
awh1-ndim = 1
awh1-dim1-coord-index = 1
awh1-dim1-start = -2.6 ; Sampling interval min value (nm)
awh1-dim1-end = 2.6 ; Sampling interval max value (nm)
awh1-dim1-force-constant = 10000 ; Force constant of the harmonic potential (kJ/(mol*nm^2))
awh1-dim1-diffusion =0.25e-3 ; Estimate of the diffusion (nm^2/ps)
awh1-dim1-coord-provider = pull

Hello everyone, I’m now facing a new problem. The simulation is working very well, but when calculating the PMF, I’m only getting positive values. Could this be due to incorrect parameters? Or is there something I need to adjust in the MDP file?

PMFs from AWH are calibrated so that the lowest value is set to 0. So, the PMF you get is the free energy difference compared to the most favourable state. You can shift the PMF yourself so that the 0 corresponds to a reference state, if that is what you are interested in.

Ah, that makes sense – thank you for the explanation!