Computational modeling of the condensed phase based on classical statistical mechanics has been rapidly developing over the last few decades and has yielded important information on various systems containing up to millions of atoms. However, if a system of interest contains important quantum effects, well-developed classical techniques cannot be used. One way of treating finite temperature quantum systems at equilibrium has been based on Feynman’s imaginary time path integral approach and the ensuing quantum-classical isomorphism. This isomorphism is exact only in the limit of infinitely many classical quasiparticles representing each physical quantum particle. In this work, we present a reductionist perspective on this problem based on the emerging methodology of coarse-graining. This perspective allows for the representations of one quantum particle with only two classical-like quasiparticles and their conjugate momenta. One of these coupled quasiparticles is the centroid particle of the quantum path integral quasiparticle distribution. Only this quasiparticle feels the potential energy function. The other quasiparticle directly provides the observable averages of quantum mechanical operators. The theory offers a simplified perspective on quantum statistical mechanics, revealing its most reductionist connection to classical statistical physics. By doing so, it can facilitate a simpler representation of certain quantum effects in complex molecular environments.

Quantum statistical mechanics1 relates macroscopic properties of quantum systems to a microscopic description of particles constituting such systems. It is vital for the correct description of hydrogen tunneling and nuclear quantum effects,2–11 electron transfer,3,12–14 infrared and Raman spectra,5,15,16 etc., in addition to the paradigmatic examples of condensed phase helium and hydrogen.17–21 In comparison to classical statistical mechanics, the quantum theory is arguably more difficult to understand, since it involves particle delocalization, therefore effectively more degrees of freedom, and uses a more advanced mathematical apparatus. From the viewpoint of numerical simulations, quantum statistical mechanics is also more complex than its classical counterpart. For these reasons, classical-like approximations of the exact quantum statistical theory, providing a reductionist picture and facilitating conceptual modeling and numerical simulation of quantum systems, are important.

Feynman’s concept of the imaginary time path integral, with the ensuing quantum-classical isomorphism, can serve as a basis for such reductionist formulations.1,22 By replacing each quantum particle with infinitely many classical quasiparticles via a discretization of the path integral, this theory expresses any observable property of a quantum system at equilibrium and at finite temperature in terms of certain properties of an isomorphic classical system.23–28 A number of computational approaches based on combining the discretized path integral representation with molecular dynamics (MD), Monte Carlo (MC), or other classical algorithms have been developed.18,25–29 These methods would provide numerically exact quantum mechanical results in the limit of an infinite number of classical quasiparticles per quantum particle, the latter being usually denoted by the integer P.

In this paper, the most essential, minimalist picture of the manifestation of quantum effects in systems at finite temperature and at equilibrium is presented, based on a coarse-grained (CG) representation of the intermediate isomorphic quasiparticles into a single CG quasiparticle. The CG methodology is commonly used for dealing with classical atomistic models of molecules (see, e.g., Ref. 30) and, to the best of our knowledge, has not yet been applied to quantum problems, except in the case of the imaginary time Feynman path integral “centroid” variable.22,31 The resulting coarse-grained path integral (CG-PI) theory provides a new reductionist perspective on quantum statistical mechanics, allowing for a simple and intuitive representation of quantum systems at equilibrium.

Starting from a standard imaginary time path integral representation of the thermal density matrix in the position basis (Fig. 1(a)), we first demonstrate that the thermal density matrix can be written as

q e β H ˆ q = m π ħ 2 β d Q e β V CG q , Q , q ,
(1)

where H ˆ is the Hamiltonian of the system, m is the mass of the quantum particle, β is 1/kBT, Q is the coordinate of the effective particle introduced by coarse-graining the (infinitely) many intermediate isomorphic classical quasiparticles, and VCG is the effective CG potential (Fig. 1(c)). It is to be noted that in this paper, the word “quasiparticle” refers only to the classical-like particles that are utilized in the discretized imaginary time Feynman path integral. We are not referring to other kinds of quasiparticles, such as excitations in a Fermi liquid. For notational convenience, Eq. (1) is written for a one-dimensional case; generalizations to more dimensions are straightforward, as shown in Sec. S7 in the supplementary material.32 It is informative to compare Eq. (1) to the following path integral approximation for the density matrix with only one intermediate imaginary time slice (Fig. 1(b)):

q e β H ˆ q m π ħ 2 β d q 2 e β V 2 q , q 2 , q ,
(2)

where V2 is Feynman’s isomorphic potential22 VP at P = 2,

V P q , q 2 , , q P , q = i = 1 P m P 2 ħ 2 β 2 q i q i + 1 2 + V q i + V q i + 1 2 P ,
(3)

where V(q) is the potential in which the quantum particle is placed. Thus, the proposed CG approach can be considered as a renormalization of Feynman’s isomorphic potential for P = 2, leading to an, in principle, exact expression for the density matrix, Eq. (1). The formal definition of the effective potential VCG is given below.

FIG. 1.

(a) In the discretized imaginary time Feynman path integral representation, any element of the quantum thermal density matrix can be computed based on the properties of a collection of quasiparticles or “beads” having classical-like properties. (b) A rough approximation to the exact Feynman’s result would be to use only one intermediate quasiparticle (in total making up three beads). (c) This paper constructs a three-quasiparticle representation of the thermal density matrix. The intermediate quasiparticle is introduced by coarse-graining all intermediate quasiparticles shown in red in (a).

FIG. 1.

(a) In the discretized imaginary time Feynman path integral representation, any element of the quantum thermal density matrix can be computed based on the properties of a collection of quasiparticles or “beads” having classical-like properties. (b) A rough approximation to the exact Feynman’s result would be to use only one intermediate quasiparticle (in total making up three beads). (c) This paper constructs a three-quasiparticle representation of the thermal density matrix. The intermediate quasiparticle is introduced by coarse-graining all intermediate quasiparticles shown in red in (a).

Close modal

This approach is based on the extension of the CG methodology to the field of quantum statistical mechanics. In the original multiscale coarse-graining (MS-CG) theory,33,34 coarse-graining is performed for systems formed by n interacting particles obeying the laws of classical mechanics, for example, for (bio)molecules described by a set of classical particles (“atoms”) governed by classical force fields. A projection from the space of the coordinates of all atoms in the system rn onto the space of N CG variables RN, where Nn, is defined by means of N linear mapping operators M R N = M R 1 r n , , M R N r n ,

M R I r n = i = 1 n c I i r i ,
(4)

where cIi are coefficients. One possible interpretation of Eq. (4) is that the system of interest is divided into N domains, and CG variables RN designate positions of the centers of mass of each domain.33 Other interpretations are also possible, for example, when only some atoms (e.g., Cα in each amino acid) are used as representative atoms, or when the CG variables are defined in terms of vibrational normal modes.35 The CG potential U R N in the MS-CG theory is defined as33 

U R N = 1 β ln d r n e β u ( r n ) δ M R N ( r n ) R N + const ,
(5)

where u(rn) is the atomistic potential.

The extension of this CG methodology to quantum systems is possible due to the quantum-classical isomorphism. Coarse-graining the intermediate beads or quasiparticles that correspond to the intermediate states in a path integral (Fig. 1(a)) allows not only for a formally exact description of a quantum system with classical-like quasiparticles but also for a major reduction in the number of quasiparticles used. The first and the last quasiparticles, however, should not be coarse-grained, since the resulting expression for the thermal density matrix must keep an explicit dependence on q and q′. This background motivates the definition of the CG variable in the present theory of CG-PI as the center of mass of all intermediate quasiparticles in the path integral,

Q = q 2 + + q P P 1 .
(6)

In the limit of P → ∞ and if q = q′, the CG variable Q coincides with the centroid defined as22,31

q c = 1 ħ β 0 ħ β d τ q τ ,
(7)

where q τ is a closed imaginary time Feynman path. Here, q τ , with τ running from 0 to ħ β, is a continuous analogue of qi, with i running from 1 to P + 1, in the limit of P → ∞.

The effective potential VCG is defined here as

V CG q , Q , q = 1 β lim P ln Z CG ( P ) q , Q , q ,
(8)

where

Z CG ( P ) q , Q , q = π ħ 2 β m m P 2 π ħ 2 β P / 2 d q 2 d q P e β V P q , q 2 , , q P , q δ Q q 2 + + q P P 1 .
(9)

With this definition, Eq. (1) is automatically satisfied, taking into consideration that1 

q e β H ˆ q = lim P m P 2 π ħ 2 β P / 2 d q 2 d q P d Q e β V P q , q 2 , , q P , q δ Q q 2 + + q P P 1 .
(10)

The proposed definition, Eqs. (8) and (9), is in general similar to the definition of the CG potential in the MS-CG theory, Eq. (5), but differs from the latter in two aspects. First, integration over variables q and q′ is not performed, and the resulting CG potential VCG depends on a set of mixed resolution variables, namely, fine-grained (q, q′) and coarse-grained (Q) ones. Second, Eqs. (8) and (9) do not include an unspecified constant; therefore, in principle, exact values of the partition function and its correct dependence on temperature can be obtained.

In the case of a free particle, the CG potential assumes the following form:

V CG q , Q , q = k Q q ̄ free particle 2 Q q ̄ 2 + k Δ q 2 free particle 2 Δ q 2 + k 0 free particle ,
(11)

where

q ̄ = q + q 2 , Δ q = q q ,
(12)

and

(13)
k Q q ̄ free particle = 12 m ħ 2 β 2 , k Δ q 2 free particle = m ħ 2 β 2 ,
k 0 free particle = ln 3 2 β .
For a derivation, see Sec. S2.2 in the supplementary material.32 

The first term on the right-hand side of Eq. (11) penalizes large deviations of the CG particle from the average position between the initial and final points, while the second term penalizes large distances between the initial and final points. In the classical or high-temperature limit, the initial, final, and CG coordinates all collapse into the same point in space, while in the low-temperature limit, corresponding to strong quantum effects, the particles tend to be equidistributed over the whole space. A number of other analytical results for the CG potential are provided in the  Appendix and the supplementary material.32 

The CG potential can be used for computing the thermal density matrix, as well as for the analysis of quantum delocalization in the system of interest at a given temperature. In particular, maximum contributions to the thermal density matrix, according to Eq. (1), come from positions of the CG particle Q close to the value Qmin defined as follows:

Q min q , q = arg min Q V CG q , Q , q .
(14)

For example, for the harmonic oscillator with the potential energy V q = m Ω 2 q 2 / 2 , Q min = tanh f / f q ̄ , where f = ħ β Ω/2. In the classical limit, this Q min q ̄ , while in the limit of strong quantum effects, Qmin → 0 regardless of the q ̄ value. This analysis can be generalized to the case of an arbitrary molecular system. For any two geometrical configurations of the system of interest, one may compute the multi-dimensional Qmin and interpret it as the most important geometry of the system as it propagates from the first given geometrical configuration to the second one.

Series expansions of VCG in terms of ℏ as a small parameter can be used to study quantum systems near the classical limit. As demonstrated in Sec. S1 in the supplementary material,32 the lowest order nontrivial terms in VCG have the following form:

V CG q , Q , q = V Q + k Q q ̄ 2 Q q ̄ 2 + k Δ q 2 2 Δ q 2 + k 0 1 + O ħ 3 β 5 / 2 m 3 / 2 V q ̄ ,
(15)

where

k Q q ̄ = k Q q ̄ free particle R Q q ̄ , k Δ q 2 = k Δ q 2 free particle R Δ q 2 ,
(16)
k 0 = 1 2 β ln f eff ( q ̄ ) tanh f eff ( q ̄ ) sinh 2 f eff ( q ̄ ) 2 f eff ( q ̄ ) 4 ,
(17)

where R Q q ̄ and RΔq2 are renormalization factors,

R Q q ̄ = f eff ( q ̄ ) 2 tanh f eff ( q ̄ ) 3 f eff ( q ̄ ) tanh f eff ( q ̄ ) ,
(18)
R Δ q 2 = f eff ( q ̄ ) tanh f eff ( q ̄ ) ,
(19)

where feff is defined as

f eff ( q ̄ ) = ħ β 2 V ( q ̄ ) m .
(20)

Equation (15) demonstrates that the most essential part of the three-body interactions between the initial, final, and CG intermediate beads reduces in the limit of ħ → 0 to a sum of three simple less-than-three-body terms.

  1. The CG particle feels the potential V in a classical way. This is a one-body potential.

  2. The CG particle is bound by a harmonic spring to the middle point between the initial and final beads. Since this middle point is, in general, “empty” which is not coinciding with q or q′, this term can be considered as a quasi-two-body potential. In the classical limit, the spring stiffness k Q q ̄ reduces to the corresponding spring constant for a free particle, while in the limit of strong quantum effects, it tends to k Q q ̄ free particle f eff ( q ̄ ) / 3 .

  3. The initial and final quasiparticles attract to each other via a two-body effective harmonic potential. In the classical limit, the spring stiffness kΔq2 is the same as that of a free particle, while in the regime of strong quantum effects, it tends to k Δ q 2 free particle f eff ( q ̄ ) .

Two general strategies of improving the accuracy of Eq. (15) are possible. First, one may analytically derive explicit expressions for higher order terms [see Eq. (A15) and Secs. S1.2-1.4 in the supplementary material32]. Second, one may choose an Ansatz expression for the CG potential, for example,

V CG trial q , Q , q = R Q V Q + R Q q ̄ k Q q ̄ free particle 2 Q q ̄ 2 + R Δ q 2 k Δ q 2 free particle 2 Δ q 2 + R 0 k 0 free particle
(21)

and find the exact values of the renormalization factors R Q , R Q q ̄ , R Δ q 2 , R 0 using one of the variational methods described in the  Appendix (see also Sec. S4 in the supplementary material32). Both strategies (the use of higher order terms or an Ansatz expression) can employ the potential energy V(q) obtained from electronic structure methods, for example ab initio molecular dynamics (AIMD), or approximated by empirical force fields.

It is valuable to further manipulate the preceding expressions to cast them in a more classical-like form. An expectation value of a quantum mechanical operator A ˆ p ˆ , q ˆ at equilibrium at finite temperature is given by

A ˆ p ˆ , q ˆ = d q d q q e β H ˆ q q A ˆ q d q q e β H ˆ q
(22)

and can be computed by expressing the thermal density matrix in terms of the CG potential, Eq. (1), and the matrix elements of the operator A ˆ p ˆ , q ˆ in terms of the corresponding Weyl map36 WA, such that

q A ˆ q = d p 2 π ħ e i p Δ q / ħ W A p , q ̄ .
(23)

(Recall that the Weyl map is a strict way to get a classical-like function corresponding to a given quantum operator. For example, if A ˆ = q ˆ n , where n is an integer, then WA = qn, and if A ˆ = p ˆ n , then WA = pn.) The final result (see Sec. S5.1 in the supplementary material for a derivation32) is given by

A ˆ p ˆ , q ˆ = d P Q d Q d p d q ̄ e β H eff P Q , p , Q , q ̄ W A p , q ̄ d P Q d Q d p d q ̄ e β H eff P Q , p , Q , q ̄ ,
(24)

where H eff P Q , p , Q , q ̄ has the form of a classical two-particle Hamiltonian. The exact expression for Heff is given by

H eff P Q , p , Q , q ̄ = P Q 2 2 M Q 1 β ln C d Δ q e i p Δ q ħ β V CG q ̄ Δ q 2 , Q , q ̄ + Δ q 2 ,
(25)

where MQ and C are arbitrary positive constants. Based on Eq. (15), the following approximation for Heff can be obtained:

H eff P Q , p , Q , q ̄ P Q 2 2 M Q + p 2 2 m eff q ̄ + V Q + k Q q ̄ q ̄ 2 Q q ̄ 2 ,
(26)

where meff is a positive-valued function of q ̄ interpreted as an effective mass. For the derivation of Eq. (26) and explicit expressions for meff, see Sec. S5.2 in the supplementary material.32 In particular, when V ( q ̄ ) > π 2 m ħ 2 β 2 , the effective mass is given by meff = mRΔq2 with RΔq2 defined by Eq. (19).

A simple classical-like picture for the most essential part of the quantum effects emerges then from Eqs. (24) and (26). A physical quantum particle is represented here by two classical quasiparticles that we will call a centroid particle and an observed particle (Fig. 2). The centroid particle, but not the observed particle, feels the physical potential. On the other hand, only the coordinates and momenta of the observed particle, and not those of the centroid particle, are used to compute the values of the function W A p , q ̄ to be averaged in order to calculate the observable A ˆ p ˆ , q ˆ . The centroid and observed particles are connected by an effective harmonic spring. An algorithmic implementation of this approach seems straightforward based on the existing codes for classical MD or MC simulations and should not add significantly to the cost of simulations of molecular and biomolecular systems, especially if only some of the atoms in the system (e.g., hydrogen) exhibit quantum behavior.

FIG. 2.

A representation of a quantum particle by only two classical quasiparticles (called the centroid and observable particles) is suggested in this paper. Only the centroid particle feels the potential, and only the coordinates and momenta of the observable particle are used to compute a desired quantum observable. In the near classical limit (ħ → 0) for an arbitrary physical potential, or exactly for a harmonic oscillator, the two particles are connected by an effective harmonic spring. An observer is represented on the left of the diagram.

FIG. 2.

A representation of a quantum particle by only two classical quasiparticles (called the centroid and observable particles) is suggested in this paper. Only the centroid particle feels the potential, and only the coordinates and momenta of the observable particle are used to compute a desired quantum observable. In the near classical limit (ħ → 0) for an arbitrary physical potential, or exactly for a harmonic oscillator, the two particles are connected by an effective harmonic spring. An observer is represented on the left of the diagram.

Close modal

In the present scheme, the spring stiffness k Q q ̄ depends on the current position of the observable particle q ̄ . Note that k Q q ̄ is well-defined (real) not only for positive, but also for negative V ( q ̄ ) . In the latter case, the value of f eff ( q ̄ ) is purely imaginary, and the right hand side of Eq. (18) is real-valued. The negative values of k Q q ̄ achieved when V ( q ̄ ) < 4 π 2 m ħ 2 β 2 [e.g., near the top of a sufficiently high and narrow barrier in the case of a double-well potential V(q)] correspond to an effective repulsion between the centroid and observed quasiparticles due to the term k Q q ̄ q ̄ Q q ̄ 2 / 2 in Eq. (26). In this case, the two quasiparticles (those with the coordinates Q and q ̄ ) tend to move away from each other, falling from the top of the barrier toward two different local minima on the potential energy surface. However, as soon as the observed quasiparticle leaves the vicinity of the top of the barrier, implying that V ( q ̄ ) increases to − 4π2−2β−2 and above, the spring stiffness becomes positive again, and the two quasiparticles start attracting each other, tending to return to the vicinity of the same local minimum. Events when the observed and centroid quasiparticles in the CG-PI model sporadically appear in basins of two different local minima on the potential energy surface reflect the ability of a quantum particle to tunnel. Therefore, possible negative values of k Q q ̄ do not contradict to the physical sense and should not lead to numerical divergences in MD or MC simulations with the use of the effective potential Heff given by Eq. (26).

The effective mass of the observed quasiparticle meff in Eq. (26) also depends on q ̄ . In this case, meff is real and positive for any value of q ̄ , as discussed in Sec. S5.2 in the supplementary material.32 This definition of the effective mass does not seem to lead to any technical complications in MC simulations. However, in the case of MD, the varying effective mass causes the dependence of the force on the momentum, making common MD integrators, such as the leapfrog algorithm, inapplicable. The problem of how to reliably and efficiently run MD simulations for particles with changing masses in a thermostat goes beyond the scope of the present paper. To deal with this issue, we introduce here one more approximation for the effective Hamiltonian Heff, namely,

H eff P Q , p , Q , q ̄ P Q 2 2 M Q + p 2 2 m eff const + V Q + k Q q ̄ q ̄ 2 Q q ̄ 2 ,
(27)

where m eff const is a positive constant. Importantly, for momentum-independent quantum mechanical operators A ˆ = A ˆ q ˆ , the error introduced by the approximation of going from Eqs. (25) to (26) and the error of going from Eqs. (26) to (27) appear to exactly cancel each other (for the proof, see Sec. S5.4 in the supplementary material32). Thus, for a wide class of practically important observables, including average potential energies, average interatomic distances and functions thereof, radial distribution functions, etc., the CG-PI model can employ quasiparticles with constant masses without introducing any errors beyond those coming from approximating the CG potential by Eq. (15) or (21) or (A15).

Equation (24) with Heff with constant meff and k Q q ̄ is exact in the case of a quantum harmonic oscillator for arbitrarily strong quantum effects (that is, to all powers of ℏ) and for arbitrary operators A ˆ p ˆ , q ˆ . Renormalization of k Q q ̄ —and, for momentum-dependent observables, of meff—accounts for the major part of quantum corrections in the case of a general potential V(q). The remaining corrections follow from higher order terms in the CG potential, Eq. (15), resulting in additional corrections to the spring stiffness k Q q ̄ , the effective mass meff and additional terms in Heff anharmonically coupling Q and q ̄ , and higher order powers of p. In particular, the O β V ( I V ) q ̄ Δ q 4 term in the CG potential [see Eq. (S79) in the supplementary material32] results in relativistic-like O p 4 and higher order in p terms in Heff. If necessary, these terms can be analytically derived in a systematic manner or accounted for by means of a variational treatment with the Ansatz representations of the CG potential, as in Eq. (21). The second approach leads to the following expression for the effective Hamiltonian Heff:

H eff P Q , p , Q , q ̄ P Q 2 2 M Q + p 2 2 m R Δ q 2 + R Q V Q + R Q q ̄ k Q q ̄ free particle 2 Q q ̄ 2 ,
(28)

where the renormalization factors RΔq2,  RQ, and R Q q ̄ are the same as in Eq. (21).

Finally, note that the CG-PI formulas discussed in this section can be generalized to multidimensional cases (see Sec. S7 in the supplementary material32). In particular, the three-dimensional generalization of Eq. (15) has the following form:

V CG q , Q , q = V Q + 1 2 Q q ̄ K Q q ̄ Q q ̄ + 1 2 Δ q K Δ q 2 Δ q + k 0 × 1 + O ħ 3 β 5 / 2 m 3 / 2 V q ̄ ,
(29)

where q , q , Q , q ̄ , and Δq are the corresponding three-dimensional vectors, K Q q ̄ and KΔq2 are 3 × 3 matrices given by

(30)
K Q q ̄ = 12 m ħ 2 β 2 U ( q ̄ ) Λ Q q ̄ ( q ̄ ) U ( q ̄ ) ,
K Δ q 2 = m ħ 2 β 2 U ( q ̄ ) Λ Δ q 2 ( q ̄ ) U ( q ̄ ) ,
where U ( q ̄ ) is the matrix diagonalizing the Hessian of the physical potential at point q ̄ , given by

V q ̄ = U q ̄ Λ H q ̄ U q ̄
(31)

and Λ Q q ̄ ( q ̄ ) , Λ Δ q 2 ( q ̄ ) and Λ H q ̄ are diagonal 3 × 3 matrices with the elements obeying the following conditions:

Λ Q q ̄ ( q ̄ ) i i = f i 2 ( q ̄ ) tanh f i ( q ̄ ) f i ( q ̄ ) tanh f i ( q ̄ ) , Λ Δ q 2 ( q ̄ ) i i = f i ( q ̄ ) tanh f i ( q ̄ ) ,
(32)
f i ( q ̄ ) = ħ β 2 Λ H q ̄ i i m ,
(33)

with k0 given by

k 0 = 1 2 β i = 1 , 2 , 3 ln sinh 2 f i ( q ̄ ) f i ( q ̄ ) tanh f i ( q ̄ ) 2 f i 4 ( q ̄ ) .
(34)

The multidimensional analogue of the expression for the effective Hamiltonian, Eq. (26), assumes the following form:

H eff P Q , p , Q , q ̄ P Q 2 2 M Q + 1 2 p M eff 1 p + V Q + 1 2 Q q ̄ K Q q ̄ Q q ̄ ,
(35)

where the effective mass matrix M eff = m U ( q ̄ ) Λ Δ q 2 ( q ̄ ) U ( q ̄ ) .

To illustrate the behavior of the CG-PI method, we consider here two model problems that represent, though in a simplified form, typical characteristics of physical potentials in realistic molecular systems. The CG-PI results are based on Eq. (27) with the physical particle mass used for m eff const . The first model problem is intended to imitate nuclear quantum effects in a condensed phase at low temperatures. In this problem, a particle is placed in an asymmetric anharmonic one-dimensional potential

V q = a q 2 + b q 3 + c q 4 ,
(36)

where the coordinate q is in angströms and V is in kcal/mol. The values of a, b, and c are 0.252, −0.570, and 0.814, respectively. This potential V(q) is essentially anharmonic, as demonstrated in Fig. 3(a) by a comparison to the corresponding harmonic potential

V HO q = a q 2 .
(37)

The temperature is taken to be 5 K, and the mass of the particle is 20.0 g/mol. (For more details on the motivation of this problem and the numerical techniques, see Sec. S6.1 in the supplementary material.32)

FIG. 3.

In the case of an anharmonic potential, the CG-PI method accounts for the major part of quantum effects. (a) The employed model anharmonic potential V(q), Eq. (36) (black solid line) and the corresponding harmonic potential VHO(q), Eq. (37) (blue dashed line) are shown. (b) The diagonal density matrix computed for the physical potential V(q) with the CG-PI method (black dotted line) is closer to the exact quantum diagonal density matrix for V(q) (thick blue line) not only in comparison to the exact classical result (thin red line) but also even in comparison to the exact quantum density matrix for the harmonic potential VHO(q) (blue dashed line).

FIG. 3.

In the case of an anharmonic potential, the CG-PI method accounts for the major part of quantum effects. (a) The employed model anharmonic potential V(q), Eq. (36) (black solid line) and the corresponding harmonic potential VHO(q), Eq. (37) (blue dashed line) are shown. (b) The diagonal density matrix computed for the physical potential V(q) with the CG-PI method (black dotted line) is closer to the exact quantum diagonal density matrix for V(q) (thick blue line) not only in comparison to the exact classical result (thin red line) but also even in comparison to the exact quantum density matrix for the harmonic potential VHO(q) (blue dashed line).

Close modal

The CG-PI approach—more specifically, the simplest version of it, with quadratic coupling between the centroid and observed quasiparticles, and in the assumption of constant masses of both quasiparticles, as in Eq. (27)—yields in this case an asymmetric non-Gaussian shape of the thermal density matrix as a function of the coordinate q (Fig. 3(b), black dotted line). This curve is much closer to the curve for the exact thermal density matrix (Fig. 3(b), thick blue line) than to the probability density for the same potential V(q) in classical statistical mechanics (thin red line), suggesting that the simplest version of the CG-PI method nearly completely accounts for the quantum effects in this model system. In particular, the CG-PI method precisely reflects broadening of the particle distribution in the quantum thermodynamic ensemble due to tunneling of the physical particle. It is important to note that the CG-PI result for this system is closer to the exact thermal density matrix for the potential V(q) than to the exact thermal density matrix for the potential VHO(q) (blue dashed line) over the most part of the q axis, except for q > 0.5 Å, confirming the statement made in Sec. II that a renormalization of the effective interaction between the centroid and observed quasiparticles allows for the major part of quantum corrections to anharmonic contributions to a thermal density matrix.

The second model problem is designed to imitate tunneling of a hydrogen atom through a potential energy barrier at temperatures typical of living organisms. A particle is placed in an asymmetric double-well one-dimensional potential

V q = c 0 + c 1 q + c 2 q 2 + c 3 q 3 + c 4 q 4 ,
(38)

where q is in angströms and V is in kcal/mol (Fig. 4(a)). The constants c0, c1, c2, c3, and c4 are 6, −0.5, −44, 2, and 88, respectively. The temperature is 310 K and the particle mass is 1.008 g/mol. (For more details, see Sec. S6.2 in the supplementary material.32)

FIG. 4.

In the case of a double-well potential with 5 kcal/mol barrier at 310 K, the CG-PI method quantitatively accounts for tunneling of a particle as light as a hydrogen atom. (a) The employed model potential V(q), Eq. (38). (b) The diagonal density matrix computed with the CG-PI method (black dotted line) is much closer to the exact quantum result (thick blue line) than to the exact classical result (thin red line).

FIG. 4.

In the case of a double-well potential with 5 kcal/mol barrier at 310 K, the CG-PI method quantitatively accounts for tunneling of a particle as light as a hydrogen atom. (a) The employed model potential V(q), Eq. (38). (b) The diagonal density matrix computed with the CG-PI method (black dotted line) is much closer to the exact quantum result (thick blue line) than to the exact classical result (thin red line).

Close modal

Under these circumstances, the particle exhibits quantum behavior, despite the high temperature and the presence of a reasonably high barrier (5 kcal/mol) to tunnel through. The diagonal values of the thermal density matrix predicted with the CG-PI method (Fig. 4(b), black dotted line) are much closer to the values of the exact thermal density matrix (Fig. 4(b), thick blue line) than to the probability density for the same potential V(q) in classical statistical mechanics (thin red line). Broadening of the particle distribution, which is particularly noticeable on the sides of the curves closer to the potential energy barrier, is quantitatively accounted for by the CG-PI approach.

The present CG-PI approach to imaginary time path integrals and thermal density matrices provides a new reductionist perspective on quantum statistical mechanics. This approach is based on a coarse-graining of Feynman’s imaginary time path integral, providing a simple way to analyze the most essential nature of quantum effects. As is demonstrated in this paper, the minimalist description employs only two quasiparticle coordinates: Q, the coordinate of the centroid particle that feels the physical potential, and q ̄ , the coordinate of a second quasiparticle that is used to evaluate the values of observables. The expectation value of a quantum operator in this representation equals an average of the classical counterpart of this operator over a classical thermodynamic ensemble. Unlike theories aiming at a representation of a quantum particle with only the centroid variable,31 the present model naturally avoids the problem of evaluating quantum operators in the centroid-only representation.

The main advantages of the CG-PI theory are the following: First, it is important from the conceptual viewpoint, since it defines a simple framework and provides a toolkit for studying quantum systems (especially near the classical limit) in a way maximally close to the well-developed and computationally efficient tools of classical statistical mechanics. For example, the use of two quasiparticles (the centroid and observed ones) to describe thermal quantum delocalization of an atom in a molecular system is a way to provide a minimalist approach for understanding quantum systems in a classical framework, including the necessary rules arising from the non-commuting position and momentum operators.

Next, using Eq. (26), one can directly run MD or MC simulations for quantum systems, employing only two quasiparticles for each physical quantum particle. Even in the simplest version of the method, namely, with the use of constant spring stiffness k Q q ̄ and constant effective mass meff, the results for a harmonic oscillator are exact at an arbitrary strength of quantum effects (that is, to all powers of ℏ). Renormalization of k Q q ̄ and, for momentum-dependent observables, of meff with changes in V ( q ̄ ) accounts in a compact, minimalistic way for the major part of quantum corrections especially near classical limit in the case of a general physical potential. The CG-PI approach also allows for adding higher order in ℏ corrections in a systematic manner by explicitly deriving higher order terms not shown explicitly in Eq. (15), as illustrated in the  Appendix, Eq. (A15), and Sec. S1 in the supplementary material.32 These corrections result in anharmonic terms in the effective potential of interaction between the centroid and observed quasiparticles, p4 and higher order terms in the effective kinetic energy in the Hamiltonian in Eq. (26), and additional terms in the expression for the effective mass meff. Another way of adding higher order correction may proceed via variational determination of the CG potential or by implementing various Ansätze to approximate and parameterize these terms. Likewise, the CG-PI theory can serve as a basis for the development of various approximate semi-phenomenological models of specific physical phenomena. In this case, an Ansatz expression for the CG potential can be chosen based on a qualitative analysis of the problem, while the specific values of the coefficients in such an expression can be related to the microscopic description of the system via a variational procedure. The computational cost of the CG-PI method is comparable to that of path integral molecular dynamics (PIMD) with only P = 2 quasiparticles per quantum particle, while the accuracy of the former is much higher, as discussed above. The use of the variable spring stiffness k Q q ̄ q ̄ , and possibly the effective mass m eff q ̄ , does not significantly increase the cost of computations, at least in the case of analytically defined potentials V(q), since the right hand sides of Eqs. (18) and (19) can be quickly computed with the use of standard computational techniques, such as look-up tables and piecewise interpolations.

Finally, the CG-PI approach can be used in a combination with the existing PIMD, path integral Monte Carlo (PIMC), etc., methods as a tool for processing and analyzing the data generated by these more exact numerical methods. Here, the advantages of coarse-graining are exhibited at the stage of the analysis and interpretation of exact numerical data. This mode of using the CG-PI theory is analogous to the use of CG models to analyze all-atom trajectory data of biomolecules (see, e.g., Ref. 37). Another potentially useful strategy of combining the CG-PI approach and the traditional PIMD and PIMC methods would be to use approximate CG-PI models, e.g., based on Eq. (15), to enhance the sampling in the PIMD or PIMC simulations.

Possible directions for the future development of the CG-PI theory will address the following questions: computations for benchmark quantum systems need to be performed with a subsequent comparison of the results with numerically exact PIMD results. Generalizations of the CG-PI approach to treat many-particle systems of bosons and fermions might also be worked out. There are also possible modifications of the CG-PI method based on various other CG approaches; the theory of ultra-coarse-graining (UCG)38,39 may prove useful for extending this approach, for example, in the description of tunneling of quantum particles between two or more deep local minima in terms of hops between different discrete states. The UCG theory is developed precisely for the situation in which the coarse-graining is very “aggressive” in going from many particles to just a few at a low CG resolution, as in the CG-PI.

To conclude, in this work, we suggest the most essential, minimalist picture of how quantum effects are manifested in a system at finite temperature and at equilibrium. The proposed CG-PI methodology may offer a new spectrum of tools for qualitative and quantitative investigation of quantum systems, especially near the classical limit such as molecular systems at room temperature.

This research was supported by the National Science Foundation (NSF Grant No. CHE-1465248). The authors are grateful to Professor David M. Ceperley, James F. Dama, Stephan Köhler, Dr. Yuxing Peng, and Dr. Frank X. Vazquez for discussing this project on different stages of its completion.

1. Exact CG potential for a harmonic oscillator

In the case of a one-dimensional harmonic oscillator with the potential energy V(q) given by

V q = m Ω 2 2 q 2 ,
(A1)

where m is the particle mass and Ω is the angular frequency of the oscillator, the exact CG potential VCG assumes the following form:

V CG q , Q , q = m Ω 2 Q 2 2 + 2 m ħ 2 β 2 f 2 tanh f f tanh f Q q ̄ 2 + m 2 ħ 2 β 2 f tanh f Δ q 2 + 1 2 β ln f tanh f sinh 2 f 2 f 4 ,
(A2)

where f = ħ βΩ/2. Note that Eq. (A2) coincides with Eq. (15) with the O term omitted, and f equals feff (q) defined by Eq. (20). For the proof of Eq. (A2), see Sec. S2.1 in the supplementary material.32 

2. General perturbative expansion for the CG potential

If the physical potential V(q) can be represented as

V q = V 0 q + λ δ V q
(A3)

such that the CG potential V CG , 0 q , Q , q for the reference potential V 0 q is known, then the CG potential for V(q) can be found in a form of a perturbative series in the small parameter λ. One variant of the perturbative series is given by

V CG q , Q , q = V CG , 0 q , Q , q + λ lim P 1 P i = 2 P δ V q i λ 2 β 2 lim P 1 P 2 i = 2 P j = 2 P δ V q i δ V q j δ V q i δ V q j + O λ 3 ,
(A4)

where the averaging is defined by

f q 2 , , q P = d q 2 d q P e β V P , 0 q , q 2 , , q P , q δ Q q 2 + + q P P 1 f q 2 , , q P d q 2 d q P e β V P , 0 q , q 2 , , q P , q δ Q q 2 + + q P P 1 .
(A5)

Alternatively, the perturbative series expansion for the CG potential can be written as

V CG q , Q , q = V CG , 0 q , Q , q + λ 0 1 d u δ V q + ( q q ) u + 2 π σ k = 1 a k k sin ( k π u ) λ 2 β 2 0 1 d u δ V 2 0 1 d u δ V 2 + O λ 3 ,
(A6)

where σ = ħ β / m , the averaging is defined by

f q , q ; a k = k = 1 d a k 2 π e k = 1 a k 2 / 2 β V ̄ 0 q , q ; a k δ Q q ̄ 2 2 π 2 σ k = 1 , 3 , 5 , a k k 2 f q , q ; a k k = 1 d a k 2 π e k = 1 a k 2 / 2 β V ̄ 0 q , q ; a k δ Q q ̄ 2 2 π 2 σ k = 1 , 3 , 5 , a k k 2
(A7)

and

V ̄ 0 q , q ; a k = 0 1 d u V 0 q + ( q q ) u + 2 π σ k = 1 a k k sin ( k π u ) .
(A8)

For the derivation of Eqs. (A4) and (A6), see Sec. S3 in the supplementary material.32 

3. Variational approach to determining the CG potential based on force-matching

In the case of an arbitrary physical potential V(q), the CG potential (up to an additive constant) can be found variationally,

V CG q , Q , q = arg min V CG trial χ 2 V CG trial + const ,
(A9)

where the variational functional χ2 is defined as

χ 2 V CG trial = lim P d q d q 2 d q P d q e β V P q , q 2 , , q P , q Σ q , q 2 , , q P , q d q d q 2 d q P d q e β V P q , q 2 , , q P , q
(A10)

and Σ is the sum of squares of deviations of the actual CG forces and the forces predicted from the trial potential V CG trial . For the exact definition of Σ and the proof of Eq. (A9), see Sec. S4.1 in the supplementary material.32 

4. Variational approach to determining the CG potential based on Gibbs-Bogolyubov-Feynman (GBF) inequality

An alternative variational approach can be based on a generalization of the GBF inequality to the context of the CG-PI theory. This generalization leads to

V CG q , Q , q = arg min V CG trial χ GBF V CG trial + g q q ,
(A11)

where g is an arbitrary function. Thus, the GBF approach can be used in the cases when the dependence of VCG on Δq is not important (e.g., for the computation of diagonal elements of the density matrix or expectation values of momentum-independent operators). The variational functional χGBF is defined by

χ GBF V CG trial = F trial + V CG q , Q , q V CG trial q , Q , q trial ,
(A12)

where Ftrial is the free energy corresponding to the trial CG potential, specifically

F trial = 1 β ln m π ħ 2 β d q d Q e β V CG trial q , Q , q ,
(A13)

and the averaging in Eq. (A12) is defined by

f q , Q , q trial = d q d Q e β V CG trial q , Q , q f q , Q , q d q d Q e β V CG trial q , Q , q .
(A14)

For the proof of Eq. (A11), see Sec. S4.2 in the supplementary material.32 

5. Third order terms in the small ħ expansion of the CG potential

The series expansion of VCG in terms of ℏ as a small parameter [see Eq. (15)] with the O(ħ3) terms explicitly shown assumes the following form:

V CG q , Q , q V Q + k Q q ̄ 2 Q q ̄ 2 + k Δ q 2 2 Δ q 2 + k 0 + κ 3 ( q ̄ ) 1 6 V q ̄ Q q ̄ 3 + φ 312 q ̄ V q ̄ Δ q 2 Q q ̄ + σ 2 φ 310 q ̄ V q ̄ Q q ̄ ,
(A15)

where σ = ħ β / m ,

κ 3 ( q ̄ ) = 5 2 f eff ( q ̄ ) 3 f eff ( q ̄ ) 2 tanh f eff ( q ̄ ) 9 15 f eff ( q ̄ ) 3 tanh 2 f eff ( q ̄ ) + 4 15 f eff ( q ̄ ) 2 tanh 3 f eff ( q ̄ ) f eff ( q ̄ ) tanh f eff ( q ̄ ) 3 ,
(A16)
φ 312 q ̄ = 1 24 f eff 3 q ̄ 3 f eff 2 q ̄ tanh f eff q ̄ + 6 f eff q ̄ 6 tanh f eff q ̄ f eff 2 q ̄ f eff q ̄ tanh f eff q ̄ ,
(A17)
φ 310 q ̄ = 3 f eff 2 q ̄ 1 + tanh 2 f eff q ̄ + f eff q ̄ 2 tanh 3 f eff q ̄ 21 tanh f eff q ̄ + 18 tanh 2 f eff q ̄ 48 f eff q ̄ tanh f eff q ̄ f eff q ̄ tanh f eff q ̄ 2 .
(A18)

Note that, in practice, the values of κ3, φ312, and φ310 can be computed as fast as the value of tanh f eff ( q ̄ ) by means of using reference tables and asymptotic expansions. For the derivation of these and some higher order corrections, see Sec. S1 in the supplementary material.32 

1.
R. P.
Feynman
,
Statistical Mechanics: A Set of Lectures
(
W. A. Benjamin
,
Reading, Massachusetts
,
1972
).
2.
P. O.
Löwdin
,
Rev. Mod. Phys.
35
,
724
(
1963
).
3.
A.
Kohen
and
J. P.
Klinman
,
Chem. Biol.
6
,
191
(
1999
).
4.
H.
Engel
,
D.
Doron
,
A.
Kohen
, and
D. T.
Major
,
J. Chem. Theory Comput.
8
,
1223
(
2012
).
5.
F.
Calvo
,
C.
Falvo
, and
P.
Parneix
,
J. Phys. Chem. A
118
,
5427
(
2014
).
6.
J.
Chen
,
X.
Ren
,
X. Z.
Li
,
D.
Alfe
, and
E.
Wang
,
J. Chem. Phys.
141
,
024501
(
2014
).
7.
P.
Durlak
and
Z.
Latajka
,
Phys. Chem. Chem. Phys.
16
,
23026
(
2014
).
8.
N.
Kungwan
,
Y.
Ogata
,
S.
Hannongbua
, and
M.
Tachikawa
,
Theor. Chem. Acc.
133
,
1553
(
2014
).
9.
O.
Marsalek
,
P.-Y.
Chen
,
R.
Dupuis
,
M.
Benoit
,
M.
Méheut
,
Z.
Bačić
, and
M. E.
Tuckerman
,
J. Chem. Theory Comput.
10
,
1440
(
2014
).
10.
S.
Wolf
,
E.
Curotto
, and
M.
Mella
,
Int. J. Quantum Chem.
114
,
611
(
2014
).
11.
K. Y.
Wong
,
Y.
Xu
, and
D. M.
York
,
J. Comput. Chem.
35
,
1302
(
2014
).
12.
A.
Kuki
and
P. G.
Wolynes
,
Science
236
,
1647
(
1987
).
13.
L.
Muhlbacher
,
J.
Ankerhold
, and
A.
Komnik
,
Phys. Rev. Lett.
95
,
220404
(
2005
).
14.
J. S.
Kretchmer
and
T. F.
Miller
3rd
,
J. Chem. Phys.
138
,
134109
(
2013
).
15.
F.
Calvo
,
P.
Parneix
, and
N. T.
Van-Oanh
,
J. Chem. Phys.
132
,
124308
(
2010
).
16.
N.
Faruk
,
M.
Schmidt
,
H.
Li
,
R. J.
Le Roy
, and
P. N.
Roy
,
J. Chem. Phys.
141
,
014310
(
2014
).
17.
P.
Sindzingre
,
D. M.
Ceperley
, and
M. L.
Klein
,
Phys. Rev. Lett.
67
,
1871
(
1991
).
18.
D. M.
Ceperley
,
Rev. Mod. Phys.
67
,
279
(
1995
).
19.
J.
Chen
,
X. Z.
Li
,
Q.
Zhang
,
M. I.
Probert
,
C. J.
Pickard
,
R. J.
Needs
,
A.
Michaelides
, and
E.
Wang
,
Nat. Commun.
4
,
2064
(
2013
).
20.
M. A.
Morales
,
J. M.
McMahon
,
C.
Pierleoni
, and
D. M.
Ceperley
,
Phys. Rev. Lett.
110
,
065702
(
2013
).
21.
C. M.
Herdman
,
A.
Rommal
, and
A.
Del Maestro
,
Phys. Rev. B
89
,
224502
(
2014
).
22.
R. P.
Feynman
and
A. R.
Hibbs
,
Quantum Mechanics and Path Integrals
(
McGraw-Hill
,
New York
,
1965
).
23.
D.
Chandler
and
P. G.
Wolynes
,
J. Chem. Phys.
74
,
4078
(
1981
).
24.
K. S.
Schweizer
,
R. M.
Stratt
,
D.
Chandler
, and
P. G.
Wolynes
,
J. Chem. Phys.
75
,
1347
(
1981
).
25.
M. F.
Herman
,
E. J.
Bruskin
, and
B. J.
Berne
,
J. Chem. Phys.
76
,
5150
(
1982
).
26.
B.
De Raedt
,
M.
Sprik
, and
M. L.
Klein
,
J. Chem. Phys.
80
,
5719
(
1984
).
27.
M.
Parrinello
and
A.
Rahman
,
J. Chem. Phys.
80
,
860
(
1984
).
28.
R. D.
Coalson
,
D. L.
Freeman
, and
J. D.
Doll
,
J. Chem. Phys.
85
,
4567
(
1986
).
29.
M. E.
Tuckerman
,
B. J.
Berne
,
G. J.
Martyna
, and
M. L.
Klein
,
J. Chem. Phys.
99
,
2796
(
1993
).
30.
G. A.
Voth
,
Coarse-Graining of Condensed Phase and Biomolecular Systems
(
CRC Press
,
Boca Raton
,
2009
).
31.
J.
Cao
and
G. A.
Voth
,
J. Chem. Phys.
100
,
5093
(
1994
).
32.
See supplementary material at http://dx.doi.org/10.1063/1.4929790 for complete derivations of the results presented in Sec. II. Technical details on the numerical simulations presented in Sec. III are given in the Supplementary Material, but not in the Appendix.
33.
W. G.
Noid
,
J.-W.
Chu
,
G. S.
Ayton
,
V.
Krishna
,
S.
Izvekov
,
G. A.
Voth
,
A.
Das
, and
H. C.
Andersen
,
J. Chem. Phys.
128
,
244114
(
2008
).
34.
S.
Izvekov
and
G. A.
Voth
,
J. Phys. Chem. B
109
,
2469
(
2005
).
35.
A. V.
Sinitskiy
and
G. A.
Voth
,
Chem. Phys.
422
,
165
(
2013
).
37.
M. G.
Saunders
and
G. A.
Voth
,
Annu. Rev. Biophys.
42
,
73
(
2013
).
38.
J. F.
Dama
,
A. V.
Sinitskiy
,
M.
McCullagh
,
J.
Weare
,
B.
Roux
,
A. R.
Dinner
, and
G. A.
Voth
,
J. Chem. Theory Comput.
9
,
2466
(
2013
).
39.
A.
Davtyan
,
J. F.
Dama
,
A. V.
Sinitskiy
, and
G. A.
Voth
,
J. Chem. Theory Comput.
10
,
5265
(
2014
).

Supplementary Material