GROMACS version: 2025.4
GROMACS modification: No
Hello fellow researchers,
I am currently figuring out wheatear the implemented AWH method is suited for the efficient calculation of the solvation free energy (SFE) simulations in our research group or not. Up to now, we are using Multistate Bennett Acceptance Ratio (MBAR) method as implemented in Alchemical Analyis and a custom tool for the optimization of amount and position of intermediate lambda states. I am looking at the performance and accuracy of reaching a reported experimental value. For this, I have chosen acetylsalicylic acid in TIP3P water as a first simple test system, using the GAFF2/RESP force field in both simulations. The initial distribution of the intermediate states is as follows:
vdw-lambdas = 0.0 0.2 0.4 0.6 0.8 1.0 1.0 1.0 1.0 1.0 1.0 1.0
coul-lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.2 0.4 0.6 0.8 1.0 1.0
From an equilibrated system onwards the MBAR approach takes 25 ns for its optimization, followed by five independent production runs of 1 ns per lambda state (13), resulting in 90 ns total simulation time. To make a fair comparison, I decided to match the length of the AWH simulationen with the 90 ns from the MBAR approach. Therefore, I performed five independent AWH simulations with 4 walkers using awh-share-multisim for 4.5 ns each. Details of the AWH settings are posted below. Most of the settings are based on the tutorial on binder. An exemplary log and the mdp files are attached to this post.
production_mbar.mdp (5.5 KB)
awh_no_covering-diameter.log (163.9 KB)
awh_no_covering-diameter.mdp (5.6 KB)
awh_no_covering-diameter_t4500.xvg (1.7 KB)
awh_with_covering-diameter.log (163.9 KB)
awh_with_covering-diameter.mdp (5.6 KB)
awh_with_covering-diameter_t4500.xvg (1.7 KB)
awh = yes
awh-potential = umbrella
awh-nstout = 50000
awh-nbias = 1
awh-nstsample = 50
awh-nsamples-update = 10
awh1-equilibrate-histogram = no
awh1-target = constant
awh1-growth = exp-linear
awh1-ndim = 1
awh1-dim1-coord-provider = fep-lambda
awh1-dim1-coord-index = 1
awh1-dim1-start = 0
awh1-dim1-end = 11
awh1-dim1-diffusion = 0.001
awh1-error-init = 10
awh-share-multisim = yes
awh1-share-group = 1
1. How do I choose an appropiate covering diameter?
At first, I did not specify a covering diameter, which resulted in the following warning:
WARNING 1 [file awh.mdp, line 113]: When simulations share an AWH bias, it is strongly recommended to use a non-zero covering diameter.
I came across a forum post, in which Magnus Lundborg suggested a value of 8 nm for a lipid membrane system. But the use case was totally different from mine. I was wondering you can suggest a quantity that I can relate this covering diameter to?
Nevertheless, I have run two sets of simulations to get a first impression of this parameter. One with the supressed warning (i.e. a covering diameter of zero) and one simulation with the suggested 8 nm. For each run I plotted the latest awh_t4500.xvg and the development of the SFE result (i.e. the difference of the first and last value of the column PMF) for each timestep.
AWH simulation without covering diameter
Regarding the first plot, I am not entirely sure what to focus on beyond the difference between the first and last lambda state. From what I recall, the friction metric can be interpreted as a measure of sampling efficiency. Does this imply that sampling at lambda state 1 is less efficient? If so, what methods are typically used to improve sampling in these regions?
The second plot shows the convergence of five independent runs toward a experimental reference and the block average (red) of the solvation free energy after 90 ns total simulation time). It seems that there is still a small drift and convergence isn’t reached yet. I am considering to extend the simulation time by 20-30 %.
| Method | Length | SFE (kJ/mol) | Est. core time |
|---|---|---|---|
| MBAR | 90 ns | -44.80 ± 0.82 | ~ 450 hours |
| AWH (covering diameter = 0 nm) | 90 ns | -44.14 ± 0.99 | ~ 890 hours |
| AWH (covering diameter = 8 nm) | 90 ns | -43.56 ± 1.23 | ~ 887 hours |
The comparison of core times suggests that the cover diameter has no noticeable impact on the required computational resources. While the MBAR approach requires roughly half the core time, this difference should be interpreted with caution, as the simulation setups are not fully identical. In particular, different pressure-coupling schemes were employed, the MBAR simulations did not utilize all available cores and were run on GROMACS 2023.0 instead of GROMACS 2025.4.
2. Unclear statement about lambda components
In each log file I find print statements of the so-called current vector of lambda components every few steps. For me it is unclear, what exactly this means.
Current vector of lambda components:[ 0.0000 0.0000 0.6000 1.0000 0.0000 0.0000 0.0000 ]
I first thought this directly relates to the intermediate states I specified in the mdp file with vdw-lambdas or coul-lambdas, but the length of the vector does not match. Additionally, only the third and fourth element of the vector are changing during the simulation.
3. Am I right to assume that this vector gives me information about the type of lambda? (e.g. vdw-lambdas, coul-lambdas, bonded-lambdas, restraint-lambdas, …) I then would interpret that the two changing elements of the vector somehow relate to the specified vdw-lambdas or coul-lambdas while the non-changing elements relate to the non-specified lambdas. What exactly do the numbers refer to?
4. Does AWH take into account the sequential activation of Van der Waals and Coulomb interactions? As you may know it is common practice that you first gradually scale vdw interactions, before touching coulomb interactions when using MBAR (and others). I was wondering if this is taken into account or if this is even relevant (due to the different nature of the method)?
5. Is there something like a ‘final lambda distribution’?
If so, how do I access that?
6. Do you have suggestions/experience on the ratio of number of walkers and OpenMP threads? For now, I only used one node with 96 physical cores and assigned 24 threads to each of the four walkers. I know this is much trial and error, I’m just curious if you have a suggestion here.
Thank you for taking the time to read this rather long post. I look forward to any comments!
