Calculating the free energy of binding using AWH approach

User discussions
, ,
1

GROMACS version: 2026.1
GROMACS modification: No

Dear Gromacs community,

I am trying to perform binding free energy calculations of a ligand binding to a protein using the AWH method. I am using a general approach described in this tutorial (https://www.alchemistry.org/wiki/Absolute_Binding_Free_Energy_-_Gromacs_2016).

image
image700×589 67.5 KB

My problem occurs in the E) to F) transition, where I am removing restraints from the complex. I am getting PMF values of zero for all lambda values when removing the restraints from the system (restraint-lambdas = 1.00 0.80 0.60 0.40 0.20 0.00), while keeping the Coulomb and vdW lambdas unchanged (vdw-lambdas = 1.00; coul-lambdas = 1.00).

Here is the gmx awh output:

@ s0 legend "PMF"

@ s1 legend "Coord bias"

@ s2 legend "Coord distr"

@ s3 legend "Ref value distr"

@ s4 legend "Target ref value distr"

@ s5 legend "sqrt(friction metric)"

# AWH metadata: target error = 0.74 kT = 1.87 kJ/mol

# AWH metadata: log sample weight = 0.49

    0.0000  0  0  0.946667  1  1  0

    1.0000  0  0  1.01333  1  1  0

    2.0000  0  0  1.14667  1  1  0

    3.0000  0  0  1.05333  1  1  0

    4.0000  0  0  0.973333  1  1  0

    5.0000  0  0  0.866667  1  1  0

I have checked the dhdl file and I can see that the system is exchanging between the states, so the system being stuck at a single lambda value is not the issue:

@ s0 legend "Thermodynamic state"
@ s1 legend "dH/d\xl\f{} coul-lambda = 1.0000"
@ s2 legend "dH/d\xl\f{} vdw-lambda = 1.0000"
@ s3 legend "dH/d\xl\f{} restraint-lambda = 1.0000"
@ s4 legend "\xD\f{}H \xl\f{} to (1.0000, 1.0000, 1.0000)"
@ s5 legend "\xD\f{}H \xl\f{} to (1.0000, 1.0000, 0.8000)"
@ s6 legend "\xD\f{}H \xl\f{} to (1.0000, 1.0000, 0.6000)"
@ s7 legend "\xD\f{}H \xl\f{} to (1.0000, 1.0000, 0.4000)"
@ s8 legend "\xD\f{}H \xl\f{} to (1.0000, 1.0000, 0.2000)"
@ s9 legend "\xD\f{}H \xl\f{} to (1.0000, 1.0000, 0.0000)"
@ s10 legend "pV (kJ/mol)"
0.0000    0 -362.77856 156.32898 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 63.345505
0.2000    3 -343.76199 57.256371 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 63.364079
0.4000    3 -362.43103 144.63574 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 63.355892
0.6000    5 -331.03439 142.25270 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 63.347473
0.8000    4 -351.15057 124.12928 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 63.360325
1.0000    2 -343.49823 107.01561 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 63.403240
1.2000    3 -369.51450 189.38162 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 63.447716
1.4000    4 -330.99670 135.78236 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 63.464073
1.6000    2 -329.36887 36.234848 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 63.500900
1.8000    1 -329.58658 79.443039 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 63.542671
2.0000    4 -343.81161 48.343582 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 63.544537

Another note, when I am simulating the transition between the ‘disappeared’ and ‘appeared’ ligand (D to E transition on the figure) with restraints applied (i.e, using a gradient of coul-lambdas and vdw-lambdas with restraint-lambdas kept at 1.00), I am getting a reasonable PMF curve.

I am using Boresch restraints specified under [ intermolecular_interactions ]. Here is the excerpt:

; restraints

[ intermolecular_interactions ]

[ bonds ]

;    ai   aj    type  bA         kA       bB    kB

  7334  6557    6     0.564      0.0    0.564   4184.00

[ angles ]

;   ai    aj    ak    type    thA      fcA      thB       fcB

  7331  7334  6557       1   109.398    0.0   109.398     41.84

  7334  6557  6567       1   87.258    0.0   87.258     41.84

[ dihedrals ]

;   ai    aj    ak    al  type    phiA      fcA    phiB      fcB

  7347  7331  7334  6557     2    106.179    0.0    106.179    41.84

  7331  7334  6557  6567     2    79.886    0.0    79.886    41.84

  7334  6557  6567  6569     2    138.827    0.0    138.827    41.84

This is the excerpt of the mdp file:

(What follows is a grossly simplified setup - a single walker, sparse lambdas, very frequent outputs, and high diffusion constant - the aim is just to check if I can even obtain the PMFs in the first place).

; Free energy parameters

free-energy              = yes

couple-moltype           = CARB        

;

sc-power                 = 1

sc-sigma                 = 0.3

sc-alpha                 = 0.5

sc_coul                  = no

;

couple-intramol          = no

couple-lambda1           = vdwq

couple-lambda0           = none

;

init-lambda-state        = 0                     

vdw-lambdas              = 1.00 1.00 1.00 1.00 1.00 1.00

coul-lambdas             = 1.00 1.00 1.00 1.00 1.00 1.00

restraint-lambdas        = 1.00 0.80 0.60 0.40 0.20 0.00

calc-lambda-neighbors    = -1        

; AWH settings

awh                        = yes

awh-potential              = umbrella   

awh-nstout                = 100    

awh-nbias                  = 1        

awh-nstsample              = 100       

awh-nsamples-update        = 10

awh1-target              = constant    

awh1-growth              = exp-linear

; multiple simulations

awh-share-multisim         = no         

awh1-share-group           = 1

awh1-error-init            = 10

awh1-equilibrate-histogram = no     

;

awh1-dim1-diffusion        = 0.01      

;

awh1-ndim                  = 1

awh1-dim1-coord-provider   = fep-lambda 

awh1-dim1-coord-index      = 1 

awh1-dim1-start            = 0

awh1-dim1-end              = 5   

awh1-dim1-cover-diameter   = 0.01

I would be very grateful for your help, because I’m stumped!

2

Is the dhdl.xvg you are showing from the same run of the gmx awh output above? That would be strange, as there are deltaH differences in there.

3

Hi,

I’ve since fiddled with the analysis so I cannot find the exact gmx awh output file, so I’ve now rerun the command (for the same mdrun, so dhdl.xvg is the same).

These are the contents of awh_t1000.xvg:

@ s0 legend "PMF"
@ s1 legend "Coord bias"
@ s2 legend "Coord distr"
@ s3 legend "Ref value distr"
@ s4 legend "Target ref value distr"
@ s5 legend "sqrt(friction metric)"
# AWH metadata: target error = 0.22 kT = 0.56 kJ/mol
# AWH metadata: log sample weight = 0.49
    0.0000  0  0  1.0056  1  1  0
    1.0000  0  0  0.9696  1  1  0
    2.0000  0  0  0.9384  1  1  0
    3.0000  0  0  1.0824  1  1  0
    4.0000  0  0  1.0164  1  1  0
    5.0000  0  0  0.9876  1  1  0

So, the trend seems to be the same.

4

Do you see non-zero restraint energies in the log/energy file?

5

dVrestraint/dl is always zero.

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
dVrestraint/dl                    0          0          0          0  (kJ/mol)
6

Thanks. Can you also see what the restraint lambda value is set to?

7

Restraint lambdas are defined as

restraint-lambdas        = 1.00 0.80 0.60 0.40 0.20 0.00

For reference, this is lambda exchange over time, just to confirm that the exchange is indeed happening:

output
output1200×800 10.1 KB