modified on 24 August 2010 at 10:33 ••• 83,126 views

QMMM Free Energy Old

From NWChem

Jump to: navigation, search

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}]:

\begin{matrix}
 \Delta W_{A\to B}^{\rm esp} = \sum_{i=0}^{n}\Delta W_{\lambda_i\to\lambda_{i+1}}^{\rm esp}
\end{matrix}

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.

In the input file coordinates of QM region (xyzi files) and ESP charges (esp files) are set using the following directives

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. First a classical MD trajectory is generated around  \,\! \lambda_i point. The parameters of this classical trajectory are set in the md block and must include

data  <integer nsteps>
record coord <integer nfreq>

directives, where nsteps is total number of MD steps and nfreq (typically set to 1) indicates the sampling frequency.

The resulting trajectory containing nsteps/nfreq frames is then processed by the QM/MM module to calculate average energy differences over all the frames. This processesing is controlled by the QM/MN block which must include the following

qmmm
  region mm
  density espfit
  trajectory <string trajectory_file> 1 nsteps/nfreq 1
  ncycles <integer n>
end

where trajectory_file indicates the MD trajectory file that should be processed. This file has an extension .trj and is typically named as <system>.trj. The ncycles setting indicates how many times the sampling process (consisting of generating if MD trajectory and its processing) should be repeated. It is possible to continue free energy sampling from the prior run by using

set qmmm:extend .true.

directive, which will rely on the availability of .thm file containing information from prior calculation.

The solvation contribution calculation to QM/MM free energy is initiated through fep task directive

 task qmmm fep