Developing a Systematic Approach for Ab Initio Path-Integral Simulations

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).

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 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, pK 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.

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

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: ˆ, mole n n n HE   where ˆm ole H is the complete (non-relativistic) molecular Hamiltonian, n  and n E 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 N n nuclei and N e electrons can fortunately be written in an analytic closed form (thanks to the inverse squaredistance proportionality in Coulomb's electrostatic force law): In Eq. (2), the units are atomic units, M j is the mass ratio of nucleus j to electron, and Z j is the atomic number of nucleus j. The Laplacian operators 2 j  and 2 i  denote the second order differentiation with respect to the coordinates of the jth nucleus and the ith 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 jth and the j'th nuclei is jj x  ; the separation between the ith and the i′th electrons is ii r  ; the distance between the jth nucleus and the ith electrons is ij r .

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 qm Q (McQuarrie 2000), by the following summation of the , B k 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 0 E , 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.

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 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: 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. www.intechopen.com

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 E n 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: where h is Planck's constant and p is the momenta associated with the nuclear coordinates x. Subsequently, the classical free energy G cl of a molecular system can be expressed in terms of the classical partition function Q cl as follows: 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 n N -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: where  is Dirac delta function, z  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 M j , 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) www.intechopen.com 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 , we obtain the average force over all ensembles or over all degrees of freedom, which is called the mean force (Kirkwood 1935 Gz , 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ńin, 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): and the unit vector for q  can be proved as (in contravariant space): 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 www.intechopen.com 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): 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 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.

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: 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 ,e x p , 2 n nn is equal to the time average in MD simulations: 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.

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 G cl corresponding to the partition function in Eq. (8) as follows: Next, we rewrite Eq. (18) (20) where 0 G is the free energy of the reference system, 0 E is the energy at a point in the phase space of the reference system, and 0  is an ensemble average for the reference system. From Eq. (20), we obtain Zwanzig's free-energy perturbation (Zwanzig 1954): 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   0 exp EE     over the ensembles sampled in the reference system. This is exactly the idea behind many enhanced sampling methods, such as the umbrella sampling technique. www.intechopen.com

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 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, pK 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: Next we expand the ensemble average in Eq. (22) and sum up the prefactors into a series of cumulants: ln exp ln exp , 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 HL G , 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: x is a position vector of 3N Cartesian coordinates of the system, D n is the nth-order tensor operator for differentiation with respect to the 3N coordinates (e.g., 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., xx 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 E HL 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)  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.
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.

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.
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: 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 Voth 1996;Ramírez, López-Ciudad et al. 1998;Ramírez and López-Ciudad 1999;Kleinert 2004;Feynman, Hibbs et al. 2005): Given the centroid potential   0 Wx , 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   0 Wx 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   0 Wx in Eq. (27) is defined as follows Voth 1996;Ramírez, López-Ciudad et al. 1998;Ramírez and López-Ciudad 1999;Kleinert 2004;Feynman, Hibbs et al. 2005 where τ is a real number and represents the component for pure imaginary time in path integral,   x  describes a path in space-time, 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' In Eq. (28), A is the quantum-statistical action: where  Vx 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).

www.intechopen.com
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 .
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 , 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 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


Fx   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 nth-order approximation, , to the centroid potential   0 Wx can be written as follows (Kleinert 2004 After using these smearing potentials given in Eq.  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   10 Wx  by virtue of the Jensen-Peierls inequality, i.e., from Eq. (32) and (35): 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   0 x n W  actually exists as a function of Ω. In this case, the least dependent Ω is obtained from the condition that the next derivative of   0 x n W  with respect to Ω is set to zero. Consequently, Ω is considered as a variational parameter in the Kleinert perturbation theory such that  opt, 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   0 x n W  with respect to Ω. Since the KP expansion is uniformly and exponentially converged, Kleinert has demonstrated that the least-dependent plateau in   0 x n W  , 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).

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 nth order KP expansion for an mth-orderpolynomial potential is denoted as KPn/Pm. For other potentials, KPn theory still involves elaborate n-dimensional space-time (2n degrees of freedom) smearing integrals in Eq. (39). The intricacy of the smearing integrals increases tremendously for multidimensional potentials, where Ω becomes a 3 3 NN  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   0 3N x q for a given configuration  3 0 N x Wong 2008;Wong and Gao 2008;Wong, Richard et al. 2009;Wong, Gu et al. 2012). Hence the multidimensional V effectively reduces to 3N 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: To obtain analytical expressions for the expectation values in Eq. (41), we use an mth order polynomial (Pm) to approximate or interpolate the potential along q i . Hereafter, an mth 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 Pm-qA. 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 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 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 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.
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 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., Wx 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 0 E 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 www.intechopen.com 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.

Molecule
Quantum We propose interpolating potential energy functions to mth order polynomials in which analytic results of path-integration can be derived Explain chemical properties in terms of frontier occupied and unoccupied molecular orbitals Quantum effects from vibration and tunneling are separated and quantified in one mathematical framework Post Hartree-Fock method to include correlation energy by systematically couple single-electron orbitals Work out a formalism to systematically couple instantaneous normal coordinates 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.

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: 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: 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.

Conclusion and outlook
In this chapter, we (wongky@biomaps.rutgers.edu; kiniu@alumni.cuhk.net) discuss developing the method to systematically generate quantum free-energy profiles at an 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 physics, we anticipate our method would be very crucial in both Life and Materials Sciences and wish that it could be used by non-specialists as a black box one day.