2019 Gromacs boresch restraints not working

User discussions
,
1

GROMACS version:2019
GROMACS modification: I don’t know

Dear gromacs community,

I am performing absolute binding free energy calculations. For the past month I had difficulty with applying boresch restraints (1 distance, 2 angles and 3 dihedrals) to my aqueous ligand-protein complex. I have tried it on couple of systems but here I will reference Absolute Binding Free Energy - Gromacs 2016 on alchemistry.org ( Absolute Binding Free Energy - Gromacs 2016 - AlchemistryWiki ) because it will be more familiar. I have run the lambda states in the tutorial without changing any of the mdp or topology files. Snippets of the files are shown below:

[ system ]
; Name
Complex

[ molecules ]
; Compound #mols
Protein_chain_A 1
ligand 1
SOL 10126
CL 8

[ intermolecular_interactions]
[ bonds ]
; ai aj type bA kA bB kB
1391 2615 6 0.654 0.0 0.654 4184.0

[ angles ]
; ai aj ak type thA fcA thB fcB
1393 1391 2615 1 88.8 0.0 88.8 41.84
1391 2615 2614 1 32.9 0.0 32.9 41.84

[ dihedrals ]
; ai aj ak al type thA fcA thB fcB
1410 1393 1391 2615 2 -159.7 0.0 -159.7 41.84
1393 1391 2615 2614 2 122.6 0.0 122.6 41.84
1391 2615 2614 2610 2 12.8 0.0 12.8 41.84

;====================================================
; Production simulation
;====================================================

;----------------------------------------------------
; RUN CONTROL
;----------------------------------------------------
integrator = sd ; stochastic leap-frog integrator
nsteps = 500000 ; 2 * 500,000 fs = 1000 ps = 1 ns
dt = 0.002 ; 2 fs
comm-mode = Linear ; remove center of mass translation
nstcomm = 100 ; frequency for center of mass motion removal

;----------------------------------------------------
; OUTPUT CONTROL
;----------------------------------------------------
nstxout = 0 ; don’t save coordinates to .trr
nstvout = 0 ; don’t save velocities to .trr
nstfout = 0 ; don’t save forces to .trr
nstxout-compressed = 1000 ; xtc compressed trajectory output every 1000 steps (2 ps)
compressed-x-precision = 1000 ; precision with which to write to the compressed trajectory file
nstlog = 1000 ; update log file every 2 ps
nstenergy = 1000 ; save energies every 2 ps
nstcalcenergy = 100 ; calculate energies every 100 steps

;----------------------------------------------------
; BONDS
;----------------------------------------------------
constraint_algorithm = lincs ; holonomic constraints
constraints = all-bonds ; hydrogens only are constrained
lincs_iter = 1 ; accuracy of LINCS (1 is default)
lincs_order = 4 ; also related to accuracy (4 is default)
lincs-warnangle = 30 ; maximum angle that a bond can rotate before LINCS will complain (30 is default)
continuation = yes ; formerly known as ‘unconstrained-start’ - useful for exact continuations and reruns

;----------------------------------------------------
; NEIGHBOR SEARCHING
;----------------------------------------------------
cutoff-scheme = Verlet
ns-type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs (default is 10)
rlist = 1.0 ; short-range neighborlist cutoff (in nm)
pbc = xyz ; 3D PBC

;----------------------------------------------------
; ELECTROSTATICS
;----------------------------------------------------
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
ewald_geometry = 3d ; Ewald sum is performed in all three dimensions
pme-order = 6 ; interpolation order for PME (default is 4)
fourierspacing = 0.10 ; grid spacing for FFT
ewald-rtol = 1e-6 ; relative strength of the Ewald-shifted direct potential at rcoulomb

;----------------------------------------------------
; VDW
;----------------------------------------------------
vdw-type = PME
rvdw = 1.0
vdw-modifier = Potential-Shift
ewald-rtol-lj = 1e-3
lj-pme-comb-rule = Geometric
DispCorr = EnerPres

;----------------------------------------------------
; TEMPERATURE & PRESSURE COUPL
;----------------------------------------------------
tc_grps = System
tau_t = 1.0
ref_t = 300
pcoupl = Parrinello-Rahman
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2 ; time constant (ps)
ref_p = 1.0 ; reference pressure (bar)
compressibility = 4.5e-05 ; isothermal compressibility of water (bar^-1)

;----------------------------------------------------
; VELOCITY GENERATION
;----------------------------------------------------
gen_vel = no ; Velocity generation is off

;----------------------------------------------------
; FREE ENERGY CALCULATIONS
;----------------------------------------------------
free-energy = yes
couple-moltype = ligand
couple-lambda0 = vdw-q
couple-lambda1 = none
couple-intramol = no
separate-dhdl-file = yes
sc-alpha = 0.5
sc-power = 1
sc-sigma = 0.3
init-lambda-state = 12
bonded-lambdas = 0.0 0.01 0.025 0.05 0.075 0.1 0.2 0.35 0.5 0.75 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
coul-lambdas = 0.0 0.00 0.000 0.00 0.000 0.0 0.0 0.00 0.0 0.00 0.0 0.25 0.5 0.75 1.0 1.00 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
vdw-lambdas = 0.0 0.00 0.000 0.00 0.000 0.0 0.0 0.00 0.0 0.00 0.0 0.00 0.0 0.00 0.0 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0
nstdhdl = 100
calc-lambda-neighbors = -1

I have gathered all the files and run a simple gmx bar -f dhdl*.xvg -o -oi command and observed delta G was 0 between the first ten restraint application lambda states as shown below:

WARNING: Some of these results violate the Second Law of Thermodynamics:
This is can be the result of severe undersampling, or (more likely)
there is something wrong with the simulations.

Final results in kJ/mol:

point 0 - 1, DG 0.00 +/- 0.00
point 1 - 2, DG -0.00 +/- 0.00
point 2 - 3, DG -0.00 +/- 0.00
point 3 - 4, DG -0.00 +/- 0.00
point 4 - 5, DG -0.00 +/- 0.00
point 5 - 6, DG 0.00 +/- 0.00
point 6 - 7, DG -0.00 +/- 0.00
point 7 - 8, DG -0.00 +/- 0.00
point 8 - 9, DG 0.00 +/- 0.00
point 9 - 10, DG -0.00 +/- 0.00
point 10 - 11, DG 9.15 +/- 0.12
point 11 - 12, DG 6.93 +/- 0.23
point 12 - 13, DG 4.24 +/- 0.30
point 13 - 14, DG 1.41 +/- 0.55
point 14 - 15, DG 1.96 +/- 0.19
point 15 - 16, DG 1.95 +/- 0.20
point 16 - 17, DG 3.82 +/- 0.24
point 17 - 18, DG 3.72 +/- 0.24
point 18 - 19, DG 4.36 +/- 0.18
point 19 - 20, DG 5.10 +/- 0.11
point 20 - 21, DG 6.03 +/- 0.02
point 21 - 22, DG 3.08 +/- 0.13
point 22 - 23, DG 3.00 +/- 0.17
point 23 - 24, DG 2.82 +/- 0.22
point 24 - 25, DG 1.20 +/- 0.32
point 25 - 26, DG 0.78 +/- 0.33
point 26 - 27, DG -0.94 +/- 0.18
point 27 - 28, DG -2.02 +/- 0.04
point 28 - 29, DG -0.61 +/- 0.02

total 0 - 29, DG 55.96 +/- 1.35

You can see how it should be when restraints are working from tutorial results:

States BAR (kcal/mol) MBAR (kcal/mol)

0 – 1 0.050 ± 0.001 0.045 ± 0.001
1 – 2 0.030 ± 0.001 0.059 ± 0.001
2 – 3 0.078 ± 0.001 0.084 ± 0.001
3 – 4 0.069 ± 0.001 0.072 ± 0.001
4 – 5 0.045 ± 0.001 0.064 ± 0.000
5 – 6 0.145 ± 0.002 0.209 ± 0.001
6 – 7 0.180 ± 0.003 0.236 ± 0.001
7 – 8 0.177 ± 0.005 0.187 ± 0.001
8 – 9 0.223 ± 0.003 0.251 ± 0.001
9 – 10 0.185 ± 0.003 0.204 ± 0.001
10 – 11 2.109 ± 0.010 2.200 ± 0.005
11 – 12 1.558 ± 0.009 1.710 ± 0.006
12 – 13 1.232 ± 0.057 1.278 ± 0.007
13 – 14 0.947 ± 0.033 0.879 ± 0.009
14 – 15 0.516 ± 0.018 0.531 ± 0.008
15 – 16 0.403 ± 0.007 0.463 ± 0.006
16 – 17 0.862 ± 0.026 0.873 ± 0.015
17 – 18 1.066 ± 0.021 1.014 ± 0.015
18 – 19 1.188 ± 0.012 1.183 ± 0.011
19 – 20 1.341 ± 0.013 1.382 ± 0.010
20 – 21 1.610 ± 0.011 1.610 ± 0.008
21 – 22 0.873 ± 0.005 0.878 ± 0.004
22 – 23 0.912 ± 0.006 0.913 ± 0.004
23 – 24 0.955 ± 0.010 0.931 ± 0.005
24 – 25 0.916 ± 0.011 0.925 ± 0.007
25 – 26 0.881 ± 0.011 0.890 ± 0.009
26 – 27 0.857 ± 0.015 0.822 ± 0.011
27 – 28 0.688 ± 0.021 0.708 ± 0.015
28 – 29 0.579 ± 0.020 0.592 ± 0.020

Coulomb: 7.026 ± 0.068 7.478 ± 0.014
vdWaals: 13.646 ± 0.058 13.714 ± 0.042
TOTAL: 20.673 ± 0.089 21.192 ± 0.044

Additionally, I measured the distance between the restrained atoms and see it drifting away after lambda 25. This, in my opinion, clearly shows that these restraints are not working in my Gromacs program. It was already installed on my university cluster so I can’t imagine this has anything to do with my program although I don’t know how to check it. On my own protein-ligand system, the ligand starts to drift away as soon as coulombic attractions are 0.8 turned off. I have tried everything that comes to mind from increasing force constants, changing bond type to 10, changing to restraint-lambdas with [angle_restraints] and [dihedral_restraints], applying restraint forces in initial and final states, etc. Nothing seems to have worked. My ligand always drifts away. Here are snippets from topology files.

[ molecules ]
; Compound #mols
Protein 1
MOL 1
SOL 21057
CL 2

[ intermolecular_interactions ]
; distance restraints
[ bonds ]
; ai aj type bA kA bB kB
3020 8308 6 0.386 0.0 0.386 4184.0

[ angles ]
; ai aj ak type thA fcA thB fcB
3020 8308 8306 1 127.76 0.0 127.76 41.84
3026 3020 8308 1 132.67 0.0 132.67 41.84

[ dihedrals ]
; ai aj ak al type thA fcA thB fcB
3026 3020 8308 8306 2 -8.13 0.0 -8.13 41.84
3020 8308 8306 8317 2 -17.88 0.0 -17.88 41.84
3027 3026 3020 8308 2 -151.14 0.0 -151.14 41.84

OR

[ molecules ]
; Compound #mols
Protein 1
MOL 1
SOL 21057
CL 2

[ intermolecular_interactions ]
; distance restraints
[ bonds ]
; i j type r0A r1A r2A fcA r0B r1B r2B fcB
3020 8308 10 0.386 0.386 10.0 0.0 0.386 0.386 10.0 4184.0

[ angle_restraints ]
; ai aj ak al type thA fcA multA thB fcB multB
3020 8308 8306 8308 1 127.76 0.0 1 127.76 41.840 1
3026 3020 8308 3020 1 132.67 0.0 1 132.67 41.840 1

[ dihedral_restraints ]
; ai aj ak al type phiA dphiA fcA phiB dphiB fcB
3026 3020 8308 8306 1 -8.13 0.0 0.0 -8.13 0.0 41.840
3020 8308 8306 8317 1 -17.88 0.0 0.0 -17.88 0.0 41.840
3027 3026 3020 8308 1 -151.14 0.0 0.0 -151.14 0.0 41.840

I tried combination of everything:

-end of my topology file where I tried everything together:

[ intermolecular_interactions ]
; distance restraints
[ bonds ]
; i j type r0A fcA r0B fcB
3020 8308 6 0.386 418400.0 0.386 418400.00

[ angle_restraints ]
; ai aj ak al type thA fcA multA thB fcB multB
3020 8308 8306 8308 1 127.76 0.0 1 127.76 41.840 1
3026 3020 8308 3020 1 132.67 0.0 1 132.67 41.840 1

[ dihedral_restraints ]
; ai aj ak al type phiA dphiA fcA phiB dphiB fcB
3026 3020 8308 8306 1 -8.13 0.0 0.0 -8.13 0.0 41.840
3020 8308 8306 8317 1 -17.88 0.0 0.0 -17.88 0.0 41.840
3027 3026 3020 8308 1 -151.14 0.0 0.0 -151.14 0.0 41.840

-my production mdp file for lambda 6:

; Production MD
;====================================================

;define = -DPOSRES_ALPHA_HELIX

; RUN CONTROL
;----------------------------------------------------
integrator = sd ; stochastic leap-frog integrator
nsteps = 50000 ; 100 ps
dt = 0.002 ; 2 fs
comm-mode = Linear ; remove center of mass translation
nstcomm = 100 ; frequency for center of mass motion removal

; OUTPUT CONTROL
;----------------------------------------------------
nstxout = 50000 ; save coordinates to .trr every 100 ps
nstvout = 0 ; don’t save velocities to .trr
nstfout = 0 ; don’t save forces to .trr

nstxout-compressed = 1000 ; xtc compressed trajectory output every 500 steps
compressed-x-precision = 1000 ; precision with which to write to the compressed trajectory file
nstlog = 1000 ; update log file every 2 ps
nstenergy = 1000 ; save energies every 2 ps
nstcalcenergy = 100 ; calculate energies every 100 steps

; BONDS
;----------------------------------------------------
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; hydrogens only are constrained
lincs_iter = 1 ; accuracy of LINCS (1 is default)
lincs_order = 6 ; also related to accuracy (4 is default)
lincs-warnangle = 30 ; maximum angle that a bond can rotate before LINCS will complain (30 is default)
continuation = yes

; NEIGHBOR SEARCHING
;----------------------------------------------------
cutoff-scheme = Verlet
ns-type = grid ; search neighboring grid cells
nstlist = 20 ; 20 fs (default is 10)
rlist = 1.2 ; short-range neighborlist cutoff (in nm)
pbc = xyz ; 3D PBC

; ELECTROSTATICS & EWALD
;----------------------------------------------------
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
ewald_geometry = 3d ; Ewald sum is performed in all three dimensions
pme-order = 6 ; interpolation order for PME (default is 4)
fourierspacing = 0.10 ; grid spacing for FFT
ewald-rtol = 1e-6 ; relative strength of the Ewald-shifted direct potential at rcoulomb

; VAN DER WAALS
;----------------------------------------------------
vdw-type = Cut-off
vdw-modifier = Potential-shift-Verlet
verlet-buffer-tolerance = 0.005
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
DispCorr = EnerPres ; apply long range dispersion corrections for Energy and Pressure

; TEMPERATURE COUPLING (SD ==> Langevin dynamics)
;----------------------------------------------------
tc_grps = Protein_MOL Water_ions ; two coupling groups - more accurate
tau_t = 1.0 1.0
ref_t = 298.15 298.15

; PRESSURE COUPLING
;----------------------------------------------------
pcoupl = Parrinello-Rahman
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2 ; time constant (ps)
ref_p = 1.0 ; reference pressure (bar)
compressibility = 4.5e-05 ; isothermal compressibility of water (bar^-1)
refcoord_scaling = com

; VELOCITY GENERATION
;----------------------------------------------------
gen_vel = no

; FREE ENERGY
;----------------------------------------------------
free-energy = yes
sc-alpha = 0.5
sc-power = 1
sc-sigma = 0.3
sc-coul = no
init-lambda-state = 6
couple-moltype = MOL
couple-lambda0 = vdw-q ; only van der Waals interactions
couple-lambda1 = none ; turn off everything, in this case only vdW
couple-intramol = no
bonded-lambdas = 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0
restraint-lambdas = 0.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0
coul-lambdas = 0.0 0.0 0.2 0.4 0.6 0.8 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0
vdw-lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.2 0.3 0.5 0.7 0.8 0.9 0.95 1.0
nstdhdl = 100
dhdl-print-energy = total
calc-lambda-neighbors = -1

-my batch script:

module purge
module load gnu/5.4.0
module load openmpi/1.10.7
module load ccs/gromacs/skylake-gpu/2019
source /opt/ohpc/pub/libs/gnu/openmpi/ccs/gromacs/2019/skylake-gpu/bin/GMXRC

FREE_ENERGY=pwd
echo “Free energy home directory set to $FREE_ENERGY”
MDP=$FREE_ENERGY/FEP_fastmdp
echo “.mdp files are stored in $MDP”

#starting from 1 because I performed 0 as a test case
for (( i=0; i<16; i++ ))
do
LAMBDA=$i

# A new directory will be created for each value of lambda and
# at each step in the workflow for maximum organization.

mkdir Lambdaabsurd_$LAMBDA
cd Lambdaabsurd_$LAMBDA

##############################
# ENERGY MINIMIZATION STEEP  #
##############################
echo "Starting minimization for lambda = $LAMBDA..."



# Iterative calls to grompp and mdrun to run the simulations

gmx grompp -f $MDP/em_$LAMBDA.mdp -c $FREE_ENERGY/md_0_10.gro -p $FREE_ENERGY/topol.top -n $FREE_ENERGY/index.ndx -o min$LAMBDA.tpr -maxwarn 1

gmx mdrun -deffnm min$LAMBDA -nb gpu

sleep 1

#####################
# NVT EQUILIBRATION #
#####################
echo "Starting constant volume equilibration..."


gmx grompp -f $MDP/nvt_$LAMBDA.mdp -c min$LAMBDA.gro -r min$LAMBDA.gro -p $FREE_ENERGY/topol.top -n $FREE_ENERGY/index.ndx -o nvt$LAMBDA.tpr -maxwarn 1

gmx mdrun -deffnm nvt$LAMBDA -nb gpu

echo "Constant volume equilibration complete."

sleep 1

#####################
# NPT EQUILIBRATION #
#####################
echo "Starting constant pressure equilibration..."


gmx grompp -f $MDP/npt_$LAMBDA.mdp -c nvt$LAMBDA.gro -r nvt$LAMBDA.gro -p $FREE_ENERGY/topol.top -t nvt$LAMBDA.cpt -n $FREE_ENERGY/index.ndx -o npt$LAMBDA.tpr -maxwarn 1


gmx mdrun -deffnm npt$LAMBDA -nb gpu
#refcood_scalin = com is inserted to all npt.mdp files. effect of this needs to be investigated

echo "Constant pressure equilibration complete."

sleep 2

#################
# PRODUCTION MD #
#################
echo "Starting production MD simulation..."


gmx grompp -f $MDP/md_$LAMBDA.mdp -c npt$LAMBDA.gro -r npt$LAMBDA.gro -p $FREE_ENERGY/topol.top -t npt$LAMBDA.cpt -n $FREE_ENERGY/index.ndx -o md$LAMBDA.tpr -maxwarn 1

gmx mdrun -deffnm md$LAMBDA -nb gpu
#maxwarn is for using parinello-rahman pressure coupling. this is OK because it is not equilibration step
echo "Production MD complete."

# End
echo "Ending. Job completed for lambda = $LAMBDA"

cd $FREE_ENERGY

done

exit;

Finally, pull-codes for distance, angle and dihedrals seems to be definitely working. I can tell it prevents my molecule from drifting away even in the last lambda. I am not sure if I can couple the pull-codes to restraint-lambdas. I am trying to do it manually by adjusting restraint force constants in the mdp files.

My main question is, why restraints of any kind using [intermolecular_interactions] seem not to be working.

Thank you very much!

2

[ intermolecular_interactions ] had been working in FEP calculations in a different lipid bilayer system built with charmm-gui. I will build my complex in solution using charmm-gui and update if it works again.

Meanwhile, I am attaching the pull code I am using below. This at least works to keep the ligand in the correct position and orientation.

; Production MD
;====================================================

; RUN CONTROL
;----------------------------------------------------
integrator = sd ; stochastic leap-frog integrator
nsteps = 2000000 ; 4 ns
dt = 0.002 ; 2 fs
comm-mode = Linear ; remove center of mass translation
nstcomm = 100 ; frequency for center of mass motion removal

; OUTPUT CONTROL
;----------------------------------------------------
nstxout = 50000 ; save coordinates to .trr every 100 ps
nstvout = 0 ; don’t save velocities to .trr
nstfout = 0 ; don’t save forces to .trr

nstxout-compressed = 1000 ; xtc compressed trajectory output every 500 steps
compressed-x-precision = 1000 ; precision with which to write to the compressed trajectory file
nstlog = 1000 ; update log file every 2 ps
nstenergy = 1000 ; save energies every 2 ps
nstcalcenergy = 100 ; calculate energies every 100 steps

; BONDS
;----------------------------------------------------
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; hydrogens only are constrained
lincs_iter = 1 ; accuracy of LINCS (1 is default)
lincs_order = 6 ; also related to accuracy (4 is default)
lincs-warnangle = 30 ; maximum angle that a bond can rotate before LINCS will complain (30 is default)
continuation = yes

; NEIGHBOR SEARCHING
;----------------------------------------------------
cutoff-scheme = Verlet
ns-type = grid ; search neighboring grid cells
nstlist = 20 ; 20 fs (default is 10)
rlist = 1.2 ; short-range neighborlist cutoff (in nm)
pbc = xyz ; 3D PBC

; ELECTROSTATICS & EWALD
;----------------------------------------------------
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
ewald_geometry = 3d ; Ewald sum is performed in all three dimensions
pme-order = 6 ; interpolation order for PME (default is 4)
fourierspacing = 0.10 ; grid spacing for FFT
ewald-rtol = 1e-6 ; relative strength of the Ewald-shifted direct potential at rcoulomb

; VAN DER WAALS
;----------------------------------------------------
vdw-type = Cut-off
vdw-modifier = Potential-shift-Verlet
verlet-buffer-tolerance = 0.005

rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
DispCorr = EnerPres ; apply long range dispersion corrections for Energy and Pressure

; TEMPERATURE COUPLING (SD ==> Langevin dynamics)
;----------------------------------------------------
tc_grps = Protein_MOL Water_ions ; two coupling groups - more accurate
tau_t = 1.0 1.0
ref_t = 298.15 298.15

; PRESSURE COUPLING
;----------------------------------------------------
pcoupl = Parrinello-Rahman
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2 ; time constant (ps)
ref_p = 1.0 ; reference pressure (bar)
compressibility = 4.5e-05 ; isothermal compressibility of water (bar^-1)
refcoord_scaling = com

; VELOCITY GENERATION
;----------------------------------------------------
gen_vel = no

; Pull code
pull = yes
pull-ngroups = 6
pull-ncoords = 6
;distance restraint
pull-group1-name = 3020
pull-group2-name = 8308
pull-coord1-type = umbrella
pull-coord1-geometry = distance ; simple distance increase
pull-coord1-groups = 1 2
pull-coord1-dim = Y Y Y
pull-coord1-rate = 0.0 ; not pull, just distance restraint at a reference
distance bw 2 groups
pull-coord1-k = 4184.0 ; kJ mol^-1 nm^-2
pull-coord1-init = 0.386 ; define initial COM distance > 0
;angle1 restraint*
pull_group3_name = 8306
pull_coord2_type = umbrella ; Type of pulling (e.g., umbrella, constant_force)
pull_coord2_geometry = angle ; Geometry of the pull coordinate
pull_coord2_dim = Y Y Y ; Dimensions for pulling (Y for active, N for inactive)
pull_coord2_groups = 1 2 2 3 ; Groups defining the angle: (group1-group2) and (group2-group3)
; This defines the angle between vector B-D and vector D-A.
pull_coord2_init = 127.76
pull_coord2_rate = 0.0 ; Pulling rate (e.g., nm/ps or degrees/ps)
pull_coord2_k = 41.84 ; Spring constant (e.g., kJ/mol/nm^2 or kJ/mol/rad^2)

;*angle2 restraint
pull_group4_name = 3026 ; Define your first group

pull_coord3_type = umbrella ; Type of pulling (e.g., umbrella, constant_force)
pull_coord3_geometry = angle ; Geometry of the pull coordinate
pull_coord3_dim = Y Y Y ; Dimensions for pulling (Y for active, N for inactive)
pull_coord3_groups = 4 1 1 2 ; Groups defining the angle: (group4-group1) and (group1-group2)
; This defines the angle between vector B-D and vector D-A.
pull_coord3_init = 132.67 ; Start pulling from the initial configuration
pull_coord3_rate = 0.0 ; Pulling rate (e.g., nm/ps or degrees/ps)
pull_coord3_k = 41.84 ; Spring constant (e.g., kJ/mol/nm^2 or kJ/mol/rad^2)
;***dihedral1 restraint
pull_coord4_type = umbrella ; Type of pulling (e.g., umbrella, constant_force)
pull_coord4_geometry = dihedral ; Geometry of the pull coordinate
pull_coord4_groups = 4 1 1 2 2 3 ; Groups defining the angle: (group4-group1) and (group1-group2)
; This defines the angle between vector B-D and vector D-A.
pull_coord4_init = -8.13 ; Start pulling from the initial configuration
pull_coord4_rate = 0.0 ; Pulling rate (e.g., nm/ps or degrees/ps)
pull_coord4_k = 41.84 ; Spring constant (e.g., kJ/mol/nm^2 or kJ/mol/rad^2)

;***dihedral2 restraint
pull_group5_name = 8317
pull_coord5_type = umbrella ; Type of pulling (e.g., umbrella, constant_force)
pull_coord5_geometry = dihedral ; Geometry of the pull coordinate
pull_coord5_groups = 1 2 2 3 3 5 ; Groups defining the angle: (group4-group1) and (group1-group2)
; This defines the angle between vector B-D and vector D-A.
pull_coord5_init = 165.47 ; Start pulling from the initial configuration
pull_coord5_rate = 0.0 ; Pulling rate (e.g., nm/ps or degrees/ps)
pull_coord5_k = 41.84 ; Spring constant (e.g., kJ/mol/nm^2 or kJ/mol/rad^2)

;***dihedral3 restraint
pull_group6_name = 3027
pull_coord6_type = umbrella ; Type of pulling (e.g., umbrella, constant_force)
pull_coord6_geometry = dihedral ; Geometry of the pull coordinate
pull_coord6_groups = 6 4 4 1 1 2 ; Groups defining the angle: (group4-group1) and (group1-group2)
; This defines the angle between vector B-D and vector D-A.
pull_coord6_init = -151.14 ; Start pulling from the initial configuration
pull_coord6_rate = 0.0 ; Pulling rate (e.g., nm/ps or degrees/ps)
pull_coord6_k = 41.84 ; Spring constant (e.g., kJ/mol/nm^2 or kJ/mol/rad^2)

; FREE ENERGY
;----------------------------------------------------
free-energy = yes
sc-alpha = 0.5
sc-power = 1
sc-sigma = 0.3
sc-coul = no
init-lambda-state = 11
couple-moltype = MOL
couple-lambda0 = vdw-q ; only van der Waals interactions
couple-lambda1 = none ; turn off everything, in this case only vdW
couple-intramol = no

coul-lambdas = 0.0 0.00 0.000 0.00 0.000 0.0 0.0 0.00 0.0 0.00 0.0 0.25 0.5 0.75 1.0 1.00 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
vdw-lambdas = 0.0 0.00 0.000 0.00 0.000 0.0 0.0 0.00 0.0 0.00 0.0 0.00 0.0 0.00 0.0 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0
nstdhdl = 100
dhdl-print-energy = total
calc-lambda-neighbors = -1

3

I realized if I don’t couple the pull-code to restraint-lambda I don’t get the delta G values for the first couple restrictions lambdas with gmx bar.

4

You don’t write which 2019 version you are using. There has been several fixes for missing intermolecular interactions with domain decomposition in 2019.3 and 2019.5.

In general we advise to use the newest version instead of a 6 years old version of GROMACS.

5

Thank you very much for your reply! I am attaching the version information configuration below. We have 2024 gmx_mpi version available to us as well but it is cumbersome to work with since I can’t use it on the login node. I will see if we can request to install 2025 version. In the meantime, 2019 version performed restrictions well with a different system before. Could you tell me if there is a systematic way to tell whether the restrictions work or don’t work except looking at the distances and delta G values between lambda states?! For example, I noticed pull code outputs pullx.xvg and pullf.xvg files automatically which are very helpful. Thank you again for your time!

gmx -version

GROMACS version: 2019
Precision: single
Memory model: 64 bit
MPI library: thread_mpi
OpenMP support: enabled (GMX_OPENMP_MAX_THREADS = 64)
GPU support: disabled
SIMD instructions: AVX_256
FFT library: fftw-3.3.8-sse2-avx-avx2-avx2_128-avx512
RDTSCP usage: enabled
TNG support: enabled
Hwloc support: disabled
Tracing support: disabled
C compiler: /opt/ohpc/pub/compiler/gcc/5.4.0/bin/gcc GNU 5.4.0
C compiler flags: -mavx -O3 -DNDEBUG -funroll-all-loops -fexcess-precision=fast
C++ compiler: /opt/ohpc/pub/compiler/gcc/5.4.0/bin/g++ GNU 5.4.0
C++ compiler flags: -mavx -std=c++11 -O3 -DNDEBUG -funroll-all-loops -fexcess-precision=fast

6

If you have potentials types that only occur in the restraints, you can check if those are non-zero in the output.

Without domain decomposition all 2019 versions should work correctly. You can prepare you files with the 2019 version and run with 2024.

7

My restrictions with [intermolecular-interactions] worked on Gromacs 2024. Thank you very much!