GROMACS 4.6 example: Ethanol solvation with expanded ensemble
| Free Energy Fundamentals |
|---|
 |
|
Methods of Free Energy Simulations
|
| Free Energy How-to's |
|---|
 |
Setting up the calculation with GROMACS
This tutorial assumes knowledge of [GROMACS]. Make sure you actually know how to use GROMACS first.
Setting up the calculations are very similar to standard free energy calculations; you can review a tutorial here for the standard way to perform the simulation. The topology and gro files are exactly the same; only the mdp files and the calls to the programs change.
In this particular case, we again set up nine intermediate states defined to perform the calculation of the absolute solvation free energy of ethanol (also known as the chemical potential at infinite dilution). However, you will only need to run one simulations, as the simulation will go back and forth between the different states.
So, start with the same (ethanol.gro) and (ethanol.top) and a new mdp file (expanded.mdp) files.
Running the calculation with GROMACS
Run grompp and mdrun as normal. Specifically:
grompp -f expanded.mdp -c ethanol.gro -p ethanol.top -o ethanol.X.tpr -maxwarn 4
There may be some warnings, and you'll need to override, hence -maxwarn.
For mdrun, we simply run:
mdrun -deffnm ethanol.X -dhdl ethanol.X.dhdl.xvg
The number of threads can be set as with any other simulation. If you don't set the dhdl file independently, it will be saved to ethanol.X.xvg,
Strictly, you will only need the dhdl files, but looking at the other output files can be useful to determine what is going on in the system if something goes wrong.
Analyze the calculation results
There are two ways to use the outputs to perform free energy calculations.
Expanded ensemble calculations build up the free energies as the simulations are performed. It does this by building up simulation weights as the simulation progresses, so each different thermodynamic state (lambda state) has a different weight. If it didn't have this weight, then the simulation would spend all the time in the lowest free energy states. When it visits each state equally, then the weights will be exactly equal to the free energies. We keep track of this in the log files. For example, after 50 ps, we have:
Step Time Lambda
25000 50.00000 0.00000
MC-lambda information
Wang-Landau incrementor is: 1
N CoulL VdwL Count G(in kT) dG(in kT)
1 0.000 0.000 38 0.00000 4.00000
2 0.200 0.000 34 4.00000 4.00000
3 0.500 0.000 30 8.00000 0.00000
4 1.000 0.000 30 8.00000 3.00000
5 1.000 0.200 27 11.00000 3.00000 <<
6 1.000 0.400 24 14.00000 1.00000
7 1.000 0.600 23 15.00000 2.00000
8 1.000 0.800 21 17.00000 -2.00000
9 1.000 1.000 23 15.00000 0.00000
the << indicate the current state, with the free energies indicated.
at 250 ps, we have:
Step Time Lambda
125000 250.00000 0.00000
MC-lambda information
Wang-Landau incrementor is: 0.0625
N CoulL VdwL Count G(in kT) dG(in kT)
1 0.000 0.000 13 0.00000 4.75000
2 0.200 0.000 25 4.75000 3.87500
3 0.500 0.000 41 8.62500 2.56250 <<
4 1.000 0.000 30 11.18750 2.12500
5 1.000 0.200 26 13.31250 1.62500
6 1.000 0.400 24 14.93750 0.93750
7 1.000 0.600 21 15.87500 -1.81250
8 1.000 0.800 6 14.06250 -6.37500
9 1.000 1.000 6 7.68750 0.00000
Note that the Wang-Landau incrementor (the amount added to the free energies of the current state) is now 0.0625
Much later, we reach:
Step Time Lambda
997000 1994.00000 0.00000
MC-lambda information
Wang-Landau incrementor is: 0.0078125
N CoulL VdwL Count G(in kT) dG(in kT)
1 0.000 0.000 299 0.00000 4.53125
2 0.200 0.000 311 4.53125 3.93750
3 0.500 0.000 335 8.46875 1.96094
4 1.000 0.000 340 10.42969 1.84375 <<
5 1.000 0.200 352 12.27344 1.64844
6 1.000 0.400 347 13.92188 0.64062
7 1.000 0.600 383 14.56250 -2.22656
8 1.000 0.800 450 12.33594 -5.54688
9 1.000 1.000 456 6.78906 0.00000
And a few steps after that, we have:
Step 998200: Weights have equilibrated, using criteria: wl-delta
So the weights now STOP equilibrating. From this point on, the simulation is an equilibrium simulation (it is not history-dependent). This is the point which we should START analyzing the dhdl file.
We can tell how close these free energies are by looking at how close the visits are to flat at the end of the simulation:
Step Time Lambda
20000000 40000.00000 0.00000
MC-lambda information
N CoulL VdwL Count G(in kT) dG(in kT)
1 0.000 0.000 24232 0.00000 4.50000
2 0.200 0.000 21616 4.50000 3.98438
3 0.500 0.000 20533 8.48438 1.96875
4 1.000 0.000 21404 10.45312 1.84375
5 1.000 0.200 20978 12.29688 1.64844
6 1.000 0.400 24308 13.94531 0.64062 <<
7 1.000 0.600 23005 14.58594 -2.22656
8 1.000 0.800 17166 12.35938 -5.54688
9 1.000 1.000 16776 6.81250 0.00000
We can see that the counts are within 2.42/1.68 = 1.44, and kBT*log(1.44) = 0.8 kJ/mol, which is not so bad since the weights were fixed after 2 ns of total simulation at all eight states.
We can now analyze the dhdl file, using pymbar with the dhdl output. Unlike last time, there is only one dhdl file, which has outputs of the form:
@ s0 legend "Thermodynamic state"
@ s1 legend "Energy (kJ/mol)"
@ s2 legend "dH/d\xl\f{} (coul -lambdas)"
@ s3 legend "dH/d\xl\f{} (vdw-lambdas)"
@ s4 legend "\xD\f{}H \xl\f{} (0,0)"
@ s5 legend "\xD\f{}H \xl\f{} (0.2,0)"
@ s6 legend "\xD\f{}H \xl\f{} (0.5,0)"
@ s7 legend "\xD\f{}H \xl\f{} (1,0)"
@ s8 legend "\xD\f{}H \xl\f{} (1,0.2)"
@ s9 legend "\xD\f{}H \xl\f{} (1,0.4)"
@ s10 legend "\xD\f{}H \xl\f{} (1,0.6)"
@ s11 legend "\xD\f{}H \xl\f{} (1,0.8)"
@ s12 legend "\xD\f{}H \xl\f{} (1,1)"
@ s13 legend "pV (kJ/mol)"
0.0000 0 -29136.337 63.503958 19.491314 0.0000000 12.700792 31.751979 63.503958 67.421648 71.375522 75.362115 79.378551 83.422410 1.6650634
0.2000 0 -29341.744 95.496064 -258.31621 0.0000000 19.099213 47.748032 95.496064 85.312750 85.791887 87.838454 90.571444 93.705941 1.6496685
0.4000 0 -28839.430 81.996287 8.7445649 0.0000000 16.399257 40.998143 81.996287 84.283621 87.281402 90.697714 94.393505 98.291373 1.6952030
The first line after the time is the thermodynamic state. We can then distinguish which state each sample is from.
Again you will need to install [[1]]. alchemical-gromacs.py will be in examples directory.
The correct invocation is:
python alchemical-gromacs.py -f directory/prefix -t 300 -p 1 -n 2500 -v > outpufile
'-t' is temperature, '-p' is pressure, '-f' is (prefix of the files, including directory), and '-v' is verbose output (not required, but helpful to understand!)
It will use all files it finds with the given prefix, and it will assume they are numbered in order. In this case, there should only be one file. If there are multiple expanded ensemble files, it will analyze all the data together. We should omit the first 2 ns, since we know that the weights were still equilibrating. Read over the usage for standard calculations.
Understanding the analysis
Let's look at the output file:
The number of files read in for processing is: 1 output is verbose Reading metadata from solvation_direct/outputs/dhdls/ethanol_expanded.dhdl.xvg... Done reading metadata from solvation_direct/outputs/dhdls/ethanol_expanded.dhdl.xvg... Reading solvation_direct/outputs/dhdls/ethanol_expanded.dhdl.xvg...
All standard stuff saying what it's doing. The next part is important. pymbar determines how many of the samples are statistically uncorrelated, and only uses every [math]\displaystyle{ 1/2\tau }[/math] samples. Although it says correlation time,s it's really the twice the correlation time, which is the length of time required between uncorrelated samples.
Now computing correlation times Correlation times: [ 1.13024145 1.53302243 1.43612724 1.08803448 1.55534344 1.37338905 1.40559608 1.08837468 1.083193 ]
number of uncorrelated samples: [24331 17939 19149 25275 17681 20024 19565 25267 25388]
The next part is initalization information. We use the (faster, ignores information Bennett's acceptance ratio to get a quick estimate of the free energies between states, and refine it with multistate Bennett's acceptance ratio. For simple calculations, probably overkill.
relative_change = 1.000 iteration 0 : DeltaF = 4.609 relative_change = 0.000 iteration 1 : DeltaF = 4.609 relative_change = 0.000 Convergence achieved. Converged to tolerance of 8.439389e-13 in 2 iterations (5 function evaluations) DeltaF = 4.609 +- 0.007
Goes on like that for a while. Then we actually start the MBAR calculation. Some diagnostic information at the beginning.
Computing free energy differences... Using embedded C++ helper code. K = 9, L = 9, N_max = 27500, total samples = 194619 There are 9 states with samples. N_k = [24331 17939 19149 25275 17681 20024 19565 25267 25388] Initial dimensionless free energies with method BAR
This is starting iteration. These f_k's are free energies times [math]\displaystyle{ \beta }[/math], so they are dimensionless. Don't get confused between these values and the final values for the free energies of each state.
f_k =
[ 0. 4.60854858 8.63821596 10.52813716 12.39909138
13.90581448 14.60066295 12.653474 7.08947627]
Determining dimensionless free energies by Newton-Raphson iteration.
self consistent iteration gradient norm is 274.37, Newton-Raphson gradient norm is 0.001627
Choosing self-consistent iteration on iteration 0
current f_k for states with samples =
[ 0. 4.60851302 8.63587897 10.52784259 12.39754897
13.90747094 14.60015956 12.65309141 7.08891004]
relative max_delta = 1.600662e-04
self consistent iteration gradient norm is 61.179, Newton-Raphson gradient norm is 1.0307e-05
Choosing self-consistent iteration for lower gradient on iteration 1
current f_k for states with samples =
[ 0. 4.60810077 8.63487764 10.52740478 12.39729092
13.90745661 14.60018318 12.65278436 7.08849961]
relative max_delta = 6.858353e-05
self consistent iteration gradient norm is 15.488, Newton-Raphson gradient norm is 4.3297e-07
Newton-Raphson used on iteration 2
And it goes on like that for a while. It switched to Newton-Raphson once it appears to be converging. This might take 1 step, it might take 5-10 steps. Once it switches to Newton-Raphson, it will go very quickly.
Soon (should just be a few seconds!) we get:
Converged to tolerance of 3.175483e-14 in 5 iterations.
Of 5 iterations, 3 were Newton-Raphson iterations and 2 were self-consistent iterations
Recomputing all free energies and log weights for storage
Final dimensionless free energies
f_k =
[ 0. 4.60763423 8.63398184 10.52688519 12.39691766
13.9073167 14.6002695 12.65255384 7.08801044]
MBAR initialization complete.
Now let's get to the the results:
First, it prints "Deltaf_ij". These are the values of [math]\displaystyle{ \beta\Delta G_{ij} }[/math]. It's a 9x9 matrix, because there are 9 states. The diagonal is zero, because the free energy difference between a state and itself is zero.
Deltaf_ij:
[[ 0. 4.60763423 8.63398184 10.52688519 12.39691766
13.9073167 14.6002695 12.65255384 7.08801044]
[ -4.60763423 0. 4.02634761 5.91925096 7.78928344
9.29968247 9.99263527 8.04491961 2.48037622]
Next, "dDeltaf_ij" entries are the uncertainties in Deltaf_ij
dDeltaf_ij:
[[ 0. 0.00647308 0.01354042 0.01840233 0.0185153 0.01912802
0.02082696 0.02560051 0.0281657 ]
[ 0.00647308 0. 0.00916897 0.01537218 0.01550672 0.01623309
0.0182042 0.02351631 0.02628569]
Finally, the full the results. Six different methods to compare! The rows only run from 0 to 7, because there are only 8 free energy DIFFERENCES -- 0 to 1, 1 to 2, etc.
TI (kJ/mol) TI-CUBIC (kJ/mol) DEXP (kJ/mol) IEXP (kJ/mol) BAR (kJ/mol) MBAR (kJ/mol)
0: 11.546 +- 0.016 11.516 +- 0.017 11.475 +- 0.031 11.588 +- 0.041 11.498 +- 0.017 11.496 +- 0.016
1: 10.326 +- 0.023 10.071 +- 0.026 10.085 +- 0.057 10.009 +- 0.088 10.054 +- 0.024 10.046 +- 0.023
2: 5.416 +- 0.027 4.706 +- 0.040 4.688 +- 0.106 4.961 +- 0.296 4.715 +- 0.027 4.723 +- 0.022
3: 4.716 +- 0.015 4.751 +- 0.015 4.922 +- 0.049 4.665 +- 0.011 4.668 +- 0.010 4.666 +- 0.007
4: 3.679 +- 0.012 3.686 +- 0.013 3.811 +- 0.029 3.746 +- 0.013 3.759 +- 0.011 3.768 +- 0.009
5: 1.419 +- 0.017 1.965 +- 0.020 1.604 +- 0.170 1.759 +- 0.022 1.734 +- 0.016 1.729 +- 0.015
6: -5.418 +- 0.023 -5.320 +- 0.026 -5.115 +- 0.207 -4.947 +- 0.062 -4.858 +- 0.031 -4.859 +- 0.031
7: -13.178 +- 0.021 -13.654 +- 0.025 -13.911 +- 0.048 -13.676 +- 0.169 -13.882 +- 0.020 -13.883 +- 0.020
------------------- ------------------- ------------------- ------------------- ------------------- -------------------
TOTAL: 18.505 +- 0.186 17.721 +- 0.198 17.559 +- 0.305 18.104 +- 0.361 17.688 +- 0.059 17.684 +- 0.070
Things to note: with 9 states, TI is not so good; it's about 1 kJ/mol away from the other values. If we use cubic interpolation to integrate, we improve the agreement with the other algorithms.
Exponential averaging in both directions, insertion and deletion (EXP and DEXP) are not that far off for this example, but the uncertainties are fairly large.
The values for BAR and MBAR are very close! But the uncertainty estimate of BAR is lower. it SEEMS like BAR is better, but that is actually because the BAR estimate is too low. If you ran the experiment lots of times, and computed the sample error you would get a result which is closer to 0.07 than 0.059. This is a known problem with BAR. The estimate of the pairwise differences (0-1, 1-2, 2-3) are all accurate, but the SUM is inaccurate, because the 0->1 and 1->2 calculations both share the same data set, data set 1.
MBAR is generally the lowest variance estimate for a given amount of sampling. In some special cases, TI can get lower, and is frequently almost as good IF enough intermediate states are included. BAR can be almost as good as MBAR, but the uncertainty estimate for the full calculation can frequently be too low.
These are very accurate calculations (in the case of MBAR). Note that they are correct to within 0.07 kJ/mol, which is 0.016 kcal/mol. This free energy can probably only be computed to within 0.05 kcal/mol experimentally.
Hamilton replica exchange
There is no point in doing expanded ensemble with Hamiltonian replica exchange -- it's ALREADY visiting all of the states.