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.
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 Zuckerman (2011); Chipot and Pohorille (2007); Frenkel and Smit (1996); Binder and Heermann (2010); Allen et al. (1987). In non equilibrium methods the work in the process of switching between the two Hamiltonians is measured Jarzynski (1997); Crooks (2000). 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 Jorgensen et al. (1987) or by normal mode classical calculations Lii and Allinger (1989) (see Entropies Section). While 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 Hayward and Go (1995).
Molecular modeling includes covalent bond, bond angle, dihedral angle, improper dihedral angle, electrostatic and vdW potentials (see Mayo et al. (1990); Phillips et al. (2005); Hess et al. (2008) 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 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 harmonic (quadratic) and by using the rigid rotor approximation and the HJR technique Boresch and Karplus (1999); Boresch et al. (2003). However, the free energies associated with dihedral potential terms in the realistic non-harmonic potentials and with non-harmonic covalent bond and bond angle terms Malde et al. (2011) 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 harmonic (see also Ref. Farhi and Singh (2015)). 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 I we present the first step of the calculation in which we transform the molecule by relaxing some potential terms of a group of atoms in that molecule and calculate the free energy difference associated with the transformation. In Section II we present the second step of the calculation in which the free energy associated with the remaining potentials (which depend on the relative spherical coordinates) of this group of atoms is calculated analytically or numerically. In Section III we suggest the potential application of free energy calculation of chemical reactions. The demonstartions of the method are explained in the Demonstration section.
Alchemical free energy transformation (as in Section I) is a standard procedure in molecular free energy calculations Hess et al. (2008); Phillips et al. (2005); Cuendet and Tuckerman (2012); Salomon-Ferrer et al. (2013); Mugnai and Elber (2012) and is usually used in the context of solvation and binding or in the context of chemical reactions in combination with QM calculations Jorgensen et al. (1987). 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 calculation with novel direct free energy calculations explained in Section II, followed by the novel application of free energy of chemical reactions presented in Section III. Here we do not need to assume harmonic dynamics.
I 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) Frenkel and Smit (1996); Kirkwood (1935); Straatsma and McCammon (1991). Alternatively, it can be performed using methods such as Exponential Averaging/ Free Energy Perturbation (FEP) Zwanzig (1954), Bennett Acceptance Ratio (BAR) Bennett (1976) and Weighted Histogram Analysis Method Kumar et al. (1992) (WHAM). We denote the Hamiltonian with the terms that are removed in the transformation by and the Hamiltonian with the other terms by . We define the Hamiltonian 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:
The free energy difference between the original and transformed systems can be written as:
where is the transformed system in which . Simulating intermediate systems at values in the range and calculating , followed by numerical integration, will enable us to calculate the free energy difference. It is noted that this calculation is valid for any environment (gas, 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 at Beutler et al. (1994).
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 includes the vdW and electrostatic interactions of the atoms O and H1 with all the atoms in the system, and the improper dihedral term .
Ii 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 , which will be chosen as the position of atoms relative to a covalently bound atom. 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, 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 , and defined with respect to the relevant atoms.
For the example of phenol we can write the partition function of transformed phenol as follows
where the degrees of freedom of atom H and O are denoted by and respectively. All the degrees of freedom of all the other atoms in the system (molecule+solvent molecules in the case of explicit solvent) and are denoted by . We define and as follows: . Our goal in this example is to calculate the free energy associated with all the degrees of freedom except for the ones denoted by . The free energy of the subsystem with the degrees of freedom denoted by 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 does depend on the coordinates of the atom which is included in and therfore 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:
The partition function of a covalent bond represented by a harmonic (quadratic) potential can be written as follows (for non-quadratic terms integrated analytically see Ref. Farhi and Singh (2015)):
where is an arbitrary length ( cancels out in comparisons). To account for the bond-dissociation energy/bond energy one has to differentiate between the potential at and (free atoms) either by introducing a factor Zavitsas (1987); Blanksby and Ellison (2003); Allinger et al. (1989) or by using potentials such as the Morse potential (numerical integration) Morse (1929).
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 harmonic (quadratic) term we can write :
This expression is real for positive and real values of and .
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 , 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
where is the modified Bessel function of the first kind. Note that this result does not depend on , which is an integer. In the context of restraints quadratic dihedral terms can be used and we give the exact expression for the corresponding partition function in the Appendix for completeness.
It can now be seen that the partition function of transformed phenol in Eq. (1) can be decomposed (exactly) into the partition function associated with the degrees of freedom and the covalent bond, bond angle and dihedral angle partition functions.
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:
where denote the bond angles of two bonds with respect to the principal bond and 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 , 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 , the integrations are decoupled. Thus, the integration for the case of one monomer that splits into two can be written as follows:
where is according to the definition in Eq. (4) and
For the general case it can be written as follows:
where can be calculated from Eq. (5) and , 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. For an explicit example with a bond junction see SM.
The free energy factors can be substituted in:
to give the free energy difference, where 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 represents the partition function of the submolecule that is fully interacting and and represent the th covalent bond and dihedral angle partition function respectively. and represent the th two bond and three or more bond junctions respectively and represents the th complex structure partition function. The arrow represents the transformation .
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 denotes the partition function of a free particle.
In Fig. 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 and degrees of freedom are associated with the submolecule A” and the degree of freedom is associated with the second system in the decomposition. The potential terms that depend on the coordinates of the and atoms in the five decomposed systems are: . The reason that the degrees of freedom and 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. (2) and the free energies of the fourth and fifth systems are calculated according to Eqs. (3) and (4) respectively.
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, for methanethiol molecule the free energies were calculated according to Eqs. (2)-(7) and were in agreement with the results obtained by performing transformations in MD simulations. The computations using these equations were faster by factors varying between and (for the same computational resources) Farhi and Singh (2015) .
Iii potential application
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 Abraham et al. (2015)) 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 one 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
we can match, if possible, molecules to molecules . Then we transform each of the molecules to and and calculate the free energy differences and respectively. Finally we calculate the free energy factors and as previously explained. Thus we get
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 . and , 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.
D.J. Bergman, A. Nitzan, D. Harries, G. Falkovich and B. Singh are acknowledged for the useful comments.
The exact expression for the partition function associated with a quadratic dihedral term is the following
- Zuckerman (2011) D. M. Zuckerman, Annual review of biophysics 40, 41 (2011).
- Chipot and Pohorille (2007) C. Chipot and A. Pohorille, Free energy calculations: theory and applications in chemistry and biology, Vol. 86 (Springer, 2007).
- Frenkel and Smit (1996) D. Frenkel and B. Smit, Understanding molecular simulation: from algorithms to applications (Academic Press, Inc. Orlando, FL, USA, 1996).
- Binder and Heermann (2010) K. Binder and D. W. Heermann, Monte Carlo simulation in statistical physics: an introduction (Springer, 2010).
- Allen et al. (1987) M. P. Allen, D. J. Tildesley, et al., Computer simulation of liquids (1987).
- Jarzynski (1997) C. Jarzynski, Physical Review Letters 78, 2690 (1997).
- Crooks (2000) G. Crooks, Phys. Rev. E 61, 2361 (2000).
- Jorgensen et al. (1987) W. L. Jorgensen, J. M. Briggs, and J. Gao, Journal of the American Chemical Society 109, 6857 (1987).
- Lii and Allinger (1989) J. H. Lii and N. L. Allinger, Journal of the American Chemical Society 111, 8566 (1989).
- Hayward and Go (1995) S. Hayward and N. Go, Annual Review of Physical Chemistry 46, 223 (1995).
- Mayo et al. (1990) S. L. Mayo, B. D. Olafson, and W. A. Goddard, Journal of Physical Chemistry 94, 8897 (1990).
- Phillips et al. (2005) J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kale, and K. Schulten, Journal of computational chemistry 26, 1781 (2005).
- Hess et al. (2008) B. Hess, C. Kutzner, D. Van Der Spoel, and E. Lindahl, Journal of chemical theory and computation 4, 435 (2008).
- Boresch and Karplus (1999) S. Boresch and M. Karplus, The Journal of Physical Chemistry A 103, 103 (1999).
- Boresch et al. (2003) S. Boresch, F. Tettinger, M. Leitgeb, and M. Karplus, The Journal of Physical Chemistry B 107, 9535 (2003).
- Malde et al. (2011) A. K. Malde, L. Zuo, M. Breeze, M. Stroet, D. Poger, P. C. Nair, C. Oostenbrink, and A. E. Mark, Journal of chemical theory and computation 7, 4026 (2011).
- Farhi and Singh (2015) A. Farhi and B. Singh, Arxiv preprint chem-ph/1501.01514 (2015).
- Cuendet and Tuckerman (2012) M. A. Cuendet and M. E. Tuckerman, Journal of Chemical Theory and Computation 8, 3504 (2012).
- Salomon-Ferrer et al. (2013) R. Salomon-Ferrer, D. A. Case, and R. C. Walker, Wiley Interdisciplinary Reviews: Computational Molecular Science 3, 198 (2013).
- Mugnai and Elber (2012) M. L. Mugnai and R. Elber, Journal of chemical theory and computation 8, 3022 (2012).
- Kirkwood (1935) J. Kirkwood, The Journal of Chemical Physics 3, 300 (1935).
- Straatsma and McCammon (1991) T. Straatsma and J. McCammon, The Journal of chemical physics 95, 1175 (1991).
- Zwanzig (1954) R. Zwanzig, The Journal of Chemical Physics 22, 1420 (1954).
- Bennett (1976) C. Bennett, Journal of Computational Physics 22, 245 (1976).
- Kumar et al. (1992) S. Kumar, J. Rosenberg, D. Bouzida, R. Swendsen, and P. Kollman, Journal of Computational Chemistry 13, 1011 (1992).
- Beutler et al. (1994) T. C. Beutler, A. E. Mark, R. C. van Schaik, P. R. Gerber, and W. F. van Gunsteren, Chemical physics letters 222, 529 (1994).
- Zavitsas (1987) A. A. Zavitsas, Journal of Physical Chemistry 91, 5573 (1987).
- Blanksby and Ellison (2003) S. J. Blanksby and G. B. Ellison, Accounts of Chemical Research 36, 255 (2003).
- Allinger et al. (1989) N. L. Allinger, Y. H. Yuh, and J. H. Lii, Journal of the American Chemical Society 111, 8551 (1989).
- Morse (1929) P. M. Morse, Physical Review 34, 57 (1929).
- Abraham et al. (2015) M. Abraham, D. Van Der Spoel, E. Lindahl, and B. Hess, GROMACS User Manual version 5 (2015).