In this repository, I've reported two methods for calculating the free energy surface and the corresponding uncertainty from alchemical metadynamics for a simple problem of argon being alchemically removed from water in GROMACS, with 6 alchemical states. For each method, I have one or two questions about the uncertainty assessment, as elaborated in the Jupyter notebooks (Method_1.ipynb and Method_2.ipynb) in folders Method_1 and Method_2.
- In the folder
Method_1, I've implemented the approach one taught in PLUMED masterclass 21-2 (Exercise 9). - In the folder
Method_2, I used a second method to calculate the uncertainty of the free energy surface, from Equation 10 in the paper Using metadynamics to explore complex free-energy landscapes. - In the folder
input_files, I stored the simulation outputs (includingCOLVARandHILLS_LAMBDA) from a 5ns 1D alchemical metadynamics, which are the inputs for our data analysis inMethod_1andMethod_2.
- As a reference, the benchmark of the free energy difference was obtained from a 5 ns expanded ensemble simulation, which was about -3.137 +/- 0.135 kT. Note that this was performed after the weights were frozen after a short initial simualtion, but that initial simulation was only 0.7 ns.
- As a result of the assessment of the influence of the block size on the uncertainty, 20 ps might be a reasonable block size, as can be seen from the jupyter notebook.
- Using 20 ps as the block size, the free energy difference estimated by Method 1 was -3.20068 +/- 0.05761 kT.
- Using 20 ps as the block size, the free energy difference estimated by Method 2 was -3.20068 +/- 0.06523 kT.
- As demonstrated in
Method_1.ipynb, Method 1 underestimated the uncertainty by a substantial amount compated to running 20 replicates of the simulation. As mentioned in the notebook, the average free energy difference calculated from 20 repetitions of alchemical metadynamics was -3.120 +/- 0.415 kT. - we expect that 1D alchemical metadynamics should have roughly the same performance as expanded ensemble simulations of the same length, which led to a much smaller uncertainty (-3.137 +/- 0.135 kT). We wonder if our way to calculate the free energy difference is correct or if our implementation of alchemical metadynamics is problematic.
- We wonder if our way of implementing Equation 10 in the paper is correct. If so, might you have any guesses for the the reason that it underestimated the true uncertainty, similarly to Method 1? The implementation is in
calculate_free_energy.pyin case that you are interested for reference. - As mentioned
Method_2.ipynb, our implementation of Method 2 was extremely computationally expensive compared to Method 1. We wonder if there is a way to accelerate the calculation of the cumulative bias B(s, t_fill + i * L_b), which was the bottleneck of the method.