Calculation of molecular free energies in classical potentials

Free energies of molecules can be calculated by quantum computations or by normal mode classical calculations. However, the first can be computationally impractical for large molecules and the second is based on the assumption of harmonic dynamics. We present a novel, accurate and complete calculation of molecular free energies in standard classical potentials. In this method we transform the molecule by relaxing potential terms which depend on the coordinates of a group of atoms in that molecule and calculate the free energy difference associated with the transformation. Then, since the transformed molecule can be treated as non interacting systems, the free energy associated with these atoms is analytically or numerically calculated. This two-step calculation can be applied to calculate free energies of molecules or free energy difference between (possibly large) molecules in a general environment. We suggest the potential application of free energy calculation of chemical reactions in classical molecular simulations.


INTRODUCTION
Free energy calculations have a variety of applications which include binding, solvation and chemical reactions. In equilibrium methods 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][2][3][4][5][6]. In non equilibrium methods the work in the process of switching between the two Hamiltonians is measured [7,8].
Molecular modeling includes covalent bond, bond angle, dihedral angle, improper dihedral angle, electrostatic and vdW potentials (see [9][10][11] and Supplementary Material). 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 works, the free energy associated with the covalent bond, bond angle and dihedral angle terms were calculated by assuming that they are harmonic and by using the rigid rotor approximation and the HJR technique [12,13]. However, the dihedral potential terms are usually not harmonic and the covalent bond and bond angle terms are not necessarily harmonic. In addition, molecular free energies of submolecules with coupled (bond angle) degrees of freedom such as the methyl group have not been calculated.
Here we present an exact calculation that does not require the potential terms to be harmonic. Moreover, we show that free energy associated with submolecules with coupled degrees of freedom can be calculated accurately. * asaffarhi@post.tau.ac.il

THE TRANSFORMATION
We now explain the first step of the calculation in which we 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) [4,14,15]. 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 as the improper dihedral, electrostatic and vdW terms that depend on the coordinates of the group of selected atoms. The λ dependent Hamiltonian that transforms between the systems can be written as follows: where λ ∈ [0, 1]. The free energy difference between the λ=1 λ=0.5 λ=0  and calculating H r , followed by numerical integration, will enable us to calculate the free energy difference.
As an example we present in Fig. 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 φ C1−C2−O−C3 .

THE REMAINING FREE ENERGY DIFFERENCE
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.
where r i ≡ r i − r i−1 , 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 standard covalent bonds, bond junctions and of more complex structures that include molecular rings. The dihedral potential term depends on the angle between two planes (see SM for details) 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. Since each of these terms depends on one independent variable, 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 these free energy factors explicitly:

Covalent bond
The partition function of a covalent bond in spherical coordinates can be written as follows: where l is an arbitrary length (l 3 cancels out in comparisons). Since at k c → 0 the atoms are free, it is clear that in order to account for the bond-dissociation energy/bond energy one has to differentiate between the potential at r = 0 and k c → 0 either by introducing a factor to this potential (multiplied by k c ) [16,17]/ [18] or by using potentials such as the Morse potential (numerical integration) [19].

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 potential term is straightforward. One of the most common bond angles potentials is the following: The integration over the corresponding degree of freedom can be written as: This expression is real for positive and real values of k θ , β and θ 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π], 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. The commonly used dihedral angle potential is the following: The integration over this degree of freedom yields the following result: where I 0 is the modified Bessel function of the first kind. Note that this result does not depend on n, which is an integer.
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 that can be written as follows: cos (θ 12 ) = cos (θ 1 ) cos (θ 2 ) + sin (θ 1 ) sin (θ 2 ) cos (∆φ) , (4) where θ 1 , θ 2 denote the bond angles of two bonds with respect to the principal bond and θ 12 , ∆φ 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 φ 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 φ 1 , the integrations are decoupled. Thus, the integration for the case of one monomer that splits into two can be written as follows: where Z d is according to the definition in Eq. (3) and For the general case it can be written as follows: where θ ij can be calculated from Eq. (4) and φ 1 , which has to be substituted in ∆φ in this equation, is an arbitrary constant. In case there are energy terms that introduce complexity they can be relaxed in the transformation. These free energy factors can be substituted in: to give the free energy difference, where A 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: where Z int represents the partition function of the submolecule that is fully interacting and Z ci and Z di represent the ith covalent bond and dihedral angle partition function respectively. Z 2bi and Z 3bi represent the ith two bond and three or more bond junctions respectively and Z complex i represents the ith complex structure partition function. The arrow represents the transformation λ = 1 → 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: where Z fp denotes the partition function of a free particle.
In Fig. 2 we present the free energy decomposition of transformed phenol (A'). We denote the degrees of freedom of atom H by r H , θ H , φ H and atom O by All the other degrees of freedom of the sys-tem (molecule+solvent molecules in the case of explicit solvent) are denoted by Ω. 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 θ O and φ O degrees of freedom are associated with the submolecule A" and the r 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 : . The reason that the degrees of freedom θ O and φ 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 Eq. (1) and the free energies of the fourth and fifth systems are calculated according to Eqs. (2) and (3)  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 Supplementary Material). In addition the results obtained by substituting in Eqs. (2)-(7) were in agreement with the results obtained by performing transformations in MD simulations and the computations using the method were faster by factors varying between 5 · 10 3 and 10 12 [20].

POTENTIAL APPLICATION
When the goal is to calculate the free energy difference of a chemical reaction, we can directly transform each of the reactant molecules into a product molecule. However, if the molecules have different number of atoms, a direct transformation is not feasible. Thus, Quantum Mechanical (QM) calculations are usually combined with molecular simulations in free energy calculations of chemical reactions. One way to calculate the free energy difference of a chemical reaction in the general case is to calculate the solvation free energies of the molecules using molecular simulations. Then, the free energy difference between the molecules in vacuum is calculated using QM methods. Thus, using the Thermodynamic Cycle the free energy difference between the molecules in liquid can be deduced [21]. Alternatively, hybrid QM/MM (Molecular mechanics) simulations in which the part of the system which is relevant for the chemical reaction has QM force fields can be performed [22]. Here we suggest that the free energy difference can be calculated by classical molecular simulations followed by analytic or numerical calculations.
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 we can match, if possible, molecules A, B to molecules C, D. Then we transform each of the molecules to A , B , C and D and calculate the free energy differences ∆F A→A , ∆F B→B , ∆F C→C and ∆F D→D respectively. Finally we calculate the free energy factors An example for such a calculation is given in Fig. 3. 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 = F C . F A and F D , which are the free energy of a free particle, are equal. 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.