Calculation of molecular free energies in classical potentials

Free energy calculations have a variety of applications which include binding, solvation and chemical reactions. For the applications of solvation and binding equilibrium and non-equilibrium methods are used. In equilibrium methods (alchemical transformations) a molecule is transformed into another through non-realistic intermediate systems. The average of the derivative of the interpolating Hamiltonian is calculated for each intermediate system and by numerically integrating, the free energy associated with the transformation is calculated. The free energies of these transformations are then used to calculate the free energy of the desired molecular process [1–5]. In non-equilibrium methods the work in the process of switching between the two Hamiltonians is measured [6, 7]. For the application of chemical reactions the free energy difference between the reactants and products needs to be calculated. Free energies of molecules can be calculated by quantum computations [8] or by normal mode classical calculations [9] (see Entropies section). While quantum mechanics (QM) calculations are highly accurate they are computationally very demanding for large molecules. Normal mode calculations can be very efficient but are based on harmonic dynamics [10].

Molecular modeling includes covalent bond, bond angle, dihedral angle, improper dihedral angle, electrostatic and van der Waals (vdW) potentials (see [11–13] and appendix A). Covalent bond, bond angle and dihedral angle terms depend on the coordinates of two, three and four nearest covalently linked atoms respectively. Electrostatic and vdW potentials are between every atom pair in the system.

In previous studies, the free energy associated with the covalent bond and bond angle terms in molecules and dihedral angle term in the context of restraints were directly calculated by assuming that they are quadratic and by using the rigid rotor approximation and the Herschbach, Johnston and Rapp technique [14–16]. However, the free energies associated with dihedral potential terms in the realistic non-quadratic potentials and with non-quadratic covalent bond and bond angle terms [17] have not been calculated. In addition, it has not been suggested to calculate molecular free energies of submolecules with coupled bond angle degrees of freedom such as the methyl group.

Here we present a direct free energy calculation which is exact and does not require the potential terms to be quadratic. We show that the partition function associated with a dihedral angle in molecules can be decoupled from the partition function of the rest of the system independently of the potential function (even though the dihedral potential term depends on coordinates of atoms in the system) and calculate analytically free energies of realistic dihedral terms. Moreover, we show that free energy associated with submolecules with coupled bond angle degrees of freedom such as the methyl group can be calculated accurately.

In section 1 we present the first step of the calculation in which we transform the molecule by relaxing all the interactions of a group of atoms in that molecule, except ones which depend on a relative spherical coordinate and calculate the free energy difference associated with the transformation. The relaxed potentials include vdW, electrostatic and all other types of potentials except bond stretching, bond angle and dihedral angle potentials. In section 2 we present the second step of the calculation in which the free energy associated with the remaining potentials (which depend on a relative spherical coordinate-bond stretching, bond angle and dihedral angle potentials) of this group of atoms is calculated analytically or numerically. This two-step calculation is applicable to any environment using the definition above, wherein for the case of a solvent in the first step we also relax the interactions of these atoms with the solvent molecules. In section 3 we demonstrate the method in free energy calculations for methanethiol and butane molecules. In section 4 we suggest the potential application of free energy calculation of chemical reactions. In section 5 we discuss the method.

Alchemical free energy transformation (as in section 1) is a standard procedure in molecular free energy calculations [12, 13, 18–20] and is usually used in the context of solvation and binding or in the context of chemical reactions in combination with QM calculations [8]. Here we apply it in the context of free energy calculations of molecules or submolecules in a general environment, where the second calculation is used to obtain the free energy difference between two molecules with submolecule in common. We combine the transformation with novel direct free energy calculations explained in section 2, followed by the novel application of free energy of chemical reactions presented in section 3. Here we do not need to assume harmonic dynamics.

We now explain the first step of the calculation in which we transform the molecule by relaxing potential terms that depend on the coordinates of a group of atoms in that molecule and calculate the free energy difference associated with the transformation using thermodynamic integration (TI) [3, 21, 22]. Alternatively, it can be performed using methods such as exponential averaging/ free energy perturbation (FEP) [23], Bennett acceptance ratio (BAR) [24] and weighted Histogram analysis method [25]. We denote the Hamiltonian with the terms that are removed in the transformation by H_r and the Hamiltonian with the other terms by H_c. We define the Hamiltonian H_r, following the definition above and the potentials in appendix A as the improper dihedral, electrostatic and vdW terms which depend on the coordinates of the group of selected atoms. The λ dependent Hamiltonian that transforms between the systems can be written as follows:

$$\begin{eqnarray*}&&H(\lambda )=\lambda {H}_{{\rm{r}}}+{H}_{{\rm{c}}},\end{eqnarray*}$$

where $\lambda \in [0,1]$. The free energy difference between the original and transformed systems can be written as:

$$\begin{eqnarray*}&&{\rm{\Delta }}{F}_{A\to {A}^{\prime }}=-{k}_{B}T{\displaystyle \int }_{0}^{1}\displaystyle \frac{{\rm{d}}\mathrm{ln}{Z}_{A}(\beta ,\lambda )}{{\rm{d}}\lambda }{\rm{d}}\lambda ={\displaystyle \int }_{0}^{1}\langle {H}_{{\rm{r}}}\rangle {\rm{d}}\lambda ,\end{eqnarray*}$$

where ${A}^{\prime }$ is the transformed system in which $\lambda =0$. Simulating intermediate systems at λ values in the range $[0,1]$ and calculating $\langle {H}_{{\rm{r}}}\rangle$, followed by numerical integration, will enable us to calculate the free energy difference. It is noted that this calculation is valid for any environment (solvent, etc.) since the interactions of the group of atoms with the environment are also relaxed in the transformation. In addition, soft core potentials are used in the vdW and electrostatic interactions to avoid divergences of $\langle {H}_{{\rm{r}}}\rangle$ at $\lambda \to 0$ [26, 27].

As an example we present in figure 1 three systems in the transformation of the molecule phenol with the selected atoms O and H1. The Hamiltonian H_r includes the vdW and electrostatic interactions of the atoms O and H1 with all the atoms in the system, and the improper dihedral term associated with ${\phi }_{{\rm{C}}1-{\rm{C}}2-{\rm{O}}-{\rm{C}}3}$.

Zoom In Zoom Out Reset image size

We now turn to explain how the free energy associated with the selected atoms can be calculated. We first switch to relative coordinates and then to spherical coordinates.

$$\begin{eqnarray*}{\rm{d}}{\boldsymbol{\Omega }} & = & \displaystyle \prod _{i=1}^{n}{\rm{d}}{{\bf{r}}}_{i}^{\prime }={\rm{d}}{{\bf{r}}}_{1}^{\prime }\displaystyle \prod _{i=2}^{k}{\rm{d}}{{\bf{r}}}_{i}\displaystyle \prod _{i=k+1}^{n}{\rm{d}}{{\bf{r}}}_{i}\\ & = & {\rm{d}}{{\bf{r}}}_{1}^{\prime }\displaystyle \prod _{i=2}^{k}{\rm{d}}{{\bf{r}}}_{i}\displaystyle \prod _{i=k+1}^{n}{r}_{i}^{2}\mathrm{sin}{\theta }_{i}{\rm{d}}{r}_{i}{\rm{d}}{\theta }_{i}{\rm{d}}{\phi }_{i},\end{eqnarray*}$$

where ${{\bf{r}}}_{i}\equiv {{\bf{r}}}_{i}^{\prime }-{{\bf{r}}}_{i-1}^{\prime }$, which will be chosen as the position of atoms relative to a covalently bound atom. $k+1$ represents the first atom in the group of selected atoms (e.g the atom O in the example of phenol).

The transformed molecule A' is first divided into elements of covalent bonds, bond junctions, dihedral terms and of more complex structures that include molecular rings. The dihedral potential term depends on the angle between two planes which is equal to the ϕ angle in spherical coordinates. This correspondence can be understood by recalling the definition of the angle between two planes as the angle between the vectors in these planes that are perpendicular to the intersection line of these planes. Hence, the covalent bond, bond angle and dihedral angle terms can be expressed in terms of the spherical variables r, θ and ϕ defined with respect to the relevant atoms.

For the example of phenol we can write the partition function of transformed phenol as follows

$$\begin{eqnarray*}&&{Z}_{{A}^{\prime }}=\int {{\rm{e}}}^{-\beta H({\boldsymbol{\Omega }})}{{\rm{e}}}^{-\beta {H}_{{r}_{{\rm{O}}}}({r}_{{\rm{O}}})}{{\rm{e}}}^{-\beta {H}_{{r}_{{\rm{H}}}}({r}_{{\rm{H}}})}{{\rm{e}}}^{-\beta {H}_{{\theta }_{{\rm{H}}}}({\theta }_{{\rm{H}}})}\mathrm{sin}{\theta }_{{\rm{H}}}{{\rm{e}}}^{-\beta {H}_{{\phi }_{{\rm{H}}}}({\phi }_{{\rm{H}}})}{\rm{d}}{\boldsymbol{\Omega }}{r}_{{\rm{O}}}^{2}{\rm{d}}{r}_{{\rm{O}}}{r}_{{\rm{H}}}^{2}{\rm{d}}{r}_{{\rm{H}}}{\rm{d}}{\theta }_{{\rm{H}}}{\rm{d}}{\phi }_{{\rm{H}}},\end{eqnarray*}$$

where the degrees of freedom of atom H and O are denoted by ${r}_{{\rm{H}}},{\theta }_{{\rm{H}}},{\phi }_{{\rm{H}}}$ and ${r}_{{\rm{O}}},{\theta }_{{\rm{O}}},{\phi }_{{\rm{O}}}$ respectively. All the degrees of freedom of all the other atoms in the system (molecule + solvent molecules in the case of explicit solvent) and ${\theta }_{{\rm{O}}},{\phi }_{{\rm{O}}}$ are denoted by ${\boldsymbol{\Omega }}$. We define ${\theta }_{{\rm{H}}}$ and ${\phi }_{{\rm{H}}}$ as follows: ${\theta }_{{\rm{H}}}\equiv {\theta }_{{\rm{H}}1-{\rm{O}}-{\rm{C}}2},{\phi }_{{\rm{H}}}\equiv {\phi }_{{\rm{H}}1-{\rm{O}}-{\rm{C}}2-{\rm{C}}1}$. Our goal in this example is to calculate the free energy associated with all the degrees of freedom except for the ones denoted by ${\boldsymbol{\Omega }}$. The free energy of the subsystem with the degrees of freedom denoted by ${\boldsymbol{\Omega }}$ can for example then cancel out with an identical subsystem when calculating the free energy difference between two molecules with a submolecule in common. Note that for example ${\phi }_{{\rm{H}}}$ does depend on the coordinates of the C1 atom which is included in ${\boldsymbol{\Omega }}$ and therefore the decomposition of the partition functions is not trivial.

As will become clear, the integration in each element is independent of the others. Thus the integrals can be performed separately and then multiplied to yield the partition function and hence the free energy difference.

We write the free energy factors explicitly:

2.1. Covalent bond

The partition function of a covalent bond represented by a quadratic potential can be written as follows (for non-quadratic terms integrated analytically see appendix E)

$$\begin{eqnarray}&&{Z}_{{\rm{c}}}=\displaystyle \frac{1}{{l}^{3}}\int {{\rm{e}}}^{-\frac{\beta {k}_{{\rm{c}}}{(r-d)}^{2}}{2}}{r}^{2}{\rm{d}}r=\displaystyle \frac{{\pi }^{2}[(2{d}^{2}\beta {k}_{{\rm{c}}}+1)(\mathrm{erf}(d\sqrt{\beta {k}_{{\rm{c}}}})+1)+2\sqrt{\pi }{{\rm{e}}}^{-\beta {k}_{{\rm{c}}}{d}^{2}}d\sqrt{\beta {k}_{{\rm{c}}}}]}{{l}^{3}{(\beta {k}_{{\rm{c}}})}^{3/2}},\end{eqnarray} \tag{ 1 }$$

where l is an arbitrary length (l³ cancels out in comparisons). To account for the bond-dissociation energy/bond energy one has to differentiate between the potential at r = 0 and ${k}_{{\rm{c}}}\to 0$ (free atoms) either by introducing a factor [28–30] or by using potentials such as the Morse potential (numerical integration) [31].

2.2. Two bonds junctions

When considering the case of a linear/bent molecular shapes, it can be seen that when varying the bond angle, the rest of the molecule moves as a rigid body. Since the spherical coordinates representation includes three independent variables, varying a bond angle is decoupled from all the other degrees of freedom of the molecule. Hence the calculation of free energy factor associated with the bond angle potential term is straightforward. For a quadratic term we can write :

$$\begin{eqnarray}{Z}_{{\rm{b}}} & = & \displaystyle \int {{\rm{e}}}^{-\beta {V}_{{\rm{b}}}}\mathrm{sin}\theta {\rm{d}}\theta\\ &&=\displaystyle \int {{\rm{e}}}^{-\frac{\beta {k}_{\theta }}{2}{(\theta -{\theta }_{0})}^{2}}\mathrm{sin}\theta {\rm{d}}\theta =\displaystyle \frac{1}{2\sqrt{\beta {k}_{\theta }}}{{\rm{e}}}^{-{\rm{i}}{\theta }_{0}-\frac{1}{2\beta {k}_{\theta }}}\times \\ & & \sqrt{\displaystyle \frac{\pi }{2}}\left\{{\rm{i}}\;\mathrm{erf}\;\left(\displaystyle \frac{{\rm{i}}-{\theta }_{0}\beta {k}_{\theta }+\beta {k}_{\theta }\pi }{\sqrt{2\beta {k}_{\theta }}}\right)+\mathrm{erfi}\;\left(\displaystyle \frac{1+{\rm{i}}{\theta }_{0}\beta {k}_{\theta }}{\sqrt{2\beta {k}_{\theta }}}\right)\right.\left.-{\mathrm{ie}}^{2{\rm{i}}{\theta }_{0}}\left[\mathrm{erf}\left(\displaystyle \frac{{\rm{i}}+{\theta }_{0}{k}_{\theta }}{\sqrt{2\beta {k}_{\theta }}}\right)-{\rm{i}}\;\mathrm{erfi}\left(\displaystyle \frac{1+{\rm{i}}\beta {k}_{\theta }(\pi -{\theta }_{0})}{\sqrt{2\beta {k}_{\theta }}}\right)\right]\right\}.\end{eqnarray} \tag{ 2 }$$

This expression is real for positive and real values of ${k}_{\theta },\beta$ and ${\theta }_{0}$.

When varying a dihedral angle, the potential term value depends on the orientation of first bond (which determines the axis from which the dihedral angle is measured). However, since the integration has to be performed over all the range $[0,2\pi ]$, integrating with respect to the ϕ angle will yield a factor which is independent of the location of the first bond. Thus, the integration does not depend on the direction of the first bond and is straightforward. For the commonly used dihedral angle potential we get

$$\begin{eqnarray}&&{Z}_{{\rm{d}}}=\int {{\rm{e}}}^{-\beta {k}_{\phi }(1+\mathrm{cos}(n\phi -{\phi }_{s}))}{\rm{d}}\phi =2\pi {{\rm{e}}}^{-\beta {k}_{\phi }}{I}_{0}(\beta {k}_{\phi }),\end{eqnarray} \tag{ 3 }$$

where I₀ is the modified Bessel function of the first kind. Note that this result does not depend on n, which is an integer.

It can now be seen that the partition function of transformed phenol can be decomposed (exactly) into the partition function associated with the ${\boldsymbol{\Omega }}$ degrees of freedom and the covalent bond, bond angle and dihedral angle partition functions.

2.3. Three or more bonds junctions

Molecule shapes can include a monomer (covalent bond) that splits into more than one monomer. Such examples are the trigonal planar, tetrahedral trigonal pyramidal etc. These cases will necessitate numerical integration which can be performed using the Spherical law of cosines:

$$\begin{eqnarray}&&\mathrm{cos}({\theta }_{12})=\mathrm{cos}({\theta }_{1})\mathrm{cos}({\theta }_{2})+\mathrm{sin}({\theta }_{1})\mathrm{sin}({\theta }_{2})\mathrm{cos}({\rm{\Delta }}\phi ),\end{eqnarray} \tag{ 4 }$$

where ${\theta }_{1},{\theta }_{2}$ denote the bond angles of two bonds with respect to the principal bond and ${\theta }_{12},{\rm{\Delta }}\phi$ denote the bond angle and the difference in ϕ angle between these two bonds respectively. Usually in these cases there is one dihedral term, which we denote by ${\phi }_{1}$, that depends on the angle defined by one of the bonds, the principal bond and a previous bond. Since the integration over the other degrees of freedom yields a factor that is independent of the value of ${\phi }_{1}$, the integrations are decoupled. Thus, the integration for the case of one monomer that splits into two can be written as follows:

$$\begin{eqnarray}&&Z=\int {{\rm{e}}}^{-\beta ({V}_{{\rm{b1}}}+{V}_{{\rm{b2}}}+{V}_{{\rm{b12}}}+{V}_{{\rm{d}}})}\mathrm{sin}{\theta }_{1}\mathrm{sin}{\theta }_{2}{\rm{d}}{\theta }_{1}{\rm{d}}{\theta }_{2}{\rm{d}}{\phi }_{1}{\rm{d}}{\phi }_{2}={Z}_{{\rm{d}}}{Z}_{3{\rm{b}}},\end{eqnarray} \tag{ 5 }$$

where Z_d is according to the definition in equation (3) and

$$\begin{eqnarray}&&{Z}_{3{\rm{b}}}=\displaystyle \int {{\rm{e}}}^{-\frac{\beta }{2}[{k}_{1}^{\theta }{({\theta }_{1}-{\theta }_{1}^{0})}^{2}+{k}_{2}^{\theta }{({\theta }_{2}-{\theta }_{2}^{0})}^{2}+{k}_{12}^{\theta }{({\theta }_{12}-{\theta }_{12}^{0})}^{2}]}\mathrm{sin}({\theta }_{1})\mathrm{sin}({\theta }_{2}){\rm{d}}{\theta }_{1}{\rm{d}}{\theta }_{2}{\rm{d}}{\phi }_{2}.\end{eqnarray} \tag{ 6 }$$

For the general case it can be written as follows:

$$\begin{eqnarray}&&{Z}_{\mathrm{nb}}=\displaystyle \int \displaystyle \prod _{i}^{n}{{\rm{e}}}^{-\frac{\beta }{2}{k}_{i}^{\theta }{({\theta }_{i}-{\theta }_{i}^{0})}^{2}}\displaystyle \prod _{i\gt j}{{\rm{e}}}^{-\frac{\beta }{2}{k}_{{ij}}^{\theta }{({\theta }_{{ij}}-{\theta }_{{ij}}^{0})}^{2}}\displaystyle \prod _{i}^{n}\mathrm{sin}{\theta }_{i}{\rm{d}}{\theta }_{i}\displaystyle \prod _{i\geqslant 2}{\rm{d}}{\phi }_{i},\end{eqnarray} \tag{ 7 }$$

where ${\theta }_{{ij}}$ can be calculated from equation (4) and ${\phi }_{1}$, which has to be substituted in ${\rm{\Delta }}\phi$ in this equation, is an arbitrary constant. In case there are energy terms that introduce complexity they can be relaxed in the transformation. For an explicit example with a bond junction see appendix B.

The free energy factors can be substituted in:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}{F}_{{A}^{\prime }\to {A}^{\prime\prime }}}{{k}_{B}T}=\left(\displaystyle \sum _{i}\mathrm{ln}{Z}_{{A}_{{c}_{i}}}+\displaystyle \sum _{i}\mathrm{ln}{Z}_{{A}_{{b}_{i}}}+\displaystyle \sum _{i}\mathrm{ln}{Z}_{{A}_{{d}_{i}}}\right),\end{eqnarray} \tag{ 8 }$$

to give the free energy difference, where ${A}^{\prime\prime }$ denotes a submolecule of the transformed molecule that is not necessarily realistic whose free energy will not be calculated.

In addition, the partition function of complex structures may be calculated numerically. In many cases the complex structures need to be compared to identical ones, eliminating the need for these calculations. Thus, we can write in terms of the partition functions:

$$\begin{eqnarray*}&&Z\to {Z}_{\mathrm{int}}\displaystyle \prod _{i=1}^{l}{Z}_{{{\rm{c}}}_{i}}\displaystyle \prod _{i=1}^{m}{Z}_{{{\rm{d}}}_{i}}\displaystyle \prod _{i=1}^{p}{Z}_{2{{\rm{b}}}_{i}}\displaystyle \prod _{i=1}^{q}{Z}_{3{{\rm{b}}}_{i}}\displaystyle \prod _{i=1}^{r}{Z}_{{\mathrm{complex}}_{i}},\end{eqnarray*}$$

where ${Z}_{\mathrm{int}}$ represents the partition function of the submolecule that is fully interacting and ${Z}_{{{\rm{c}}}_{i}}$ and ${Z}_{{{\rm{d}}}_{i}}$ represent the ith covalent bond and dihedral angle partition function respectively. ${Z}_{2{{\rm{b}}}_{i}}$ and ${Z}_{3{{\rm{b}}}_{i}}$ represent the ith two bond and three or more bond junctions respectively and ${Z}_{{\mathrm{complex}}_{i}}$ represents the ith complex structure partition function. The arrow represents the transformation $\lambda =1\to 0$.

We can also choose the group of selected atoms as all the atoms of the molecule, and calculate the free energy in a similar manner. The partition functions in this case can be written as follows:

$$\begin{eqnarray*}&&Z\to {Z}_{\mathrm{fp}}\displaystyle \prod _{i=1}^{l}{Z}_{{{\rm{c}}}_{i}}\displaystyle \prod _{i=1}^{m}{Z}_{{{\rm{d}}}_{i}}\displaystyle \prod _{i=1}^{p}{Z}_{2{{\rm{b}}}_{i}}\displaystyle \prod _{i=1}^{q}{Z}_{3{{\rm{b}}}_{i}}\displaystyle \prod _{i=1}^{r}{Z}_{{\mathrm{complex}}_{i}},\end{eqnarray*}$$

where ${Z}_{\mathrm{fp}}$ denotes the partition function of a free particle.

In figure 2 we present the free energy decomposition of transformed phenol (A'). The degrees of freedom associated with each system are written on top. F denotes the free energy and A'' which is the reference (possibly imaginary) submolecule, is presented in the first term in the decomposition. It should be noted that the ${\theta }_{{\rm{O}}}$ and ${\phi }_{{\rm{O}}}$ degrees of freedom are associated with the submolecule A'' and the ${r}_{{\rm{O}}}$ degree of freedom is associated with the second system in the decomposition. The potential terms that depend on the coordinates of the O and H atoms in the five decomposed systems are: $1\;:{V}_{{\rm{b}}}({\theta }_{{\rm{C}}1-{\rm{C}}2-{\rm{O}}}),{V}_{{\rm{b}}}({\theta }_{{\rm{C}}3-{\rm{C}}2-{\rm{O}}}),{V}_{{\rm{d}}}({\phi }_{{\rm{C}}4-{\rm{C}}1-{\rm{C}}2-{\rm{O}}})$ $2\;:{V}_{{\rm{c}}}({r}_{{\rm{O}}})$, $3\;:{V}_{{\rm{c}}}({r}_{{\rm{H}}})$, $4\;:{V}_{{\rm{b}}}({\theta }_{{\rm{C}}2-{\rm{O}}-{\rm{H}}1})\;\mathrm{and}\;$ $5\;:{V}_{{\rm{d}}}({\phi }_{{\rm{C}}1-{\rm{C}}2-{\rm{O}}-{\rm{H}}1})$. The reason that the degrees of freedom ${\theta }_{{\rm{O}}}$ and ${\phi }_{{\rm{O}}}$ are associated with the first system in the decomposition is that the corresponding potential terms usually do not depend on the atom (e.g O). The free energies of the second and third systems are calculated according to equation (1) and the free energies of the fourth and fifth systems are calculated according to equations (2) and (3) respectively.

Zoom In Zoom Out Reset image size

The covalent bond free energy calculation has been demonstrated and compared with numerical integration for two molecules of two atoms in a spherical potential (see appendix C). In addition, the free energies for methanethiol (${\mathrm{CH}}_{3}\mathrm{SH}$) and butane (${{\rm{C}}}_{4}{{\rm{H}}}_{10}$) molecules were calculated in vacuum and water environments. In the first step we removed the vdW and Coulomb interactions in a transformation and calculated the free energy difference. We then performed the second step of the calculation using the method and compared to the results obtained by MD simulations. It should be noted that ${{\rm{C}}}_{4}{{\rm{H}}}_{10}$ has a very low solubility in water and our goal was to calculate free energies when the molecule is in an environment (other than vacuum). In the second step of the calculation Helmholtz free energy (denoted by F) is calculated. Since we expect ${\rm{\Delta }}{F}_{\mathrm{method}}={\rm{\Delta }}{F}_{\mathrm{MD}}$ and usually ${\rm{\Delta }}F\backsimeq {\rm{\Delta }}G$ [32] (G is Gibbs free energy) we compared ${\rm{\Delta }}{F}_{\mathrm{method}}$ to ${\rm{\Delta }}{G}_{\mathrm{MD}}$ obtained in constant pressure and temperature (NPT) MD simulations. In case that the difference between the results was larger than 0.15 kcal mol⁻¹ we also compared the results to ${\rm{\Delta }}{F}_{\mathrm{MD}}$ obtained by performing constant temperature and volume (NVT) MD simulations to show that it originates from a difference between F and G. This occurred only once (0.485 kcal mol⁻¹) and when we compared ${\rm{\Delta }}{F}_{\mathrm{method}}$ to ${\rm{\Delta }}{F}_{\mathrm{MD}}$ the free energy difference was 0.038 kcal mol⁻¹.

3.1. Free energy calculations for methanethiol (${\mathrm{CH}}_{3}\mathrm{SH}$)

In the calculations in vacuum due to the fact that the only existing potentials are bond stretching, bond angle and dihedral angle we performed only the second step of the calculation (in solvent both steps were performed). We compared the results obtained according to the equations in section 2 to these obtained by MD simulations in vacuum and water environments and compared the running times.

We performed the calculations at $T=300\;{\rm{K}}$ and used the GROMOS ATB force field in which each molecule is individually parametrized [17] (see figure 3). The force field parametrization was according to bond stretching and bond angle potentials which are slightly different from the standard ones and we performed our calculations accordingly. We give the equations for the potentials here for completeness:

$$\begin{eqnarray*}&&{V}_{{\rm{c}}}=\frac{1}{2}{k}_{{\rm{c}}}{({r}^{2}-{r}_{0}^{2})}^{2},{V}_{{\rm{b}}}=\frac{1}{2}{k}_{\theta }{(\mathrm{cos}\theta -\mathrm{cos}{\theta }_{0})}^{2},{V}_{{\rm{d}}}({\phi }_{{ijkl}})={k}_{\phi }(1+\mathrm{cos}(n\phi -{\phi }_{s})).\end{eqnarray*}$$

Zoom In Zoom Out Reset image size

Our goal is to demonstrate the free energy calculations associated with the potentials which depend on the relative spherical coordinates presented in section 2. This calculation is for atoms which do not interact via vdW and electrostatic interactions, after that these interactions have been removed in a transformation. In vacuum state since only atoms that are distanced four covalent bonds or more interact via vdW and electrostatic interactions, for methanethiol these interactions do not exist. Thus the free energy of the molecule can be decomposed into the free energies associated with the bond stretching terms, the bond angle term, the dihedral angle term and the CH₃ bond junction. When the molecule is in a solvent we first removed the vdW and Coulomb interactions of the atoms involved in a transformation (first step of the calculation) and then calculated the free energies associated with the bond and dihedral angle terms (second step of the calculation). We compared the results and the running times in the second step of the calculation to the ones obtained by performing transformations in which the corresponding potentials were relaxed (MD simulations of the intermediate systems). In these transformations we also treated the molecule as non-interacting systems and in that sense they are also novel at least for the bond junctions.

3.1.1. Results

We calculated the free energies associated with removing the the vdW and Coulomb interactions of atom 1 and atoms 4–6 in water environment and obtained

$$\begin{eqnarray*}&&{\rm{\Delta }}{G}_{\mathrm{vdW}-\mathrm{Coulomb}1\mathrm{MD}\mathrm{water}}=-7.266\pm 0.019\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\\ &&{\rm{\Delta }}{G}_{\mathrm{vdW}-\mathrm{Coulomb}4-6\mathrm{MD}\mathrm{water}}=-14.91\pm 0.019\;\mathrm{kcal}\;{\mathrm{mol}}^{-1}.\end{eqnarray*}$$

After performing these transformations we could calculate the free energy associated with ${\phi }_{1234},{\theta }_{123}$ and the CH₃ bond junction respectively. There was no need to perform these transformations in vacuum $({\rm{\Delta }}G=0).$

The free energy contributions in the second step were calculated by substituting the ATB bond stretching and bond angle potentials in the corresponding partition functions presented in section 2 and integrating (see appendix E) and by using the partition functions for the dihedral potential derived in equation (3). The results for the free energies associated with all the components of the molecule defined in section 2 are presented below:

$$\begin{eqnarray*}&&{F}_{{\rm{c}}12}=5.33\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},{F}_{c23}=5\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\\ &&{F}_{{\rm{c}}34}={F}_{c35}={F}_{c36}=5.54\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},{F}_{{\rm{b}}123}=0.97\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\\ &&{F}_{{\rm{d}}1234}=-0.58\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},{F}_{{\mathrm{CH}}_{3}}=5.39\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\end{eqnarray*}$$

where we have used $F=-{k}_{B}T\mathrm{ln}Z$ and ${F}_{{\mathrm{CH}}_{3}}$ is the free energy associated with the CH₃ bond junction. These results are valid both for vacuum and solvent.

The free energy associated with the bond and dihedral angle terms can be calculated using MD simulations by removing terms in a transformation and subtracting the free energy of the corresponding element ${\text{}}{without}$ the potentials (this is simply $F=-{k}_{B}T\mathrm{ln}{\rm{\Omega }}$, where Ω is phase space volume). The free energy associated with the bond stretching terms on the other hand are less practical to calculate by MD simulations since removing these terms will result in atoms which are unbound to the molecule.

To verify the analysis and the calculations we first performed 21 MD simulations of intermediate systems which interpolate between the molecule with the dihedral term and without the dihedral term associated with ${\phi }_{1234}$ both in vacuum/dilute gas and in solvent. The calculated free energy differences associated with these transformations were

$$\begin{eqnarray*}&&{\rm{\Delta }}{G}_{{\rm{d}}\mathrm{MD}\mathrm{vacuum}}=-0.5095\pm 0.0023\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\\ &&{\rm{\Delta }}{G}_{{\rm{d}}\mathrm{MD}\mathrm{water}}=-0.5095\pm 0.0023\;\mathrm{kcal}\;{\mathrm{mol}}^{-1}.\end{eqnarray*}$$

These results can be compared to

$$\begin{eqnarray*}{\rm{\Delta }}{F}_{{\rm{d}}\mathrm{equation}(3)} & = & -{k}_{B}T\mathrm{ln}{{\rm{\Omega }}}_{{\rm{d}}}-{F}_{{\rm{d}}}=-{k}_{B}T\mathrm{ln}2\pi +{k}_{B}T\mathrm{ln}2\pi {{\rm{e}}}^{-\beta {k}_{\phi }}{I}_{0}(\beta {k}_{\phi })\\ & = & {k}_{B}T\mathrm{ln}[{{\rm{e}}}^{-\beta {k}_{\phi }}{I}_{0}(\beta {k}_{\phi })]=-0.507982\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\end{eqnarray*}$$

obtained by simply substituting ${k}_{\phi }$ and β in this equation which is based on equation (3). This validates both the free energy calculation and the partition function decomposition in equation (5).

We then performed transformations in which the bond angle term associated with ${\theta }_{123}$ is relaxed. We performed the transformation both using the bond angle potential of the force field and by treating the parameters as if they were given for the standard quadratic bond angle term. The results obtained from the MD simulations in vacuum and water for the quadratic potential are the following:

$$\begin{eqnarray*}&&{\rm{\Delta }}{G}_{{\rm{b}}\mathrm{MD}\mathrm{vacuum}\mathrm{quadratic}}=-1.39\pm 0.0023\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\\ &&{\rm{\Delta }}{G}_{{\rm{b}}\mathrm{MD}\mathrm{water}\mathrm{quadratic}}=-1.39\pm 0.01\;\mathrm{kcal}\;{\mathrm{mol}}^{-1}.\end{eqnarray*}$$

These results are in agreement with the calculation according to equation (2):

$$\begin{eqnarray*}&&{\rm{\Delta }}{F}_{{\rm{b}}\mathrm{Eq}.(2)}=-{k}_{B}T\mathrm{ln}2-{F}_{{\rm{b}}\mathrm{quadratic}}=-1.38152\;\mathrm{kcal}\;{\mathrm{mol}}^{-1}.\end{eqnarray*}$$

The results obtained from MD simulations in vacuum and water for the ATB bond angle potential are the following:

$$\begin{eqnarray*}&&{\rm{\Delta }}{G}_{{\rm{b}}\mathrm{MD}\mathrm{vacuum}\mathrm{cos}}=-1.388\pm 0.0023\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},{\rm{\Delta }}{G}_{{\rm{b}}\mathrm{MD}\mathrm{water}\mathrm{cos}}=-1.421\pm 0.024\;\mathrm{kcal}\;{\mathrm{mol}}^{-1}.\end{eqnarray*}$$

These results are in agreement with the calculation according to the same potential:

$$\begin{eqnarray*}&&{\rm{\Delta }}{F}_{{\rm{b}}\mathrm{cos}}=-{k}_{B}T\mathrm{ln}2-{F}_{{\rm{b}}\mathrm{cos}}=-1.38152\;\mathrm{kcal}\;{\mathrm{mol}}^{-1}.\end{eqnarray*}$$

We then performed the transformation of the CH₃ bond junction in which the bond angle terms associated with the angles ${\theta }_{234},{\theta }_{235},{\theta }_{236},{\theta }_{435},{\theta }_{436}$ and ${\theta }_{536}$ were relaxed. The results obtained from MD simulations for the quadratic potential are the following:

$$\begin{eqnarray*}&&{\rm{\Delta }}{G}_{{\mathrm{CH}}_{3}\mathrm{MD}\mathrm{vacuum}\mathrm{quadratic}}=-9.019\pm 0.007\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\\ &&{\rm{\Delta }}{G}_{{\mathrm{CH}}_{3}\mathrm{MD}\mathrm{water}\mathrm{quadratic}}=-9.038\pm 0.012\;\mathrm{kcal}\;{\mathrm{mol}}^{-1}.\end{eqnarray*}$$

These results are in agreement with result obtained by numerically integrating according to equation (7)

$$\begin{eqnarray*}&&{\rm{\Delta }}{F}_{{\mathrm{CH}}_{3}\mathrm{equation}(7)}=-{k}_{B}T[\mathrm{ln}(32{\pi }^{2})-\mathrm{ln}({Z}_{{\mathrm{CH}}_{3}\mathrm{quadratic}})]=-9.02383\pm 0.0021\;\mathrm{kcal}\;{\mathrm{mol}}^{-1}.\end{eqnarray*}$$

We also calculated the free energy difference associated with the CH₃ transformation using the ATB force field bond angle potential. The results obtained from MD simulations are the following:

$$\begin{eqnarray*}&&{\rm{\Delta }}{G}_{{\mathrm{CH}}_{3}\mathrm{MD}\mathrm{vacuum}\mathrm{cos}}=-8.83\pm 0.0047\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\\ &&{\rm{\Delta }}{G}_{{\mathrm{CH}}_{3}\mathrm{MD}\mathrm{water}\mathrm{cos}}=-8.82\pm 0.014\;\mathrm{kcal}\;{\mathrm{mol}}^{-1}.\end{eqnarray*}$$

These results are in agreement with the result obtained by the numerical integration

$$\begin{eqnarray*}&&{\rm{\Delta }}{F}_{{\mathrm{CH}}_{3}\mathrm{cos}}=-8.804\pm 0.00365\;\mathrm{kcal}\;{\mathrm{mol}}^{-1}.\end{eqnarray*}$$

The running times of the calculations in the first step were 0 s for vacuum (no transformation) and $1.27\times 21=26.7\;{\rm{h}}$ and $1.38\times 21=29\;{\rm{h}}$ for removing vdW and Coulomb terms of atom 1 and atoms 4–6, respectively, in water environment (Six-Core AMD Opteron(tm) processor 2427, 2.19 Ghz).

The calculations of the free energies associated with the bond stretching, bond angle and dihedral angle terms according to section 2 were immediate ∼10⁻⁸–10⁻⁷ with eight digits of precision as they are merely substitutions in equations (1)–(3) (for the ATB force field potentials equations (1) and (2) were slightly different). MD simulations in a total running time of $3.5\times 21=73.5$ min and $1\times 21=21$ h (Six-Core AMD Opteron(tm) processor 2427, 2.19 Ghz) for vacuum and water environments respectively, yielded a precision of 2–3 digits, that may be lower in practice due to the spacing between intermediates.

The running time of the numerical integration in the free energy calculation of the CH₃ bond junction according to equation (7) was $2.34\;\mathrm{sec}$ on a 2 GHz dual core intel processor. The total running time of the MD simulations was $3.5\times 21=73.5$ min and $9\times 21=189$ h (Six-Core AMD Opteron(tm) processor 2427, 2.19 Ghz) for vacuum and water environments respectively.

In conclusion, the calculations according to equations (2)–(7) were in agreement with the results obtained by performing transformations using MD simulations. The computations using these equations (second step of the calculation) were much faster for the same computational resources. This originates from the fact that using the method we perform calculations on subsystems with a small number of degrees of freedom, compared with the MD calculations which include all the degrees of freedom.

When using the method the dominant running times are of the first step of the calculation, which are reported above. Note that calculating free energies of a covalent bond can be performed by a direct calculation (see for example equations (1) and (E1)) and removing a covalent bond term in MD transformation will involve breaking the covalent bond and difficult sampling. The same applies to calculating free energies of molecular rings, which may be performed only by a direct calculation. Thus, the complete free energy calculation is difficult or impractical to perform using only MD transformations. In addition, the running times of the method were compared in some cases (e.g. bond junctions) to the running times of novel MD transformations.

3.2. Free energy calculations for butane (${{\rm{C}}}_{4}{{\rm{H}}}_{10})$

We performed free energy calculations for butane in vacuum and water at T = 300 K using the method and by performing MD simulations. In figure 4 we present butane molecule with the ATB bond and dihedral angle potential (force field) constants [17]. The dihedral terms in this molecule can have high energy barriers, which lead to difficult sampling. This can be overcome by using techniques such as parallel tempering/ replica exchange [33–35], umbrella sampling [25] or transforming variables without introducing a bias [36, 37]. In our simulations and for the purpose of our calculations we could perform the sampling without using these techniques.

Zoom In Zoom Out Reset image size

We first removed the potentials which do not depend on the spherical coordinates of atoms 1–4 in vacuum and in water environments. The Gibbs free energies associated with these transformations were the following:

$$\begin{eqnarray*}&&{\rm{\Delta }}{G}_{\mathrm{vdW}-\mathrm{Coulomb}\mathrm{MD}\mathrm{vacuum}}=-3.368\pm 0.0024\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\\ &&{\rm{\Delta }}{G}_{\mathrm{vdW}-\mathrm{Coulomb}\mathrm{MD}\mathrm{water}}=-3.348\pm 0.0024\;\mathrm{kcal}\;{\mathrm{mol}}^{-1}.\end{eqnarray*}$$

We then calculated the free energies associated with the four bonds junction which involves atoms 1–5 and with the dihedral angle ${\phi }_{3-2-5-8}$. These free energies were calculated using the method and by transforming the molecule and performing MD simulations in vacuum and water environments.

We obtained the following results

$$\begin{eqnarray*}&&{\rm{\Delta }}{G}_{{{\rm{C}}}_{2}{{\rm{H}}}_{3}\mathrm{MD}\mathrm{vacuum}}=-8.442\pm 0.0024\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\\ &&{\rm{\Delta }}{G}_{{{\rm{C}}}_{2}{{\rm{H}}}_{3}\mathrm{MD}\mathrm{water}}=-8.535\pm 0.0335\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\\ &&{\rm{\Delta }}{F}_{{{\rm{C}}}_{2}{{\rm{H}}}_{3}\mathrm{method}}=-8.4089\pm 0.0051\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\\ &&{\rm{\Delta }}{G}_{{\rm{d}}3-2-5-8\mathrm{MD}\mathrm{vacuum}}=-0.777\pm 0.00\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\\ &&{\rm{\Delta }}{G}_{{\rm{d}}3-2-5-8\mathrm{MD}\mathrm{water}}=-0.837\pm 0.0024\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\\ &&{\rm{\Delta }}{F}_{{\rm{d}}3-2-5-8\mathrm{method}}=-0.762\;\mathrm{kcal}\;{\mathrm{mol}}^{-1}.\end{eqnarray*}$$

We then removed in a transformation the potentials which do not depend on the spherical coordinates of atoms 5–7 in vacuum and water environments and obtained

$$\begin{eqnarray*}&&{\rm{\Delta }}{G}_{\mathrm{vdW}-\mathrm{Coulomb}\mathrm{MD}\mathrm{vacuum}}=7.311\pm 0.0072\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\\ &&{\rm{\Delta }}{G}_{\mathrm{vdW}-\mathrm{Coulomb}\mathrm{MD}\mathrm{water}}=14.166\pm 0.012\;\mathrm{kcal}\;{\mathrm{mol}}^{-1}.\end{eqnarray*}$$

We calculated the free energies associated with the four bonds junction which involves atoms 2, 5–8 and with the dihedral angle ${\phi }_{2-5-8-11}$. These free energies were calculated using the method and by transforming the molecule and performing MD simulations in vacuum and water. We obtained the following results

$$\begin{eqnarray*}&&{\rm{\Delta }}G{}_{{{\rm{C}}}_{3}{{\rm{H}}}_{2}\mathrm{MD}\mathrm{vacuum}}=-9.582\pm 0.0096\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\\ &&{\rm{\Delta }}G{}_{{{\rm{C}}}_{3}{{\rm{H}}}_{2}\mathrm{MD}\mathrm{water}}=-9.981\pm 0.0143\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\\ &&{\rm{\Delta }}F{}_{{{\rm{C}}}_{3}{{\rm{H}}}_{2}\mathrm{MD}\mathrm{water}}=-9.534\pm 0.0143\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\\ &&{\rm{\Delta }}F{}_{{{\rm{C}}}_{3}{{\rm{H}}}_{2}\mathrm{method}}=-9.4958\pm 0.0064\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\\ &&{\rm{\Delta }}G{}_{{\rm{d}}2-5-8-11\mathrm{MD}\mathrm{vacuum}}=-0.762\pm 0.0024\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\\ &&{\rm{\Delta }}G{}_{{\rm{d}}2-5-8-11\mathrm{MD}\mathrm{water}}=-0.770\pm 0.00\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\\ &&{\rm{\Delta }}F{}_{{\rm{d}}2-5-8-11\mathrm{metod}}=-0.762\;\mathrm{kcal}\;{\mathrm{mol}}^{-1},\end{eqnarray*}$$

where we have calculated ${\rm{\Delta }}F{}_{{{\rm{C}}}_{3}{{\rm{H}}}_{2}\mathrm{MD}\mathrm{water}}$ to show that the difference between the results originates from the difference between F and G.

Here we suggest the potential application of calculating free energies of chemical reactions using classical molecular simulations followed by analytic or numerical calculations. We will calculate the free energies of only small parts of the molecules, possibly in solvent, to obtain the free energy of the chemical reaction which can involve large molecules.

To obtain the equilibrium constant of a chemical reaction the standard Gibbs free energies (which are similar to Helmholtz free energies in many cases [32]) are usually calculated. The standard state is the hypothetical state with the standard state concentration but exhibiting infinite-dilution behavior (the interactions between e.g the solute molecules are negligible). Hence, when we are interested in the standard free energy of a substance a single copy of the molecule of interest can be simulated either in vacuum or solvent environments.

The idea is to transform the reactants and the products, between which the free energy difference is calculated, into molecules that have the same partition functions up to factors that can be calculated. First, we match reactant molecules with product molecules that have submolecules in common if possible so that the free energy of these submolecules will cancel out. Then, we transform the molecules by relaxing potential terms of the atoms that are different between the molecules and calculate the free energy differences associated with each transformation. Finally the free energy factors associated with the different atoms are calculated and we can deduce the free energy of the chemical reaction. For example, given the chemical reaction

$$\begin{eqnarray}&&A+B\to C+D,\end{eqnarray} \tag{ 9 }$$

we can match, if possible, molecules $A,B$ to molecules $C,D$. Then we transform each of the molecules to ${A}^{\prime },{B}^{\prime },{C}^{\prime }$ and ${D}^{\prime }$ and calculate the free energy differences ${\rm{\Delta }}{F}_{A\to {A}^{\prime }},{\rm{\Delta }}{F}_{B\to {B}^{\prime }},{\rm{\Delta }}{F}_{C\to {C}^{\prime }}$ and ${\rm{\Delta }}{F}_{D\to {D}^{\prime }}$ respectively. Finally, we calculate the free energy factors ${\rm{\Delta }}{F}_{{A}^{\prime }\to {A}^{\prime\prime }},{\rm{\Delta }}{F}_{{B}^{\prime }\to {B}^{\prime\prime }},{\rm{\Delta }}{F}_{{C}^{\prime }\to {C}^{\prime\prime }}$ and ${\rm{\Delta }}{F}_{{D}^{\prime }\to {D}^{\prime\prime }}$ as previously explained. Thus, for the standard state free energy we get

$$\begin{eqnarray}{\rm{\Delta }}{F}_{A+B\to C+D} & = & {\rm{\Delta }}{F}_{A\to {A}^{\prime }}+{\rm{\Delta }}{F}_{B\to {B}^{\prime }}-{\rm{\Delta }}{F}_{C\to {C}^{\prime }}-{\rm{\Delta }}{F}_{D\to {D}^{\prime }}\\ & & +{\rm{\Delta }}{F}_{{A}^{\prime }\to {A}^{\prime\prime }}+{\rm{\Delta }}{F}_{{B}^{\prime }\to {B}^{\prime\prime }}-{\rm{\Delta }}{F}_{{C}^{\prime }\to {C}^{\prime\prime }}-{\rm{\Delta }}{F}_{{D}^{\prime }\to {D}^{\prime\prime }}.\end{eqnarray} \tag{ 10 }$$

An example for such a calculation is given in figure 5. The molecules that participate in the chemical reaction and the transformed molecules are presented on top and bottom, respectively. In this case molecule B is matched with C and ${F}_{{B}^{\prime\prime }}={F}_{{C}^{\prime\prime }}$. ${F}_{{A}^{\prime\prime }}$ and ${F}_{{D}^{\prime\prime }}$, which are the free energy of a free particle, are equal.

Zoom In Zoom Out Reset image size

To apply the method, the force field should contain potentials that depend on a relative spherical coordinate (bond stretching, bond angle and dihedral angle), which is the case in most force fields. In addition in order for the free energies associated with bond stretching potentials to fully match experimental results, the force field should account for bond energies. An example for a force field which satsifies both requirements is MM3 [30]. Free energies associated with bond stretching, bond angle, dihedral angle and bond junctions can be readily calculated. It is assumed that free energy associated with an aromatic ring can be calculated. The free energy calculation of more complicated molecular rings needs to be assessed and treating them as a non-interacting system is expected to simplify their free energy calculations (also in other direct calculation methods). Fortunately, molecular rings are in many cases stable so the free energy associated with them often does not need to be calculated. Currently, the software which supports MM3 force field enables to perform alchemical transformations (first step of the calculation) only for single atoms. We intend to extend the code and demonstrate the method in a free energy calculation of a chemical reaction in our next project.

We presented an accurate and complete method for calculating molecular free energies in classical potentials. We suggested the potential application of free energy calculation of chemical reactions in classical potentials. This application can find use in organic chemistry, biochemistry and drug discovery.

D J Bergman, A Nitzan, O Hod, D Harries and G Falkovich are acknowledged for the useful comments. H Suchowski is acknowledged for funding the publication.

Here we briefly describe the molecular potentials.

The covalent bond between two atoms is modeled with the following potential:

$$\begin{eqnarray}&&{V}_{{\rm{c}}}={k}_{{\rm{c}}}{(r-d)}^{2},\end{eqnarray} \tag{ A1 }$$

where r is the distance between the atoms. The bond angle term between three atoms is usually modeled with the following potential:

$$\begin{eqnarray}&&{V}_{{\rm{b}}}(\theta )=\displaystyle \frac{1}{2}{k}_{\theta }{(\theta -{\theta }_{0})}^{2}.\end{eqnarray} \tag{ A2 }$$

The commonly used dihedral angles potential is of the following type:

$$\begin{eqnarray}&&{V}_{{\rm{d}}}({\phi }_{{ijkl}})={k}_{\phi }(1+\mathrm{cos}(n{\phi }_{{ijkl}}-{\phi }_{s})),\end{eqnarray} \tag{ A3 }$$

where ${\phi }_{{ijkl}}$ is defined by the angle between the plane formed by the ijk atoms and the plane formed by the jkl atoms. The improper dihedral term that is used to enforce planarity is defined as follows:

$$\begin{eqnarray}&&V({\phi }_{{ijkm}})=\frac{1}{2}k{({\phi }_{{ijkm}}-{\phi }_{s})}^{2},\end{eqnarray} \tag{ A4 }$$

where ${\phi }_{{ijkm}}$ is defined by the angle between the plane formed by the ijm atoms and the plane formed by the imk atoms. The Coulomb and vdW terms are defined by the $1/r$ and ${r}^{-12}-{r}^{-6}$ potentials respectively. The coefficients for these interactions are determined by the atom charge and vdW coefficients of both atoms.

To illustrate a case with a bond junction we write the partition function of transformed Toluene (see figure 7), where as in the case of phenol, the Hamiltonian H_r includes the improper dihedral, electrostatic and vdW potential terms which depend on the coordinates of the group of selected atoms

$$\begin{eqnarray*}{Z}_{{\mathrm{Toluene}}^{\prime }} & = & \displaystyle \int {{\rm{e}}}^{-\beta H({\boldsymbol{\Omega }})}\displaystyle \prod _{i=1}^{3}{{\rm{e}}}^{-\beta {H}_{{{\rm{r}}}_{{Hi}}}({r}_{{Hi}})}{{\rm{e}}}^{-\beta {H}_{{\theta }_{{Hi}-C2}}({\theta }_{{Hi}})}{{\rm{e}}}^{-\beta {H}_{{\theta }_{H1-H2}}({\theta }_{H1-H2})-\beta {H}_{{\theta }_{H2-H3}}({\theta }_{H2-H3})-\beta {H}_{{\theta }_{H3-H1}}({\theta }_{H3-H1})}\\ & & \times {{\rm{e}}}^{-\beta {H}_{{\phi }_{H1}}({\phi }_{H1})}{\rm{d}}{\boldsymbol{\Omega }}{r}_{{Hi}}^{2}{\rm{d}}{r}_{{Hi}}\mathrm{sin}{\theta }_{{Hi}}{\rm{d}}{\theta }_{{Hi}}{\rm{d}}{\phi }_{Hi}.\end{eqnarray*}$$

It can be readily seen that the partition function can be decomposed into the one associated with the ${\boldsymbol{\Omega }}$ degrees of freedom and the covalent bond, dihedral angle and three bonds junction partition functions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to demonstrate the method it was used with all its ingredients in MC simulation to calculate the free energy difference between the systems A and B and this value was compared to the one calculated by numerical integration. Then, in order to assess the efficiency of the method, the free energy difference between the systems was calculated using TI combined with Hamiltonian parallel tempering (for sampling rugged energy landscape).

The compared systems are composed of a molecule of two atoms in which one atom is fixed to the origin and the second one is bound to the first by a covalent bond. The second atom in each system is in a θ dependent potential (in spherical coordinates), containing ${\theta }^{-12}$ term to represent the vdW repulsive term used in molecular modeling. The potential barrier was chosen to be of typical value of systems with tens of atoms, having rugged energy landscape. The covalent bond length difference was chosen to represent systems with few different atom lengths (the values of the pairs of spring constant and covalent bond length were taken from molecular simulation software). The following potential and parameters were used:

$$\begin{eqnarray*}&&{V}_{A}=\displaystyle \frac{1}{2}{k}_{{cA}}{({r}_{12}-{r}_{{\mathrm{eq}}_{A}})}^{2}-5.5{k}_{B}T\cdot \left[\mathrm{sin}(4\theta )-\displaystyle \frac{{10}^{-8}}{{\left(\theta -\frac{3\pi }{8}\right)}^{12}}\right],\end{eqnarray*}$$

$$\begin{eqnarray*}&&{V}_{B}=\displaystyle \frac{1}{2}{k}_{{cB}}{({r}_{12}-{r}_{{\mathrm{eq}}_{B}})}^{2}+5{k}_{B}T\cdot \left[\mathrm{sin}(4\theta )+\displaystyle \frac{{10}^{-8}}{{\left(\theta -\frac{5\pi }{8}\right)}^{12}}\right],\end{eqnarray*}$$

$$\begin{eqnarray*}&&{r}_{{\mathrm{eq}}_{A}}=2.1\;\mathring{\rm A} ,{r}_{{\mathrm{eq}}_{B}}=1.3\;\mathring{\rm A} ,{k}_{{cA}}=123\;\mathrm{kcal}\;{\mathrm{mol}}^{-1}\;{\mathring{\rm A} }^{-2},\end{eqnarray*}$$

$$\begin{eqnarray*}&&{k}_{{cB}}=428\;\mathrm{kcal}\;{\mathrm{mol}}^{-1}\;{\mathring{\rm A} }^{-2}.\end{eqnarray*}$$

Here we used the partition function of two atoms with a covalent bond term represented by a quadratic potential that can be written as follows:

$$\begin{eqnarray*}&&Z(\beta )=\displaystyle \frac{1}{{l}^{3}}{\displaystyle \int }_{-\infty }^{\infty }{{\rm{e}}}^{-\beta {k}_{{\rm{c}}}{(r-d)}^{2}}{{\rm{d}}}^{3}x=\displaystyle \frac{4\pi }{{l}^{3}}{\displaystyle \int }_{0}^{\infty }{{\rm{e}}}^{-\beta {k}_{{\rm{c}}}{(r-d)}^{2}}{r}^{2}{\rm{d}}r,\end{eqnarray*}$$

$$\begin{eqnarray}&&=\displaystyle \frac{\pi [(2{d}^{2}\beta {k}_{{\rm{c}}}+1)(\mathrm{erf}(d\sqrt{\beta {k}_{{\rm{c}}}})+1)+2\sqrt{\pi }{{\rm{e}}}^{-\beta {k}_{{\rm{c}}}{d}^{2}}d\sqrt{\beta {k}_{{\rm{c}}}}]}{4\;{l}^{3}{(\beta {k}_{{\rm{c}}})}^{3/2}},\end{eqnarray} \tag{ C1 }$$

where l is an arbitrary length. This was used to calculate the free energy difference between the systems with the coupling terms relaxed.

The comparison of the method to the numerical integration yielded exactly the same results.

We have used Gromacs simulation package [38] and Gromos53a7 force field parameters (from ATB server [17]) along with spc water model for the simulations. For all the simulations a time step of 2 fs was used. The cubic or dodecahedron box with a minimum distance of 1 nm between the solute and the box edge was considered during the simulations. After minimization, the systems were equilibrated (only for simulations in water) in NVT and NPT ensembles for 100 ps. All the production simulations were performed under NPT ensemble, except the bond angle transformation for ${{\rm{C}}}_{4}{{\rm{H}}}_{10}$ in water, where NVT ensemble was also used due to the difference between ${\rm{\Delta }}{G}_{\mathrm{MD}}$ and ${\rm{\Delta }}{F}_{\mathrm{method}}$. In both vacuum and water, the systems were simulated at each of the 21 equi-spaced intermediate states including the initial ($\lambda =0$) and final states ($\lambda =1$). The length of simulation was 20 ns in vacuum, 1 ns for ${\mathrm{CH}}_{3}\mathrm{SH}$ in water and 10 ns for ${{\rm{C}}}_{4}{{\rm{H}}}_{10}$ in water. However, in the case of the CH₃ bond angle transformation in water (in ${\mathrm{CH}}_{3}\mathrm{SH}$), the length of the simulation was 20 ns at each intermediate state. We have computed the free energy change using the Bennett's acceptance ratio method in which a series of individual free energies is combined into a free energy estimate [24].

We substituted the (ATB force field) bond stretching and bond angle potentials in the corresponding equations in section 2 and integrated. We obtained the following expressions

$$\begin{eqnarray}{Z}_{{\rm{c}}} & = & \displaystyle \frac{\pi {{\rm{e}}}^{-\frac{1}{4}a}}{8\beta {k}_{{\rm{c}}}{r}_{0}}\left\{a\left[{I}_{\frac{3}{4}}(a/4)+{I}_{\frac{5}{4}}(a/4)\right]+{{aI}}_{-\frac{1}{4}}(a/4)+(a+2){I}_{\frac{1}{4}}(a/4)\right\},\end{eqnarray} \tag{ E1 }$$

where $a\equiv \beta {k}_{{\rm{c}}}{r}_{0}^{4}$ and ${I}_{n}(x)$ is the modified Bessel function of the first kind

$$\begin{eqnarray}&&{Z}_{{\rm{b}}}=\sqrt{\displaystyle \frac{\pi }{2\beta {k}_{\theta }}}\left\{{\rm{erf}}\;\left[\sqrt{2\beta {k}_{\theta }}{\mathrm{sin}}^{2}\left(\displaystyle \frac{{\theta }_{0}}{2}\right)\right]+{\rm{erf}}\;\left[\sqrt{\displaystyle \frac{\beta {k}_{\theta }}{2}}(\mathrm{cos}({\theta }_{0})+1)\right]\right\},\end{eqnarray} \tag{ E2 }$$

$$\begin{eqnarray}&&{Z}_{n{\rm{b}}}=\displaystyle \int \displaystyle \prod _{i}^{n}{{\rm{e}}}^{-\frac{\beta }{2}{k}_{i}^{\theta }{(\mathrm{cos}{\theta }_{i}-\mathrm{cos}{\theta }_{i}^{0})}^{2}}\displaystyle \prod _{i\gt j}{{\rm{e}}}^{-\frac{\beta }{2}{k}_{{ij}}^{\theta }{(\mathrm{cos}{\theta }_{{ij}}-\mathrm{cos}{\theta }_{{ij}}^{0})}^{2}}\displaystyle \prod _{i}^{n}\mathrm{sin}{\theta }_{i}{\rm{d}}{\theta }_{i}\displaystyle \prod _{i\geqslant 2}{\rm{d}}{\phi }_{i}.\end{eqnarray} \tag{ E3 }$$