FEP Tutorial
In this tutorial you will learn how to perform a Free Energy Perturbation (FEP) calculation in CAST and how to deal with the different options for this task. Before doing this you should make yourself familiar with the theory behind it with the help of some literature as it will not be explained here. It is also recommended to do the MD tutorial before as FEP is based on MD simulations.
- First example: Zero-Transformation
- Bennets Acceptance Ratio (BAR)
- Graphical Analysis (Python required!)
- Heating
- Softcore Potentials
- Conclusion
For the sake of demonstration we will first calculate the free energy difference for a transformation from ethane to ethane, so the free energy should be zero.
Files needed for this example:
- structure file: ethan_FEP.arc
- forcefield parameter file: oplsaa.prm (download from tinker homepage)
- CAST inputfile: CAST.txt
First take a look at the structure file. It looks like this:
13 (oplsaa.prm param.)
1 C -1.022410 1.992230 0.125010 80 2 6 7 8 12 13
2 C 0.259490 1.183960 0.016170 80 1 3 4 5 IN
3 H 0.267580 0.359880 0.736780 85 2 IN
4 H 1.127190 1.820060 0.217110 85 2 IN
5 H 0.372470 0.763140 -0.988080 85 2 IN
6 C -2.251410 1.127300 -0.098230 80 1 9 10 11 OUT
7 H -1.006100 2.801940 -0.613070 85 1
8 H -1.078980 2.455910 1.115950 85 1
9 H -2.225100 0.653100 -1.084410 85 6 OUT
10 H -3.159380 1.735350 -0.036950 85 6 OUT
11 H -2.318590 0.340110 0.659750 85 6 OUT
12 H 0.259490 1.183960 0.016170 85 1 OUT
13 H -2.251410 1.127300 -0.098230 85 1 IN
You see that it contains more atoms than a "normal" ethane molecule (that would be 8 whereas we have 13 here). Some of the atoms are marked with "IN" or "OUT". Those are the atoms that appear (IN) or disappear (OUT) during the simulation. In this case a methyl group (consisting of atoms 6, 9, 10, and 11) disappears and is replaced by H atom 13 while H atom 12 is replaced by a methyl group formed by atoms 2 to 5.
This is your CAST.txt file:
#### GENERAL INFORMATION ####
verbosity 3
name ethan_FEP.arc
outname output
inputtype TINKER
#### TASK ####
task FEP
#### OPTIONS FOR ENERGY INTERFACE ###
interface OPLSAA
paramfile oplsaa.prm
cutoff 12.0
switchdist 10.0
#### OPTIONS FOR MD ####
MDsteps 1000
MDintegrator 0
MDheat 0 300
MDtimestep 0.001
MDtrack 1
MDsnap 10
#### OPTIONS FOR FEP TASK ####
FEPlambda 1.0
FEPdlambda 0.1
FEPvdwcouple 0.8
FEPeleccouple 0.2
FEPvshift 0.5
FEPcshift 0.1
FEPequil 20000
FEPsteps 50000
FEPfreq 1000
The options for general stuff and MD simulations are already explained in the MD tutorial so only the FEP options will be described here:
The transformation is done stepwise in a number of windows that differ in their value λ which normally goes from 0 (educt) to 1 (product). If you don't want to go until the product you can change the value for FEPlambda to a value smaller than 1 (normally you don't do this). The value FEPdlambda defines the size of the window Δλ. In this case it is 0.1 so you have 10 windows.
FEPvdwcouple, FEPeleccouple, FEPvshift and FEPcshift are parameters for softcore-potentials that are used in order to avoid numerical instabilities (details see manual and literature). FEPvshift and FEPcshift control the strength of the softcore-potential. FEPvdwcouple and FEPeleccouple are there because the contributions by Van der Waals-interactions and electrostatic interactions are scaled independently. FEPvdwcouple is the lambda value where the vdw-interactions for appearing atoms reach their maximum, FEPeleccouple is the lambda value where the coulomb-interactions to appearing atoms start. The interactions for disappearing atoms are scaled down in a similar way. After running the program you can see how van der Waals interactions and electrostatic interactions are scaled down in a table at the beginning of the console output:
FEP Coupling Parameters:
Van-der-Waals Coupling:
1 1 1 0 0 0.125
1 1 1 0 0.125 0.25
1 1 0.875 0.125 0.25 0.375
1 0.875 0.75 0.25 0.375 0.5
0.875 0.75 0.625 0.375 0.5 0.625
0.75 0.625 0.5 0.5 0.625 0.75
0.625 0.5 0.375 0.625 0.75 0.875
0.5 0.375 0.25 0.75 0.875 1
0.375 0.25 0.125 0.875 1 1
0.25 0.125 0 1 1 1
0.125 0 0 1 1 1
Electrostatic Coupling:
1 1 0.875 0 0 0
1 0.875 0.75 0 0 0
0.875 0.75 0.625 0 0 0.125
0.75 0.625 0.5 0 0.125 0.25
0.625 0.5 0.375 0.125 0.25 0.375
0.5 0.375 0.25 0.25 0.375 0.5
0.375 0.25 0.125 0.375 0.5 0.625
0.25 0.125 0 0.5 0.625 0.75
0.125 0 0 0.625 0.75 0.875
0 0 0 0.75 0.875 1
0 0 0 0.875 1 1
The first three columns are for the disappearing atoms (first for λ - Δλ, then for λ, then for λ + Δλ), the last three columns for the appearing atoms (in the same order).
The calculation consists of one initialization run whose length is determined by MDsteps (1000) and a number of windows (see above). In each of these windows an equilibration is done whose length is determined by FEPequil (20000) and after that the production with length FEPsteps (50000). During the production run all the conformations are saved to calculate the free energy of each window. FEPfreq controls the output frequency of alchemical.txt and does not affect the calculation.
Now we can start the calculation... We see that several MDs are running for different values of λ.
As outputfiles we get the "normal" MD output output_MD_TRACE.csv and output_MD_SNAP.arc that is quite similar to the one described in the MD tutorial. Furthermore there is a file called alchemical.txt and one with the name FEP_Results.txt.
This last file contains the most important output so we open it first. The first column shows the lambda value for the current window so it goes from 0 to 1. The remaining columns contain the free energy change until this window, each calculated with a slightly different method. The numbers in the second column are calculated by a "normal" forward transformation, the energies in the third column are calculated by backward transformation and those in the fourth column are calculated by Simple Overlap Sampling (SOS) which combines forward and backward transformation and should thus give the best results (for explanation of these methods see literature). The fourth column contains currently only zeros, we will see later what it is for. So what you are probably interested in is the total free energy change of the reaction which can be found in the last line (columns 2 to 4). The numbers there should be something close to zero as we have a zero transformation (in my case 0.0006, -0.0008 and -0.0001 but they can vary a bit as there is some randomness in the MD simulations).
The file alchemical.txt contains some more detailed information about the simulations. For every FEPfreq'th step of every equilibration and production run the following things are given:
- electrostatic interactions for λ - Δλ
- electrostatic interactions for λ
- electrostatic interactions for λ + Δλ
- Van der Waals interactions for λ - Δλ
- Van der Waals interactions for λ
- Van der Waals interactions for λ + Δλ
- temperature
- ΔEpot = U(λ + Δλ) - U(λ) [U = electrostatic + vdW interactions]
- Free Energy difference ΔG calculated from a forward transformation out of ΔEpot
- ΔEpot, back = U(λ) - U(λ - Δλ)
- Free Energy difference ΔGback calculated from a backward transformation out of ΔEpot, back
After each production run the free energy change for the complete window is given and the total free energy change until the current window (sum of all free energy changes until this point). Those values are given only for the forward transformation so they are identical to the values of the second column in the file FEP_Results.txt.
Now you probably want to know what the fifth column in the FEP_Results.txt is for. In order to find out we start the same calculation again but add the following line:
FEPbar 1
The resulting files and their content is quite the same as before but the fifth column is now filled. The values there are also free energy changes but again calculated with a different method, namely with Bennets Acceptance Ratio (BAR) which is a more sophisticated method to combine forward and backward transformation so the results should be even better than with SOS. But as our system is so small this must not be the case. When I tried this calculation I got 0.0000 for the forward transformation, -0.0013 for the backward transformation, -0.0006 for SOS and 0.0018 for BAR so BAR gave the value with the biggest difference from the "real" value of 0.0. Your results might again be different.
If you don't have the possibility to compile CAST with python just jump to the next section!
Now we want to see a few nice images in order to check the validity of our calculation. The production of graphics in CAST is done by the python module matplotlib so you have to use the configuration Python_Release. Furthermore Python 2.7 with matplotlib has to be installed on your system and the file FEP_analysis.py has to be provided. So in the folder where you want to run the program create a subfolder python_modules and copy the python script into it. You can download it from here.
Then take your inputfile CAST.txt (with or without BAR activated) and add the following line:
FEPanalyze 1
Now run CAST (with Python)...
We again get the quite the same output as before but there are 10 .png files which show the distributions of the potential energy difference for the forward and backward transformation for each window. The better they overlap the better is the calculation. As we have quite a small system the overlap should be quite good in our case. There is also a file overlap.txt which gives the overlap for each window in percent. In my case the numbers there ranged from 87% to 96%.
In order to see what happens if we apply a temperature gradient we replace the heating option (MDheat) in the inputfile above by the following lines:
MDheat 0 0.1
MDheat 500 300
As you can see in the outputfile output_MD_TRACE.csv this temperature gradient only applies to the initialization run. After that the temperature is just kept constant over all equilibration and production runs.
In the explanation of the very first example softcore potentials were mentioned. In this section you will see what they are for and get recommendations how to use them. For these examples you need another structurefile ubq_FEP_met2cys.arc and a modified OPLSAA parameterfile oplsaa_mod.prm. Notice that this is not a zero transformation so don't expect the energy difference to be zero.
This is our inputfile CAST.txt:
#### GENERAL INFORMATION ####
verbosity 3
name ubq_FEP_met2cys.arc
outname softcore
inputtype TINKER
#### TASK ####
task FEP
#### OPTIONS FOR ENERGY INTERFACE ####
interface OPLSAA
paramfile oplsaa_mod.prm
cutoff 12.0
switchdist 10.0
#### OPTIONS FOR MD ####
MDsteps 100
MDintegrator 0
MDheat 0 1
MDheat 100 300
MDtimestep 0.001
MDtrack 1
MDsnap 10
#### OPTIONS FOR FEP TASK ####
FEPlambda 1.0
FEPdlambda 0.1
FEPvdwcouple 0.8
FEPeleccouple 0.2
FEPvshift 0.0
FEPcshift 0.0
FEPequil 1000
FEPsteps 5000
FEPfreq 100
As you see we have switched off all softcore potentials by setting FEPvshift and FEPcshift to zero.
If you look at the file alchemical.txt you see that for the runs where lambda = 0 and dlamda = 0.1 the energies of the forward transformation ΔEpot fluctuate very strongly. The reason is that the conformation is found by simulating the molecule with λ = 0 which means that the atoms that are appearing during the mutation are not taken into account. So there is nothing to prevent them from crashing with other atoms. However for the calculation of the free energies the potential energy for the conformation with λ = 0.1 where these atoms are (partly) present is also taken into account. Since the potential energy function has a singularity at r=0 it gets very big if atoms are crashing and that's the reason for the big fluctuations that have a bad influence on the calculation of the free energy.
Another thing that might happen is the explosion of the molecule when the appearing atoms are finally switched on for the conformational search (i. e. in the equilibration of the next step).
To avoid these things it is recommended to use a softcore potential that has no singularity at r = 0 instead of the normal Lennard-Jones Potential. It is not so important to use a softcore potential for the electrostatic interactions because they are normally switched on later (by FEPeleccouple so in our case in window λ = 0.2) when the van der Waals interactions already prevent the atoms from crashing.
In the last example you saw that it is a good idea to use a softcore potential in order to avoid simulation instabilities. This example now shows you that too much softcore is also not a good idea.
So take the inputfile of the last section and replace the FEPvshift option by
FEPvshift 100.0
This causes the simulation to break. (see alchemical.txt, FEP_Results.txt, softcore_MD_TRACE.txt or the console output) The reason for this can be understood if you look at these softcore potentials for different FEPvshift values:
For small values (green and blue) the potential keeps more or less the original form of the Lennard Jones potential but just without the singularity. For bigger values (red and orange) there is no minimum but the lowest energy distance is zero. So as a consequence it is not only not avoided that atoms crash but they are even forced to do so. And in the next window when the softcore potential approximates the "normal" potential better this can be the reason for the program to break.
So you should always choose a sensible value for FEPvshift. If you see the warning "WARNING! You should choose a smaller value for vshift!" at the beginning of the console output (see screenshot below) this is a sign that the value is too large and the softcore potential for your smallest lambda has no minimum anymore.
For some more information about the FEP task also look at our manual (LaTeX file).
Date: 13 Sep 2019