QMMM Free Energy¶
Overview¶
Free energy capabilities of QM/MM module are at this point restricted to calculations of free energy differences between two fixed configurations of the QM region.
Users must be warned that this the least automated QM/MM functionality containing several calculation stages. Solid understanding of free energy calculations is required to achieve a meaningful calculation.
Description of the implemented methodology can be found in the following paper. In this approach the free energy difference between the two configurations of the QM region (e.g. A and B):
\(\Delta W_{A\to B}=-1/\beta \ln \langle (e^{-\beta (E_B-E_A)}) \rangle_{A}\)
is approximated as a sum of internal QM contribution and solvation:
\(\Delta W_{A\to B}\approx\Delta W_{A\to B}^{int}+\Delta W_{A\to B}^{solv}\)
It is presumed that structures of A and B configurations are available as restart files sharing common topology file.
Internal contribution¶
The internal QM contribution is given by the differences in the internal QM energies evaluated at the optimized MM environment:
\(\Delta W_{A\to B}^{int}=E_{B}^{int}-E_{A}^{int}\)
The internal QM energy is nothing more but a gas phase expression total energy but evaluated with the wavefunction obtained in the presence of the environment. To calculate internal QM contribution to free energy difference one has to
- Optimize MM environment for each configuration
- Perform total energy calculation for each configuration
- Calculate internal energy difference
Note that internal QM energy can be found in the QM/MM output file under “quantum energy internal” name.
Solvation contribution¶
The solvation contribution is evaluated by averaging energy difference between A and B configurations of the QM system represented by a set of ESP charges.
\(\Delta W_{A\to B}^{solv}=-1/\beta \ln \langle (e^{-\beta (E_B^{ESP}-E_A^{ESP})}) \rangle_{A}\)
where \(E_A^{ESP}\) is the total energy of the system where QM region is replaced by a set of fixed point ESP charges.
In majority of cases the A and B configuration are “too far apart” and one step free energy calculation as shown above will not lead to meaningful results. One solution is to introduce intermediate points that bridge A and B configurations by linear interpolation
\(R_{\lambda_i} = (1-\lambda_i)R_A + \lambda_i R_B\)
\(Q_{\lambda_i} = (1-\lambda_i)Q_A + \lambda_i Q_B\)
where
\(\lambda_i = \frac {i}{n}, \quad i=0,..,n \,\!\)
The solvation free energy difference can be then written as sum of differences for the subintervals \(\,\! [\lambda_i\to\lambda_{i+1}]\) :
\( \Delta W_{A\to B}^{\rm esp} = \sum_{i=0}^{n}\Delta W_{\lambda_i\to\lambda_{i+1}}^{\rm esp} \)
To expedite the calculation it is convenient to use a double wide sampling strategy where the free energy differences for the intervals \(\,\! [\lambda_{i-1}\to\lambda_{i}]\) and \(\,\! [\lambda_i\to\lambda_{i+1}]\) are calculated simultaneously by sampling around \(\,\!\lambda_{i}\) point. In the simplest case where we use two subintervals (n=2)
\(\Delta W_{A\to B}^{solv} \equiv \Delta W_{0\to 1}= \Delta W_{0\to 0.5}^{solv}+\Delta W_{0.5\to 1}^{solv}\)
or
\(\Delta W_{A\to B}^{\rm solv} = -\Delta W_{0.5\to 0}^{\rm solv}+\Delta W_{0.5 \to 1}^{\rm solv}\)
The following items are necessary:
- Restart file corresponding to either A or B configuration of the QM region (sharing the same topology file)
- ESP charges for QM region in .esp format for both configurations
- Coordinates for QM regions in .xyzi format for both configurations
Both .esp and .xyzi files would be typically obtained during the calculation of internal free energy (see above). ESP charges would be generated in the perm directory during optimization of the MM region. The xyzi is basically xyz structure file with an extra column that allows to map coordinates of QM atoms to the overall system. The xyzi file can also be obtained as part of calculation of internal free energy by inserting
set qmmm:region_print .true.
anywhere in the input file during energy calculation. Both xyzi and esp files should be placed into the perm directory!!!
In the input file the restart file is specified in the MD block following the standard notation
md
system < name of rst file without extension>
...
end
while coordinates of QM region (xyzi files) and ESP charges (esp files) are set using the following directives (at the top level outside of any input blocks)
set qmmm:fep_geom xxx_A.xyzi xxx_B.xyzi
set qmmm:fep_esp xxx_A.esp xxx_B.esp
The current interpolation interval \(\,\! [\lambda_i\to\lambda_{i+1}]\)
for which free energy difference is calculated is defined as
set qmmm:fep_lambda lambda_i lambda_i+1
To enable double wide sampling use the following directive
set qmmm:fep_deriv .true.
If set, the above directive will perform both \(\,\! [\lambda_i\to\lambda_{i+1}]\) and \(\,\! [\lambda_i\to\lambda_{i-1}]\) calculations, where
\(\,\! \lambda_{i-1}=\lambda_i - (\lambda_{i+1}-\lambda_i)\)
The calculation proceeds in cycles, each cycle consisting of two phases. First phase is generation of classical MD trajectory around \(\,\! \lambda_i\) point. Second phase is processing of the generated trajectory to calculate averages of relevant energy differences. The number of MD steps in the first phase is controlled by the QM/MM directive
This is a required directive for QM/MM free energy calculations.
Number of overall cycles is defined by the QM/MM directive
In most cases explicit definition of QM/MM density and region should not be required. The QM/MM density will automatically default to espfit and region to mm.
Prior to data collection for free energy calculations user may want to prequilibrate the system, which can be achieved by equil keyword in the MD block:
md
...
equil <number of equilibration steps>
end
Other parameters (e.g. temperature and pressure can be also set in the MD block.
The actual QM/MM solvation free energy calculation is invoked through the following task directive
task qmmm fep
The current value of solvation free energy differences may be tracked though
grep free <name of the output file>
The first number is a forward (\(\,\! [\lambda_i\to\lambda_{i+1}]\))
free energy difference and second number is backward
(\(\,\! [\lambda_i\to\lambda_{i-1}]\)) free energy difference, both in
kcal/mol. The same numbers can also be found in the 4th and 6th columns
of
The same
set qmmm:extend .true.
directive, the
Here is an example of the input file for QM/MM solvation free energy calculation.