# 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

1. Optimize MM environment for each configuration
2. Perform total energy calculation for each configuration
3. 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:

1. Restart file corresponding to either A or B configuration of the QM region (sharing the same topology file)
2. ESP charges for QM region in .esp format for both configurations
3. 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 .thm file but this time in atomic units.

The same .thm file can also be used to continue from the prior calculation. This will require the presence of

 set qmmm:extend .true.


directive, the .thm file, and the appropriate rst file.