Quantum monte carlo methods for constrained systems

REVIEW

WWW.Q-CHEM.ORG

Quantum Monte Carlo Methods for Constrained Systems Sarah Wolf,[a] Emanuele Curotto,*[a] and Massimo Mella[b] The torsional ground state for ethane, the torsional, rotational, and mixed torsional and rotational ground state of propane are computed with a version of diffusion Monte Carlo adapted to handle the geometric complexity of curved spaces such as the Ramachandra space. The quantum NVT ensemble average for the mixed torsional and rotational degrees of freedom of propane is computed, using a version of Monte Carlo path integral, also adapted to handle curved spaces. These three problems are selected to demonstrate the generality and the applicability of the approaches described. The spaces of coordinates can be

Introduction The time scale problem is a familiar one to quantum chemists. Rooted into our common basic training, the Born Oppenheimer approximation is the essential tool needed to unravel fundamental problems that would otherwise be intractable. In its most familiar form, it allows one to separate the electronic degrees of freedom, associated with a set of much lighter bodies, from the heavier nuclear degrees of freedom. With some notable exceptions of marked nonadiabatic behavior, the unraveling of which remains an very active field of research, the Born Oppenheimer approximation yields very accurate results, and it is often used to build accurate models for the phenomenological potential energy surface felt by the nuclei. The latter set of bodies is traditionally treated by classical mechanics. Generally, at sufficiently elevated temperatures, for sufficiently heavy nuclei, and for sufficiently harmonic interactions, classical mechanics provides reasonably accurate answers. In this article, however, we depart from the traditional theme and explore the time scale problem one encounters when by necessity the quantum laws of motion are applied to the nuclei themselves. This implies the nuclei are light, the temperature is too low, and the interactions are highly anharmonic. There is a vast number of fundamental problems in chemical physics, where all these conditions take place at the same time, and insight into these is fundamental to a myriad of disciplines ranging from astrophysics to computational biology. More specifically, the set of problems that have preoccupied our two groups intensely for a number of years are the theoretical estimation of physical properties of molecular clusters. Clusters in general and molecular clusters in particular are models of condensed matter that can be studied both theoretically and experimentally. Insight gained from these investigations has already created a vast improvement in our understanding of the complicated phenomena that take place in the assembly process, thermodynamic stability as function of size of condensed matter, microsolvation,

best constructed from the parameters of continuous Lie groups, and alternative methods based on vector spaces, where extended Lagrangian terms would be too cumbersome to implement. We note that the geometric coupling between the torsions and the rotations of propane produces a substantial effect on the ground state energy of propane, and that the quantum effects on the energy of propane are quite large even C 2014 Wiley Periodicals, Inc. well above room temperature. V DOI: 10.1002/qua.24647

adsorption versus absorption, the effects of surface tension just to name a few. At the temperature and pressure conditions that produce stable molecular aggregates, nuclei of elements in the first three periods are sufficiently light to produce significant quantum effects, whereas at the same time, the weak interaction between the molecules in the aggregate are far from harmonic. More importantly for the discussion in this article, most intramolecular degrees of freedom are associated with relatively deeper dissociation energies, as in the case of stretching modes, and relatively stiff force constants compared to the intermolecular degrees of freedom. Atomistic quantum simulation of molecular aggregates is rendered either particularly challenging or practically intractable by the large difference in the time scales. The classical equations of motion satisfied by a typical molecular aggregate are stiff. A small time step is required to sample properly the high frequency dynamics, whereas a long time scale is required to capture the effects on the system from the set of lower frequency events. Additionally, this problem exacerbates the lack of ergodicity and the occurrence of rare events when exploring thermodynamic properties in the classical limit with random walks. Therefore, judicious use of holonomic constraints is routine in classical simulations of molecular clusters. The typical outcome for thermodynamic properties is a significant reduction in the statistical error,

[a] S. Wolf, E. Curotto Department of Chemistry and Physics, Arcadia University, Glenside, Pennsylvania, 19038-3295 E-mail [email protected] [b] Massimo Mella Dipartimento di Scienza ed Alta Tecnologia, via Valleggio 11, Universit a degli studi dell’Insubria, 22100, Como, Italy Contract grant sponsor: ACS (Petroleum Research Fund; E.C.); contract grant number: 48146-B6. Contract grant sponsor: Universita degli Studi dell’Insubria via the fund Fondi di Ateneo per la Ricerca (FAR; M.M). C 2014 Wiley Periodicals, Inc. V

International Journal of Quantum Chemistry 2014, 114, 611–625

611

REVIEW

WWW.Q-CHEM.ORG

Figure 1. Ground state energy of a relatively stiff monodimensional harmonic oscillator for various values of the DMC step size. The population size is maintained around 105 replicas. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

making simulations more efficient. In quantum simulations, there is an added benefit gained by constraining stiff modes, at least for the two stochastic approaches we discuss in this article. The convergence properties of the diffusion Monte Carlo[1–12] (DMC) and the imaginary time path integral[13–31] (MCPI) are greatly enhanced by constraining high frequency modes. In a recent article,[25] we investigate a simple onedimensional (1D) harmonic chain with 1000 particles. All particles have the same mass, and every particle is connected to two neighboring particles on a line with one stiff harmonic spring on one side and a one soft one on the other side. The simple harmonic chain model constructed this way mimics a set of condensed molecules. Using analytical solutions of the imaginary time path integral at finite Trotter number, we compare the convergence of the analytical solutions of the path integral expression for the heat capacity as a function of temperature and Trotter number. Our results show that the adiabatic approximation is accurate for temperatures below 1 pffiffiffi ; 10 2bhxh

(1)

where xh is the smallest of the high frequency set. Below this temperature, a fully flexible simulation is highly inefficient as it requires many hundreds of time slices to converge, where the constrained simulations require less than 20 slices to converge for most of the temperature range. Consequently, the efficiency gains produced by performing the MCPI simulation with constraints are massive. Similar gains are quantifiable for DMC.[1] In Figure 1, we show the estimate of the ground state energy of a particle of unit mass in a harmonic potential V5k x 2 =2 with k510 a.u. as a function of simulation time for 612

International Journal of Quantum Chemistry 2014, 114, 611–625

various values of the step size Ds. This is the parameter that has to be systematically reduced until convergence is achieved in DMC simulations. To properly interpret the result in Figure 1, it is important to keep in mind that a similar computation with k51 a.u. is fully converged with a Ds50:02 a.u. (the largest value of Ds in the set of simulations represented). It is evident that an increase in frequency by a mere factor of 3.2 requires a DMC step smaller than a factor of 10. As is also evidenced in Figure 1, DMC simulations must also converge with respect to simulation pseudotime (i.e., the number of DMC steps, on the x-axis in Fig. 1) and a smaller Ds increases proportionally the total number of DMC steps required to reach the asymptotic distribution, that is, the ground state wavefunction in this case. When DMC is used to simulate systems with multiple time scales, the stiffer mode demands a small Ds, while the softer one demands a relative longer time scale to reach convergence, and consequently a much greater number of steps. Therefore, when in a system, the subdivision of stiff and soft modes is clear, and the adiabatic approximation is expected to work well, it is highly advantageous to use the proper holonomic constraints, and the efficiency gained can be on the same order of magnitude as what we measure in MCPI simulations. However, constraining intermolecular degrees of freedom, like bond stretches, bond angles, and the like, create curved spaces and special techniques have to be developed to handle the geometric complexity of these. The purpose of this article is to review some of the recent advances in the form of algorithms specifically designed to carry out quantum Monte Carlo simulations in curved spaces, and to demonstrate their powers. In the methodology section, we review briefly the basic objects needed for the geometry and the dynamics in general nonWWW.CHEMISTRYVIEWS.ORG

REVIEW

WWW.Q-CHEM.ORG

Euclidean spaces. Then, we develop the general theory for DMC and MCPI in non-Euclidean spaces. Our results section contains selected numerical examples and in the conclusions section, we discuss a set of directions that our group is currently undertaking to further develop our tools and continue to explore the rich and important field of condensed molecular matter.

Methodology The methods discussed in this article are based on four fields of modern applied mathematics and computation theory. These are differential geometry,[32,33] Lie algebra and Lie groups,[34] the Feynman path integral,[35–38] and Stochastic Calculus.[39] In this article, combinations of these theories are applied to quantum Monte Carlo simulations.[40–44] The books cited are excellent sources for those readers that wish to deepen their own understanding in these rich theoretical areas. Before working through the main theory for diffusion and path integral simulations in nonEuclidean spaces, we briefly introduce some of the fundamentals of differential geometry to clarify the notation used in the article. It is important to keep in mind that even Euclidean spaces can be mapped with curvilinear coordinates, and these give rise to the same tensorial objects introduced in the next section. Curvilinear coordinates are often chosen with the eventual purpose of exploiting symmetries (as e.g., an isotropic potential in R3 ), that eventually are used to identify cyclic coordinates, and simplify the treatment of systems. As the formalism of Euclidean spaces mapped with curvilinear coordinates takes us more than half way there, we start our discussion from curved spaces. However, there are numerous kinds of curved spaces, and the methods in this article have been developed to work only in orientable Riemann spaces with a positive definite signature, no torsion (vide infra), and that can be represented with a single coordinate map except for a set of points with zero measure.[39] In this article, we use the symbol Rn to notate an Euclidean space mapped with an orthogonal set of Cartesian coordinates. Non-Euclidean spaces Most n-dimensional spaces used in molecular physics can be derived with the aid of coordinate maps, U : x0

x

(2)

a set of n algebraic (typically nonlinear) equations that transform a set of coordinates x to a set x 0 in one or two main ways: i. They are subspaces of Rn1c with the introduction of c holonomic constraints, such as, for example, the twosphere (S2 ) mapped from R3 with polar angles h; /,

ii. They are the spaces of parameters of some continuous Lie group. We have specifically chosen examples that are best derived by the second method to demonstrate the power of the approaches contained in this article, and these are presented in detail in the results section. In all cases, the spaces must be sufficiently smooth and analytical, meaning the set of coordinates that make up the map U are everywhere continuous, have a finite first derivative and a symmetric Hessian derivative @ @ U5@l0 @m0 U5@m0 @l0 U; (5) @x l0 @x m0 at every point in the space. These features allow one to expand the coordinate map locally, that is, in an infinitesimal neighborhood of a mapped point, into a linear “Cartesian-like” map. Locally flat spaces of this nature are known as Riemann spaces and the local Cartesian-like coordinates are known as Riemann coordinates. The machinery of the chain rule is implemented in the analysis of the geometric features of the space, and a powerful notation using Greek subscripts and superscripts is commonly introduced to handle tensorial objects and their transformation properties. For example, common vectors are represented by their components. For example, V a is the component along the unit vector ea that form the basis for the vector space. The placement of the Greek index indicates the contravariant manner by which compo0 nents of vectors transform upon coordinate changes x l ! x l , 0

0

Vl 5

@x l l V : @x l

(6)

The notation in the last equation also includes an implied sum from one to the number of dimensions of the space every time the same Greek letter occurs in an expression in both the lower and upper position. One Greek index in the lower position indicates the covariant way in which components of one-form (such as derivatives of scalar functions @r f ) transform. The inverse of the Jacobian matrix elements in Eq. (6) are used instead. The process of categorizing covariant and contravariant behavior is generalized to ðn; n0 Þ tensors with n indexes in the contravariant position and n0 in the covariant position. A particularly useful (0,2) tensor, or two-form, is the metric tensor, gl0 ;m0 5

@x m @x l glm : @x m0 @x l0

(7)

This entity contains all the geometric information about the space. For example, the dot product in a non-Euclidean space is defined using the metric tensor

1=2

ðx 2 1y 2 Þ z 21 y ; /5tan x

h5tan 21

;

subject to the constraint x 2 1y 2 1z 2 5r 2 5 constant.

(3)

A  B5glm Al Bm 5Am Bm ;

(4)

where the right-hand most expression demonstrates how the metric tensor is used to lower the index of a vectorial quantity A and produce its one-form equivalent.

(8)

International Journal of Quantum Chemistry 2014, 114, 611–625

613

REVIEW

WWW.Q-CHEM.ORG

Am 5glm Al :

dV5

In a similar manner, the inverse of the metric tensor glm can be used to raise indexes. Other useful information such as connections between the degrees of freedom of the space, and the curvature of the space can be extracted from the derivatives of the metric tensor. A space with a symmetric tensor glm 5gml is a torsion free space, and a space that produces a positive definite metric tensor is a space with positive definite signature. In the molecular physics literature, the Hessian metric tensor (always symmetric and positive definite) is commonly known as the mass matrix, and this object additionally contains information about the effective mass of the system. The Hessian metric tensor is normally generated with Eq. (7) by transforming from a Cartesian space, where it is represented as a diagonal matrix with masses along the diagonal. The Hessian metric tensor is an essential object for the development of dynamic theories in manifolds. The classical Lagrangian of a system in non-Euclidean space, for example, is written as follows 1 L5 glm x_ l x_ m 2V: 2

(9)

The derivation of the laws of motion follows the same variational homotopy calculus as in Rn . The end result is the celebrated Euler–Lagrange equations, and these have exactly the same expression in all Riemann spaces, Euclidean or nonEuclidean alike. One interesting feature of the metric tensor is that the coordinates chosen to map a given space may not be orthogonal, and this is reflected by the coupling matrix element between the velocities x_ l and x_ m . Therefore, the laws of motion, classical or otherwise, couple the two degrees of freedom. To distinguish the latter types of couplings from those that occur, when the potential energy term contains product terms like x l x m we call them “geometric couplings,” while the traditional “dynamic couplings” are those produced by the potential energy. In most systems of interest, degrees of freedom are coupled in both manners, making it difficult to extract the subtle effects that one kind of coupling has on simulated physical properties relative to the other. Both types of coupling terms can be safely ignored whenever the difference in time scale between the two degrees of freedom is such that an adiabatic approximation is accurate. However, when the time scales are similar, significant contributions to the physical properties of the system can be expected.

pffiffiffiffiffiffi 1 jgjdx Ù dx 2   Ù dx n ;

(11)

pffiffiffiffiffiffi where jgj is the square root of the determinant of the metric tensor, and the symbol Ù notates the wedge product. The square root of the determinant of the metric tensor makes the totally antisymmetric density dV invariant under change of variables and is known as the Jacobian. The reader will have noticed that Fick’s first law produces a one form for the flux. Fick’s second law is a statement regarding the conservation of flux across an oriented boundary in M. In integral form, for a source and sink free conditions it translates to a statement of conservation of probability, ð

ð @P dm50: div J dm1 M M @t

(12)

The proper divergence expression for a (1,0) vector that satisfies the conservation laws in integral form is pffiffiffiffiffiffi  1 div X5 pffiffiffiffiffiffi @m jgjX m : jgj

(13)

The divergence takes this form because the derivative of nonscalar quantities in Riemann spaces does not transform covariantly. Rather, the covariant derivative of a (1,0) vector is given by, @l X m 1 Cmlk X k :

(14)

One defines the Christoffell connections of the second kind for Riemann spaces (Cmlk ), by forcing the expression in Eq. (14) to transform as a (1,1) tensor. Furthermore, it can be shown that the Christoffell connections can be evaluated using the following symmetric derivative of the metric tensor  1  Cmlb 5 gmq @l gbq 1@b gql 2@q glb : 2

(15)

Equation (15) is obtained by taking the derivative of glm and forcing it to transform as a (0,3) tensor under change of coordinates.[32,33] The expression in Eq. (13) is the trace of the object in Eq. (14), as notated by the repeated index in the upper and lower position. Additionally, as the flux is a (0,1) vector, its index is raised using the inverse of the metric tensor, pffiffiffiffiffiffi  1 div J5 pffiffiffiffiffiffi @m jgjglm Jl : jgj

(16)

Finally, inserting Fick’s first law, Eq. (10), on the right to eliminate J gives,

DMC in non-Euclidean spaces We now proceed to the task of deriving the diffusion equation in general Riemann spaces. Fick’s first law of diffusion defines the flux vector J as the derivative of Pðx; tÞ, the scalar probability of finding a diffusion particle between x and x1dV,

pffiffiffiffiffiffi  1 div J52D pffiffiffiffiffiffi @m jgjglm @l P : jgj

(17)

The operator on the right-hand side of Eq. (17) Jm 52D@m P;

(10)

where D is the diffusion coefficient, and x represents the set of generalized coordinates of the n-dimensional Riemann space M. dV is the differential volume element in M, 614

International Journal of Quantum Chemistry 2014, 114, 611–625

pffiffiffiffiffiffi 1 pffiffiffiffiffiffi @m jgjglm @l 5DLB ; jgj is known as the Laplace–Beltrami operator. WWW.CHEMISTRYVIEWS.ORG

(18)

REVIEW

WWW.Q-CHEM.ORG

€dinger equation, in imaginary The time-dependent Schro time (t5isÞ is isomorphic to the diffusion equation, with a source-sink term V2Vref where Vref is simply a guess for the ground state energy of the system.[1–12] In M it reads, 2

2 h @w : DLB w1ðV2Vref Þw5h @s 2

(19)

This differential equation can be simulated in general Riemann spaces [6–11] as well as its counterpart in Rn . The proof of this statement begins with the local expansion of the map U at a point x0 up the linear term. The resulting map, valid in the neighborhood of x0, is well defined at every point in M by the requirements for U that have been elaborated earlier. Then, the partial derivatives in Eq. (7) are constant (i.e., evaluated at x0). Consequently, the metric tensor is constant in the infinitesimal neighborhood of x0 as well, and the Laplace–Beltrami operator takes a relatively simple form, DLB ! glm ðx0 Þ@m @l :

(20)

From the local diffusion equation h

@w h2 glm ðx0 Þ !2 @m @l w1ðV2Vref Þw; @s 2

(21)

one can write a solution which is useful for the DMC approach. Branching is carried out the same way as in Rn , whereas the solution of the source-sink free diffusion equation can be used to sample random numbers that have hDsglm ðxk Þ as the correlation matrix. rffiffiffiffiffiffiffiffiffiffiffiffiffi   jgj 1 exp 2 glm ðxk ÞDx l Dx m : 2phDs 2hDs

(22)

The random numbers (Dx) in Eq. (22), are used to simulate diffusion for a particle from xk to xk11 , where xk is the point from which the map is expanded, rather than x0. We have explored several methods to sample the distribution in Eq. (22). The most general approach is to use the Cholesky decomposition of the correlation matrix, and the Box–Muller algorithm[45] together. The Box–Muller algorithm is used to generate n values for the set Dy, distributed as Gaussian variables with zero mean, no correlation, and with unit variance. Then, the Cholesky decomposition of the correlation matrix hDsglm ðxk Þ5LLT ;

(23)

produces the proper transformation from the Dy uncorrelated set to the correlated one,[46] Dx5LDy:

(24)

case where a non-Euclidean space has a constant metric tensor, as a modified branching procedure can be used to accelerate the convergence of the DMC method up to a third-order behavior.[5] This approach relies on the possibility of decomposing a diffusion step in two consecutive ones with half Ds either exactly (e.g., in Rn ), or within a well-specified order of error in Ds. In the case of the diffusion equation on the two-sphere, S2 , a particular choice of coordinate involving the geodesic distance allows one to implement a perturbation treatment[47] leading to a short time approximation of the diffusion Green’s function that is a third-order in Ds and fourth in the ratio between the diffusion displacement along a great circle and the sphere radius. Such approximation has been tested, and confirmed to have a cumulative second-order error in Ds, when computing dynamical observables such as the average diffused angle along a great circle after a chosen elapsed time nstep Ds (nstep  102100).[11] The need for testing the diffusive dynamics directly, rather than the steady state estimate of a few observables as commonly done while testing DMC, is rooted into the fact that the s ! 1 distribution for the diffusion in a curved space is a constant. Coupled with the commonly used symmetric splitting of the branching step based on the average of prediffusion and postdiffusion potential energy, the approximated Green’s function delivers a robust total second-order in the energy observables, which can be improved by means of an a posteriori extrapolation.[4] Additional usage of the algorithm just discussed, beside the initial application to O2Hen systems,[11] has indicated that its performance does not deteriorate even when dealing with highly quantum objects such as para H2 adsorbed on ammonia clusters.[48] Albeit not strictly necessary in the applications so far presented, importance sampling introduced in the branching step, or simulated with a drifting term via optimized guiding functions is possible in non-Euclidean spaces as well.[8] The gain in efficiency is comparable[2,3] to that obtained in Rn . Feynman path integral in non-Euclidean spaces The Feynman quantization for Riemann spaces was first formulated by DeWitt.[49] Its formalism is based on the following construct for the matrix element hxi ; ti jxi21 ; ti21 i of the time evolution operator,

1 d=2 21=4 1=2 21=4 hxi ; ti jxi21 ; ti21 i5 gi D ðxi ; ti jxi21 ; ti21 Þgi21 2pih (25)  i 3 exp Sðxi ; ti jxi21 ; ti21 Þ ; h where, gi is the determinant of the metric tensor evaluated at xi 2 M. The action Sðxi ; ti jxi21 ; ti21 Þ is defined as the time integral of the Lagrangian, and Dðxi ; ti jxi21 ; ti21 Þ is the Van Vleck determinant, Dlm 52

Second-order DMC on the two-sphere The diffusion procedure simulated using the aforementioned process is convergent to first-order in Ds, the time step, for the most general case. This situation can be improved in the special

@2S m @xil @xi21

:

(26)

Equations (25) and (26) pertain to one of the N time intervals used to subdivide the path t0 ; x0 ! tN ; xN . These equations can be used to formulate the regular time sliced approach. International Journal of Quantum Chemistry 2014, 114, 611–625

615

REVIEW

WWW.Q-CHEM.ORG

If the action S, and the Van Vleck determinant is expanded about xi21 ; ti21 (the so called prepoint expansion), up to firstorder in ðti 2ti21 Þ one obtains an approximate expression for hxi ; ti jxi21 ; ti21 i, 

 3exp

n=2 N 1=2 gi hxi ; ti jxi21 ; ti21 i  2piðti 2ti21 Þh

   i i l  m m xi 2xi21 2 ðti 2ti21 ÞVeff : glm xil 2xi21 2hðti 2ti21 Þ h (27)

The effective potential energy in this expression contains a quantum correction,[49] Veff 5V2

2 ðti 2ti21 Þ h R; 6

r5hb1=2 ;

(29)

(30)

International Journal of Quantum Chemistry 2014, 114, 611–625

~ k ðuÞ; alk K

(32)

t : bh

(33)

k5km 11

b5ðkB T Þ21 ;

u5

~ k ðuÞ so that the partial averagand constructs the functions K ing expansion about the core path derived from Eq. (32) is equal to the same derived with the infinite series. The require~ k ðuÞ is easily derived.[30] ment for K 0

km X

(31)

The resolution of the identity in curved spaces can be rigorously applied to derive the path integral expression using the same limiting procedure as in Euclidean spaces mapped by Cartesian coordinates. However, the short time approximation has to be obtained by expanding the propagator into a power series and keeping all the terms that are first-order in Dt. This produces the correction term in Eq. (28). The term proportional to the Riemann curvature scalar results from the proper transformation of the path integral measure.[39] Unless the Riemann curvature scalar is constant, it cannot be omitted from the expression of the path integral or the latter will not converge. Problems that require curved spaces with a nonuniform Riemann curvature can be quantized in a number of equally valid ways.[17,49] However, for the ellipsoids of inertia, and the Ramachandra spaces we explore in this work, the Riemann curvature is constant and these additional difficulties can be safely ignored. Furthermore, it is important to realize that as the Riemann curvature is a scalar quantity, it is inherently independent of the coordinate map chosen to reach points in these spaces. Path integral simulations of molecular aggregates do not abound in the chemical physics literature.[13–27] It is clear that the primitive approach is rather limited and that methods with accelerated convergence would be necessary.[28–31] Accelerating the convergence of the traditional time-sliced path integral simulations require either a first- or a second-order derivative of the potential. The gains obtained by reducing the number of slices N, must be weighted against the additional cost for 616

0

km X

alk Kk ðuÞ1r

where

and the Riemann–Cartan curvature scalar is obtained by contracting the Ricci tensor with the inverse of the metric tensor, R5glm Rlm :

km X k51

The Riemann curvature tensor contracts to the Ricci tensor, Rlm 5Rklkm ;

x l ðuÞ5x0l 1r

(28)

where R is the Riemann–Cartan curvature scalar a contraction[32,33] of the Riemann curvature tensor, Rqrlm 5@l Cqmr 2@m Cqlr 1Cqlk Ckmr 2Cqmk Cklr :

computing the derivatives. In multidimensional applications, especially when the potential energy models are sophisticated, the computation of the gradient, or the Hessian can in fact offset the efficiency gained using a second-order method. It is the combination of these factors that made the reweighted random series, coupled with the finite difference estimators, such powerful contributions[30] even in Rn . These approaches have been recently extended to simulations in Riemann spaces.[23–27] With the reweighted random series, one gains quadratic convergence and statistically stable estimators without the need to evaluate any derivatives of V. For the reweighted Fourier–Wiener path integral in imaginary time, one begins by 0 redefining the random path with km > km terms,

~ 2 ðuÞ5uð12uÞ2 K k

k5km 11

km X

K2k ðuÞ:

(34)

k51

For the Fourier–Wiener path integral, the path functions are rffiffiffiffiffi 2 sin ðkpuÞ ; Kk ðuÞ5 p2 k rffiffiffiffiffi ~ k ðuÞ5f ðuÞ 2 sin ðkpuÞ ; K p2 k

(35) (36)

and vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Pm 2 u 2 uuð12uÞ2 kk51 p2 k2 sin ðkpuÞ t f ðuÞ5 : Pk0m 2 2 k5km p2 k2 sin ðkpuÞ

(37)

Equation (32) is inserted into Eq. (27) with iðti 2ti21 Þ5

bh : N

(38)

Then, for a N point quadrature, and an n-dimensional Riemann space the density matrix qRW becomes, Nn=2  2 2n=2 1 h b JK q ðx; x ; bÞ5 2p   ð ð 1 pffiffiffiffiffiffi d ½ar jgjexp 2b duUðx ðuÞÞ ; RW

0

0

where JK is the Jacobian of the transformation xi ! ai , and WWW.CHEMISTRYVIEWS.ORG

(39)

REVIEW

WWW.Q-CHEM.ORG

1 Uðx ðuÞÞ5 glm x_ l x_ m 1Veff : 2

(40)

This algorithm, complete with the numerical difference estimators for the energy and heat capacity

 ð1  n @ 1 b duU ; 2b @b 0

 ð1  2 CV n n @ b duU 1 5 1 1nb @b kB 2 4 0 *  ð 2 +

2  ð1  1 @ @ b duU b duU 2 2b2 b2 @b @b2 0 0 

 ð1 2 n @ 2 2b b duU ; 2 @b 0 hEib 5

2. Carbon atom 2 is placed along the z-axis, its coordinates ðCÞ are given by the vector x2 5rC ez 5ð0; 0; rC Þ where ez is the unit vector along z, and the carbon–carbon average distance is rC 52:9101 bohr. 3. The position of the other carbon atoms is obtained by rotating and translating the coordinates of carbon 2,

(41) ðCÞ

(42)

pffiffiffi where, k varies from 1 to n 22, the angle c52tan 21 2 is the tetrahedral angle, and /k are the dihedral angles with /n 50. 4. The coordinates of the hydrogen atoms are constructed by performing a similar set of rotations and translations to the vector rH ez, where the carbon–hydrogen average distance is rH 52:0598 bohr. The first hydrogen atom on carbon 1 depends on /1 ,

produces the desired convergence properties of the energy and heat capacity. The derivatives,  ð1 @ b duU ; @b 0

ðCÞ

xk12 5rC Rz ð/1 ÞRx ðp2cÞ    Rz ð/k ÞRx ðp2cÞez 1xk11 ;

(43)

are evaluated numerically.

ðH Þ

x1 1 5rH Rz ð/1 ÞRx ðcÞez : The other two hydrogen atoms are found simply by changing the angle of rotation about the z-axis to /1 12p=3 and /1 14p=3. For the hydrogen atoms on carbon atom k52 to n22, the first set of rotation angles about z are augmented by 2p=3 and 4p=3, whereas the angle about x is always p2c.

2p ðCÞ Rx ðp2cÞez 1xk ; 5rH Rz ð/1 ÞRx ðp2cÞ    Rz /k21 1 3

4p ðH Þ ðCÞ Rx ðp2cÞez 1xk : xk 2 5rH Rz ð/1 ÞRx ðp2cÞ    Rz /k21 1 3 ðH1 Þ

Results

xk

Ramachandra spaces We have selected two molecular systems that are small enough to generate reproducible results in a short amount of time, but complicated enough to demonstrate the power and the applicability of the methods reviewed in the article. These are ethane and propane. In both the cases, it is more difficult to define a map using the first approach discussed earlier, that is, explicitly write down the equations of constraints. Rather, the use of Lie groups of continuous transformations creates a systematic approach to construct Cartesian coordinates from a well defined and physically meaningful parameter space, and this can be generalized and automated to produce linear chains of any dimension. The process we use to construct the coordinate map for a linear saturated hydrocarbon CnH2n12 , U : R3ð3n12Þ ! Tn21  I3 for ethane (n 5 2) and propane (n 5 3) is an adaptation of the algorithm outlined in Patriciu et al.[50] The vector space Tn21  I3 is the Cartesian product of the Ramachandra space of torsions and the inertia ellipsoid for the rotational degrees of freedom. The latter is typically mapped with three Eulerian angles using the Rz ð/ÞRx ðhÞRz ðwÞ convention, whereas Ramachandra’s space is typically mapped with n 2 1 dihedral angles. We first build a body-fixed Cartesian coordinate xðBF Þ representation, and then rotate it into a space-fixed representation, using the following procedure: ðCÞ

1. Carbon atom 1 is placed at the origin, x1 5ð0; 0; 0Þ.

For carbon n 2 1, the last rotation about z is by 12p=3 for one hydrogen atom and 14p=3 for the second hydrogen atom. For carbon n, the last rotation about z is by /n21 for one hydrogen atom, by /n21 12p=3, and by /n21 14p=3, respectively, for the other two. 5. The coordinates constructed in the previous four steps are translated to the center of mass xðBF Þ 5x2xCM : For ethane and propane, the center of mass is independent of the dihedral angles, and this fact simplifies the analysis. 6. The coordinate of each atom are operated upon by the rotations about the three Euler angles, xðSF Þ 5Rz ð/ÞRx ðhÞRz ðwÞxðBF Þ ; where, clearly, the body-fixed configuration only depends on the torsion angle(s), whereas the space-fixed coordinate set is a function of possibly all n 1 2 variables. The expressions for the body-fixed coordinates contain several fixed parameters, the dihedral angle c, the carbon–hydrogen distance rH, the carbon–carbon distance rC, and a Lie group parameter space of all permissible values of the dihedral angles /k . Equation (7) is used to transform a 3ð3n12Þ33 ð3n12Þ diagonal matrix with the proper masses along the diagonal that correspond to the element associated with the International Journal of Quantum Chemistry 2014, 114, 611–625

617

REVIEW

WWW.Q-CHEM.ORG

particular body-fixed Cartesian coordinate. If we arrange the variables so that x 1 5/1 ; x 2 5/2 ; . . ., up to n 2 1 and the last three variables for /; h, and w, then the lowest 333 block of the metric tensor has a standard general form[44] that depend on trigonometric functions of the Euler angles and elements of the inertia tensor, n X   3ði21Þ1j0 3ði21Þ1j00 glm 5 Clm j0 j00 mi xðBF Þ xðBF Þ : (44) i51

n  l; m  n12. The symbol Clm represents a set of nine distinct two-forms that operate in R3 associated with the bodyfixed Cartesian coordinates for atom i, 

Clm



j0 j00

 j 000 5djj000 @l R j0 ½@m Rjj00 :

(45)

Expressions for these 333 tensorial elements can be obtained readily.[44]

DMC simulations of ethane and propane For ethane, it is possible to obtain relatively simple analytical expressions for the body-fixed frame Cartesian coordinates as the center of mass obtained by the procedure in the previous section is simply rC =2. From these, it is straightforward to obtain analytical expressions for the inertia and the metric tensor. Ethane is a prolate top, and its inertia tensor is independent of /. ( )

2 1 1 Ixx 5Iyy 5 mc rc2 1mH 6 rH cos c2 rC 13rH2 sin 2 c ; 2 2

(46)

Izz 52lR2e :

(47)

lR2e 53mH rH2 sin 2 c:

(48)

where

The nonzero matrix elements of the metric are shown later. Only the elements on the main diagonal and above are displayed as the metric tensor is symmetric. g11 5lR2e

(49)

g12 5lR2e cos h

(50)

g14 5lR2e

(51)

g22 5Ixx sin 2 h12lR2e cos 2 h;

(52)

g24 52lR2e cos h;

(53)

g33 5Ixx ;

(54)

g44 52lR2e :

(55)

The symmetry Ixx 5Iyy has been used. We verify these expressions with mathematica,[51] and we obtain the same result, although the full simplification algorithm does not yield these expressions directly, rather we obtain expressions containing the masses of the two elements and the parameters rH and rC. 618

International Journal of Quantum Chemistry 2014, 114, 611–625

Nevertheless, it is straightforward to insert Eqs. (46) and (47) into Eq. (49) through (55) and verify the equivalences. For propane, the center of mass is independent of the two torsion angles /1 ; /2 , however, the expressions are rather complex, and inspection of analytic expressions is substantially less insightful. A similar observation is made for the inertia tensor. Propane is an asymmetric rotor, and the elements of the inertia tensor are constants, but in the body-fixed frame we have defined earlier the inertia tensor is not diagonal, rather a relatively small nonzero Iyz term in present. Therefore, the expressions contained in Eq. (44) do not simplify significantly. It is instructive to inspect the nonzero structure of the 535 matrix representation of the metric tensor nonetheless. 0 2 I3 @ gT 5 lm

2

gT lm

gClm

gClm

gIlm

3

1 A

(56)

The top 232 block related to the 2D Ramachandra space, 2 gT lm , is diagonal, and the two diagonal elements have the same expression, g11 5g22 5lR2e ;

(57)

where the right-hand side is the same as in ethane, given in Eq. (48). This fact implies that the two torsional degrees of freedom are not coupled with one another geometrically. As in ethane, the rotational degrees of freedom couple geometrically with the torsional ones. In the coupling block gClm , the expressions connecting /1 and the two Euler angles / and w, are identical to those in Eqs. (50) and (51), and the element connecting /1 with h vanishes. However, for the secondtorsional degree of freedom the coupling elements are more complex and all three Euler angles contribute nonvanishing coupling elements with /2 . 1 g23 5 lR2e (58) 3 pffiffiffi 2 2 2 lRe g24 5 (59) 3 n o pffiffiffi 1 g25 5 lR2e cos h22 2cos /sin h : (60) 3  3 The lower 333 block gIlm contains functions of the two Euler angles / and h, as in the most general case, and is independent of the two torsion angles. Therefore, from all these expressions, we note that Christoffel connection coefficients between the torsions and the rotations vanish. Nonvanishing Christoffel connection coefficients are another potential source of geometric couplings for more subtle but potentially significant effects. The Christoffel connection coefficients are present in the classical equations of motion derived from the general Lagrangian in Eq. (9). For the ensuing analysis, it is critical to realize that for propane, none of the degrees of freedom are dynamically coupled, as the potential energy for propane is simply a Fourier expansions[52] for the torsions and an external field, WWW.CHEMISTRYVIEWS.ORG

REVIEW

WWW.Q-CHEM.ORG

Figure 2. Convergence of the ground state energy for propane and ethane is shown by graphing the ground state energy estimate versus the step size used in the DMC simulation. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]



1 1 V5V3 11 cos 3/1 1 cos 3/2 1cos h ; 2 2

(61)

with V3 54:78083 mhartree as for ethane. The last term on the right-hand side is introduced to mimic the typical level of hindrance to rotation a propane molecule experiences on a surface. As a result, it is possible to gauge directly the effects that geometrical couplings have on the ground state energy and other important physical properties by considering both the sets of degrees of freedom separately, and then, compare the results with a simulation that includes all five degrees of freedom at once. We could have used ethane for this investigation, but we selected propane because the analysis we present here would be very cumbersome to carry out with the method of extended Lagrangian,[53] or with vector spaces.[54] For ethane, the simulation in the torsional space alone provides valuable insight and a convenient mean of comparing our DMC simulations results with a vector space computation of the ground state energy and wavefunction. For ethane, we use the following potential energy model[52] 1 V5 V3 ð11cos 3/Þ; 2

(62)

with V3 54:78083 mhartree. The results of four DMC simulations are graphed in Figure 2. We use the basic DMC algorithm,[1] amended only by the procedure for the diffusion steps, as elaborated by the discussion surrounding Eqs. (23) and (24). The population size is controlled by the usual feedback mechanism,[1] and maintained around 105 replica. At each value of Ds, we run 23 independent

simulations allowing 105 steps to reach equilibrium, and collecting averages for another 105 steps. Using the standard deviation in the mean from the 23 independent simulations, we construct the error bars shown in Figure 2. For ethane [panel(d)], we compare the ground state energy obtained from the DMC simulations over the potential surface in Eq. (62) against the lowest eigenvalue obtained by a vector space diagonalization. To carry out the latter, we construct a representation of the Hamiltonian in the vector space 1 wm 5 pffiffiffiffiffiffi exp ðim/Þ; 2p

m50; 61; 62; . . .

and use diagonalization to find the energy and eigenvectors of the system. The ground state is triply degenerate and with lR2e 5g11 520; 787:24 a.u. we obtain 0.4947397(7) mhartree with 601 basis sets as its energy. The error indicated in the energy is from the basis set convergence, and we estimated it from a separate computation using 801 basis sets. The ethane molecule in T1 is isomorphic to a particle in a ring with a mass of 3mH and a radius equal to rH sin cT . A step size of Ds  60 a.u. is sufficient to converge to the correct answer. The ground state energy for the ethane molecule in the space fixed frame is 0.494(7) mhartrees in excellent agreement with the vector space result. The statistical error is also relatively small compared to typical unguided DMC simulations. The graph in Figure 2a contains the outcome of the same simulation for the propane molecule in the body-fixed frame, namely in the 2D toroid space T2 . The potential energy model used is in Eq. (61) but for this simulation the h dependent term is left out. A step size of 50 a.u. produces a ground state International Journal of Quantum Chemistry 2014, 114, 611–625

619

REVIEW

WWW.Q-CHEM.ORG

energy of ET2 50:989ð2Þ mhartree. The figure in parenthesis is uncertain as the result of statistical fluctuations. The diagonalization results is obtained simply by doubling the ground state energy of ethane (0.989479402 mhartree), and the two estimates agree quantitatively with one another. We then perform two other sets of simulations for propane, using the same parameters, and varying only the step size. One simulation is performed by holding the two torsions rigidly at their minimum energy values (staggered), and using only the V3 cos h term from Eq. (61). The other, is performed in the 5D space T2  I3 . The convergence profiles with respect to Ds are in Figure 2b and 2c, respectively. From all the data in Figure 2, we make several observations. First, we note a convergence behavior for ethane in T1 and propane in T2 to an order in Ds grater than linear. The other two simulations of propane converge linearly instead. The reasons for the observed enhanced convergence behavior are as follows. For all our simulations, we make use of a secondorder branching expression,[4]    1   wi 5exp 2 V x 0i 1V ðxi Þ22Vref Ds : 2

(63)

The metric tensor in the torsional subspaces for ethane and propane is diagonal and independent of the values of the dihedral angles, therefore, a number of higher order terms that enter generally in the expression of the Green function propagator are absent. The homotopy of the space (the ring is a multiply connected set of spaces) and the boundary conditions that the propagator must satisfy, may contribute additional linear terms in general, but these are clearly too small to be observed in ethane and propane. All these conditions create the observed nonlinear behavior of the ground state energy with respect to Ds. Furthermore, systems with larger number of torsions create more complicated expressions for the metric tensor in Tn , and this object need not be constant. Therefore, in general the procedure outlined in this article is expected to converge linearly in a Ramachandra space of arbitrary dimension. €dinger equation by vector Second, the solution of the Schro space and diagonalization is considerably more involved for the propane molecule in both the I3 and the T2  I3 spaces. For the 3D case, the asymmetry of the inertia tensor alone requires a perturbation approach to handle the free rotation part, together with the Lanczos algorithm, and other sparse matrix strategies to handle angular momentum states greater than l > 14. The coupling block gClm has to be included in the pentadimensional treatment, and its contribution can only be handled by perturbation methods. Unfortunately, the relatively large size of the perturbation created by the gClm block is such that the pentadimensional propane problem is particularly challenging from the vector space prospective, whereas its DMC and MCPI simulations are not any more involved than the 3D and 2D counterparts. Last, we are able to quantify the excess energy caused by the coupling block gClm from the three ground state energies of propane, 620

International Journal of Quantum Chemistry 2014, 114, 611–625

Eexcess 5ET2 I3 2ET2 2EI3

(64)

Using ET2 I3 51:269ð9Þ, EI3 50:118ð6Þ, and ET2 50:989ð2Þ mhartree, we obtain an excess energy equal to Eexcess 50:162ð 1Þ mhartree. This value is clearly statistically significant and quite substantial (12.8% of ET2 I3 ). Remarkably, the geometric coupling affect the estimate of the ground state energy without contributing any explicit term in the estimator of the property as this is obtained from the population average potential energy, and the degrees of freedom in question are not dynamically coupled. Rather, the geometric coupling terms in this example affect the equilibrium population distribution via the correlation terms included in the sampling of the diffusion steps according to the procedure in Eqs. (22–24). MCPI simulations of propane For the stochastic evaluation of the Feynman path integral, there are difficulties in using angular variables. For instance, the random series expansion of the Brownian bridge becomes much more complicated to implement into an algorithm when periodic boundary conditions are used on the values of the angles. In trying to use the Euler angle h in conjunction with Eq. (32), the random variables (the coefficients of the path) must be constrained themselves to values that produce the correct range for the angle h. Fortunately, we have found a set of coordinates that can map spaces like the ellipsoid of inertia I3 , or the space T2  I3 for the example problem. These are the stereographic projections. The toroid spaces created by the dihedral angles can be easily remapped by defining nl, as follows, nl 5

2 cos /l 12sin /l

(65)

for all the n 2 1 dihedral angles. The transformation of variables from a dihedral angle /l to the projection nl is identical to that used for the particle in a ring of unit radius,[44] and for propane, with a mass equal to 3mH rH sin 2 c. The expressions in Eq. (65) follow from straightforward trigonometric identities. The transformation of variables from Euler angles to projections is slightly more involved.[44] One begins by defining a 4D space of quaternions, constrained to the surface of a three2 2 2 2 sphere with unit radius, that is, ðq1 Þ 1ðq2 Þ 1ðq3 Þ 1ðq4 Þ 51. The conversion map from Euler angles to quaternion coordinates that satisfies this equation of constraint is,

h w1/ q 5cos cos ; 2 2

h w2/ ; q2 5sin cos 2 2

h w2/ ; q3 5sin sin 2 2

h w1/ : q4 5cos sin 2 2 1

(66) (67) (68) (69)

Congruently with this result, we have shown that the onesphere result in Eq. (65) can be generalized to a n-sphere WWW.CHEMISTRYVIEWS.ORG

REVIEW

WWW.Q-CHEM.ORG

Figure 3. Reweighted random series simulation of propane in T2  I3 mapped with stereographic projection coordinates. The average energy as a function of b51=kB T is graphed for several values of km [cf. Eqs. (32–42)].

space or unit radius, defined as dlm ql qm 51, where dlm is the usual Kronecker delta and the sum is from 1 to n 1 1. The n stereographic projections are computed directly from, nl 52

ql 12qn11

l51; 2; . . . ; n

(70)

If we define hll0 as the partial derivatives involved in the map of I3, dh5h11 dn1 1h12 dn2 1h13 dn3 ;

(71)

d/5h21 dn1 1h22 dn2 1h23 dn3 ;

(72)

dw5h31 dn1 1h32 dn2 1h33 dn3 :

(73)

and we compute the equivalent quantities from Eq. (65), we have all the elements for the transformation of the metric tensor in Eq. (7). Explicit expressions for the nine hll0 , we define here can be found,[44] h11 52

16d 2 n1 ; d 1 ðr14Þ

h12 5

2d5 n2 d1 d2 ðr14Þ

h13 5

;

2d 5 n3 d 1 d2 ðr14Þ

;

(74) h21 5h31 5

5

d ; d4 h32 5

1 2

h22 5

3

8n n n 1 3; d4 d

8n1 n2 n3 2 3; d4 d

h33 5

1 3

h23 5

2

8n n n 2 3; d4 d

8n1 n3 n2 1 3: d4 d

(75) (76)

where the auxiliary expressions d1 through d5 are defined as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 d 1 5 16 n1 1ðr24Þ2 ;

(77)

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 d2 5 r2 n1 ;

(78)

 2 d3 5 d2 ;

(79)

 2  2 d 4 516 n1 1ðr24Þ2 5 d 1 : h  2 i d 5 54 2 n1 2ðr24Þ ;

(80) (81)

All these equations can be coded readily, and the computation of the metric tensor of propane in T2  I3 is fast. Therefore, the technical problems surrounding Eq. (32) for the MCPI simulations of propane in T2  I3 are eliminated, if we perform our random walk with the stereographic projections and the auxiliary random path coefficients. The results of several MCPI simulations on the pentadimensional propane molecule in an external field are presented in Figure 3. The quantum NVT ensemble average of the total energy in hartree [cf. Eq. (41)] is graphed as a function of b51=kB T for several values of the core coefficients km defined in Eq. (32). The parameter k0m regulates the length of the tail used to accelerate the convergence of the energy and heat capacity estimators to third-order. The nonlinear convergence pattern for hEib is evidenced by the data graphed in Figure 3, when comparing the relative distance between the tree sets of simulations. More importantly, we find that the quantum effects on the equilibrium thermodynamic properties of propane are quite large over a vast range of temperatures. At room temperature, for example, b  1050 hartree21, the km 5 10 simulation estimates an energy (17.39 mhartree) three times larger than the classical one (5.48 mhartree). The data in Figure 3 demonstrates that the classical limit is being reached asymptotically as b ! 0, and that the km 5 4 simulation is essentially converged for b < 250 hartree21, while the km 5 10 simulation International Journal of Quantum Chemistry 2014, 114, 611–625

621

REVIEW

WWW.Q-CHEM.ORG

(white squares) is converged for b > 2000 a.u. In fact, a km 5 14 simulation (black triangles) is statistically indistinguishable from the km 5 10 simulation (white squares) in the entire range of b in Figure 3. However, the quantum effects at b550 hartree21 ( 6300 K) are still visible on the scale of the graph. This observation alone demonstrates the importance of including quantum effects even when studying the gaseous state of propane.

Conclusions In this article, we review the details of two powerful quantum Monte Carlo techniques (DMC and MCPI) that our groups have adapted to generic manifolds. These adapted techniques find important applications in molecular physics, whenever it is sensible to separate time scales. The advantage of constraining high frequency degrees of freedom manifests itself in the convergence properties with respect to the imaginary time step of the methods, and the gains in efficiency can be massive. We have chosen two examples from molecular physics where the need for a curved space is evident. Both ethane and propane have stretching and bending modes with frequencies that are 4–5 times larger than those associated with the torsional degree of freedom, therefore, isolating the internal and external rotations from the rest of the normal modes is a natural choice. Assuming the potential energy of a system under investigation is known accurately, DMC and MCPI permit to introduce controllable approximations. For example, the finiteness in the time step in DMC, related to the Trotter number for MCPI, can be systematically varied until convergence is confirmed. A similar statement applies for DMC and MCPI in manifolds, whereby one essentially introduces one additional approximation that can be carefully controlled, namely the adiabatic approximation. For instance, in the range of temperatures considered in Figure 3, there are intervals where the bending and the stretching modes are no longer predominantly in the ground state[26] and these will contribute to the energy profile in Figure 3 on the left side. However, for these higher frequency degrees of freedom, a harmonic approximation[55,56] is valid for a good range of values of b, and this approximation too can be systematically corrected. Alternatively, and for even higher temperatures, path integrals treated using molecular dynamics and multiple time stepping techniques[57] can be used to simulate with greater efficiency. Notice that the range of temperatures where the approaches proposed in this work and in Ref. 57] may be both efficiently applied is expected to be rather narrow, as the adiabatic approximation is valid only when temperatures are low with respect to bending and stretching frequencies. The latter situation requires many “slices” in the ring polymers to correctly describe the associated degrees of freedom with the approach in Ref. 57], possibly making any comparison based on relative efficiency unfair. Also, the range of temperatures where the adiabatic approximation is valid, if at all, depends on the system, however, the number of problems in molecular physics that can be addressed with the two methods we have developed is truly vast given the large number of small molecules like ethane and propane, and the plethora of possible applications in academic 622

International Journal of Quantum Chemistry 2014, 114, 611–625

and industrial endeavors. For instance, let us mention the simulation of condensation/evaporation for gases used in the refrigeration technology, as well as in the study of stability and structure of aerosols. Although we have not “canned” our approaches into black box methods, it seems feasible at the moment. Besides, our experience has contributed to develop the idea that facing the grand challenges of introducing quantum effects in material science simulations will require juxtaposing the use of holonomic constraints and multiple time stepping techniques. This is the task, we are currently facing while working on simulating ammonia clusters[10] and droplets. The molecule of propane is simulated using a combination of curved spaces, the ellipsoid of inertia to treat its rotations relative to the laboratory frame, and the 2D toroid to handle the two torsions. The internal coordinate map is best generated from the parameter set of the continuous Lie groups of rotations and translations that are used to generate the configuration in Cartesian coordinates as elaborated in the results section. Using holonomic constraints is significantly more challenging for this particular problem given the number of constraining conditions. For propane, for example, one would have to write a total of 28 equations, and many of these are not trivial to express, especially those involving the angles. Therefore, the method of extended Lagrangian[53] is more challenging to implement. Also challenging for propane is the treatment of the pentadimensional space by vector space methods.[54] The impact that the geometric couplings between the torsions and the rotations has on the ground state energy is substantial for propane, and these terms cannot be ignored when constructing a vector space representation of the Hamiltonian operator. However, projecting the Laplace–Beltrami operator in a suitably chosen vector space is a nontrivial task for the asymmetric ellipsoid of inertia space alone, and the aforementioned coupling terms complicate the treatment substantially. We have demonstrated that DMC and MCPI can be carried out in Riemann spaces with positive signature, by simply introducing the proper expressions of the metric tensor, which is a unit matrix only in flat spaces mapped by Cartesian coordinates. However, simulations in manifolds force one to introduce the Jacobian factors in Eq. (11) as well. In MCPI, the Jacobian has to be included in the move strategy to maintain detailed balance. Additionally, the transformation of the Wiener measure[49] introduces a quantum correction to the potential energy in the form of the curvature scalar [cf. Eqs. (27) and (28)]. In DMC, the diffusion and branching process do not require explicitly Jacobian factors. Nevertheless, the interpretation of the population distribution for the unguided algorithm, for example, is that it approaches asymptotically the ground state wavefunction w0 , and in turn this object is related to a probability distribution P ðx Þ to locate a particle between x and x1dV, that is, P ðx Þ5 Ð

g1=2 jw0 j2

Mg

1=2 jw j2 dV 0

:

(82)

The Jacobian factor in Ramachandra’s spaces is still a matter of controversy, and it is this topic that motivated the work of Patriciu et al.[50] It is worth discussing the issue at this point. WWW.CHEMISTRYVIEWS.ORG

REVIEW

WWW.Q-CHEM.ORG

In their book,[43] Frenkel and Smit argue that the outcome of a simulation carried out with rigid modes should be the same as the result of a simulations where the high frequency degrees of freedom are treated with stiff springs. There is a fundamental difference between the spaces these two systems require. The former is typically a curved Riemann manifold M with all the geometric intricacies that have been discussed in this article. The latter is a curvilinear remapping of R3n ! M  N for n atoms with 3n-Cartesian coordinates. The metric tensors and     the Jacobi factors gM , and gMN associated with the two maps are inherently different. For DMC simulations, the difference between the two types of spaces only affects properties derived from the integral of the wavefunction, it does not impact, for example, the ground state energy. In MCPI simulations, however, the resolution of the identity is needed to derive the density matrix expression and this gives rise to the Jacobian factors in Eq. (11) and the Riemann curvature scalar[44] in Eq. (28). We have recently shown[44] that the ratio 

det gMN

However, the second-order methods for DMC with respect to the time interval Ds are only available for Rn and for the twosphere S2 . By simply expanding the map to second-order one generates all the terms in the Laplace–Beltrami operator that are missing in Eq. (20), and these can be incorporated into a second-order diffusion branching scheme in two ways: (a) by treating the terms neglected in the first-order scheme as drifting terms and (b) by making use of Ito–Taylor expansions for the metric tensor and its Cholesky decomposition. Another area that is under development currently in our laboratories is the adaptation of the ring polymer dynamics approach[58–66] to curved spaces. Briefly, in curved spaces one can write an expression for the Ring-polymer Hamiltonian for a n-bead system, where n is now the Trotter number, n n    X 1 lm 1 X l m Hn 5 g plj pmj 1 2 2 glm xjl 2xj21 xjm 2xj21 2 2bn h j51 j51 n X   1 V xj ;

1=2

ðdet gM Þ1=2

(84)

j51

(83)

and the curvature is generally a constant for all the types of ellipsoids of inertia, for the n2 sphere and for the monodimensional toroid. Therefore, the simulations in these types of spaces are in fact equivalent to those predicted by simulations with stiff modes in place of rigid constraints.[43] For the Ramachandra spaces the issue remains, even though we do calculate the Riemann curvature scalar for both ethane and propane and we find that it is constant. Patriciu et al. have argued in their work that simulations in the Ramachandra space should have a unit Jacobian factor.[50] They arrive at such conclusion by comparing the statistics of classical torsional Monte Carlo simulations with those obtained by dynamic simulations where the constraints are replaced by stiff bond lengths and bond angles. However, this argument does not take into consideration how quantum simulations differ in the two types of spaces, and the controversy is far from being put to rest. It is also clear that it would be difficult to resolve the matter experimentally as the zero-point energy contributed by the five degrees of freedom in propane is the same regardless of the space chosen M or M  N. Furthermore, we have made the general assertion that geometric and dynamic couplings can be neglected, whenever the adiabatic approximation is accurate between the stiff and the soft modes. With the tools discussed in this article one can, at least in principle, carefully verify such statements. In fact, a systematic study of this nature to provide general guidelines useful for the community is a goal for our group in the near future. The tools we have been adapting for applications in condensed molecular matter with the adiabatic approximations will continue to need refinement, and we take this opportunity to describe some of the future directions for further developments in this area. These endeavors are occupying our groups at the time of this writing. Higher order convergence for MCPI in manifolds have been developed, and we have implemented them for the propane problem in this article.

and where the second term on the right is the Matzubara harmonic term. The subscript j on each coordinate is used to identify each “bead.” We have recently developed an integrator in M using the variational principle but again avoiding the explicit use of holonomic constraints expressions, eliminating the need to use extended Lagrangian.[67,68] Rather one begins by setting to zero the variation of the action, ð d Ldt

(85)

where L is as in Eq. (9). When the integral is replaced by its discretized numerical equivalent, we obtain the discretized form of the equation of motion, m

m

1 x 2x m 1 x 2x m plk 5 glmk k11 k 1 glmk11 k11 k 2 2 Dt Dt r

m

r m Dt x 2x xk11 2xk Dt 2 @l grmk k11 k 1 @ l Vk 4 2 Dt Dt

(86)

and m

m

x 2x m x 2x m 1 1 plk11 5 glmk k11 k 1 glmk11 k11 k 2 2 Dt Dt r

m

r m x 2x xk11 2xk Dt Dt 1 @l grmk11 k11 k 2 @l Vk11 ; 4 2 Dt Dt

(87)

where in Eqs. (86) and (87) the notation xkl and glmk is an abbreviation of the respective quantities at time t5tk . Equation (86) represents a set of n independent coupled nonlinear l algebraic equations. The root xk11 depends, in general, on the values of all the other roots, and the solution for such nonlinear coupled system can only be found iteratively by evaluating l the metric tensor at an initial guess for xk11 and then solving for the latter from Eq. (86). Once the procedure is iterated to self consistency, the solution is used to compute the updates in momenta in Eq. (87). We are currently testing this numerical International Journal of Quantum Chemistry 2014, 114, 611–625

623

REVIEW

WWW.Q-CHEM.ORG

integrator in ellipsoids of inertia and we measure second-order convergence in the position as well as the energy. More importantly, we observe no significant drift in the energy over very long time scales. The details of this work will be published elsewhere. Last, this work on the smallest examples of Ramachandra’s spaces leaves us with a multitude of interesting questions and the answers to which can be potentially transforming for the community interested in protein modeling. Specifically, how do the quantum effects and the coupling between hindered rotation and torsion depend on the size of the molecule? One would expect the effects to get smaller as the mass associated with individual torsions grows. However, in the Ramachandra spaces, concerted motion of sets of torsional degrees of freedom take the system from one set of minima to another, and the effective mass involved does not necessarily need to be large. Quantum methods like DMC adapted for torsional spaces can yield not only the correct ground state but can also alleviate some of the difficulties regarding the lack of convergence to the ergodic limit for the state of the art stochastic methods available today.[69–89] This problem has plagued the work in the Ramachandra spaces from the outset. We have found that DMC alleviates quasiergodicity, and that for sufficiently large systems DMC alone can become trapped in local minima and become nonergodic. However, our groups have already found ways to improve the convergence of DMC toward the ergodic sampling limit.[12] Keywords: quantum Monte Carlo  Holonomic constraints  Ramachandra space  diffusion Monte Carlo  path integral Monte Carlo

How to cite this article: S. Wolf, E. Curotto, M. Mella. Int. J. Quantum Chem. 2014, 114, 611–625. DOI: 10.1002/qua.24647

[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]

624

J. B. Anderson, J. Chem. Phys. 1975, 63, 1499. J. B. Anderson, J. Chem. Phys. 1980, 73, 3897. M. H. Kalos, D. Levesque, L. Verlet, Phys. Rev. A 1974, 9, 2178. P. Ha˚kansson, M. Mella, D. Bressanini, G. Morosi, M. Patrone, J. Chem. Phys. 2006, 125, 184106. P. Ha˚kansson, M. Mella, J. Chem. Phys. 2007, 126, 104106. S. Chiesa, M. Mella, G. Morosi, D. Bressanini, J. Chem. Phys. 2003, 119, 5601. M. W. Aviles, E. Curotto, J. Phys. Chem. A 2007, 111, 2610. T. Luan, E. Curotto, M. Mella, J. Chem. Phys. 2008, 128, 164102. E. Asare, A.-R. Musah, E. Curotto, D. L. Freeman, J. D. Doll, J. Chem. Phys. 2009, 131, 184508. E. Curotto, M. Mella, J. Chem. Phys. 2010, 133, 214301. M. Mella, J. Chem. Phys. 2011, 135, 114504. K. Roberts, R. Sebsebie, E. Curotto, J. Chem. Phys. 2012, 136, 074104. R. A. Kuharski, P. J. Rossky, J. Chem. Phys. 1985, 82, 5164. D. Marx, M. H. Muser, J. Phys.: Condens. Matter 1999, 11, R117. D. Marx, P. Nielaba, Phys. Rev. A 1992, 45, 8968. T. F. Miller, III, D. C. Clary, J. Chem. Phys. 2003, 119, 68. L. Kaplan, N. T. Maitra, E. J. Heller, Phys. Rev. A 1997, 56, 2592. G. Guillon, T. Zeng, P.-N. Roy, J. Chem. Phys. 2013, 138, 184101. G. Guillon, T. Zeng, P. N. Roy, J. Chem. Phys. 2013 139, 184115. M. F. Russo, Jr., E. Curotto, J. Chem. Phys. 2003, 118, 6806. M. F. Russo, Jr., E. Curotto, J. Chem. Phys. 2005, 120, 2110.

International Journal of Quantum Chemistry 2014, 114, 611–625

[22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] [64] [65] [66] [67] [68] [69] [70] [71]

M. W. Avil es, E. Curotto, J. Chem. Phys. 2005, 122, 164109. E. Curotto, J. Chem. Phys. 2005, 123, 134102. M. W. Avil es, P. T. Gray, E. Curotto, J. Chem. Phys. 2006, 124, 174305. S. F. Langley, E. Curotto, D. L. Freeman, J. D. Doll, J. Chem. Phys. 2007, 126, 084506. E. Curotto, D. L. Freeman, J. D. Doll, J. Chem. Phys. 2008, 128, 204107. M. W. Avil es, M. L. McCandless, E. Curotto, J. Chem. Phys. 2008, 128, 124517. J. D. Doll, R. D. Coalson, D. L. Freeman, Phys. Rev. Lett. 1985, 55, 1. R. D. Coalson, D. L. Freeman, J. D. Doll, J. Chem. Phys. 1986, 85, 4567. C. Predescu, D. Sabo, J. D. Doll, D. L. Freeman, J. Chem. Phys. 2003, 119, 12119. S. Kunikeev, D. L. Freeman, J.D. Doll, Int. J. Quantum Chem. 2009, 109, 2916. W. D. Curtis, F. R. Miller, Differential Manifolds and Theoretical Physics; Academic Press: San Diego, 1985. B. F. Schutz, Geometrical Methods of Mathematical Physics; Cambridge University Press: Cambridge 1980. R. Gilmore, Lie Groups, Lie Algebras, and Some of Their Applications; Dover: New York, 2006. R. P. Feynman, Rev. Mod. Phys. 1948, 20, 367. R. P. Feynman, H. R. Hibbs, Quantum Mechanics and Path Integrals; McGraw-Hill: New York, 1965. L. S. Schulman, Techniques and Applications of Path Integration; Wiley: New York, 1981. H. Kleinert, Path integrals in Quantum Mechanics, Statistics and Polymer Physics; World Scientific: Singapore, 1990. K. Ito, Lectures on Stochastic Processes; Tata Institute of Fundamental Research: Bombay, India, 1960. N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. M. Teller, E. Teller, J. Chem. Phys. 1953, 21, 1087. M. H. Kalos, P. A. Whitlock, Monte Carlo Methods; Wiley: New York, 1986. M. P. Allen, D. J. Tildesley, Computer Simulations of Liquids; Claredon Press: Oxford, 1987. D. Frenkel, B. Smit, Understanding Molecular Simulations, 2nd ed.; Academic Press: New York, 2002. E. Curotto, Stochastic Simulations of Clusters: Quantum Methods in Flat and Curved Spaces; CRC: Boca Raton, FL, 2010. W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipes in FORTRAN, 2nd ed.; Cambridge Press: Cambridge, 1992. S. T. Li, J. L. Hammond, IEEE Tran. Syst. Man Cybern. 1975, SMC-5, 557. J. Faraudo, J. Chem. Phys. 2002, 116, 5831. M. Mella, E. Curotto, J. Chem. Phys. 2013, 139, 124319. B. S. DeWitt, Rev. Mod. Phys. 1957, 29, 377. A. Patriciu, G. S. Chirikjian, R. V. Pappu, J. Chem. Phys. 2004, 121, 12708. Wolfram Research, Inc., Mathematica, Version 8.0; Champaign, IL, 2010. G. Kamiski, E. M. Duffy, T. Matsui, W. L. Jorgensen, J. Phys. Chem. 1994, 98, 13077. A. Sarsa, K. E. Schmidt, J. W. Moskowitz, J. Chem. Phys. 2000, 113, 44. V. Buch, J. Chem. Phys. 1992, 97, 726. F. Calvo, C. Falvo, P. Parneix, J. Chem. Phys. 2013, 138, 034305. E. Balog, A. L. Hughes, G. J. Martyna, J. Chem. Phys. 2000, 112, 870. M. E. Tuckerman, B. J. Berne, G. J. Martyna, M. L. Klein, J. Chem. Phys. 1993, 99, 2796. I. R. Craig, D. E. Manolopoulos, J. Chem. Phys. 2004, 121, 3368. A. P erez, M. E. Tuckerman, M. H. Muser, J. Chem. Phys. 2009, 130, 184105. A. Yoshimori, J. Chem. Phys. 2008, 128, 234105. S. Habershon, B. J. Braams, D. E. Manolopoulos, J. Chem. Phys. 2007, 127, 174108. B. J. Braams, D. E. Manolopoulos, J. Chem. Phys. 2006, 125, 124105. T. D. Hone, P. J. Rossky, G. A. Voth, J. Chem. Phys. 2006, 124, 154103. R. CraigI, D. E. Manolopoulos, J. Chem. Phys. 2005, 123, 034102. A. Horikoshi, K. Kinugawa, J. Chem. Phys. 2005, 122, 174104. S. Wolf, E. Curotto, J. Chem. Phys. 2012, 137, 014109. J. P. Ryckaert, G. Ciccotti, H. J. C. Berendsen, J. Comput. Phys. 1977, 23, 327. H. C. Andersen, J. Comput. Phys. 1983, 52, 24. D. D. Frantz, J. Chem. Phys. 2001, 115, 6136. I. Andricioaei, J. E. Straub, A. F. Voter, J. Chem. Phys. 2001, 114, 6994. P. Nigra, D. L. Freeman, J. D. Doll, J. Chem. Phys. 2005, 122, 114113.

WWW.CHEMISTRYVIEWS.ORG

REVIEW

WWW.Q-CHEM.ORG

[72] V. Sharapov, D. Meluzzi, V. Mandelshtam, Phys. Rev. Lett. 2007, 98, 105701. [73] I. R. McDonald, K. Singer, Discuss. Faraday Soc. 1967, 43, 40. [74] C. Bichara, J. Gaspard, J. C. Mathieu, Phys. Lett. A 1987, 119, 462. [75] A. M. Ferrenberg, R. H. Swendsen, Phys. Rev. Lett. 1988, 61, 2635. [76] A. M. Ferrenberg, R. H. Swendsen, Phys. Rev. Lett. 1989, 63, 1195. [77] P. Labastie, R. L. Whetten, Phys. Rev. Lett. 1990, 65, 1567. [78] D. D. Frantz, D. L. Freeman, J. D. Doll, J. Chem. Phys. 1990, 93, 2769. [79] G. M. Torrie, J. P. Valleau, Chem. Phys. Lett. 1974, 28, 578. [80] I. Andricioaei, J. E. Straub, J. Chem. Phys. 1997, 107, 9117. [81] Y. Pak, M. S. Wang, J. Chem. Phys. 1999, 111, 4359. [82] R. Zhou, B. Berne, J. Chem. Phys. 1997, 107, 9185. [83] M. Falcioni, M. W. Deem, J. Chem. Phys. 1999, 110, 1754. [84] J. P. Neirotti, F. Calvo, D. L. Freeman, J. D. Doll, J. Chem. Phys. 2000, 112, 10340.

[85] F. Calvo, J. P. Neirotti, D. L. Freeman, J. D. Doll, J. Chem. Phys. 2000, 112, 10350. [86] J. P. Neirotti, D. L. Freeman, J. D. Doll, Phys. Rev. E 2000, 62, 7445. [87] C. Predescu, M. Predescu, C. V. Ciobanu, J. Chem. Phys. 2004, 120, 4119. [88] D. Sabo, M. Meuwly, D. L. Freeman, J. D. Doll, J. Chem. Phys. 2008, 128, 174109. [89] N. Plattner, J. D. Doll, P. Dupuis, H. Wang, Y. Liu, J. E. Gubernatis, J. Chem. Phys. 2011, 135, 134111.

Received: 19 November 2013 Revised: 31 January 2014 Accepted: 4 February 2014 Published online 28 February 2014

International Journal of Quantum Chemistry 2014, 114, 611–625

625