Free energy calculations

The computation of free energy can be applied to a variety of chemical and biochemical phenomena [192]. Solvation processes, molecular stability, molecular association such as ligand binding ...etc. We are interested in free energy differences between reactants, products and transition states. While in gas phase and small molecular systems free energy differences can be computed by means of analytical expressions, in condensed phase we have to compute the Potential of Mean Force (PMF) along a distinguished coordinate using sampling techniques.

Free energy calculations for small molecules:
Before describing the PMF calculation in condensed phase systems, it can be useful to indicate the particular case of free energy calculations in small molecular systems. In the canonical ensemble (NVT) the Helmholtz free energy $F$ is computed from the canonical partition function $Z$

$$F(NVT)=-k_BT \ln(Z)$$ (2.114)

For $N$ independent and identical particles in the classical limit, $Z$ is computed from the partition function of a single molecule $z$.

$$Z=\frac{z^N}{N!}\quad ; \quad z=\sum_i^{states} e^{E_{tot}i/k_BT}$$ (2.115)

The total energy of every molecular state $E_{tot}i$ can be approximated as a sum involving translational, rotational, vibrational and electronic states. This assumption implies that the molecule partition function $z_{tot}$ can be written as a product of terms

$$E_{tot} = E_{trans} + E_{rot} + E_{vib} + E_{elec}$$

$$z_{tot} = z_{trans} \cdot z_{rot} \cdot z_{vib}\cdot z_{elec}$$ (2.116)

In small molecular systems the partition functions can be approximated by analytical expressions. The term $z_{trans}$ is computed with the free particle model, $z_{rot}$ as the rigid rotor and the $z_{vib}$ is described as a factorization of normal modes of vibration within the harmonic oscillator. Some improvements exist for anharmonic and hindered rotor [193] models but in any case an exploration of the phase space is never performed.

In non-rigid condensed phase systems, the independent particle assumption $Z=z^N/N!$ is not valid, and the translational partition function $z_{trans}$ cannot be modeled as a free particle. Therefore, we will have to compute the whole partition function through the phase space integral which in the NVT ensemble takes into account all the degrees of freedom of the configuration space.

$$Z=\frac{1}{h^{3N}}\int_{\bf q}\int_{\bf p} d{\bf q}d{\bf p} ^{- \frac{ H( {\bf q},{\bf p} ) }{k_BT} }$$ (2.117)

Note that if we consider that our particle system is made up of atoms (nuclei), the vibrational and translational will be considered altogether.

The usual Hamiltonian expressions permits to separate the atomic momenta from the potential energy.

$$H({\bf p},{\bf q}) = \sum_i^N \frac{{\bf p}_i^2}{2m_i} + V({\bf q})$$ (2.118)

Therefore, the integral over the phase space $(p,q)$ becomes the configurational integral multiplied by a constant

$$Z\propto \int_{\bf q} d{\bf q} ^{- \frac{ V( {\bf q} ) }{k_BT} }$$ (2.119)

We could consider our condensed phase system (solution or biomolecule) as a supermolecule and apply the same formula than equations 1.116. However, when computing the vibrational motion of this supermolecule the local and harmonic approximation is not valid since there are many minima separated by barriers lower than the k$_B$T factor. In consequence we still need to compute the internal nuclear motion by MD or MC techniques until convergence of the integral of configuration.

After the assumption of Born-Oppenheimer approximation the nuclei have a classical behavior and it can be entirely reproduced by classical statistical mechanics. However, it has been seen that the absence of the zero point energy and the effects of quantized vibrational motion not contemplated in classical statistical mechanics can be a source of error, mainly when computing activation free energies in reactions involving hydrogen transfer [194].