Kin-Yiu Wong

*BioMaPS Institute for Quantitative Biology Rutgers, The State University of New Jersey USA* 

#### **1. Introduction**

An ultimate level of theory in molecular simulations [e.g., molecular dynamics (MD) and Monte Carlo (MC) simulations], which can accurately reproduce or even predict many experimental values, should be *ab initio* path integral. In *ab initio* path-integral simulations, both electrons and nuclei are treated quantum mechanically and adiabatically. No empirical parameter is involved, other than those fundamental physical constants (e.g., electronic mass and Planck's constant). The only *inherent* approximations are the Born-Oppenheimer approximation (to decouple internuclear dynamics from electronic motions) and the ergodicity in MD simulations or the importance samplings in MC simulations (to partly integrate the entire phase space). Consequently, correlation energy among electrons, anharmonic zero-point motions and tunnelling effects in nuclei, and isotope effects can all be incorporated in the simulations. Proper consideration of the electronic and internuclear quantum effects, even just partially, can be critical to compare computed values with state-of-the-art experiments, e.g., (I) hydrogen adsorption in carbon nanotechnology (Tanaka, Kanoh et al. 2005; Kowalczyk, Gauden et al. 2007; Kowalczyk, Gauden et al. 2008); (II) electronic redistributions and isotope effects (Wong and Gao 2007; Wong and Gao 2008; Wong, Richard et al. 2009; Gao and Wong 2008) on biochemical reactions in protein (Wong and Gao 2007; Wong and Gao 2011; Wu and Wong 2009; Warshel, Olsson et al. 2006; Gao, Major et al. 2008; Major, Heroux et al. 2009) and RNA enzymes (Wong, Lee et al. 2011; Wong, Gu et al. 2012).

However, owing to the extraordinarily high computational cost, *ab initio* path-integral simulations are thus far not practical even for modest size molecules, and are limited to only some relatively simpler or smaller molecular systems, e.g., thirty-two water molecules, and malonaldehyde [CH2(CHO)2]. Nevertheless, the unique information and invaluable insight for a molecular system, which can be provided perhaps only from *ab initio* path-integral simulations, have already been recognized in a number of pure computational publications in some high-profile journals, e.g., *Nature*, *Science*, and *Physical Review Letters*, etc (Marx and Parrinello 1995; Tuckerman, Marx et al. 1997; Marx, Tuckerman et al. 1999; Tuckerman and Marx 2001; Tuckerman, Marx et al. 2002; Ohta, Ohta et al. 2004; Hayashi, Shiga et al. 2006; Paesani, Iuchi et al. 2007).

In this chapter, after quickly going over the fundamental physical laws tailoring MD simulations, we (wongky@biomaps.rutgers.edu; kiniu@alumni.cuhk.net) discuss a new theoretical method that combines our novel systematic free-energy expansion approach, based on Zwanzig's free-energy perturbation theory, with our recently developed automated integration-free path-integral method, based on Kleinert's variational perturbation theory, (Wong and Gao 2007; Wong and Gao 2008; Wong, Richard et al. 2009; Wong, Gu et al. 2012) to perform *ab initio* path-integral simulations for realistic macromolecules at an affordable computational cost. Since in this new method, we can progressively choose computationally affordable levels of theory, now important physical quantities, e.g., free-energy barrier, change of binding energy, p*K*a value, and isotope effect, can all be computed at an *ab initio* path-integral level. Therefore, we anticipate this new systematic approach will become an essential computational tool to catch up with or even predict experimental results for breaking down subtle mechanisms underlying a variety of molecular systems in Life and Materials Sciences.

#### **2. Fundamental physical laws governing molecular dynamics simulations**

In this section, we lay the theoretical foundation for molecular dynamics (MD) simulations.

#### **2.1 Molecular Schrödinger equation**

Ever since quantum mechanics was constructed in the 1920s, solving the non-relativistic time-independent Schrödinger equation for a system of nuclei and electrons has become an essential step to understand *every single* detail of atomic or molecular properties (Kleppner and Jackiw 2000). The non-relativistic time-independent Schrödinger equation for a molecular system (the molecular Schrödinger equation) is:

$$
\hat{H}\_{mole}\Psi\_n = E\_n\Psi\_{n\prime} \tag{1}
$$

Developing a Systematic Approach for *Ab Initio* Path-Integral Simulations 109

Once the energy eigenvalues or the quantized energy spectrum in Eq. (1) are calculated, it is straightforward to obtain a central physical quantity in thermodynamics, i.e., the quantum canonical partition function *Qqm* (McQuarrie 2000), by the following summation of the

> *qm n* exp , *n Q E*

thermodynamic quantities for a system of nuclei and electrons, e.g., free energy, internal energy, entropy, pressure, etc., can be derived from it. In Eq. (3), the lowest energy level *E*<sup>0</sup> , which is often called the ground state energy or zero-point energy (ZPE), is usually the dominant energy level contributing to the partition function. Further, by virtue of Heisenberg's uncertainty principle, the ZPE is always larger than the minimum value of potential energy because a particle can never be at rest anywhere in a given potential or a

Unfortunately, even though *all* physics and chemistry of a (time-independent) molecular system is essentially in the molecular Schrödinger equation [Eq. (1)], it can be exactly solved only for simplest one-electron atoms or ions. For other systems, approximations must be introduced to calculate numerical solutions with the aid of computers. The most common and perhaps the mildest approximation is the Born-Oppenheimer approximation (Born and Oppenheimer 1927; Hirschfelder and Meath 1967; Kolos 1970; Ballhausen and Hansen 1972; Hehre, Radom et al. 1986; Szabo and Ostlund 1996; Helgaker, Jørgensen et al. 2000; Springborg 2000; Mielke, Peterson et al. 2003). It decouples internuclear motions from electrons so that nuclei effectively move on a potential energy surface (PES) obtained by

This approximation is based on the fact that an electron is much lighter than any nucleus (e.g., a proton, the lightest nucleus, is about 1840 times heavier than an electron). Nuclei move, consequently, much slowlier. As a result, from the electronic perspective, for a given set of nuclear positions, electrons adjust their positions 'instantly' before nuclei have a chance to move. On the other hand, from the standpoint of nuclei, electrons are moving so fast that their effects on nuclei are averaged out over the electronic wave functions. Mathematically, to simplify the molecular Hamiltonian, we first solve the electronic part of the Schrödinger equation for a particular set of nuclear configurations *xj* . The electronic part of the complete molecular Hamiltonian [Eq. (2)] is called electronic Hamiltonian:

> <sup>ˆ</sup> 1 1 <sup>2</sup> . <sup>2</sup> *N NN N e ne e <sup>j</sup>*

With this electronic Hamiltonian, we can obtain the electronic energy *Eelec* from the

*i j i ii ij ii*

*Z*

*r r*

(4)

*elec i*

*H*

corresponding electronic Schrödinger equation:

particle with a particular momentum can be everywhere in a given potential.

**2.3 Origin of potential energy surface: Born-Oppenheimer approximation** 

solving the electronic part of Schrödinger equation.

*k T* , *Bk* is Boltzmann's constant, and *T* is temperature. All standard

(3)

**2.2 Central quantity in quantum thermodynamics: Quantum partition function** 

Boltzmann energy distribution:

where 1 *<sup>B</sup>* 

where ˆ *Hmole* is the complete (non-relativistic) molecular Hamiltonian, *n* and *En* are an energy eigenfunction (or wave function) and an energy eigenvalue at an eigenstate *n*, respectively. In contrast to the (intra)nuclear or nucleon Hamiltonian (Dean 2007), the complete molecular Hamiltonian (Hehre, Radom et al. 1986; Szabo and Ostlund 1996; Kohn 1999; Pople 1999; Helgaker, Jørgensen et al. 2000; Springborg 2000) for *Nn* nuclei and *Ne* electrons can fortunately be written in an analytic closed form (thanks to the inverse squaredistance proportionality in Coulomb's electrostatic force law):

$$\hat{H}\_{mloc} = \sum\_{j}^{N\_{\text{tr}}} -\frac{1}{2M\_{j}}\nabla\_{j}^{2} + \sum\_{j$$

In Eq. (2), the units are atomic units, *Mj* is the mass ratio of nucleus *j* to electron, and *Zj* is the atomic number of nucleus *j*. The Laplacian operators <sup>2</sup> *j* and <sup>2</sup> *i* denote the second order differentiation with respect to the coordinates of the *j*th nucleus and the *i*th electron, respectively. The first term in Eq. (2) represents the kinetic energy operator for nuclei; the second term is the Coulomb repulsion between nuclei; the third term is the operator for the kinetic energy of electrons; the fourth and fifth terms indicate the Coulomb attraction between electrons and nuclei, and the repulsion between electrons, respectively. The distance between the *j*th and the *j'*th nuclei is *jj x* ; the separation between the *i*th and the *i′*th electrons is *ii r* ; the distance between the *j*th nucleus and the *i*th electrons is *ij r* .

theoretical method that combines our novel systematic free-energy expansion approach, based on Zwanzig's free-energy perturbation theory, with our recently developed automated integration-free path-integral method, based on Kleinert's variational perturbation theory, (Wong and Gao 2007; Wong and Gao 2008; Wong, Richard et al. 2009; Wong, Gu et al. 2012) to perform *ab initio* path-integral simulations for realistic macromolecules at an affordable computational cost. Since in this new method, we can progressively choose computationally affordable levels of theory, now important physical quantities, e.g., free-energy barrier, change of binding energy, p*K*a value, and isotope effect, can all be computed at an *ab initio* path-integral level. Therefore, we anticipate this new systematic approach will become an essential computational tool to catch up with or even predict experimental results for breaking down subtle mechanisms underlying a variety of

**2. Fundamental physical laws governing molecular dynamics simulations**  In this section, we lay the theoretical foundation for molecular dynamics (MD) simulations.

Ever since quantum mechanics was constructed in the 1920s, solving the non-relativistic time-independent Schrödinger equation for a system of nuclei and electrons has become an essential step to understand *every single* detail of atomic or molecular properties (Kleppner and Jackiw 2000). The non-relativistic time-independent Schrödinger equation for

*Hmole* is the complete (non-relativistic) molecular Hamiltonian, *n* and *En* are an energy eigenfunction (or wave function) and an energy eigenvalue at an eigenstate *n*, respectively. In contrast to the (intra)nuclear or nucleon Hamiltonian (Dean 2007), the complete molecular Hamiltonian (Hehre, Radom et al. 1986; Szabo and Ostlund 1996; Kohn 1999; Pople 1999; Helgaker, Jørgensen et al. 2000; Springborg 2000) for *Nn* nuclei and *Ne* electrons can fortunately be written in an analytic closed form (thanks to the inverse square-

> <sup>ˆ</sup> 11 1 2 2 . 2 2 *Nn Nn N NN N e ne e j j <sup>j</sup>*

In Eq. (2), the units are atomic units, *Mj* is the mass ratio of nucleus *j* to electron, and *Zj* is the atomic number of nucleus *j*. The Laplacian operators <sup>2</sup> *j* and <sup>2</sup> *i* denote the second order differentiation with respect to the coordinates of the *j*th nucleus and the *i*th electron, respectively. The first term in Eq. (2) represents the kinetic energy operator for nuclei; the second term is the Coulomb repulsion between nuclei; the third term is the operator for the kinetic energy of electrons; the fourth and fifth terms indicate the Coulomb attraction between electrons and nuclei, and the repulsion between electrons, respectively. The distance between the *j*th and the *j'*th nuclei is *jj x* ; the separation between the *i*th and the

*i′*th electrons is *ii r* ; the distance between the *j*th nucleus and the *i*th electrons is *ij r* .

*j j j j jj i j i ii ij ii*

*M x r r*

*Z Z Z*

 (2)

*H E mole n n n* , (1)

molecular systems in Life and Materials Sciences.

a molecular system (the molecular Schrödinger equation) is:

distance proportionality in Coulomb's electrostatic force law):

ˆ

*mole j i*

**2.1 Molecular Schrödinger equation** 

*H*

where ˆ

#### **2.2 Central quantity in quantum thermodynamics: Quantum partition function**

Once the energy eigenvalues or the quantized energy spectrum in Eq. (1) are calculated, it is straightforward to obtain a central physical quantity in thermodynamics, i.e., the quantum canonical partition function *Qqm* (McQuarrie 2000), by the following summation of the Boltzmann energy distribution:

$$Q\_{qm} = \sum\_{n} \exp\left(-\beta E\_n\right),\tag{3}$$

where 1 *<sup>B</sup> k T* , *Bk* is Boltzmann's constant, and *T* is temperature. All standard thermodynamic quantities for a system of nuclei and electrons, e.g., free energy, internal energy, entropy, pressure, etc., can be derived from it. In Eq. (3), the lowest energy level *E*<sup>0</sup> , which is often called the ground state energy or zero-point energy (ZPE), is usually the dominant energy level contributing to the partition function. Further, by virtue of Heisenberg's uncertainty principle, the ZPE is always larger than the minimum value of potential energy because a particle can never be at rest anywhere in a given potential or a particle with a particular momentum can be everywhere in a given potential.

#### **2.3 Origin of potential energy surface: Born-Oppenheimer approximation**

Unfortunately, even though *all* physics and chemistry of a (time-independent) molecular system is essentially in the molecular Schrödinger equation [Eq. (1)], it can be exactly solved only for simplest one-electron atoms or ions. For other systems, approximations must be introduced to calculate numerical solutions with the aid of computers. The most common and perhaps the mildest approximation is the Born-Oppenheimer approximation (Born and Oppenheimer 1927; Hirschfelder and Meath 1967; Kolos 1970; Ballhausen and Hansen 1972; Hehre, Radom et al. 1986; Szabo and Ostlund 1996; Helgaker, Jørgensen et al. 2000; Springborg 2000; Mielke, Peterson et al. 2003). It decouples internuclear motions from electrons so that nuclei effectively move on a potential energy surface (PES) obtained by solving the electronic part of Schrödinger equation.

This approximation is based on the fact that an electron is much lighter than any nucleus (e.g., a proton, the lightest nucleus, is about 1840 times heavier than an electron). Nuclei move, consequently, much slowlier. As a result, from the electronic perspective, for a given set of nuclear positions, electrons adjust their positions 'instantly' before nuclei have a chance to move. On the other hand, from the standpoint of nuclei, electrons are moving so fast that their effects on nuclei are averaged out over the electronic wave functions. Mathematically, to simplify the molecular Hamiltonian, we first solve the electronic part of the Schrödinger equation for a particular set of nuclear configurations *xj* . The electronic part of the complete molecular Hamiltonian [Eq. (2)] is called electronic Hamiltonian:

$$\hat{H}\_{\text{elec}} = -\frac{1}{2} \sum\_{i}^{N\_c} \nabla\_i^2 - \sum\_{j}^{N\_d} \sum\_{i}^{N\_c} \frac{Z\_j}{r\_{ij}} + \sum\_{i$$

With this electronic Hamiltonian, we can obtain the electronic energy *Eelec* from the corresponding electronic Schrödinger equation:

$$
\hat{H}\_{elec}\psi\_{elec} = E\_{elec}\left(\left\{\mathbf{x}\_{j}\right\}\right)\psi\_{elec} \tag{5}
$$

Developing a Systematic Approach for *Ab Initio* Path-Integral Simulations 111

In practice, quantum effects on internuclear motions are much smaller than those on the electronic part. In many applications, the internuclear quantum effects are insignificant and could even be neglected. Thus, the eigenenergy spectrum *En* in Eq. (1) would become continuous. Given an internuclear potential *V*, the quantum canonical partition function in

3 3 3 2

*cl N j*

where *h* is Planck's constant and *p* is the momenta associated with the nuclear coordinates *x*. Subsequently, the classical free energy *Gcl* of a molecular system can be expressed in terms

*n*

Note that the partition function and the free energy defined above are 'state' functions, which is independent of any nuclear coordinate and momentum (as we integrate out the entire phase space). Given a particular 3 *Nn* -degree-of-freedom molecular system described by a particular potential energy function *V* at particular temperature, the partition function

On the other hand, of significant interest in simulating a many-body biochemical or physical event is to examine how the free energy of a molecular system *varies* during the event. Conventionally, we first predetermine a coordinate which should be able to describe the event of interest from the start to the end. Next, we generate a free energy profile, which is an energy function of that predetermined coordinate, to investigate how the profile changes during the event. In fact, such a kind of free-energy profile can also be termed as potential energy of ensemble-average or mean force (Kirkwood 1935). Reasons are given below.

The free energy profile of a molecular system as a function of a predetermined coordinate of

freedom along *z*-direction, and *C* is a normalization factor dependent on the inverse of the thermal de Broglie wavelengths for all degrees of freedom (the wavelength is a function of the nuclear mass *Mj*, and temperature *T*). *C* should be a constant during the biochemical or physical event of our interest. The integrand in the final configurational integral of Eq. (10)

ln exp <sup>2</sup>

*dx dp <sup>p</sup> G z kT x z V x*

 

*z B N z z j*

3 3 3 2

*h M*

*n n n*

*N N N*

3

*k T dx Vx x z C*

ln exp , ,

*n*

3 1

*N*

*n*

*B j z*

*cl B cl B N j*

*dx dp <sup>p</sup> G kT Q kT V x*

*n n n*

*N N N*

*dx dp <sup>p</sup> <sup>Q</sup> V x h M* 

<sup>3</sup> exp , <sup>2</sup>

*j*

(8)

*j j*

3 3 3 2

*h M* 

*N N N*

<sup>3</sup> ln ln exp . <sup>2</sup> *n n n*

*j*

*j*

(10)

*j j*

*<sup>z</sup>* is the thermal de Broglie wavelength for the degree of

*j j*

(9)

**2.4 Classical free-energy profile vs classical potential of mean force** 

Eq. (3) consequently reduces to the classical canonical partition function as:

*n*

of the classical partition function *Qcl* as follows:

and the free energy are *constants*.

interest *z* can be written as follows:

is Dirac delta function,

where  where *elec* is the electronic wave function. Note that the electronic energy *E x elec <sup>j</sup>* depends parametrically on the nuclear positions *xj* . With this electronic energy, the molecular Hamiltonian in Eq. (2) can be simplified as follows:

$$\begin{split} \hat{H}\_{mule} &\approx \sum\_{j}^{N\_{\text{e}}} -\frac{1}{2M\_{j}} \nabla\_{j}^{2} + \sum\_{j$$

where signifies the average over electronic wave functions or the expectation value. In Eq. (6), *V* is defined as the sum of the nuclear repulsion energy and electronic energy, which effectively turns out to be the internuclear potential energy function as a consequence of the Born-Oppenheimer approximation:

$$V\left(\left\{\mathbf{x}\_{j}\right\}\right) \equiv \sum\_{j$$

There are many systematic and rigorous theories in electronic structure calculations to derive the internuclear potential energy from first principles (i.e., besides the universal fundamental constants in physics, there is no other empirical parameter involved in the calculations), e.g., Hartree-Fock theory, configuration interaction method, Møller-Plesset perturbation theory, coupled cluster approach, and Kohn-Sham density functional theory. All these quantum mechanical (QM) approaches for electronic structure calculations are often known as *ab initio* methods (Hehre, Radom et al. 1986; Szabo and Ostlund 1996; Kohn 1999; Pople 1999; Helgaker, Jørgensen et al. 2000; Springborg 2000).

In contrast, a complete empirical method to determine an internuclear potential energy surface is to parameterize the potential energy as an analytic function without treating electronic degrees of freedom. This type of approach is referred to as molecular mechanical (MM) method and the empirical potential energy is called force-field energy. Comparing to *ab initio* approach, MM methods are computationally much less expensive and can be applied to describe equilibrium properties in macromolecular systems involving over tens of thousands of heavy atoms (Hagler, Huler et al. 1974; Brooks, Bruccoleri et al. 1983; Weiner, Kollman et al. 1984; Jorgensen and Tirado-Rives 1988; Mayo, Olafson et al. 1990). But for the process involving electronic redistributions (e.g., electronic transfer, chemical bond breaking or forming, etc.), MM force field is often unable to describe it. Later, a hybrid approach called combined QM/MM method has emerged to synthesize the efficiency of MM force field with the accuracy of QM calculations (Field, Bash et al. 1990; Gao and Truhlar 2002). For the rest of this chapter, discussions are limited to the Born-Oppenheimer approximation, which adiabatically decouples nuclear and electronic degrees of freedom.

<sup>ˆ</sup> , *H Ex elec elec elec*

2 2

*n n e ne e*

where signifies the average over electronic wave functions or the expectation value. In Eq. (6), *V* is defined as the sum of the nuclear repulsion energy and electronic energy, which effectively turns out to be the internuclear potential energy function as a consequence of the

> . *Nn j j j elec j j j jj*

There are many systematic and rigorous theories in electronic structure calculations to derive the internuclear potential energy from first principles (i.e., besides the universal fundamental constants in physics, there is no other empirical parameter involved in the calculations), e.g., Hartree-Fock theory, configuration interaction method, Møller-Plesset perturbation theory, coupled cluster approach, and Kohn-Sham density functional theory. All these quantum mechanical (QM) approaches for electronic structure calculations are often known as *ab initio* methods (Hehre, Radom et al. 1986; Szabo and Ostlund 1996; Kohn

In contrast, a complete empirical method to determine an internuclear potential energy surface is to parameterize the potential energy as an analytic function without treating electronic degrees of freedom. This type of approach is referred to as molecular mechanical (MM) method and the empirical potential energy is called force-field energy. Comparing to *ab initio* approach, MM methods are computationally much less expensive and can be applied to describe equilibrium properties in macromolecular systems involving over tens of thousands of heavy atoms (Hagler, Huler et al. 1974; Brooks, Bruccoleri et al. 1983; Weiner, Kollman et al. 1984; Jorgensen and Tirado-Rives 1988; Mayo, Olafson et al. 1990). But for the process involving electronic redistributions (e.g., electronic transfer, chemical bond breaking or forming, etc.), MM force field is often unable to describe it. Later, a hybrid approach called combined QM/MM method has emerged to synthesize the efficiency of MM force field with the accuracy of QM calculations (Field, Bash et al. 1990; Gao and Truhlar 2002). For the rest of this chapter, discussions are limited to the Born-Oppenheimer approximation, which adiabatically decouples nuclear and electronic degrees of freedom.

*Z Z V x E x x* 

*j j j j jj i j i ii ij ii*

*M x r r*

*E x*

*N N N NN N*

11 1 ˆ

*j j j elec j*

*Z Z*

 

 *elec* is the electronic wave function. Note that the electronic energy *E x elec <sup>j</sup>* depends parametrically on the nuclear positions *xj* . With this electronic energy, the

 

*j j j*

*Z Z Z*

*<sup>j</sup> elec* (5)

(7)

(6)

2 2

*mole j i*

molecular Hamiltonian in Eq. (2) can be simplified as follows:

1 2

*n*

*j j*

*M*

*N*

*H*

Born-Oppenheimer approximation:

2

*j j j j jj*

*n n*

*N N*

2

1999; Pople 1999; Helgaker, Jørgensen et al. 2000; Springborg 2000).

<sup>1</sup> , <sup>2</sup>

*M x*

*j j*

*V x*

where

#### **2.4 Classical free-energy profile vs classical potential of mean force**

In practice, quantum effects on internuclear motions are much smaller than those on the electronic part. In many applications, the internuclear quantum effects are insignificant and could even be neglected. Thus, the eigenenergy spectrum *En* in Eq. (1) would become continuous. Given an internuclear potential *V*, the quantum canonical partition function in Eq. (3) consequently reduces to the classical canonical partition function as:

$$Q\_{cl} = \bigcap\_{-\alpha \to -\infty}^{\alpha} \frac{dx^{3N\_n} dp^{3N\_n}}{h^{3N\_n}} \exp\left\{-\beta \left[ \left(\sum\_{j}^{3N\_n} \frac{p\_j^2}{2M\_j}\right) + V\left(\left\{x\_j\right\}\right) \right] \right\},\tag{8}$$

where *h* is Planck's constant and *p* is the momenta associated with the nuclear coordinates *x*. Subsequently, the classical free energy *Gcl* of a molecular system can be expressed in terms of the classical partition function *Qcl* as follows:

$$\mathbf{G}\_{cl} = -k\_B T \ln Q\_{cl} = -k\_B T \ln \int\_{-\infty}^{\infty} \int \frac{d\mathbf{x}^{3N\_n} dp^{3N\_n}}{h^{3N\_n}} \exp\left\{-\beta \left[ \left(\sum\_{j}^{3N\_p} \frac{p\_j^2}{2M\_j}\right) + V\left(\left\{x\_j\right\}\right) \right] \right\}.\tag{9}$$

Note that the partition function and the free energy defined above are 'state' functions, which is independent of any nuclear coordinate and momentum (as we integrate out the entire phase space). Given a particular 3 *Nn* -degree-of-freedom molecular system described by a particular potential energy function *V* at particular temperature, the partition function and the free energy are *constants*.

On the other hand, of significant interest in simulating a many-body biochemical or physical event is to examine how the free energy of a molecular system *varies* during the event. Conventionally, we first predetermine a coordinate which should be able to describe the event of interest from the start to the end. Next, we generate a free energy profile, which is an energy function of that predetermined coordinate, to investigate how the profile changes during the event. In fact, such a kind of free-energy profile can also be termed as potential energy of ensemble-average or mean force (Kirkwood 1935). Reasons are given below.

The free energy profile of a molecular system as a function of a predetermined coordinate of interest *z* can be written as follows:

$$\begin{split} \mathcal{G}\_{z}\{\boldsymbol{z}\} &= -k\_{B}T \ln \left[ \int\_{-\alpha-\alpha}^{\alpha} \int \frac{d\boldsymbol{x}^{\boldsymbol{3N}\_{n}} dp^{\boldsymbol{3N}\_{n}}}{h^{\boldsymbol{3N}\_{n}}} \mathcal{A}\_{z} \mathcal{S}\{\boldsymbol{x}\_{z} - \boldsymbol{z}\} \exp\left\{-\beta \left[ \left(\sum\_{j}^{\boldsymbol{3N}\_{n}} \frac{p\_{j}^{2}}{2\mathcal{M}\_{j}}\right) + V\left(\left\{\boldsymbol{x}\_{j}\right\}\right) \right] \right\} \right] \\ &= -k\_{B}T \ln \left[ \int\_{-\alpha}^{\alpha} d\boldsymbol{x}^{\boldsymbol{3N}\_{n} - 1} \exp\left\{-\beta \left[ V\left(\left\{\boldsymbol{x}\_{j}\right\}, \boldsymbol{x}\_{z} = \boldsymbol{z}\right) \right] \right\} \right] + \mathcal{C}, \end{split} \tag{10}$$

where is Dirac delta function, *<sup>z</sup>* is the thermal de Broglie wavelength for the degree of freedom along *z*-direction, and *C* is a normalization factor dependent on the inverse of the thermal de Broglie wavelengths for all degrees of freedom (the wavelength is a function of the nuclear mass *Mj*, and temperature *T*). *C* should be a constant during the biochemical or physical event of our interest. The integrand in the final configurational integral of Eq. (10)

Developing a Systematic Approach for *Ab Initio* Path-Integral Simulations 113

rest of coordinates, then the PMF will be unique and its relation with the free-energy profile

In Eq. (14), the Leibnizian contribution is nil, the first term on RHS is the mean force for *q*

Finally, the Fixman potential (Fixman 1974), which corrects the velocity-bias in constrained MD, will also be presented with correct dependence on mass in our forthcoming paper.

By assuming the molecular system of our interest is ergodic, molecular dynamics (MD) simulation techniques can be employed to compute the ensemble average of a physical quantity. In essence, MD simulations is numerically solving, integrating or propagating the Newtonian equations of motion, one-time-step by one-time-step. Given an internuclear potential *V* (regardless of using QM, MM, or hybrid QM/MM to construct), the motion or

0 0

 

(14)

2 2 . *j*

*d x*

(15)

, exp , <sup>2</sup>

*<sup>t</sup>* (17)

*x p V x*

(16)

*j*

*dt*

 

> 

is the ensemble average over all

,

*q q*

 

2 2

*B*

can be proved as follows (den Otter 2000), after using Eq. (12) and Eq. (13):

*dG <sup>q</sup> <sup>q</sup> <sup>V</sup> k T d q q* 

**2.5 Simulating classical thermodynamics: Molecular dynamics simulations** 

trajectory of a nucleus *j* as a function of time *t* is governed by Newton's second law:

can be either a scalar or a vector) over the *entire* phase space, i.e.,

3

*n*

*f*

*t f*

1

is equal to the time average in MD simulations:

**Extended Forces**

*j j j*

Note the extended forces in Eq. (15) are *essential* for having *canonical* ensemble (constant temperature) instead of *microcanonical* ensemble (constant energy) in MD simulations (Hünenberger 2005). In the ergodic hypothesis (Lebowitz and Penrose 1973; Cogswell 1999) [the dynamical version of ergodic theory was first proposed by Birkhoff (Birkhoff 1931), in which Liouville's theorem was applied to ensure the ensemble distribution in phase-space is invariant with time], if the simulation time for propagating the trajectory *<sup>j</sup> x t* of the nucleus *j* is *infinitely long*, the ensemble average of a physical quantity **f***x p* , (which

> **f f** 3 3 3 2

> > **f f** 0

In other words, longer MD simulation time allows us to sample more phase space for computing the corresponding ensemble average, which in turn could be in higher accuracy.

*f*

*t*

<sup>1</sup> , . lim

*cl j j*

*Q M h*

*N N N*

*dx dp p*

*n n n*

*N j*

*x t p t dt*

*Vxt M*

0

 .

configurations with 0 *q*

0

the second term is the Jacobian contribution, and 0 *<sup>q</sup>*

is basically the probability density of the molecular system as a function of *z*. In practice, it is rare to determine the value of *C* because what we often care about is the free-energy difference at various values of *z*.

Notably, by taking the negative derivative of *G z <sup>z</sup>* , i.e., *dG z dz <sup>z</sup>* , we obtain the average force over all ensembles or over all degrees of freedom, which is called the mean force (Kirkwood 1935), based on the ensemble average definition in Eq. (16):

$$\begin{aligned} & \left\{ \left[ \mathbf{F} \left( \mathbf{z} \right) \right] \right\} \\ &= \left[ \int\_{-\infty}^{\infty} \frac{dx^{3N\_s} dp^{3N\_s}}{h^{3N\_s}} \lambda\_s \delta \left( x\_z - z \right) \right] - \bar{\nabla}\_{x\_j} V \left( \left\{ x\_j \right\} \right) \right]\_{x\_j = z} \exp\left\{ -\rho \left[ \left( \sum\_{j}^{3N\_s} \frac{p\_j^2}{2M\_j} \right) + V \left( \left\{ x\_j \right\} \right) \right] \right\} \\ &= \frac{\stackrel{\circ}{\mathbb{1}}}{\stackrel{\circ}{\to}} \frac{\stackrel{\circ}{\text{e}}}{h^{3N\_s}} \frac{dx^{3N\_s} dp^{3N\_s}}{h^{3N\_s}} \lambda\_z \delta \left( x\_z - z \right) \exp\left\{ -\rho \left[ \left( \sum\_{j}^{3N\_s} \frac{p\_j^2}{2M\_j} \right) + V \left( \left\{ x\_j \right\} \right) \right] \right\}} \\ &= \frac{\stackrel{\circ}{\text{e}}}{\stackrel{\circ}{\to}} \frac{dx^{3N\_s - 1} \left[ -\bar{\nabla}\_z V \left( \left\{ x\_j \right\} , x\_z = z \right) \right] \exp\left\{ -\rho \left[ \sum V \left( \left\{ x\_j \right\} , x\_z = z \right) \right] \right\}}{\stackrel{\circ}{\to}} = \frac{-dG\_z \left( z \right)}{dz} \hat{\mathbf{z}}. \end{aligned} \tag{11}$$

Thus *G z <sup>z</sup>* , the free energy profile as a function of a predetermined coordinate, is also called the potential of mean force (PMF) (Kirkwood 1935).

However, please note that if the predetermined coordinate of interest is *not* a linear combination of rectilinear coordinates, or in other words, if it is a curvilinear coordinate, then PMF is oftentimes *not* exactly equal to free-energy profile. Not only the Jacobiandeterminant contribution makes their difference (Ruiz-Montero, Frenkel et al. 1997; Hénin, Fiorin et al. 2010), but also in a forthcoming paper, we will show that actually change of domains with respect to the coordinate of interest can also contribute to the free-energy profile, i.e., the Leibnizian contribution (Flanders 1973).

In addition, we will also show that according to differential geometry and general relativity, once we realize the equivalence between orthogonal covariant and contravariant vectors (Arfken and Weber 2001), then the Jacobian scale factor for a predetermined curvilinear coordinate of interest, *q*, can be proved to be (in contravariant space):

$$\mathcal{H}\_{q\_{\vec{\xi}}} = \left| \vec{\nabla} q\_{\vec{\xi}} \right|^{-1} \tag{12}$$

and the unit vector for *q*can be proved as (in contravariant space):

$$
\hat{q}\_{\xi} = \bar{\nabla} q\_{\xi} \not\!\Big| \bar{\nabla} q\_{\xi} \bigg| \tag{13}
$$

In Eq. (12) and (13), *q* must belong to at least one *complete* set of curvilinear coordinates, hypothetically. In general, unless we explicitly define the rest of the *complete* curvilinear coordinates, the sole definition of *q* is *not* sufficient to make the PMF be unique. But, we will show the free-energy profile does *not* suffer from this uniqueness problem. In fact, if we restrict ourselves to a *complete* set of curvilinear coordinates in which *q*is orthogonal to the

is basically the probability density of the molecular system as a function of *z*. In practice, it is rare to determine the value of *C* because what we often care about is the free-energy

Notably, by taking the negative derivative of *G z <sup>z</sup>* , i.e., *dG z dz <sup>z</sup>* , we obtain the average force over all ensembles or over all degrees of freedom, which is called the mean force

3 2 3 3

*N N <sup>N</sup> <sup>j</sup>*

*n n n*

3 2 3 3

*N N <sup>N</sup> <sup>j</sup>*

*n n n*

*j z*

Thus *G z <sup>z</sup>* , the free energy profile as a function of a predetermined coordinate, is also

However, please note that if the predetermined coordinate of interest is *not* a linear combination of rectilinear coordinates, or in other words, if it is a curvilinear coordinate, then PMF is oftentimes *not* exactly equal to free-energy profile. Not only the Jacobiandeterminant contribution makes their difference (Ruiz-Montero, Frenkel et al. 1997; Hénin, Fiorin et al. 2010), but also in a forthcoming paper, we will show that actually change of domains with respect to the coordinate of interest can also contribute to the free-energy

In addition, we will also show that according to differential geometry and general relativity, once we realize the equivalence between orthogonal covariant and contravariant vectors (Arfken and Weber 2001), then the Jacobian scale factor for a predetermined curvilinear

, can be proved to be (in contravariant space):

can be proved as (in contravariant space):

*q qq* ˆ

 

hypothetically. In general, unless we explicitly define the rest of the *complete* curvilinear

will show the free-energy profile does *not* suffer from this uniqueness problem. In fact, if we

restrict ourselves to a *complete* set of curvilinear coordinates in which *q*

*<sup>q</sup> h q* 

1

 

must belong to at least one *complete* set of curvilinear coordinates,

 

exp <sup>2</sup>

 

*z*

*dG z dz*

(12)

(13)

is orthogonal to the

is *not* sufficient to make the PMF be unique. But, we

ˆ.

(11)

*z*

*j j*

*x z*

,

 

*j z*

exp <sup>2</sup>

*<sup>N</sup> zz x j <sup>j</sup> x z <sup>j</sup> <sup>j</sup>*

*dx dp <sup>p</sup> x z Vx V x h M*

> *dx dp <sup>p</sup> x z V x h M*

*N z z j*

(Kirkwood 1935), based on the ensemble average definition in Eq. (16):

*dx V x x z V x*

, exp

exp , *<sup>n</sup>*

*dx Vx x z*

*<sup>z</sup> <sup>n</sup> <sup>z</sup>*

 

3

 

*z jz*

3 1

called the potential of mean force (PMF) (Kirkwood 1935).

profile, i.e., the Leibnizian contribution (Flanders 1973).

coordinates, the sole definition of *q*

*N*

 

*n*

difference at various values of *z*.

*z*

3

3 1

*N*

*n*

**F**

coordinate of interest, *q*

and the unit vector for *q*

In Eq. (12) and (13), *q*

rest of coordinates, then the PMF will be unique and its relation with the free-energy profile can be proved as follows (den Otter 2000), after using Eq. (12) and Eq. (13):

$$\frac{d\mathbf{G}\_{\xi}\left(\xi\_{0}\right)}{d\xi\_{0}} = \left\langle \bar{\nabla}V \cdot \left(\frac{\bar{\nabla}q\_{\xi}}{\left|\bar{\nabla}q\_{\xi}\right|^{2}}\right) \right\rangle\_{q\_{\xi}=\xi\_{0}} - k\_{B}T \left\langle \bar{\nabla} \cdot \left(\frac{\bar{\nabla}q\_{\xi}}{\left|\bar{\nabla}q\_{\xi}\right|^{2}}\right) \right\rangle\_{q\_{\xi}=\xi\_{0}} \tag{14}$$

In Eq. (14), the Leibnizian contribution is nil, the first term on RHS is the mean force for *q* , the second term is the Jacobian contribution, and 0 *<sup>q</sup>* is the ensemble average over all configurations with 0 *q* .

Finally, the Fixman potential (Fixman 1974), which corrects the velocity-bias in constrained MD, will also be presented with correct dependence on mass in our forthcoming paper.

#### **2.5 Simulating classical thermodynamics: Molecular dynamics simulations**

By assuming the molecular system of our interest is ergodic, molecular dynamics (MD) simulation techniques can be employed to compute the ensemble average of a physical quantity. In essence, MD simulations is numerically solving, integrating or propagating the Newtonian equations of motion, one-time-step by one-time-step. Given an internuclear potential *V* (regardless of using QM, MM, or hybrid QM/MM to construct), the motion or trajectory of a nucleus *j* as a function of time *t* is governed by Newton's second law:

$$-\bar{\nabla}\_{\dot{\boldsymbol{\beta}}}V\left(\left\{\boldsymbol{x}\_{\dot{\boldsymbol{\beta}}}(t)\right\}\right) + \left(\mathbf{Extended}\ \mathbf{Fores}\right) = M\_{\dot{\boldsymbol{\beta}}}\frac{d^2\bar{\boldsymbol{x}}\_{\dot{\boldsymbol{\beta}}}}{dt^2}.\tag{15}$$

Note the extended forces in Eq. (15) are *essential* for having *canonical* ensemble (constant temperature) instead of *microcanonical* ensemble (constant energy) in MD simulations (Hünenberger 2005). In the ergodic hypothesis (Lebowitz and Penrose 1973; Cogswell 1999) [the dynamical version of ergodic theory was first proposed by Birkhoff (Birkhoff 1931), in which Liouville's theorem was applied to ensure the ensemble distribution in phase-space is invariant with time], if the simulation time for propagating the trajectory *<sup>j</sup> x t* of the nucleus *j* is *infinitely long*, the ensemble average of a physical quantity **f***x p* , (which can be either a scalar or a vector) over the *entire* phase space, i.e.,

$$\left\langle \mathbf{f} \right\rangle = \frac{1}{Q\_{cl}} \int\_{-\infty}^{\infty} \int \frac{d\mathbf{x}^{3N\_{\text{s}}} dp^{3N\_{\text{s}}}}{h^{3N\_{\text{s}}}} \, \mathbf{f} \left( \left\langle \mathbf{x} \right\rangle, \left\langle p \right\rangle \right) \exp\left[ -\beta \left[ \left( \sum\_{j}^{3N\_{\text{s}}} \frac{p\_{j}^{2}}{2M\_{j}} \right) + V \left( \left\langle \mathbf{x}\_{j} \right\rangle \right) \right] \right],\tag{16}$$

is equal to the time average in MD simulations:

$$\langle \mathbf{f} \rangle = \lim\_{t\_f \to \infty} \frac{1}{t\_f} \Big| \mathbf{f} \left( \{ \mathbf{x}(t) \} , \{ p(t) \} \right) dt. \tag{17}$$

In other words, longer MD simulation time allows us to sample more phase space for computing the corresponding ensemble average, which in turn could be in higher accuracy.

#### **3. Zwanzig's free-energy perturbation theory**

Owing to the Boltzmann exponential energy distribution, one of the major difficulties in computing a converged free-energy profile or potential of mean force [Eq. (10)] via MD and MC sampling techniques is that it takes longer simulation time or runs more MC steps to have enough higher-energy samples. Yet, many interesting biochemical or physical molecular properties could be in higher-energy regions, e.g., the transition state during protein folding or biochemical reaction.

In practice, in order for having effective samplings on both the lower-energy (e.g., reactant state) and higher-energy regions (e.g., transition state), Zwanzig's free-energy perturbation (Zwanzig 1954) [which is also referred to as statistical-mechanical perturbation theory (McQuarrie 2000)] has been extensively applied. The feature of the perturbation is relating the change of free energy between *two* systems (both have the same number of degrees of freedom) by an ensemble average taken in only *one* of the two systems. This can be illustrated by first writing the classical free energy *Gcl* corresponding to the partition function in Eq. (8) as follows:

$$\mathcal{G}\_{cl} = -k\_B T \ln Q\_{cl} = \int\_0^\alpha \int \frac{d\mathbf{x}^{3N\_n} dp^{3N\_n}}{h^{3N\_n}} \exp(-\beta E),\tag{18}$$

Developing a Systematic Approach for *Ab Initio* Path-Integral Simulations 115

**4. Systematic** *ab initio* **molecular dynamics approach: Free-energy expansion** 

A fundamental key to have successful molecular simulations is the accuracy of internuclear potential for describing atomic motions during biochemical or physical events. By exploiting Zwanzig's free-energy perturbation (FEP) theory, we are developing a new rigorous method to systematically obtain accurate free-energy profiles, in which the internuclear potential energy is effectively computed at a high-level *ab initio* theory. Our new method is a systematic free-energy expansion (FEE) in terms of a series of covariance tensors. The new expansion will enable us to have a free-energy profile at a level as high as the coupled cluster theory at an affordable computational cost, which is currently known as the gold standard but unreachable level of theory for free-energy simulations. The focus of our FEE method will be on the difference of free energy calculated by two different internuclear potential. Furthermore, in contrast to Car-Parrinello MD (CPMD) which is limited to potential energy derived from DFT (Car and Parrinello 1985), our method is independent of how the potential energy functions being constructed. Therefore, by combining it with our novel automated integration-free path-integral (AIF-PI) method together (See Section 5; Wong and Gao 2007; Wong and Gao 2008; Wong, Richard et al. 2009; Wong, Gu et al. 2012), we will also be able to compute free-energy barriers, changes of binding energy, p*K*a values,

Let's begin with the FEP theory. From Eq. (21), the free energy difference between using

*HL LL* ln exp . *<sup>B</sup> HL LL LL G G G kT*

Next we expand the ensemble average in Eq. (22) and sum up the prefactors into a series of

*HL LL LL LL c <sup>n</sup>*

 *LL c*, is a cumulant, and *n* is the order of a cumulant. In his original 1954 paper (Zwanzig 1954), Zwanzig showed that the cumulant expansion is fast converging when the change of

of computational cost, this cumulant expansion does not provide an advantage for correcting lower-level free energy. This is because the time required for calculating the cumulant average *LL c*, with computer is basically as much as the time needed to directly compute the higher-level free energy *GHL* , regardless of whether the perturbation Δ*E* is big. In order to ease up this situation, in a forthcoming paper we will prove that each cumulant can be further expanded as a Taylor series expansion fluctuating about the ensemble

energy Δ*E* in the ensemble is reasonably small relative to the inverse of

 

(23)

  , <sup>1</sup> 1 1 ln exp ln exp , *<sup>n</sup>*

*E E*

*E E* (22)

. However, in terms

 

*EE E HL LL* , (24)

**method as a series of covariance tensors** 

and isotope effects at an *ab initio* path integral level (see Section 6).

cumulants:

where

*GG G*

average position **x***LL* in the form:

lower-level (LL) and higher-level (HL) *ab initio* methods can be expressed as:

where *E* is the energy at a point *pj j* ,*x* in the phase space, i.e.,

$$E = E\left(\left\{p\_j\right\}, \left\{\mathbf{x}\_j\right\}\right) = \left(\sum\_{j}^{3N\_v} \frac{p\_j^2}{2M\_j}\right) + V\left(\left\{\mathbf{x}\_j\right\}\right). \tag{19}$$

Next, we rewrite Eq. (18) as:

$$\mathbf{G}\_{cl} = -k\_B T \ln \left[ \frac{\int^{\alpha} \int \frac{\alpha}{h^{3N\_v}} d\boldsymbol{p}^{3N\_v}}{h^{3N\_v}} \exp\left[ -\beta (\mathcal{E} - \mathcal{E}\_0) \right] \mathbf{e}^{-\beta \mathcal{E}\_0}}{\int^{\alpha} \int \frac{\alpha^{3N\_v} d\boldsymbol{p}^{3N\_v}}{h^{3N\_v}} \mathbf{e}^{-\beta \mathcal{E}\_0}} \right] - k\_B T \ln \int^{\alpha} \int \frac{\alpha^{3N\_v} d\boldsymbol{p}^{3N\_v}}{h^{3N\_v}} \mathbf{e}^{-\beta \mathcal{E}\_0} \tag{20}$$
 
$$= -k\_B T \ln \left\{ \exp\left[ -\beta (\mathcal{E} - \mathcal{E}\_0) \right] \right\}\_0 + \mathbf{G}\_{0'}$$

where *G*0 is the free energy of the reference system, *E*0 is the energy at a point in the phase space of the reference system, and <sup>0</sup> is an ensemble average for the reference system. From Eq. (20), we obtain Zwanzig's free-energy perturbation (Zwanzig 1954):

$$G - G\_0 = -k\_B T \ln \left\langle \exp \left[ -\mathcal{J} \left( E - E\_0 \right) \right] \right\rangle\_0. \tag{21}$$

As a result, by taking the advantage of the perturbation [Eq. (21)], we can readily have enough samples in higher-energy regions in a reference frame where their original high potential energy values intentionally get lowered. Afterwards, the corrected free energy can straightforwardly be recovered by taking the average of the exponential factor exp *E E*<sup>0</sup> over the ensembles sampled in the reference system. This is exactly the idea behind many enhanced sampling methods, such as the umbrella sampling technique.

#### **4. Systematic** *ab initio* **molecular dynamics approach: Free-energy expansion method as a series of covariance tensors**

A fundamental key to have successful molecular simulations is the accuracy of internuclear potential for describing atomic motions during biochemical or physical events. By exploiting Zwanzig's free-energy perturbation (FEP) theory, we are developing a new rigorous method to systematically obtain accurate free-energy profiles, in which the internuclear potential energy is effectively computed at a high-level *ab initio* theory. Our new method is a systematic free-energy expansion (FEE) in terms of a series of covariance tensors. The new expansion will enable us to have a free-energy profile at a level as high as the coupled cluster theory at an affordable computational cost, which is currently known as the gold standard but unreachable level of theory for free-energy simulations. The focus of our FEE method will be on the difference of free energy calculated by two different internuclear potential. Furthermore, in contrast to Car-Parrinello MD (CPMD) which is limited to potential energy derived from DFT (Car and Parrinello 1985), our method is independent of how the potential energy functions being constructed. Therefore, by combining it with our novel automated integration-free path-integral (AIF-PI) method together (See Section 5; Wong and Gao 2007; Wong and Gao 2008; Wong, Richard et al. 2009; Wong, Gu et al. 2012), we will also be able to compute free-energy barriers, changes of binding energy, p*K*a values, and isotope effects at an *ab initio* path integral level (see Section 6).

Let's begin with the FEP theory. From Eq. (21), the free energy difference between using lower-level (LL) and higher-level (HL) *ab initio* methods can be expressed as:

$$\mathbf{G}\_{HL} - \mathbf{G}\_{LL} = \Delta \mathbf{G} = -k\_B T \ln \left\langle \exp \left[ -\mathcal{J} \left( E\_{HL} - E\_{LL} \right) \right] \right\rangle\_{LL} \tag{22}$$

Next we expand the ensemble average in Eq. (22) and sum up the prefactors into a series of cumulants:

$$\left\|\mathbf{G}\_{HL} - \mathbf{G}\_{LL} = \Delta \mathbf{G} = -\frac{1}{\mathcal{\beta}} \ln \left\langle \exp \left[ -\boldsymbol{\mathcal{\beta}} \,\Delta \mathbf{E} \right] \right\rangle\_{LL} = -\frac{1}{\mathcal{\beta}} \ln \left\langle \exp \left[ \sum\_{n=1}^{x} \left\langle \left( -\boldsymbol{\mathcal{\beta}} \,\Delta \mathbf{E} \right)^{n} \right\rangle\_{LL,\mathcal{\varepsilon}} \right] \right\rangle\_{\mathcal{\varepsilon}} \tag{23}$$

where

114 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy

Owing to the Boltzmann exponential energy distribution, one of the major difficulties in computing a converged free-energy profile or potential of mean force [Eq. (10)] via MD and MC sampling techniques is that it takes longer simulation time or runs more MC steps to have enough higher-energy samples. Yet, many interesting biochemical or physical molecular properties could be in higher-energy regions, e.g., the transition state during

In practice, in order for having effective samplings on both the lower-energy (e.g., reactant state) and higher-energy regions (e.g., transition state), Zwanzig's free-energy perturbation (Zwanzig 1954) [which is also referred to as statistical-mechanical perturbation theory (McQuarrie 2000)] has been extensively applied. The feature of the perturbation is relating the change of free energy between *two* systems (both have the same number of degrees of freedom) by an ensemble average taken in only *one* of the two systems. This can be illustrated by first writing the classical free energy *Gcl* corresponding to the partition

*cl B cl N*

where *E* is the energy at a point *pj j* ,*x* in the phase space, i.e.,

ln exp ,

*k T EE G*

3

*h*

*N*

*n*

From Eq. (20), we obtain Zwanzig's free-energy perturbation (Zwanzig 1954):

exp

*dx dp EE e*

 

ln ln

*cl B N N B N*

*<sup>h</sup> dx dp G kT k T <sup>e</sup>*

*e*

0 0 <sup>0</sup>

3 3

*N N*

*n n*

*dx dp G kT Q <sup>E</sup>*

*<sup>p</sup> EEp x V x*

3 3 <sup>3</sup> ln exp , *n n n N N*

*h*

 3 2 , . <sup>2</sup> *Nn <sup>j</sup> j j j j j*

0

*E*

where *G*0 is the free energy of the reference system, *E*0 is the energy at a point in the phase space of the reference system, and <sup>0</sup> is an ensemble average for the reference system.

0 0 <sup>0</sup> ln exp . *G G kT <sup>B</sup>*

As a result, by taking the advantage of the perturbation [Eq. (21)], we can readily have enough samples in higher-energy regions in a reference frame where their original high potential energy values intentionally get lowered. Afterwards, the corrected free energy can straightforwardly be recovered by taking the average of the exponential factor

 *E E*<sup>0</sup> over the ensembles sampled in the reference system. This is exactly the idea behind many enhanced sampling methods, such as the umbrella sampling technique.

 

*M*

0

*E N N N*

*n n n n n n*

(20)

3 0 3 3 3 3 3

*dx dp h*

(18)

(19)

*E E* (21)

0

*E*

**3. Zwanzig's free-energy perturbation theory** 

protein folding or biochemical reaction.

function in Eq. (8) as follows:

Next, we rewrite Eq. (18) as:

*B*

exp 

$$
\Delta E = E\_{\rm HL} - E\_{\rm LL} \tag{24}
$$

 *LL c*, is a cumulant, and *n* is the order of a cumulant. In his original 1954 paper (Zwanzig 1954), Zwanzig showed that the cumulant expansion is fast converging when the change of energy Δ*E* in the ensemble is reasonably small relative to the inverse of . However, in terms of computational cost, this cumulant expansion does not provide an advantage for correcting lower-level free energy. This is because the time required for calculating the cumulant average *LL c*, with computer is basically as much as the time needed to directly compute the higher-level free energy *GHL* , regardless of whether the perturbation Δ*E* is big.

In order to ease up this situation, in a forthcoming paper we will prove that each cumulant can be further expanded as a Taylor series expansion fluctuating about the ensemble average position **x***LL* in the form:

$$\left\langle \left( -\mathcal{J}\,\Delta E \right)^{n} \right\rangle\_{\rm LL} = \left\langle f\_{n}\left( \mathbf{x} \right) \right\rangle\_{\rm LL} = f\_{n}\left( \overline{\mathbf{x}}\_{\rm LL} \right) + \frac{1}{2!} \left[ \hat{\mathbf{D}}^{2} f\_{n}\left( \overline{\mathbf{x}}\_{\rm LL} \right) \right] \text{cov}\left( \mathbf{x}^{\rm T\_{r}}, \mathbf{x} \right) + \dotsb \tag{25}$$

Developing a Systematic Approach for *Ab Initio* Path-Integral Simulations 117

The preliminary results using this new systematic FEE method, i.e., Eq. (23), are very encouraging. Table 1 shows the free energy correction *G* for a single water molecule from HF/6-31G(d) to MP2/6-311G(d,p). Even just up to the first cumulant at the zeroth order correction, the computed error is in the order of magnitude ~0.1 kcal/mol. The first

All the above discussions on simulating internuclear thermodynamics are limited to classical mechanics (regardless of using QM, MM, hybrid QM/MM to construct potential energy). However, the real world is described by quantum mechanics, including nuclei. In some important applications of Life and Materials Sciences, such as hydrogen adsorption in carbon nanotechnology, the transport mechanism of hydrated hydroxide ions in aqueous solution, and kinetic isotope effects on a proton-transfer reaction, actually internuclear quantum-statistical effects (e.g., quantization of vibration and quantum tunneling) are not negligible. A popular choice for incorporating such internuclear quantum-statistical effects in the conventional molecular dynamics (MD) or Monte Carlo (MC) simulations (Tanaka, Kanoh et al. 2005; Warshel, Olsson et al. 2006; Kowalczyk, Gauden et al. 2007; Gao, Major et al. 2008; Kowalczyk, Gauden et al. 2008; Major, Heroux et al. 2009; Wong, Gu et al. 2012) is using Feynman's path integral (Feynman 1948; Feynman 1966; Kleinert 2004; Brown 2005;

The essence of Feynman's path integral is to transform the Schrödinger *differential* equation to become an *integral* equation. As a result, the many-body path integrations can be carried out by the conventional MD or MC sampling techniques. In addition, the quantum canonical partition function can be directly obtained with no need to compute individual

**5.1 Kleinert's variational perturbation theory for centroid density of path integrals** 

simulations which have been widely used in condensed phases.

Kleinert's variational perturbation (KP) theory (Kleinert 2004) for the centroid density (Gillan 1987; Gillan 1987; Voth 1996; Ramírez, López-Ciudad et al. 1998; Ramírez and López-Ciudad 1999; Feynman, Hibbs et al. 2005) of Feynman path integrals (Feynman 1948; Feynman 1966; Kleinert 2004; Brown 2005; Feynman, Hibbs et al. 2005) provides a complete theoretical foundation for developing non-stochastic methods to systematically incorporate internuclear quantum-statistical effects in condensed phase systems. Similar to the complementary interplay between the rapidly growing quantum Monte Carlo simulations (Anderson 1975; Grossman and Mitas 2005; Lester and Salomon-Ferrer 2006; Wagner, Bajdich et al. 2009) and the well-established *ab initio* or density-functional theories (DFT) for electronic structure calculations (Hehre, Radom et al. 1986; Szabo and Ostlund 1996; Kohn 1999; Pople 1999; Helgaker, Jørgensen et al. 2000; Springborg 2000), non-stochastic pathintegral methods can complement the conventional Fourier or discretized path-integral Monte-Carlo (PIMC) (MacKeown 1985; Coalson 1986; Ceperley 1995; Mielke and Truhlar 2001; Sauer 2001) and molecular dynamics (PIMD) (Cao and Voth 1994; Voth 1996)

To simplify the illustration of the essence of Kleinert's variational perturbation theory, we now consider a one-particle one-dimensional system. For a one-particle one-dimensional

system, the classical canonical partition function in Eq. (8) reduces to become:

cumulant is basically converged as soon as the second order correction is included.

**5. Simulating quantum thermodynamics: Feynman's path integral** 

Feynman, Hibbs et al. 2005).

energy eigenvalues.

where *Tp* is transpose, **<sup>x</sup>***<sup>n</sup> nf E* , **x** is a position vector of 3*N* Cartesian coordinates of the system, **D**ˆ *<sup>n</sup>* is the *n*th-order tensor operator for differentiation with respect to the 3*N* coordinates (e.g., **D x** <sup>ˆ</sup> *n LL <sup>f</sup>* is the gradient and **<sup>2</sup> D x** <sup>ˆ</sup> *n LL <sup>f</sup>* is the Hessian matrix), and cov , **x x** *Tp* is the covariance matrix. The higher order terms in Eq. (25) involve higher order covariance tensors. Note that the term associated with the gradient is not shown in Eq. (25) because the first order central moment, i.e., **x x***LL LL* , is always zero by definition.

By combining Eq. (25) with Eq. (23), we have enough equations to systematically approach the exact value of high-level free energy at a reduced computational cost. The number of calculations involving *EHL* is now considerably decreased to only a single-point energy calculation at **x***LL* for the zeroth order correction, and merely a normal-mode frequency analysis at **x***LL* for the second order correction.

To increase the converging property for the expansion in Eq. (25) as well as to overcome the problem of multi-model probability distribution, we can further generalize the FEE method by considering a decomposition of the ensemble average into subgroups by clustering methods. The clustering scheme will be determined in a way such that the FEE expansion is converged up to the second order correction in each group or each cluster. Please note that, in the limit that the number of clusters becomes as many as the number of ensembles, the formalism reduces back to the original ensemble average, and inclusion of only the zeroth order term in Eq. (25) is able to return us back the exact result of Eq. (23).


Table 1. Free-Energy correction *G* for H2O from HF/6-31G(d) to MP2/6-311G(d,p).

Since single-point energy calculations and a normal-mode frequency analysis at high-level electronic structure calculations are actually very common in literature (which are often used for minimized structures, though), we anticipate this new free-energy expansion method would be particularly useful for coupling accurate results from high-level *ab initio* theory with computational efficiency of lower-level samplings in free-energy calculations.

*Ef f f* (25)

*E* , **x** is a position vector of 3*N* Cartesian coordinates

<sup>ˆ</sup> *n LL <sup>f</sup>* is the Hessian matrix), and

*G* (kcal/mol) Error

 **2 x x D x xx** <sup>1</sup> <sup>ˆ</sup> cov , 2! *n Tp <sup>n</sup> n LL n LL LL LL*

of the system, **D**ˆ *<sup>n</sup>* is the *n*th-order tensor operator for differentiation with respect to the 3*N*

cov , **x x** *Tp* is the covariance matrix. The higher order terms in Eq. (25) involve higher order covariance tensors. Note that the term associated with the gradient is not shown in Eq. (25)

By combining Eq. (25) with Eq. (23), we have enough equations to systematically approach the exact value of high-level free energy at a reduced computational cost. The number of calculations involving *EHL* is now considerably decreased to only a single-point energy calculation at **x***LL* for the zeroth order correction, and merely a normal-mode frequency

To increase the converging property for the expansion in Eq. (25) as well as to overcome the problem of multi-model probability distribution, we can further generalize the FEE method by considering a decomposition of the ensemble average into subgroups by clustering methods. The clustering scheme will be determined in a way such that the FEE expansion is converged up to the second order correction in each group or each cluster. Please note that, in the limit that the number of clusters becomes as many as the number of ensembles, the formalism reduces back to the original ensemble average, and inclusion of only the zeroth

6th ∞ −170.917 0.000

Since single-point energy calculations and a normal-mode frequency analysis at high-level electronic structure calculations are actually very common in literature (which are often used for minimized structures, though), we anticipate this new free-energy expansion method would be particularly useful for coupling accurate results from high-level *ab initio* theory with computational efficiency of lower-level samplings in free-energy calculations.

Table 1. Free-Energy correction *G* for H2O from HF/6-31G(d) to MP2/6-311G(d,p).

0th −170.601 0.316 1st −170.601 0.316 2nd −170.680 0.237 3rd −170.680 0.237 ∞ −170.680 0.237

0th −170.601 0.316 1st −170.601 0.316 2nd −170.875 0.042 3rd −170.877 0.040 ∞ −170.883 0.034

because the first order central moment, i.e., **x x***LL LL* , is always zero by definition.

order term in Eq. (25) is able to return us back the exact result of Eq. (23).

*nf* 

coordinates (e.g., **D x** <sup>ˆ</sup> *n LL <sup>f</sup>* is the gradient and **<sup>2</sup> D x**

where *Tp* is transpose, **<sup>x</sup>***<sup>n</sup>*

analysis at **x***LL* for the second order correction.

Cumulant Tensor

1st

2nd

The preliminary results using this new systematic FEE method, i.e., Eq. (23), are very encouraging. Table 1 shows the free energy correction *G* for a single water molecule from HF/6-31G(d) to MP2/6-311G(d,p). Even just up to the first cumulant at the zeroth order correction, the computed error is in the order of magnitude ~0.1 kcal/mol. The first cumulant is basically converged as soon as the second order correction is included.
