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

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; Feynman, Hibbs et al. 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 energy eigenvalues.

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

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) simulations which have been widely used in condensed phases.

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:

$$Q\_{cl} = \sqrt{\frac{\mathcal{M}k\_B T}{2\pi\hbar^2}} \int\_{-\infty}^{\infty} e^{-\beta V(x\_0)} d\mathbf{x}\_0. \tag{26}$$

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

A number of non-stochastic approaches have been developed to approximately estimate the centroid potential. For example, Feynman and Hibbs described a first-order cumulant expansion by introducing a Gaussian smearing function in a free-particle reference frame to yield an upper bound on the centroid potential (Feynman, Hibbs et al. 2005). This was subsequently modified by Doll and Myers (DM) by using a Gaussian width associated with the angular frequency at the minimum of the original potential (Doll and Myers 1979). Mielke and Truhlar employed a free-particle reference state and approximated the sum over paths by a minimal set of paths constrained for a harmonic oscillator. The action integral is obtained by using the three-point trapezoidal rule for the potential to yield the displaced-

A closely related theoretical approach to the KP theory is the variational method independently introduced by Giachetti and Tognetti (Giachetti and Tognetti 1985), and by Feynman and Kleinert (hereafter labeled as GTFK) (Feynman and Kleinert 1986), which formally corresponds to the first order approximation in the KP theory, i.e., KP1. The GTFK approach is a variational method that adopts a harmonic reference state by variationally optimizing the angular frequency. This variational method has been applied to a variety of systems, including quantum dynamic processes in condensed phases (e.g., water and helium). Although the original GTFK approach is among the most accurate approximate methods for estimating the path-integral centroid potential in many applications (Mielke and Truhlar 2001), significant errors can exist in situations in which quantum effects are dominant, especially at low temperatures. Higher order perturbations of KP theory can significantly and systematically improve computational accuracy over the KP1 results.(Kleinert 2004; Wong and Gao 2007; Wong and Gao 2008; Wong, Richard et al. 2009;

In essence, what Kleinert's variational perturbation (KP) theory does is to systematically builds up anharmonic corrections to the harmonic centroid potential calculated in a harmonic reference state characterized by a trial angular frequency Ω (Kleinert 2004). Given

**.** <sup>0</sup> <sup>2</sup> <sup>2</sup> <sup>2</sup>

the centroid potential *W x* <sup>0</sup> in Eq. (28) can be expressed as a path integral of the harmonic

<sup>0</sup> 0 0 <sup>0</sup> <sup>0</sup> <sup>0</sup>

0 <sup>2</sup> / 2 , sinh / 2

<sup>2</sup> ,

2 2 *<sup>x</sup> <sup>M</sup> d x Mx x*

0

*W x x*

*<sup>e</sup> x xxe e Q e Mk T*

<sup>0</sup> <sup>0</sup>

 

and 0 *<sup>x</sup>* is the expectation value over all closed paths of the action in Eq. (31):

*<sup>x</sup> x*

action which is perturbed by the anharmonicity of the original potential:

 

where 0 *<sup>x</sup> Q* is the local harmonic partition function given as follows:

*Q x xxe*

2

*B*

*Mk T* 

1

0

*x x <sup>x</sup>*

 

 <sup>D</sup> *<sup>A</sup>* (33)

<sup>D</sup> *<sup>A</sup> A A A A* (32)

*x*

 

*A* (31)

point path integral (DPPI) centroid potential (Mielke and Truhlar 2001).

Wong, Gu et al. 2012)

the reference, or trial harmonic action:

0

2

*B*

The traditional way to obtain the quantum canonical partition function, i.e., Eq. (3), is to solve the internuclear Schrödinger equation to get the individual energy eigenvalues. But in the path-integral (PI) formulation, we do not know the individual energy eigenvalues for obtaining the quantum partition function. This is because the PI representation of the quantum partition function can be written in terms of the centroid effective potential *W* as a classical configuration integral (Gillan 1987; Gillan 1987; Voth 1996; Ramírez, López-Ciudad et al. 1998; Ramírez and López-Ciudad 1999; Kleinert 2004; Feynman, Hibbs et al. 2005):

$$Q\_{qm} = \sum\_{n} \exp\left(-\beta E\_n\right) = \sqrt{\frac{Mk\_BT}{2\pi\hbar^2}} \int\_{-\infty}^{\infty} e^{-\beta W(\mathbf{x}\_0)} d\mathbf{x}\_0. \tag{27}$$

Given the centroid potential *W x* <sup>0</sup> , thermodynamic and quantum dynamic quantities can be accurately determined, including molecular spectroscopy of quantum fluids and the rate constant of chemical and enzymatic reactions. The mass-dependent nature of *W x* <sup>0</sup> is also of particular interest because isotope effects can be obtained, and it has been applied to carbon nanotubes (Tanaka, Kanoh et al. 2005; Kowalczyk, Gauden et al. 2007; Kowalczyk, Gauden et al. 2008), and biochemical reactions in protein (Warshel, Olsson et al. 2006; Gao, Major et al. 2008; Major, Heroux et al. 2009) and RNA enzymes (Wong, Gu et al. 2012).

The centroid potential *W x* <sup>0</sup> in Eq. (27) is defined as follows (Gillan 1987; Gillan 1987; Voth 1996; Ramírez, López-Ciudad et al. 1998; Ramírez and López-Ciudad 1999; Kleinert 2004; Feynman, Hibbs et al. 2005):

$$\mathcal{W}\left(\mathbf{x}\_{0}\right) = -k\_{B}T\ln\left[\sqrt{\frac{2\pi\hbar^{2}}{Mk\_{B}T}}\oint \mathcal{D}\left[\mathbf{x}\left(\tau\right)\right]\delta\left(\overline{\mathbf{x}} - \mathbf{x}\_{0}\right)\exp\left\{-\mathcal{A}\left[\mathbf{x}\left(\tau\right)\right]\right\}
\hbar\right\}\right],\tag{28}$$

where *τ* is a real number and represents the component for pure imaginary time in path integral, *x* describes a path in space-time, *x xx* <sup>0</sup> <sup>D</sup> denotes a summation over *all* possible closed paths in which *x* is equal to 0 *x* (i.e., a functional integration), and *x* is the time-average position, called 'centroid'

$$\overline{\mathfrak{X}} \equiv \frac{1}{\beta \hbar} \int\_0^{\beta \hbar} \mathfrak{x}(\tau) d\tau \,. \tag{29}$$

In Eq. (28), *A* is the quantum-statistical action:

$$\mathcal{A}\left[\mathbf{x}\left(\tau\right)\right] = \int\_0^{\mu\_0} d\tau \left\{\frac{M}{2}\dot{\mathbf{x}}\left(\tau\right)^2 + V\left[\mathbf{x}\left(\tau\right)\right]\right\},\tag{30}$$

where *V x* is the original potential energy of the system. Generalization of Eq. (28) to a multi-dimensional system is straightforward (Kleinert 2004; Feynman, Hibbs et al. 2005).

*Mk T <sup>Q</sup> e dx*

*cl*

*qm n n*

2004; Feynman, Hibbs et al. 2005):

the time-average position, called 'centroid'

In Eq. (28), *A* is the quantum-statistical action:

integral, *x*

 <sup>0</sup> <sup>2</sup> 0. <sup>2</sup> *B V x*

(26)

 <sup>0</sup> <sup>2</sup> <sup>0</sup> exp . <sup>2</sup>

*Mk T <sup>Q</sup> <sup>E</sup> e dx*

Given the centroid potential *W x* <sup>0</sup> , thermodynamic and quantum dynamic quantities can be accurately determined, including molecular spectroscopy of quantum fluids and the rate constant of chemical and enzymatic reactions. The mass-dependent nature of *W x* <sup>0</sup> is also of particular interest because isotope effects can be obtained, and it has been applied to carbon nanotubes (Tanaka, Kanoh et al. 2005; Kowalczyk, Gauden et al. 2007; Kowalczyk, Gauden et al. 2008), and biochemical reactions in protein (Warshel, Olsson et al. 2006; Gao, Major et al. 2008; Major, Heroux et al. 2009) and RNA enzymes (Wong, Gu et al. 2012).

The centroid potential *W x* <sup>0</sup> in Eq. (27) is defined as follows (Gillan 1987; Gillan 1987; Voth 1996; Ramírez, López-Ciudad et al. 1998; Ramírez and López-Ciudad 1999; Kleinert

 

where *τ* is a real number and represents the component for pure imaginary time in path

*all* possible closed paths in which *x* is equal to 0 *x* (i.e., a functional integration), and *x* is

0 <sup>1</sup> *x xd* 

 

<sup>0</sup> 2 *<sup>M</sup> x d x Vx*

 

<sup>2</sup> ln exp , *<sup>B</sup>*

 <sup>D</sup> *<sup>A</sup>* (28)

> 

**.**

 

*A* (30)

 

**,** <sup>2</sup>

 

where *V x* is the original potential energy of the system. Generalization of Eq. (28) to a multi-dimensional system is straightforward (Kleinert 2004; Feynman, Hibbs et al. 2005).

2 0 0

*W x kT x xx x*

*B*

describes a path in space-time, *x xx*

*Mk T* 

*B W x*

(27)

<sup>0</sup> <sup>D</sup> denotes a summation over

(29)

The traditional way to obtain the quantum canonical partition function, i.e., Eq. (3), is to solve the internuclear Schrödinger equation to get the individual energy eigenvalues. But in the path-integral (PI) formulation, we do not know the individual energy eigenvalues for obtaining the quantum partition function. This is because the PI representation of the quantum partition function can be written in terms of the centroid effective potential *W* as a classical configuration integral (Gillan 1987; Gillan 1987; Voth 1996; Ramírez, López-Ciudad et al. 1998; Ramírez and López-Ciudad 1999; Kleinert 2004; Feynman, Hibbs et al. 2005):

A number of non-stochastic approaches have been developed to approximately estimate the centroid potential. For example, Feynman and Hibbs described a first-order cumulant expansion by introducing a Gaussian smearing function in a free-particle reference frame to yield an upper bound on the centroid potential (Feynman, Hibbs et al. 2005). This was subsequently modified by Doll and Myers (DM) by using a Gaussian width associated with the angular frequency at the minimum of the original potential (Doll and Myers 1979). Mielke and Truhlar employed a free-particle reference state and approximated the sum over paths by a minimal set of paths constrained for a harmonic oscillator. The action integral is obtained by using the three-point trapezoidal rule for the potential to yield the displacedpoint path integral (DPPI) centroid potential (Mielke and Truhlar 2001).

A closely related theoretical approach to the KP theory is the variational method independently introduced by Giachetti and Tognetti (Giachetti and Tognetti 1985), and by Feynman and Kleinert (hereafter labeled as GTFK) (Feynman and Kleinert 1986), which formally corresponds to the first order approximation in the KP theory, i.e., KP1. The GTFK approach is a variational method that adopts a harmonic reference state by variationally optimizing the angular frequency. This variational method has been applied to a variety of systems, including quantum dynamic processes in condensed phases (e.g., water and helium). Although the original GTFK approach is among the most accurate approximate methods for estimating the path-integral centroid potential in many applications (Mielke and Truhlar 2001), significant errors can exist in situations in which quantum effects are dominant, especially at low temperatures. Higher order perturbations of KP theory can significantly and systematically improve computational accuracy over the KP1 results.(Kleinert 2004; Wong and Gao 2007; Wong and Gao 2008; Wong, Richard et al. 2009; Wong, Gu et al. 2012)

In essence, what Kleinert's variational perturbation (KP) theory does is to systematically builds up anharmonic corrections to the harmonic centroid potential calculated in a harmonic reference state characterized by a trial angular frequency Ω (Kleinert 2004). Given the reference, or trial harmonic action:

$$\mathbf{A}\_{\Omega}^{x\_0} = \int\_0^{\beta \hbar} d\tau \left\{ \frac{M}{2} \dot{\mathbf{x}} \left( \tau \right)^2 + \frac{1}{2} M \Omega^2 \left[ \mathbf{x} \left( \tau \right) - \mathbf{x}\_0 \right]^2 \right\}. \tag{31}$$

the centroid potential *W x* <sup>0</sup> in Eq. (28) can be expressed as a path integral of the harmonic action which is perturbed by the anharmonicity of the original potential:

$$e^{-\beta \mathcal{W}(\mathbf{x}\_0)} = \sqrt{\frac{2\pi \hbar^2}{Mk\_B T}} \oint \mathcal{D}\mathbf{x}\left(\tau\right) \mathcal{S}\left(\overline{\mathbf{x}} - \mathbf{x}\_0\right) e^{-\mathcal{A}\_{\Omega}^{\text{up}} \hbar \frac{1}{\epsilon} \left(\mathsf{A} - \mathcal{A}\_{\Omega}^{\text{up}}\right) \hbar} = Q\_{\Omega}^{\text{x\_0}} \left\langle e^{-\left(\mathsf{A} - \mathcal{A}\_{\Omega}^{\text{up}}\right) \hbar \right\rangle^{\text{x\_0}}}\right\rangle\_{\Omega},\tag{32}$$

where 0 *<sup>x</sup> Q* is the local harmonic partition function given as follows:

$$Q^{x\_0}\_{\Omega} = \sqrt{\frac{2\pi\hbar^2}{Mk\_BT}} \oint \text{D}x\left(\tau\right) \delta\left(\overline{\mathbf{x}} - \mathbf{x}\_0\right) e^{-\mathcal{A}^{x\_0}\_{\Omega}/\hbar} = \frac{\beta\hbar\Omega\left/2}{\sinh\left(\beta\hbar\Omega/2\right)},\tag{33}$$

and 0 *<sup>x</sup>* is the expectation value over all closed paths of the action in Eq. (31):

$$\left\langle e^{-F\left[\mathbf{x}(\tau)\right]} \right\rangle\_{\Omega}^{\hbar} = \frac{1}{Q\_{\Omega}^{x\_0}} \sqrt{\frac{2\pi\hbar^2}{M\mathbf{k}\_B T}} \oint \text{D}\mathbf{x}\left(\tau\right) \delta\left(\overline{\mathbf{x}} - \mathbf{x}\_0\right) e^{-F\left[\mathbf{x}(\tau)\right] \hbar} e^{-\mathbf{A}\_{\Omega}^{\alpha} \hbar} \,,\tag{34}$$

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

is an element of the inverse matrix of <sup>2</sup>

0

*k T k T kT Q d V x d dVx V x*

int 1 1 0 ,

 

<sup>1</sup> , !

approach to determine the optimal value of Ω, **opt,***<sup>n</sup> x*<sup>0</sup> (Kleinert 2004).

*B jk <sup>n</sup> j k <sup>c</sup>*

*<sup>n</sup> <sup>x</sup> n n <sup>x</sup>*

 

2

now be written in terms of *ordinary* integrals as follows (Kleinert 2004):

*M*

*k T d Vx*

1

where <sup>0</sup> <sup>2</sup> <sup>2</sup> int 0 1 2 *<sup>x</sup> V x Vx M x x*

Jensen-Peierls inequality, i.e., from Eq. (32) and (35):

Kleinert perturbation theory such that **opt,** <sup>0</sup>

*n*

Eq. (31) cancel each other out).

is the determinant of the *n n* **-matrix** consisting of the Gaussian

<sup>1</sup> cosh / 2 1 . 2 sinh / 2

<sup>0</sup> <sup>0</sup> 0 0 0 0

*x x x x B B x x <sup>B</sup> <sup>c</sup>*

0

(the kinetic energy terms in Eq. (30) and

int <sup>1</sup> <sup>2</sup> 12 1 2 int int , <sup>0</sup> 0 0

 

 

After using these smearing potentials given in Eq. (39), the *n*th-order Kleinert variational perturbation (KP*n*) approximation, *W x <sup>n</sup>* <sup>0</sup> , shown in Eq. (35) as *functional* integrals, can

2!

 

As *n* tends to infinity, *W x <sup>n</sup>* <sup>0</sup> approaches the exact value of the centroid potential *W x* <sup>0</sup> in Eq. (28), which is independent of the trial Ω. But the truncated sum in Eq. (41) does depend on Ω, and the optimal choice of this trial frequency at a given order of KP expansion and at a particular centroid position 0 *<sup>x</sup>* is determined by the least-dependence of <sup>0</sup> *<sup>x</sup> Wn* on Ω itself. This is the so-called frequency of least dependence, which provides a variational

Of particular interest is the special case when 1 *n* , which turns out to be identical to the original GTFK variational approach. An important property of KP1 or the GTFK variational approach is that there is a definite upper bound for the computed *W x* 1 0 by virtue of the

> <sup>0</sup> <sup>0</sup> 0 0 <sup>0</sup> 0 0 1 0 exp exp . *<sup>x</sup> <sup>x</sup> x x W x x x W x e Q <sup>Q</sup> <sup>e</sup>*

0

*<sup>x</sup> W x n n* is least-dependent on Ω.

*A A A A* (42)

 

Note that by choosing 0 (i.e., the reference state is for a free particle), KP1 or GTFK (Giachetti and Tognetti 1985; Feynman and Kleinert 1986) reduces to the Feynman-Hibbs approach (Feynman, Hibbs et al. 2005). For higher orders of *n*, unfortunately, it is not guaranteed that a minimum of <sup>0</sup> *<sup>x</sup> Wn* actually exists as a function of Ω. In this case, the least dependent Ω is obtained from the condition that the next derivative of <sup>0</sup> *<sup>x</sup> Wn* with respect to Ω is set to zero. Consequently, Ω is considered as a variational parameter in the

  *k k a* 

> 

, and the Gaussian

  (41)

(40)

where **Det** <sup>2</sup> *k k a* 

0

ln

( )

*n*

*W x*

 , <sup>2</sup> *k k a* 

width is a function of the trial frequency Ω:

2

*a*

width <sup>2</sup> *k k a* 

In Eq. (34), *F x* denotes an arbitrary functional. It is of interest to note that Eq. (32) is similar to the starting point of Zwanzig's free-energy perturbation (Section 3), which has been extensively used in free-energy calculations through Monte Carlo and molecular dynamics simulations. Their difference is one is for ordinary ensemble average, while another one is for closed-path average, i.e., *functional* average.

If we expand the exponential functional in Eq. (32) and sum up the prefactors into an exponential series of cumulants, then the *n*th-order approximation, *W x <sup>n</sup>* <sup>0</sup> , to the centroid potential *W x* <sup>0</sup> can be written as follows (Kleinert 2004):

$$\begin{split} e^{-\beta W\_n^{\Omega}(\mathbf{x}\_0)} = & \mathbf{Q}\_{\Omega}^{\mathbf{x}\_0} \exp\left\{-\frac{1}{\hbar} \int\_0^{\beta \hbar} d\tau \left\langle \mathbf{A}\_{\mathrm{int}}^{\mathbf{x}\_0} \right\rangle\_{\Omega,\boldsymbol{\varepsilon}}^{\times\_0} + \frac{1}{2!\hbar^2} \int\_0^{\beta \hbar} d\tau\_1 \int\_0^{\beta \hbar} d\tau\_2 \left\langle \mathbf{A}\_{\mathrm{int}}^{\mathbf{x}\_0} \right[ \mathbf{x}(\tau\_1) \Big] \mathbf{A}\_{\mathrm{int}}^{\mathbf{x}\_0} \left[ \mathbf{x}(\tau\_2) \right] \right\rangle\_{\Omega,\boldsymbol{\varepsilon}}^{\times\_0} \\ & + \dots + \left\{ \prod\_{j=1}^n \int\_0^{\beta \hbar} d\tau\_j \right\} \frac{(-1)^n}{n!\hbar^n} \left\langle \prod\_{k=1}^n \mathbf{A}\_{\mathrm{int}}^{\mathbf{x}\_0} \right[ \mathbf{x}(\tau\_k) \Big] \right\rangle\_{\Omega,\boldsymbol{\varepsilon}}^{\times\_0}, \end{split} \tag{35}$$

where 0 0 int *x x A AA* is the so-called inter-*action*, representing the perturbation to the harmonic reference state, 0 , *<sup>x</sup> <sup>c</sup>* is a cumulant which can be written in terms of expectation values 0 *<sup>x</sup>* by the cumulant expansion (Zwanzig 1954; Kubo 1962; Kleinert 2004), e.g.,

$$\left\langle \mathcal{A}\_{\text{int}}^{x\_0} \left[ \boldsymbol{\chi} \left( \boldsymbol{\tau} \right) \right] \right\rangle\_{\Omega, c}^{x\_0} \equiv \left\langle \mathcal{A}\_{\text{int}}^{x\_0} \left[ \boldsymbol{\chi} \left( \boldsymbol{\tau} \right) \right] \right\rangle\_{\Omega}^{x\_0} \,\tag{36}$$

$$
\left\langle \mathsf{A}\_{\mathrm{init}}^{\mathrm{x\_{0}}} \left[ \mathbbm{x} \left( \tau\_{1} \right) \right] \mathsf{A}\_{\mathrm{init}}^{\mathrm{x\_{0}}} \left[ \mathbbm{x} \left( \tau\_{2} \right) \right] \right\rangle\_{\Omega, \varepsilon}^{\mathrm{x\_{0}}} \equiv \left\langle \mathsf{A}\_{\mathrm{init}}^{\mathrm{x\_{0}}} \left[ \mathbbm{x} \left( \tau\_{1} \right) \right] \mathsf{A}\_{\mathrm{init}}^{\mathrm{x\_{0}}} \left[ \mathbbm{x} \left( \tau\_{2} \right) \right] \right\rangle\_{\Omega}^{\mathrm{x\_{0}}} - \left\langle \left\langle \mathsf{A}\_{\mathrm{init}}^{\mathrm{x\_{0}}} \left[ \mathbbm{x} \left( \tau \right) \right] \right\rangle\_{\Omega}^{\mathrm{x\_{0}}} \right\rangle,\tag{37}
$$

$$
\begin{split}
\left\langle \mathbf{A}\_{\text{int}}^{\text{x}\_{0}} \left[ \mathbf{x} \left( \tau\_{1} \right) \right] \mathbf{A}\_{\text{int}}^{\text{x}\_{0}} \left[ \mathbf{x} \left( \tau\_{2} \right) \right] \mathbf{A}\_{\text{int}}^{\text{x}\_{0}} \left[ \mathbf{x} \left( \tau\_{3} \right) \right] \right\rangle\_{\Omega, \varepsilon}^{\text{x}\_{0}} &= \left\langle \mathbf{A}\_{\text{int}}^{\text{x}\_{0}} \left[ \mathbf{x} \left( \tau\_{1} \right) \right] \mathbf{A}\_{\text{int}}^{\text{x}\_{0}} \left[ \mathbf{x} \left( \tau\_{2} \right) \right] \mathbf{A}\_{\text{int}}^{\text{x}\_{0}} \left[ \mathbf{x} \left( \tau\_{3} \right) \right] \right\rangle\_{\Omega}^{\text{x}\_{0}} \\ &- 3 \left\langle \mathbf{A}\_{\text{int}}^{\text{x}\_{0}} \left[ \mathbf{x} \left( \tau\_{1} \right) \right] \mathbf{A}\_{\text{int}}^{\text{x}\_{0}} \left[ \mathbf{x} \left( \tau\_{2} \right) \right] \right\rangle\_{\Omega}^{\text{x}\_{0}} \left\langle \mathbf{A}\_{\text{int}}^{\text{x}\_{0}} \left[ \mathbf{x} \left( \tau \right) \right] \right\rangle\_{\Omega}^{\text{x}\_{0}} + 2 \left\langle \left[ \mathbf{A}\_{\text{int}}^{\text{x}\_{0}} \left[ \mathbf{x} \left( \tau \right) \right] \right] \right\rangle\_{\Omega}^{\text{x}\_{0}} \right\rangle \end{split} \tag{38}
$$

More importantly, Kleinert and co-workers derived a math equation for expressing the expectation value <sup>0</sup> 1 1 0 *<sup>x</sup> n n j kk j k d Fx* from the *functional-integral* form to be in terms of Gaussian smearing convolution integrals (*ordinary* integrals) (Kleinert 2004):

$$\begin{aligned} \left\langle \prod\_{j=1}^{n} \prod\_{0}^{\theta \hbar} d\pi\_{j} \right\rangle \left\langle \prod\_{k=1}^{n} F\_{k} \left[ \mathbf{x} \left( \mathbf{r}\_{k} \right) \right] \right\rangle\_{\Omega}^{\mathbf{x}\_{0}} &= \left\langle \prod\_{j=1}^{n} \prod\_{0}^{\theta \hbar} d\pi\_{j} \right\rangle \left\langle \prod\_{k=1}^{n} \prod\_{-\alpha}^{\alpha} d\mathbf{x}\_{k} F\_{k} \left( \mathbf{x}\_{k} \right) \right\rangle \\ &\times \frac{1}{\sqrt{\left( 2\pi \right)^{n} \mathbf{D} \mathbf{t} \left[ \mathbf{a}\_{\tau\_{k}\tau\_{\mathbf{t}}}^{2} \left( \mathbf{Q} \right) \right]}} \exp\left\{ -\frac{1}{2} \sum\_{k=1}^{n} (\mathbf{x}\_{k} - \mathbf{x}\_{0}) a\_{\tau\_{k}\tau\_{\mathbf{t}}}^{-2} \left( \mathbf{Q} \right) (\mathbf{x}\_{k'} - \mathbf{x}\_{0}) \right\}, \end{aligned} \tag{39}$$

2

*e x xxe e*

<sup>0</sup> <sup>0</sup>

 

similar to the starting point of Zwanzig's free-energy perturbation (Section 3), which has been extensively used in free-energy calculations through Monte Carlo and molecular dynamics simulations. Their difference is one is for ordinary ensemble average, while

If we expand the exponential functional in Eq. (32) and sum up the prefactors into an exponential series of cumulants, then the *n*th-order approximation, *W x <sup>n</sup>* <sup>0</sup> , to the centroid

<sup>0</sup> <sup>0</sup> <sup>0</sup> <sup>0</sup> <sup>0</sup> 0 0

*<sup>n</sup> <sup>x</sup> n n <sup>x</sup> j k <sup>n</sup> j k <sup>c</sup>*

*d x*

 

*x x W x x x x x*

1 1 0 ,

*x x A AA* is the so-called inter-*action*, representing the perturbation to the

*AA AA A* (37)

 

0

<sup>1</sup> , !

 

 

*<sup>x</sup> <sup>c</sup>* is a cumulant which can be written in terms of expectation

*A A* (36)

int

*A*

 

*A A A*

int <sup>2</sup> 12 1 2 int int , , <sup>0</sup> 0 0

 

denotes an arbitrary functional. It is of interest to note that Eq. (32) is

1 2 , *<sup>x</sup> <sup>x</sup> F x F x*

0

 

0

 

*c c*

<sup>D</sup> *<sup>A</sup>* (34)

 

2

 

0 0 <sup>2</sup> <sup>1</sup> 1

*n k k <sup>n</sup> <sup>k</sup> k*

 

 

from the *functional-integral* form to be in

2

 

*x xa x x*

3

(35)

(38)

,

(39)

0

another one is for closed-path average, i.e., *functional* average.

potential *W x* <sup>0</sup> can be written as follows (Kleinert 2004):

,

In Eq. (34), *F x*

( )

*n*

where 0 0 int

harmonic reference state, 0

expectation value <sup>0</sup>

11 11 0 0

*j k j k*

*<sup>x</sup> n n n n*

*<sup>x</sup> n n j kk*

 

*d Fx*

 

*j kk j kk k*

*d Fx d dx F x*

 

1 1 0

*j k*

0

*<sup>x</sup> <sup>B</sup>*

*eQ d dd x x*

*n*

values 0 *<sup>x</sup>* by the cumulant expansion (Zwanzig 1954; Kubo 1962; Kleinert 2004), e.g.,

 0 0 0 0 int int , , *x x x x <sup>c</sup> x x*

<sup>0</sup> 0 0 0 0 0 0 <sup>0</sup>

int 123 123 int int int int int ,

*xx x x*

*xxx xxx*

int 1 2 int int int 3 2

 

 

**Det**

terms of Gaussian smearing convolution integrals (*ordinary* integrals) (Kleinert 2004):

*a*

 

*k k*

1 1

exp <sup>2</sup> <sup>2</sup> *k k*

*AA A A*

*AAA AAA*

int 12 12 int int int int , , *<sup>x</sup> x x xx xx x <sup>c</sup> xx xx x*

<sup>0</sup> <sup>0</sup> <sup>000</sup> <sup>000</sup>

*x x xxx xxx c xx x xx x x*

00 0 00 0 <sup>0</sup>

More importantly, Kleinert and co-workers derived a math equation for expressing the

1 1 exp 2!

 

*Q Mk T*

where **Det** <sup>2</sup> *k k a* is the determinant of the *n n* **-matrix** consisting of the Gaussian width <sup>2</sup> *k k a* , <sup>2</sup> *k k a* is an element of the inverse matrix of <sup>2</sup> *k k a* , and the Gaussian width is a function of the trial frequency Ω:

$$a\_{rr'}^2(\Omega) = \frac{1}{\beta M \Omega^2} \left\{ \frac{\beta \hbar \Omega}{2} \frac{\cosh\left[ \left( \left| \tau - \tau' \right| - \beta \hbar \left< \ 2 \right.\right.\right) \Omega}{\sinh\left( \beta \hbar \Omega \left. \ / \ 2 \right.\right)} - 1 \right].\tag{40}$$

After using these smearing potentials given in Eq. (39), the *n*th-order Kleinert variational perturbation (KP*n*) approximation, *W x <sup>n</sup>* <sup>0</sup> , shown in Eq. (35) as *functional* integrals, can now be written in terms of *ordinary* integrals as follows (Kleinert 2004):

$$\begin{split} & \quad \mathcal{W}\_{\boldsymbol{\pi}}^{\Omega}(\boldsymbol{x}\_{0}) \\ &= -k\_{B}T \ln Q\_{\Omega}^{\boldsymbol{\pi}\_{0}} + \frac{k\_{B}T}{\hbar} \int\_{0}^{\beta\hbar} d\tau \Big\langle V\_{\mathrm{int}}^{\boldsymbol{x}\_{0}} \Big[\boldsymbol{x}\left(\boldsymbol{\tau}\_{1}\right)\Big] \Big\rangle\_{\Omega}^{\boldsymbol{x}\_{0}} - \frac{k\_{B}T}{2!\hbar^{2}} \int\_{0}^{\beta\hbar} d\tau\_{1} \int\_{0}^{\beta\hbar} d\tau\_{2} \Big\langle V\_{\mathrm{int}}^{\boldsymbol{x}\_{0}} \Big[\boldsymbol{x}\left(\boldsymbol{\tau}\_{1}\right)\Big] V\_{\mathrm{int}}^{\boldsymbol{x}\_{0}} \Big\rangle\_{\Omega}^{\boldsymbol{x}\_{0}} \\ & + \cdots + k\_{B}T \frac{\left(-1\right)^{n+1}}{n!\hbar^{n}} \Big\langle \prod\_{j=1}^{n} \int\_{0}^{\beta\hbar} d\tau\_{j} \Big\rangle \Big\langle \prod\_{k=1}^{n} V\_{\mathrm{int}}^{\boldsymbol{x}\_{0}} \Big[\boldsymbol{x}\left(\boldsymbol{\tau}\_{k}\right)\Big] \Big\rangle\_{\Omega,c}^{\boldsymbol{x}\_{0}} \end{split} \tag{41}$$

where <sup>0</sup> <sup>2</sup> <sup>2</sup> int 0 1 2 *<sup>x</sup> V x Vx M x x* (the kinetic energy terms in Eq. (30) and Eq. (31) cancel each other out).

As *n* tends to infinity, *W x <sup>n</sup>* <sup>0</sup> approaches the exact value of the centroid potential *W x* <sup>0</sup> in Eq. (28), which is independent of the trial Ω. But the truncated sum in Eq. (41) does depend on Ω, and the optimal choice of this trial frequency at a given order of KP expansion and at a particular centroid position 0 *<sup>x</sup>* is determined by the least-dependence of <sup>0</sup> *<sup>x</sup> Wn* on Ω itself. This is the so-called frequency of least dependence, which provides a variational approach to determine the optimal value of Ω, **opt,***<sup>n</sup> x*<sup>0</sup> (Kleinert 2004).

Of particular interest is the special case when 1 *n* , which turns out to be identical to the original GTFK variational approach. An important property of KP1 or the GTFK variational approach is that there is a definite upper bound for the computed *W x* 1 0 by virtue of the

Jensen-Peierls inequality, i.e., from Eq. (32) and (35):

$$e^{-\beta \mathcal{W}(\mathbf{x}\_0)} = Q\_{\Omega}^{\mathbf{x}\_0} \left\langle \exp \left( -\frac{\mathsf{A} - \mathsf{A}\_{\Omega}^{\mathbf{x}\_0}}{\hbar} \right) \right\rangle\_{\Omega}^{\mathbf{x}\_0} \geq Q\_{\Omega}^{\mathbf{x}\_0} \exp \left\langle -\frac{\mathsf{A} - \mathsf{A}\_{\Omega}^{\mathbf{x}\_0}}{\hbar} \right\rangle\_{\Omega}^{\mathbf{x}\_0} = e^{-\beta \mathcal{W}\_1^{\Omega}(\mathbf{x}\_0)}. \tag{42}$$

Note that by choosing 0 (i.e., the reference state is for a free particle), KP1 or GTFK (Giachetti and Tognetti 1985; Feynman and Kleinert 1986) reduces to the Feynman-Hibbs approach (Feynman, Hibbs et al. 2005). For higher orders of *n*, unfortunately, it is not guaranteed that a minimum of <sup>0</sup> *<sup>x</sup> Wn* actually exists as a function of Ω. In this case, the least dependent Ω is obtained from the condition that the next derivative of <sup>0</sup> *<sup>x</sup> Wn* with respect to Ω is set to zero. Consequently, Ω is considered as a variational parameter in the Kleinert perturbation theory such that **opt,** <sup>0</sup> 0 *<sup>x</sup> W x n n* is least-dependent on Ω.

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

approximation sacrifices some accuracy, in exchange, it allows analyses of quantum mechanical vibration and tunneling, and their separate contributions to the *W*. Positive and negative values of *wi* raise (vibration) and lower (tunneling) the original potential, respectively. In practice, real frequencies from the INM analysis often yields positive *wi* 's in Eq. (43) with dominant contributions from zero-point-energy effects. For imaginary frequencies in the INM, the values of *wi* are often negative, due to tunneling contributions. To obtain analytical expressions for the expectation values in Eq. (41), we use an *m*th order polynomial (P*m*) to approximate or interpolate the potential along *qi*. Hereafter, an *m*th order polynomial representation of the original potential energy function obtained with an interpolating step size *q* Å both in the forward and backward directions along the normal mode coordinate at *x*0 is denoted as P*m*-*q*A. Note that analytical results for P4 have been used by Kleinert for a quadratic-quartic anharmonic potential and a double-well potential (Kleinert 2004); however, higher order polynomials are needed to achieve the desired accuracy in real systems. We have thus derived the analytical closed forms of Eq. (41) up to P20 (Wong and Gao 2007; Wong 2008; Wong and Gao 2008; Wong, Richard et al. 2009; Wong, Gu et al. 2012). Consequently, the *W* as a function of an arbitrary Ω can be promptly obtained. This provides a convenient way to determine the least dependent Ω value without computing the complicated smearing integrals [Eq. (39)] iteratively for different trial values of Ω by Monte Carlo multi-dimensional numerical integrations. In fact, after the interpolating potential along each instantaneous normal-mode coordinate is determined, there is little computational cost for obtaining the *W*. Thereby, high level *ab initio* or densityfunctional (DFT) methods can be used to evaluate the potential energy function for *ab initio* 

path-integral calculations (Wong, Richard et al. 2009; Wong, Gu et al. 2012).

2. The original potential *V* is scanned from the configuration <sup>3</sup>

for these polynomials have been analytically integrated.

derivative for KP1 and usually second derivative for KP2).

3. After the P20-0.1A interpolations, each 0

, ( ) *<sup>x</sup> w q i n <sup>i</sup>*

instantaneous normal-mode coordinates.

et al. 2009; Wong, Gu et al. 2012):

0 *N*

4. The values of 0

0 , ( ) *<sup>x</sup> w q i n <sup>i</sup>*

1. For each <sup>3</sup>

The computational procedure for obtaining the first and second order KP approximations to the centroid potential using our automated integration-free path-integral (AIF-PI) method is summarized below (Wong and Gao 2007; Wong 2008; Wong and Gao 2008; Wong, Richard

*x* , the mass-scaled Hessian matrix is diagonalized to obtain <sup>0</sup>

points both in the forward and backward directions to interpolate *V* as P20-0.1A. A step size of 0.1 Å should be a reasonable choice to yield *W* in a few per cent of the exact.

, ( ) *<sup>x</sup> w q i n <sup>i</sup>*

using the analytical expressions of KP1/P20 or KP2/P20. Note that the path integrals

The procedure presented above is integration-free and essentially automated (Wong and Gao 2007; Wong 2008; Wong and Gao 2008; Wong, Richard et al. 2009; Wong, Gu et al. 2012). We hope it could be used by non-path-integral experts or experimentalists as a "black-box" for any given system. We are currently developing a formalism to systematically couple

on Ω, i.e., zeroing the lowest order derivative of 0

are determined by numerically locating the least dependence of

0 *N*

as a function of Ω is readily obtained

, ( ) *<sup>x</sup> w q i n <sup>i</sup>*

*x* along each 0 *<sup>x</sup>*

<sup>3</sup>*<sup>N</sup> <sup>x</sup> <sup>q</sup>* .

w.r.t. Ω (first

*<sup>i</sup> q* for 10

is the centroid potential for normal mode *i*. Although the INM

where <sup>0</sup> , *<sup>x</sup> w q i n <sup>i</sup>*

This variational criterion relies on the uniformly and exponentially convergent property demonstrated from the KP theory. Kleinert and coworkers proved that his theory exhibits this property in several strong anharmonic-coupling systems. More importantly, this remarkably fast convergent property can also be observed even for computing the *electronic* ground state energy of a hydrogen atom (3 degrees of freedom). The ground state energy was determined by calculating the electronic centroid potential at the zero-temperature limit. The accuracies of the first three orders of the KP theory for a hydrogen atom are 85%, 95%, and 98%, respectively (Kleinert 2004).

In practice, for odd *n*, there is typically a minimum point in Ω, but due to the alternating sign of the cumulants in Eq. (41), there is usually *no* minimum in Ω for even *n*. Nevertheless, the frequency of least-dependence for an even order perturbation in *n* can be determined by locating the inflexion point, i.e., the zero-value of the second derivative of <sup>0</sup> *<sup>x</sup> Wn* with respect to Ω. Since the KP expansion is uniformly and exponentially converged, Kleinert has demonstrated that the least-dependent plateau in <sup>0</sup> *<sup>x</sup> Wn* , which is characterized by a minimum point for odd *n* or by an inflexion point for even *n*, grows larger and larger with increasing orders of *n* (Kleinert 2004)*.*

#### **5.2 Automated integration-free path-integral method**

An especially attractive feature of Eq. (41) is that the if the real system potential is expressed as a series of polynomials or Gaussians, then analytic expressions of Eq. (41) can be obtained, making the computation extremely efficient because the time-demanding Monte Carlo samplings for multi-dimensional numerical integrations could be avoided. Hereafter, the level of calculations up to *n*th order KP expansion for an *m*th-orderpolynomial potential is denoted as KP*n*/P*m*. For other potentials, KP*n* theory still involves elaborate *n*-dimensional space-time (2*n* degrees of freedom) smearing integrals in Eq. (39). The intricacy of the smearing integrals increases tremendously for multidimensional potentials, where Ω becomes a 3 3 *N N* matrix Ω*ij* for *N* nuclei. This complexity is a major factor limiting applications of the KP theory beyond KP1, the original FK approach.

To render the KP theory feasible for many-body systems with *N* particles, we decouple the instantaneous normal mode (INM) coordinates <sup>0</sup> <sup>3</sup>*<sup>N</sup> <sup>x</sup> <sup>q</sup>* for a given configuration <sup>3</sup> 0 *N x*

(Wong and Gao 2007; Wong 2008; Wong and Gao 2008; Wong, Richard et al. 2009; Wong, Gu et al. 2012). Hence the multidimensional *V* effectively reduces to 3*N* one-dimensional potentials along each normal mode coordinate. Note that INM are naturally decoupled through the second order Taylor expansion of *V*. The approximation of decoupling the INM coordinates has also been used elsewhere (Stratt 1995; Deng, Ladanyi et al. 2002). This approximation is particularly suited for the KP theory because of the exponential decaying property of the Gaussian convolution integrals in Eq. (39). In the decoupling INM approximation, the total effective centroid potential for *N* nuclei can be simplified as:

$$\mathcal{W}\_n^{\Omega} \left( \{ \mathbf{x}\_0 \} ^{\text{3N}} \right) \approx V \left( \{ \mathbf{x}\_0 \} ^{\text{3N}} \right) + \sum\_{i=1}^{\text{3N}} w\_{i,n}^{\Omega} \left( q\_i^{\mathbf{x}\_0} \right) \tag{43}$$

This variational criterion relies on the uniformly and exponentially convergent property demonstrated from the KP theory. Kleinert and coworkers proved that his theory exhibits this property in several strong anharmonic-coupling systems. More importantly, this remarkably fast convergent property can also be observed even for computing the *electronic* ground state energy of a hydrogen atom (3 degrees of freedom). The ground state energy was determined by calculating the electronic centroid potential at the zero-temperature limit. The accuracies of the first three orders of the KP theory for a hydrogen atom are 85%,

In practice, for odd *n*, there is typically a minimum point in Ω, but due to the alternating sign of the cumulants in Eq. (41), there is usually *no* minimum in Ω for even *n*. Nevertheless, the frequency of least-dependence for an even order perturbation in *n* can be determined by locating the inflexion point, i.e., the zero-value of the second derivative of <sup>0</sup> *<sup>x</sup> Wn* with respect to Ω. Since the KP expansion is uniformly and exponentially converged, Kleinert has demonstrated that the least-dependent plateau in <sup>0</sup> *<sup>x</sup> Wn* , which is characterized by a minimum point for odd *n* or by an inflexion point for even *n*, grows larger and larger with

An especially attractive feature of Eq. (41) is that the if the real system potential is expressed as a series of polynomials or Gaussians, then analytic expressions of Eq. (41) can be obtained, making the computation extremely efficient because the time-demanding Monte Carlo samplings for multi-dimensional numerical integrations could be avoided. Hereafter, the level of calculations up to *n*th order KP expansion for an *m*th-orderpolynomial potential is denoted as KP*n*/P*m*. For other potentials, KP*n* theory still involves elaborate *n*-dimensional space-time (2*n* degrees of freedom) smearing integrals in Eq. (39). The intricacy of the smearing integrals increases tremendously for multidimensional potentials, where Ω becomes a 3 3 *N N* matrix Ω*ij* for *N* nuclei. This complexity is a major factor limiting applications of the KP theory beyond KP1, the

To render the KP theory feasible for many-body systems with *N* particles, we decouple the

(Wong and Gao 2007; Wong 2008; Wong and Gao 2008; Wong, Richard et al. 2009; Wong, Gu et al. 2012). Hence the multidimensional *V* effectively reduces to 3*N* one-dimensional potentials along each normal mode coordinate. Note that INM are naturally decoupled through the second order Taylor expansion of *V*. The approximation of decoupling the INM coordinates has also been used elsewhere (Stratt 1995; Deng, Ladanyi et al. 2002). This approximation is particularly suited for the KP theory because of the exponential decaying property of the Gaussian convolution integrals in Eq. (39). In the decoupling INM

approximation, the total effective centroid potential for *N* nuclei can be simplified as:

 <sup>0</sup> <sup>3</sup> 3 3 00,

*n i n i*

*W x Vx w q*

*<sup>N</sup> N N <sup>x</sup>*

1

*i*

,

(43)

<sup>3</sup>*<sup>N</sup> <sup>x</sup> <sup>q</sup>* for a given configuration <sup>3</sup>

0 *N x*

95%, and 98%, respectively (Kleinert 2004).

increasing orders of *n* (Kleinert 2004)*.*

original FK approach.

**5.2 Automated integration-free path-integral method** 

instantaneous normal mode (INM) coordinates <sup>0</sup>

where <sup>0</sup> , *<sup>x</sup> w q i n <sup>i</sup>* is the centroid potential for normal mode *i*. Although the INM approximation sacrifices some accuracy, in exchange, it allows analyses of quantum mechanical vibration and tunneling, and their separate contributions to the *W*. Positive and negative values of *wi* raise (vibration) and lower (tunneling) the original potential, respectively. In practice, real frequencies from the INM analysis often yields positive *wi* 's in Eq. (43) with dominant contributions from zero-point-energy effects. For imaginary frequencies in the INM, the values of *wi* are often negative, due to tunneling contributions.

To obtain analytical expressions for the expectation values in Eq. (41), we use an *m*th order polynomial (P*m*) to approximate or interpolate the potential along *qi*. Hereafter, an *m*th order polynomial representation of the original potential energy function obtained with an interpolating step size *q* Å both in the forward and backward directions along the normal mode coordinate at *x*0 is denoted as P*m*-*q*A. Note that analytical results for P4 have been used by Kleinert for a quadratic-quartic anharmonic potential and a double-well potential (Kleinert 2004); however, higher order polynomials are needed to achieve the desired accuracy in real systems. We have thus derived the analytical closed forms of Eq. (41) up to P20 (Wong and Gao 2007; Wong 2008; Wong and Gao 2008; Wong, Richard et al. 2009; Wong, Gu et al. 2012). Consequently, the *W* as a function of an arbitrary Ω can be promptly obtained. This provides a convenient way to determine the least dependent Ω value without computing the complicated smearing integrals [Eq. (39)] iteratively for different trial values of Ω by Monte Carlo multi-dimensional numerical integrations. In fact, after the interpolating potential along each instantaneous normal-mode coordinate is determined, there is little computational cost for obtaining the *W*. Thereby, high level *ab initio* or densityfunctional (DFT) methods can be used to evaluate the potential energy function for *ab initio*  path-integral calculations (Wong, Richard et al. 2009; Wong, Gu et al. 2012).

The computational procedure for obtaining the first and second order KP approximations to the centroid potential using our automated integration-free path-integral (AIF-PI) method is summarized below (Wong and Gao 2007; Wong 2008; Wong and Gao 2008; Wong, Richard et al. 2009; Wong, Gu et al. 2012):


The procedure presented above is integration-free and essentially automated (Wong and Gao 2007; Wong 2008; Wong and Gao 2008; Wong, Richard et al. 2009; Wong, Gu et al. 2012). We hope it could be used by non-path-integral experts or experimentalists as a "black-box" for any given system. We are currently developing a formalism to systematically couple instantaneous normal-mode coordinates.

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

equation (Wong 2008; Wong and Gao 2008), e.g., see Table 3. Together with the accurate low-lying excitation energies (Ramírez and López-Ciudad 2001) which could be obtained by the frequency analysis of the Hessian matrix at the sole minimum point at absolute zero temperature (including tunneling splitting), potentially one day our AIF-PI method could replace MC or MD *simulations* to have highly *reproducible* and *precise* free-energy *calculations*

> Molecule Quantum Harmonic KP1 KP2 HCl 4.231 4.274 4.253 4.234 HF 5.732 5.793 5.762 5.736 H2 6.193 6.284 6.238 6.202

Table 3. Ground state energy values (kcal/mol) for hydrogen chloride, hydrogen fluoride,

approximation, and our AIF-PI method based on first and second orders of the Kleinert's

**Born-Oppenheimer Approximation** 

**Hartree-Fock** (HF) theory **Kleinert's variational perturbation** 

Table 4. Comparison (1) between Kleinert's variational perturbation (KP) theory and Hartree-Fock (HF) theory, (2) between our decoupled instantaneous normal coordinate approximation and independent electron approximation, and (3) between our integrationfree path-integral results for polynomials in the KP theory and Roothann-Hall basis function

**Internuclear** Schrödinger equation Systematic internuclear thermodynamics theory

All thermodynamic properties virtually can be derived from **quantum partition functions**

**theory** for centroid effective potential

**Decoupled instantaneous normal coordinate** approximation (DINCA)

We propose interpolating potential energy functions to *m*th order **polynomials** in which analytic results of path-integration can be derived

Quantum effects from **vibration** and **tunneling** are separated and quantified in one mathematical framework

Work out a formalism to **systematically couple instantaneous normal coordinates**

and hydrogen molecules from the Morse potential using the harmonic-oscillator

for many-body systems.

variational perturbation theory (KP1 and KP2).

**Electronic** Schrödinger equation *Ab initio* molecular orbital theory

Most molecular properties of interest are at **low lying electronic energy states**

**Independent electron** (single-electron) approximation

**Roothaan and Hall** expressed the Fock operator in terms of **basis functions** for solving HF equations in matrix algebra self-consistently (SCF)

Explain chemical properties in terms of **frontier occupied** and **unoccupied**  molecular orbitals

Post Hartree-Fock method to include **correlation energy by systematically couple single-electron orbitals**

approach for HF theory.

Due to the integration-free feature, our AIF-PI method is computationally efficient such that the potential energy can be evaluated using *ab initio* or density-functional theory (DFT) for performing the so-called *ab initio* path-integral calculations. Consequently, we used DFT to construct the internuclear potential energy function for computing kinetic isotope effects (KIE) on several series of proton transfer reactions in water with the AIF-PI method. These reactions are relevant to biosynthesis of cholesterol. The computed KIE results at the KP2 level are in good agreement with experiment (Wong, Richard et al. 2009). Recently, we also employed the same computational technique to perform *ab initio* path-integral calculations of KIE on some RNA model reactions. Again, as shown in Table 2, the calculated values are in good agreement with experiments (Wong, Gu et al. 2012).


Table 2. Calculated primary kinetic isotope effects (KIEs) on 2' nucleophile (18*k*Nu) and 5' leaving (18*k*Lg or 34*k*Lg) oxygens for RNA-model reactions using our AIF-PI method based on second order of Kleinert's variational perturbation theory (KP2), along with the most relevant available experimental (Expt) results for comparison. Experimental errors in the last decimal place are given in parenthesis.

Another compelling feature of the AIF-PI method is that it does not suffer the convergence difficulties of PIMC or PIMD simulations at the zero-temperature limit, i.e., absolute zero temperature. At the zero-temperature limit (*T =* 0 K), in principle, minimizing the centroid effective potential with respect to the nuclear positions can give us two important physical quantities: the exact value of the eigenenergy for zero-point motion (i.e., the zero-point energy ZPE or the ground state energy) and the exact expectation values of the nuclear positions at the ground state (Ramírez, López-Ciudad et al. 1998; Ramírez and López-Ciudad 1999), i.e.,

$$\lim\_{T \to 0} \mathcal{W}\_{\text{min}} \left( \mathfrak{x}\_{\text{min}} \right) = E\_{0\prime} \tag{44}$$

and

$$\mathbf{x}\_{\min} = \langle \boldsymbol{\nu}\_{0} | \boldsymbol{x} | \boldsymbol{\nu}\_{0} \rangle\_{\prime} \tag{45}$$

where *x* is the position operator, and min *x* and *W x* min min are, respectively, the coordinate and value at the (global) minimum of the centroid potential. In Eq. (44) and (45), 0 is the nuclear ground state wave function and *E*0 is the lowest eigenvalue of the Hamiltonian, i.e., the zero-point energy. In a forthcoming paper, we will have a rigorous proof showing that in fact at absolute zero temperature, there is *only one* stationary and minimum point in centroid potential, which is true even for any many-body systems. Hence, our recently derived analytical zero-temperature-limit results provide a convenient way to compute these two important physical quantities without solving the Schrödinger

Due to the integration-free feature, our AIF-PI method is computationally efficient such that the potential energy can be evaluated using *ab initio* or density-functional theory (DFT) for performing the so-called *ab initio* path-integral calculations. Consequently, we used DFT to construct the internuclear potential energy function for computing kinetic isotope effects (KIE) on several series of proton transfer reactions in water with the AIF-PI method. These reactions are relevant to biosynthesis of cholesterol. The computed KIE results at the KP2 level are in good agreement with experiment (Wong, Richard et al. 2009). Recently, we also employed the same computational technique to perform *ab initio* path-integral calculations of KIE on some RNA model reactions. Again, as shown in Table 2, the calculated values are

Reaction KP2 Expt

Native 0.968 1.059 0.981(3) 1.034(4) S3′ 1.043 1.008 1.119(6) 1.0118(3) S5′ 1.042 1.002 1.025(5) 1.0009(1)

relevant available experimental (Expt) results for comparison. Experimental errors in the last

Another compelling feature of the AIF-PI method is that it does not suffer the convergence difficulties of PIMC or PIMD simulations at the zero-temperature limit, i.e., absolute zero temperature. At the zero-temperature limit (*T =* 0 K), in principle, minimizing the centroid effective potential with respect to the nuclear positions can give us two important physical quantities: the exact value of the eigenenergy for zero-point motion (i.e., the zero-point energy ZPE or the ground state energy) and the exact expectation values of the nuclear positions at the ground state (Ramírez, López-Ciudad et al. 1998; Ramírez and López-

> min min 0 <sup>0</sup> lim , *<sup>T</sup> Wx E*

min 0 0 *x x* 

where *x* is the position operator, and min *x* and *W x* min min are, respectively, the coordinate and value at the (global) minimum of the centroid potential. In Eq. (44) and (45),

 0 is the nuclear ground state wave function and *E*0 is the lowest eigenvalue of the Hamiltonian, i.e., the zero-point energy. In a forthcoming paper, we will have a rigorous proof showing that in fact at absolute zero temperature, there is *only one* stationary and minimum point in centroid potential, which is true even for any many-body systems. Hence, our recently derived analytical zero-temperature-limit results provide a convenient way to compute these two important physical quantities without solving the Schrödinger

Table 2. Calculated primary kinetic isotope effects (KIEs) on 2' nucleophile (18*k*Nu) and 5' leaving (18*k*Lg or 34*k*Lg) oxygens for RNA-model reactions using our AIF-PI method based on second order of Kleinert's variational perturbation theory (KP2), along with the most

<sup>18</sup>*k*Nu 18,34*k*Lg

(44)

, (45)

<sup>18</sup>*k*Nu 18,34*k*Lg

in good agreement with experiments (Wong, Gu et al. 2012).

decimal place are given in parenthesis.

Ciudad 1999), i.e.,

and

equation (Wong 2008; Wong and Gao 2008), e.g., see Table 3. Together with the accurate low-lying excitation energies (Ramírez and López-Ciudad 2001) which could be obtained by the frequency analysis of the Hessian matrix at the sole minimum point at absolute zero temperature (including tunneling splitting), potentially one day our AIF-PI method could replace MC or MD *simulations* to have highly *reproducible* and *precise* free-energy *calculations* for many-body systems.


Table 3. Ground state energy values (kcal/mol) for hydrogen chloride, hydrogen fluoride, and hydrogen molecules from the Morse potential using the harmonic-oscillator approximation, and our AIF-PI method based on first and second orders of the Kleinert's variational perturbation theory (KP1 and KP2).


#### **Born-Oppenheimer Approximation**

Table 4. Comparison (1) between Kleinert's variational perturbation (KP) theory and Hartree-Fock (HF) theory, (2) between our decoupled instantaneous normal coordinate approximation and independent electron approximation, and (3) between our integrationfree path-integral results for polynomials in the KP theory and Roothann-Hall basis function approach for HF theory.

Finally, we make a quite interesting table (Table 4) to compare the traditional *ab initio*  molecular orbital theory for electronic structure calculations with our systematic approach for computing internuclear quantum effects. In short, the rigor and the spirit of both types of methods is the same. We first breakdown or dissect a complicated many-body problem into many one-body problems. Then we identify which one bodies are more important. Next we couple back those important one bodies to systematically approach the exact.

#### **6. Systematic** *ab initio* **path-integral free-energy expansion approach**

In order to systematically refine a classical free-energy profile to become ultimate quantum free-energy profile, in which both electrons and nuclei are treated quantum mechanically and adiabatically, we are developing a systematic *ab initio* path-integral free-energy expansion (SAI-PI-FEE; ) approach. In this approach, we combine our novel freeenergy expansion (FEE) method (Section 4) with our automated integration-free pathintegral (AIF-PI) method (Section 5.2) such that we can perform *ab initio* path-integral simulations for realistic molecular systems. The key of this combination is that first we realize the quantum partition function can be computed as a classical configuration shown in Eq. (27), then now in Eq. (23), we treat the *E* as:

$$
\Delta E = \mathcal{W} - \mathcal{V}\_{\prime} \tag{46}
$$

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

This UPDF can be used to determine the change of free energy *G* in Eq. (23), by simply locating the intersection point of two probability density functions (Nanda, Lu et al. 2005; Chipot and Pohorille 2007). In Eq. (47), *E* is a variable for the difference of the Hamiltonian or energy between two levels of theory, while *K*, *a*, *b*, and *s* are the fitting parameters. In Figure 1, we demonstrate the simultaneous fitting to the UPDF to determine the change of free-energy for a water molecule from HF/6-31G(d) to MP2/6-311G(d,p). The intersection point of the two probability functions at −170.917 kcal/mol is the best estimate value for the

In this chapter, we (wongky@biomaps.rutgers.edu; kiniu@alumni.cuhk.net) discuss

*ab initio* path-integral level in molecular simulations. Since quantum free energy or partition function is a universal central quantity in thermodynamics of biology, chemistry, and

Anderson, J. B. (1975). Random-walk simulation of the Schroedinger equation. Hydrogen

Arfken, G. B. and H.-J. Weber (2001). *Mathematical methods for physicists*. San Diego,

Ballhausen, C. J. and A. E. Hansen (1972). Electronic spectra. *Annual Review of Physical* 

Birkhoff, G. D. (1931). Proof of the Ergodic Theorem. *Proceedings of the National Academy of* 

Born, M. and J. R. Oppenheimer (1927). Zur Quantentheorie der Molekeln (On the Quantum

Brooks, B. R., R. E. Bruccoleri, et al. (1983). CHARMM: a program for macromolecular

Brown, L. M. (2005). *Feynman's thesis : a new approach to quantum theory*. Singapore ;

Cao, J. and G. A. Voth (1994). The formulation of quantum statistical mechanics based on the

Car, R. and M. Parrinello (1985). Unified approach for molecular dynamics and densityfunctional theory. *Physical Review Letters*, Vol. 55, No. 22, pp. 2471-2474. Ceperley, D. M. (1995). Path integrals in the theory of condensed helium. *Reviews of Modern* 

*Journal of Chemical Physics*, Vol. 101, No. 7, pp. 6168-6183.

Academic Press. See exercise 2.2.3 for the fundamental definition of an orthogonal unit vector; See section 2.10 for the discussion of the covariant and contravariant

Theory of Molecules). *Annalen der Physik*. 84,pp. 457-484. English translation: *Quantum chemistry : classic scientific papers,* Hettema, H., Ed.; World Scientific:

energy, minimization, and dynamics calculations. *Journal of Computational* 

Feynman path centroid density. IV. Algorithms for centroid molecular dynamics.

Sciences and wish that it could be used by non-specialists as a black box one day.

ion (H3+). *Journal of Chemical Physics*, Vol. 63, No. 4, pp. 1499-1503.

*Sciences of the United States of America*, Vol. 17, No. 12, pp. 656-660.

method to systematically generate quantum free-energy profiles at an

method would be very crucial in both Life and Materials

*G* in Table 1 above.

developing the

**8. References** 

**7. Conclusion and outlook** 

physics, we anticipate our

vectors.

*Chemistry*, Vol. 23, pp. 15-38.

Singapore; London, 2000, pp 1-24.

*Chemistry*, Vol. 4, No. 2, pp. 187-217.

Hackensack, NJ, World Scientific.

*Physics*, Vol. 67, No. 2, pp. 279-355.

where *V* is the original internuclear potential and *W* is the centroid potential. So once we get the accurate value of *W* using our AIF-PI method, we can go ahead using our FEE method to systematically upgrade the level of our classical free-energy profile to an *ab initio* pathintegral level, in which zero-point energy and tunnelling effects in nuclei, and isotope effects could all be incorporated.

In order to rigorously validate our method in a more effective way, the free-energy perturbation (FEP) in the Hamiltonian space will be performed, using the recently derived "universal" probability density function (UPDF), which is defined as follows:

exp . *b Es P E K ab E s e* (47)

Fig. 1. Free energy perturbation for a water molecule in the Hamiltonian space using the universal probability density function (UPDF).

This UPDF can be used to determine the change of free energy *G* in Eq. (23), by simply locating the intersection point of two probability density functions (Nanda, Lu et al. 2005; Chipot and Pohorille 2007). In Eq. (47), *E* is a variable for the difference of the Hamiltonian or energy between two levels of theory, while *K*, *a*, *b*, and *s* are the fitting parameters. In Figure 1, we demonstrate the simultaneous fitting to the UPDF to determine the change of free-energy for a water molecule from HF/6-31G(d) to MP2/6-311G(d,p). The intersection point of the two probability functions at −170.917 kcal/mol is the best estimate value for the *G* in Table 1 above.
