Enhancing Quantum Monte Carlo sampling - Methods ...

Enhancing Quantum Monte Carlo sampling Methods in manifolds, and a Smart Darting approach

E. Curotto, S. E. Wolf,1 1 Arcadia 2 Universita‘

and M. Mella2

University

Dell’ Insubria, Como

Pacifichem, 2015

E. Curotto, S. E. Wolf,, and M. Mella

Enhancing Quantum Monte Carlo sampling

Pacifichem, 2015

1 / 32

Outline

1

RPMD in manifolds A Leap Frog - like integrator A rigid ammonia molecule in an external field

2

DMC in manifolds The time - dependent Schr¨ odinger equation The Langevin equation DMC simulations of Ethane, and propane

3

Smart Darting and DMC Some facts about Smart darting A validation test: Multiple double wells



Pacifichem, 2015

2 / 32

Outline

1


2


3




Pacifichem, 2015

3 / 32

A Leap Frog - like integrator

Using the discrete version of the Euler-Lagrange equations, D2 Ld (qk−1 , qk , ∆t) + D1 Ld (qk , qk+1 , ∆t) = 0,

(1)

For Langrangians of the type 1 L(q, q, ˙ t) = gµν q˙ µ q˙ ν − V 2 Where,

0

gµν =

(2)

0

∂x µ ∂x ν gµ0 ν 0 . ∂q µ ∂q ν

(3)

is the metric tensor.



Pacifichem, 2015

4 / 32


The algorithm becomes, pµ k = pµ k+1 =

1 ∆t ∆t (gµν k + gµν k+1 ) γkν − ∂µ gσν k γkσ γkν + ∂µ Vk , 2 4 2

(4)

1 ∆t ∆t (gµν k + gµν k+1 ) γkν + ∂µ gσν k+1 γkσ γkν − ∂µ Vk+1 , (5) 2 4 2

where,

µ qk+1 − qkµ ∆t is the average velocity vector at time tk .

γkµ =



Pacifichem, 2015

5 / 32


µ Solve for qk+1 from eq. (4) and insert into eq. (5) to update the momentum pµ k+1



Pacifichem, 2015

6 / 32



Integrates using the generalized coordinates directly.



Pacifichem, 2015

6 / 32



Integrates using the generalized coordinates directly. Makes no use of Lagrange multipliers.



Pacifichem, 2015

6 / 32



Integrates using the generalized coordinates directly. Makes no use of Lagrange multipliers. Is quadratically convergent in ∆t



Pacifichem, 2015

6 / 32



Integrates using the generalized coordinates directly. Makes no use of Lagrange multipliers. Is quadratically convergent in ∆t Reduces to the momentum version of the Leap Frog algorithm if gµν = mδµν



Pacifichem, 2015

6 / 32



Integrates using the generalized coordinates directly. Makes no use of Lagrange multipliers. Is quadratically convergent in ∆t Reduces to the momentum version of the Leap Frog algorithm if gµν = mδµν An Iterative solution is necessary if gµν is configuration dependent.



Pacifichem, 2015

6 / 32



Integrates using the generalized coordinates directly. Makes no use of Lagrange multipliers. Is quadratically convergent in ∆t Reduces to the momentum version of the Leap Frog algorithm if gµν = mδµν An Iterative solution is necessary if gµν is configuration dependent. However, the system to be solved can be substantially smaller than SHAKE and RATTLE.



Pacifichem, 2015

6 / 32

Outline

1


2


3




Pacifichem, 2015

7 / 32

A rigid ammonia molecule in an external field

10

n=1 n=3 n=8 n = 10

2

(bohr )

8

~ C (t)

6

4

2

0

5000


10000

15000 t (a.u.)

20000


25000

30000

Pacifichem, 2015

8 / 32

A rigid ammonia molecule in an external field

n=1

0.06 0.04 0.02 0

0.06 n=3

~ f( ω )

0.04 0.02 0.2

0 n = 10

0.15 0.1 0.05 0

0

1000


2000

~ ω

3000

4000

5000

-1

(cm )


Pacifichem, 2015

9 / 32

Outline

1


2


3




Pacifichem, 2015

10 / 32

Diffusion in manifolds 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 x + dV, Jν = −D∂ν P, (6) Jν is the flux vector, P(x, t), the probability of finding a diffusion particle between x, and x + dV Fick’s second law on manifolds reads p 1 |g |g µν ∂µ P . (7) div J = −D p ∂ν |g | The operator on the right hand side of eq. (7) p 1 p ∂ν |g |g µν ∂µ = ∆LB , |g |

(8)

is known as the Laplace-Beltrami operator. E. Curotto, S. E. Wolf,, and M. Mella


Pacifichem, 2015

11 / 32

The time - dependent Schr¨odinger equation In M the TDSE is, −

~2 ∂ψ ∆LB ψ + (V − Vref ) ψ = ~ . 2 ∂τ

(9)

It can also be written as an advection - diffusion equation, ~2 µν ∂ψ g ∂µ ∂ν ψ + F µ ∂µ ψ − (V − E0 ) ψ = −~ . 2 ∂t

(10)

The first term on the left and the advection vector F µ is attainable directly by differentiating eq. (8) once through, Fµ =


i ~2 h (∂ν g µν ) + g µν Γλλν . 2


(11)

Pacifichem, 2015

12 / 32

Outline

1


2


3




Pacifichem, 2015

13 / 32

The Langevin equation

The diffusion part of the algorithm is derived from the theory of Stochastic Differential equations, χµn+1 = χµn + σνµ Wnν + F µ ∆t +

1 λ σκ ∂λ σνµ (Wnκ Wnν + Ωκν n ) 2

where, Ωλν represents a set of two-point distributed random variables  Ωλν = ∆t if ν < λ  1/2 λν P Ωn = (12)  λν 1/2 Ω = −∆t if ν > λ,



Pacifichem, 2015

14 / 32


χµn+1 = χµn + σνµ Wnν + F µ ∆t +


Wnν is a Gaussian variate with zero mean and ∆t variance



Pacifichem, 2015

15 / 32




Wnν is a Gaussian variate with zero mean and ∆t variance σνµ is the “square root” of the inverse of the metric tensor



Pacifichem, 2015

16 / 32




Wnν is a Gaussian variate with zero mean and ∆t variance σνµ is the “square root” of the inverse of the metric tensor The “drift term” F µ ∆t is important only if the mass is sufficiently light and if the potential energy is sufficiently shallow for the wavefunction to “feel” the boundaries of the manifold.



Pacifichem, 2015

17 / 32




Wnν is a Gaussian variate with zero mean and ∆t variance σνµ is the “square root” of the inverse of the metric tensor The “drift term” F µ ∆t is important only if the mass is sufficiently light and if the potential energy is sufficiently shallow for the wavefunction to “feel” the boundaries of the manifold. The terms

1 λ σκ ∂λ σνµ (Wnκ Wnν + Ωκν n ) 2 are necessary for the diffusion part of the algorithm to converge quadratically in ∆t. E. Curotto, S. E. Wolf,, and M. Mella


Pacifichem, 2015

18 / 32


The method was shown to be quadratically convergent in general if second order branching is used This approach has been successfully tested only on bosonic ground states so far, but the group is trying to find a fermion - like manifold to gain insight. Methods for DMC using extended Lagrangians have been developed [A. Sarsa, K. E. Schmidt, and J. W. J. Moskowitz, J. Chem. Phys. 113, 44 (2000).] V. Buch pioneered the first order version of this approach in 1992 on the H2 -H2 O clusters [ V. Buch, J. Chem. Phys. 97, 726 (1992).]



Pacifichem, 2015

19 / 32

Outline

1


2


3




Pacifichem, 2015

20 / 32

DMC simulations of Ethane, and propane

For ethane,

1 V = V3 (1 + cos 3φ) , 2

(13)

1 1 V = V3 1 + cos 3φ1 + cos 3φ2 + cos θ , 2 2

(14)

For propane,

with V3 = 4.78083 mhartree.



Pacifichem, 2015

21 / 32


< V > hartree

0.00115

0.00019

(a) propane (torsions)

0.0011

0.00018

0.00105

0.00017

Diagonalization

0.001

0.00016

0.00095 0.0009

0

500

1000

1500

2000

0.00129

< V > (hartree)

(b) propane (rotations)

0.00015 0.00014

0

100

200

300

400

0.0005

0.00126

Diagonalization

0.00048 0.00123 0.00046 0.0012 0.00044

(c) propane (rotations + torsions)

0.00117 0

100

200 300 ∆τ (a.u.)


400

500

0.00042

(d) ethane (torsion) 0


500

∆τ (a.u)

1000

1500

Pacifichem, 2015

22 / 32




Pacifichem, 2015

23 / 32

Outline

1


2


3




Pacifichem, 2015

24 / 32

Some facts about Smart darting

The original approach is called Smart Walking. [R. Zhou, and B. J. Berne, J. Chem. Phys. 107,9185 (1997).]



Pacifichem, 2015

25 / 32


The original approach is called Smart Walking. [R. Zhou, and B. J. Berne, J. Chem. Phys. 107,9185 (1997).] Smart Walking is modified to Smart Darting (SD) so that detailed balance is satisfied for the moves [I. Andricioaei, J. E. Straub, and A. F. Voter, J. Chem. Phys. 114, 6994 (2001).]



Pacifichem, 2015

25 / 32


The original approach is called Smart Walking. [R. Zhou, and B. J. Berne, J. Chem. Phys. 107,9185 (1997).] Smart Walking is modified to Smart Darting (SD) so that detailed balance is satisfied for the moves [I. Andricioaei, J. E. Straub, and A. F. Voter, J. Chem. Phys. 114, 6994 (2001).] A Jacobian for SD moves is necessary for Boltzmann integrals. [ P. Nigra, D. L. Freeman, and J. D. Doll, J. Chem. Phys. 122, 114113 (2005).]



Pacifichem, 2015

25 / 32


The original approach is called Smart Walking. [R. Zhou, and B. J. Berne, J. Chem. Phys. 107,9185 (1997).] Smart Walking is modified to Smart Darting (SD) so that detailed balance is satisfied for the moves [I. Andricioaei, J. E. Straub, and A. F. Voter, J. Chem. Phys. 114, 6994 (2001).] A Jacobian for SD moves is necessary for Boltzmann integrals. [ P. Nigra, D. L. Freeman, and J. D. Doll, J. Chem. Phys. 122, 114113 (2005).] Roberts et. al. First introduced into DMC. No Jacobian is needed. SD moves attempted with frequency αSD create a ground state energy bias sd proportional to αSD .



Pacifichem, 2015

25 / 32

Outline

1


2


3




Pacifichem, 2015

26 / 32

Multiple double wells

V (x) =

n X

Vi (xi ) ,

(15)

i=1

where, Vi (xi ) = axi4 + bxi3 + cx 2 + 1,

(16)

and the parameters a, b, and c are chosen to produce two minima (xi = 1, V = 0), (xi = −0.968, V = 0.10 hartree) and a barrier (xi = 0, V = 1 hartree) ∀i.



Pacifichem, 2015

27 / 32


V (x) =

n X

Vi (xi ) ,

(15)

i=1


(16)

and the parameters a, b, and c are chosen to produce two minima (xi = 1, V = 0), (xi = −0.968, V = 0.10 hartree) and a barrier (xi = 0, V = 1 hartree) ∀i. Separable with ground state energy is nE0 , where is the ground state energy of the n = 1 system computed with vector spaces. E.g. for a mass of 200 a.u. E0 = 0.099029hartree



Pacifichem, 2015

27 / 32


V (x) =

n X

Vi (xi ) ,

(15)

i=1


(16)

and the parameters a, b, and c are chosen to produce two minima (xi = 1, V = 0), (xi = −0.968, V = 0.10 hartree) and a barrier (xi = 0, V = 1 hartree) ∀i. Separable with ground state energy is nE0 , where is the ground state energy of the n = 1 system computed with vector spaces. E.g. for a mass of 200 a.u. E0 = 0.099029hartree The system has 2n minima and DMC becomes trapped with as few as 30 dimensions. E. Curotto, S. E. Wolf,, and M. Mella


Pacifichem, 2015

27 / 32

100 double wells

12.5

Ground State Energy (hartree)

Darting stops 12

11.5 IS-DMC (no smart darting) 11

10.5

10

Exact Ground state

0

50000

1e+05 Number of time steps

1.5e+05

2e+05

Figure: n = 100, m = 200 a.u. E0 = -0.099029 hartree, ≈ 1030 minima. E. Curotto, S. E. Wolf,, and M. Mella


Pacifichem, 2015

28 / 32

100 double wells with a lighter mass

24 IS - DMC SD + IS DMC Ground State Energy (hartree)

23

22

21

20 Exact Ground State 19

0

50000

1e+05 Number of time steps

1.5e+05

2e+05

Figure: n = 100, m = 50 a.u. E0 = 0.1950(2) hartree. E. Curotto, S. E. Wolf,, and M. Mella


Pacifichem, 2015

29 / 32

About the Smart Darting bias sd

It originates from the fact that SD moves do not satisfy the correct detailed balance. Moves are accepted or rejected based on some assumed wavefunction (The same used to provide IS). It’s a linear function of the frequency of attempted SD moves αsd , the size of the system n, and the reciprocal of the mass, sd = K m−1/2 αsd n

(17)

It can be eliminated by turning off SD moves after a sufficiently long walk. (At least for the multiple double well system)



Pacifichem, 2015

30 / 32

Summary A Leap Frog - like algorithm for MD and RPMD is developed and tested. The Langevin equation for second order DMC in manifolds is derived. A simple method promises to eliminate the Smart Darting bias in DMC RPMD in manifolds needs to be developed further. A contraction scheme similar to the one Markland and Manoulopolos derive in Cartesian coordinates is in the early stages of development [T. E. Markland, D. E. Manolopoulos, J. Chem. Physics. 129, 024105 (2008).] Currently RPMD converges linearly with respect to the Trotter number N. Can this convergence be accelerated? And is that even useful?

All three methods have shown promising results, but need additional work to become firmly established. E. Curotto, S. E. Wolf,, and M. Mella


Pacifichem, 2015

31 / 32

Aknowledgments

Student collaborators K. Roberts M. Aviles T. Luan H. Christensen S. E. Wolf M. Hayers A. Fodor L. Jake D. Bierswish


Senior Collaborators D. L. Freeman M. Mella J. D. Doll K. Jordan


Pacifichem, 2015

32 / 32