Transcript Document

Universita’ degli Studi dell’Insubria
Quantum Monte Carlo
Simulations of Mixed
3He/4He Clusters
Dario Bressanini
[email protected]
http://www.unico.it/~dario
Overview

Introduction to quantum monte carlo methods


Helium clusters simulations


VMC, QMC, advantages and drawbacks
Problems, solutions
Mixed 3He/4He clusters

Trimers
 3

He4HeN
Future directions
© Dario Bressanini
2
Monte Carlo Methods

How to solve a deterministic problem using a Monte
Carlo method?

Rephrase the problem using a probability distribution
A   P(R) f (R)dR

R  N
“Measure” A by sampling the probability distribution
1
A 
N
© Dario Bressanini
N
 f (R )
i 1
i
R i ~ P(R )
3
Monte Carlo Methods

The points Ri are generated using random numbers
This is why the methods are
called Monte Carlo methods

Metropolis, Ulam, Fermi, Von Neumann (-1945)

We introduce noise into the problem!!

Our results have error bars...

... Nevertheless it might be a good way to proceed
© Dario Bressanini
4
Quantum Mechanics

We wish to solve H Y = E Y to high accuracy


The solution usually involves computing integrals in
high dimensions: 3-30000
The “classic” approach (from 1929):

Find approximate Y ( ... but good ...)

... whose integrals are analitically computable (gaussians)

Compute the approximate energy
chemical accuracy ~ 0.001 hartree ~ 0.027 eV
© Dario Bressanini
6
VMC: Variational Monte Carlo

Start from the Variational Principle
H

Y (R ) HY (R)dR


E
 Y (R)dR
0
2
Translate it into Monte Carlo language
H   P(R ) EL (R )dR
HY (R )
EL (R ) 
Y (R )
© Dario Bressanini
P(R ) 
Y (R )
2
2
Y
 (R)dR
7
VMC: Variational Monte Carlo
E  H   P(R ) EL (R )dR

E is a statistical average of the local energy EL over P(R)
1
E H 
N

N
E
i 1
L
(R i )
R i ~ P(R )
Recipe:
 take an appropriate trial wave function
 distribute N points according to P(R)
 compute the average of the local energy
© Dario Bressanini
8
VMC: Variational Monte Carlo

There is no need to analytically compute integrals, so there is
complete freedom in the choice of the trial wave function

Quantum chemistry uses a
product of single particle
functions. With VMC we can
use any function: explicitly
correlated wave functions
can be used
r12
r1
r2
He atom
Ye
© Dario Bressanini
a r1 b r2  g r12
9
The Metropolis Algorithm
P(R ) 
Y (R )
2
?

How do we sample

Use the Metropolis algorithm (M(RT)2 1953) ...
... and a powerful computer

The algorithm is a random
walk (markov chain) in
configuration space
2
Y
 (R)dR
Anyone who consider
arithmetical methods of
producing random digits
is, of course, in a state of sin.
John Von Neumann
© Dario Bressanini
10
The Metropolis Algorithm
move
Ri
Call the Oracle
Rtry
reject
accept
Ri+1=Ri
Ri+1=Rtry
Compute
averages
© Dario Bressanini
12
The Metropolis Algorithm
The Oracle
 Y (new) 

p  
 Y (old ) 
2
if p  1
/* accept always */
accept move
If 0  p  1 /* accept with probability p
*/
if p > rnd()
accept move
else
reject move
© Dario Bressanini
13
VMC advantages

No need to make the single-particle approximation

Can use  for which no analytical integrals exist

Use explicitly correlated wave functions

Can satisfy the cusp conditions
He atom ground state
E19 terms = -2.9037245 a.u.
Exact = -2.90372437 a.u.
© Dario Bressanini
Y   ci e
ai r1  bi r2  g i r12
i
14
VMC advantages

Can compute difficult quantities, e.g.
1
H 
N
  H
2

N
E
L
(R i ) 
i
2
 H
2
H2

N
H
2
1

N
N
E
2
L
(R i )
i
Can compute lower bounds
H  E0  H  
© Dario Bressanini
15
VMC advantages

Can easily go beyond the Born-Oppenheimer
approximation.
H2+ molecule ground state
Y   ci e
2
ai rA  bi rB  g i rAB  f i rAB
i
E1 term = -0.596235(9)a.u.
E10 terms = -0.597136(3)a.u.
Exact = -0.597139
a.u.
© Dario Bressanini
16
VMC advantages

Can work with ANY potential, in ANY number of
dimensions.
Ps2 molecule (e+e+e-e-) in 2D and 3D

Optimization of nonlinear parameters

min  ( H )  min H
2
2
 H
2


Numerically stable
 Minimum known in advance (0)
 Can be used for excited states with same symmetry too
© Dario Bressanini
17
VMC drawbacks

Error bar goes down as N-1/2

It is computationally demanding

The optimization of Y becomes difficult as the
number of nonlinear parameters increases

It depends critically on our skill to invent a good Y

There exist exact, automatic ways to get better wave
functions.
Let the computer do the work ...
© Dario Bressanini
19
Diffusion Monte Carlo

VMC is a “classical” simulation method
Nature is not classical, dammit, and if you
want to make a simulation of nature, you'd
better make it quantum mechanical, and by
golly it's a wonderful problem, because it
doesn't look so easy.
Richard P. Feynman

Suggested by Fermi in 1945, but implemented only in
the 70’s
© Dario Bressanini
20
Diffusion equation analogy
 The time dependent
Schrödinger equation
is similar to a diffusion
equation
 The diffusion
equation can be
“solved” by directly
simulating the system
Y
2 2
i

 Y  VY
t
2m
C
2
 D C  kC
t
Time
evolution
Diffusion
Branch
Can we do the same with the Schrödinger equation ?
© Dario Bressanini
21
Imaginary Time Sch. Equation

The analogy is only formal

Y is a complex quantity, while C is real and positive
Y (R, t )  e iEnt /   n (R)

If we let the time t be imaginary, then Y can be real!
Y
 D 2 Y  VY

Imaginary time Schrödinger equation
© Dario Bressanini
22
Y as a concentration

Y is interpreted as a concentration of fictitious
particles, called walkers

The schrödinger equation
is simulated by a process
of diffusion, growth and
disappearance of walkers
Y
2
 D Y  VY

Y (R, )   ai  i (R)e  ( Ei  ER )
i
Y(R,  )   0 (R)e ( E0  ER )
Ground State
© Dario Bressanini
23
Diffusion Monte Carlo
SIMULATION: discretize time
•Diffusion process
Y
 D 2 Y

Y(R,  )  e
 ( R  R 0 ) 2 / 4 D
•Kinetic process (branching)
Y
 (V (R )  ER )Y

Y(R,  )  e
© Dario Bressanini
 (V ( R )  ER ) 
Y(R,0)
24
The DMC algorithm
© Dario Bressanini
25
The Fermion Problem

Wave functions for fermions have nodes.

Diffusion equation analogy is lost. Need to introduce
positive and negative walkers.
The (In)famous Sign Problem

Restrict random walk to a positive region bounded by nodes.
Unfortunately, the exact nodes are unknown.

Use approximate nodes from
a trial Y. Kill the walkers if
they cross a node.
© Dario Bressanini
+
-
26
Helium

A helium atom is an
elementary particle. A weakly
interacting hard sphere.

Interatomic potential is known
more accurately than any
other atom.
Two isotopes:
• 3He (fermion: antisymmetric trial function, spin 1/2)
• 4He (boson: symmetric trial function, spin zero)
• The interaction potential is the same
© Dario Bressanini
27
Helium Clusters

They are not rigid, and show large-amplitude motion

Normal mode analysis useless

“Equilibrium structure” is ill-defined

Stochastic methods well suited to study helium
clusters, both pure or with impurities

VMC, DMC, GFMC, PIMC etc...
© Dario Bressanini
28
Helium Clusters
1. Small mass of helium atom
2. Very weak He-He interaction
0.02 Kcal/mol
0.9 * 10-3 cm-1
0.4 * 10-8 hartree
10-7 eV
Highly non-classical systems
Superfluidity
High resolution
spectroscopy
Low temperature
chemistry
© Dario Bressanini
30
The Simulations

Both VMC and DMC simulations

Potential = sum of two-body TTY pair-potential
V (R )  VHe He (rij )
i j


Three-body terms not important for small clusters
Standard Y 
N
  (r
4
i j
N
He  4 He
)  (r4 He3 He )
k
p5 p2
 (r )   (r )  exp(  5  2  p0 ln( r )  p1r )
r
r
© Dario Bressanini
31
Pure
4He
n
Clusters
0
DMC
VMC
Energy cm-1
-1
-2
-3
-4
The quality of
the VMC
simulations
decreases as the
cluster increases
-5
-6
-7
-8
2 3 4 5 6 7 8 9 10 11 12
n
© Dario Bressanini
32
Y for

4He
n
Clusters
Wave function quality decreases as N increases

It was optimized to get minimum (H), not minimum <H>

Are three- and many-body terms in Y important ?

Very difficult to optimize. Unstable process
especially for the trimers.
Can we improve Y ?
© Dario Bressanini
33
Mixed
0.0
(0,2)
Energy cm-1
-0.5
3He/4He
Clusters
(1,2) (2,2)
(0,3)
(m,n) = 3Hem4Hen
(1,3)
(0,4)
(1,4)
Bressanini et. al.
J.Chem.Phys.
112, 717 (2000)
(0,5)
-1.0
(1,5)
-1.5
(0,6)
-2.0
(1,6)
(0,7)
-2.5
1
© Dario Bressanini
2
3
4
5
Number of atoms
6
7
34
Helium Clusters: energy (cm-1)
N DMC 4HeN VMC 4HeN-13He DMC 4HeN-13He
2
3
4
5
6
7
11
20
-0.00089(1)
-0.08784(7)
-0.3886(1)
-0.9015(3)
-1.6077(4)
-2.4805(7)
-7.286(1)
-23.04(1)
© Dario Bressanini
-0.00666(2)
-0.19199(2)
-0.57484(6)
-1.1505(2)
-1.8595(2)
-5.975(3)
-19.98(1)
-0.00984(5)
-0.2062(1)
-0.6326(2)
-1.2626(4)
-2.0718(5)
-6.679(4)
-22.234(9)
35
Helium Clusters: stability
 4HeN

is destabilized by substituting a 4He with a 3He
The structure is only weakly perturbed.
Dimers
4He4He
Bound
Trimers
4He
3
Bound
Tetramers
4He
4
Bound
© Dario Bressanini
4He3He
3He3He
Unbound
Unbound
4He 3He
2
4He3He
2
Bound
Unbound
4He 3He
3
Bound
4He 3He
2
2
Bound
36
Trimers and Tetramers Stability
4He
3
E = -0.08784(7) cm-1
4He 3He
2
E = -0.00984(5) cm-1
Bonding interaction
Non-bonding interaction
4He
4
E = -0.3886(1) cm-1
4He 3He
3
E = -0.2062(1) cm-1
4He 3He
2
2
E = -0.071(1) cm-1
Five out of six unbound pairs!
© Dario Bressanini
37
3He/4He
Distribution Functions
Pair distribution functions
0.016
4
4
He- He
3
4
He- He
0.012
g(r)
3He(4He)
5
0.008
0.004
0.000
0
© Dario Bressanini
10
20
r (u.a.)
30
40
38
3He/4He
Distribution Functions
Distributions with respect to the center of mass
0.016
3He(4He)
5
0.012
g(r)
4
He
3
He
0.008
c.o.m
0.004
0.000
0
© Dario Bressanini
10
20
r (u.a.)
30
40
39
4He 3He
N
Distribution Functions in
r(4He-4He)
r(3He-4He)
0.015
0.3
N=4
N=19
N=5
N=6
N=10
N=3
0.010
N=19
P(r)
P(r)
0.2
N=2
0.005
0.1
N=2
0.0
0
0.000
10
20
r (bohr)
© Dario Bressanini
30
0
10
20
30
r (bohr)
40
Distribution Functions in 4HeN3He
c.o.m. = center of mass
r(4He-C.O.M.)
r(3He-C.O.M.)
0.015
0.015
N=3
0.010
3He4He
N=19
r (r) (bohr-3)
r (r)
-3
(bohr )
N=19
2
0.005
0.000
0.010
3He4He
2
N=3
0.005
0.000
0
10
r (bohr)
Similar to pure clusters
© Dario Bressanini
20
0
10
20
30
r (bohr)
Fermion is pushed away
41
What is the shape of
4He

Some people say is an equilateral triangle ...

... some say it is linear (almost) ...

... some say it is both.
3
?
Pair distribution function
© Dario Bressanini
42
What is the shape of
© Dario Bressanini
4He
3
?
43
The Shape of the Trimers
Ne trimer
r(Ne-center of mass)
He trimer
r(4He-center of mass)
© Dario Bressanini
44
Ne3 Angular Distributions
a
b
Ne trimer
b
b
a
© Dario Bressanini
a
45
4He
3
a
b
Angular Distributions
b
b
a
a
© Dario Bressanini
b
a
46
3He4He
a
2
Angular Distributions
b
b
b
a
© Dario Bressanini
a
47
Different wave function form
p5 p2
( b  rij )
 (r )  exp(  5  2  p0 ln( rij )  p1rij  a e
)
rij rij
Spline
0.60
0.60
0.40
0.40
f(r)
f(r)
0.20
0.20
0.00
0.00
0.00
2.00
4.00 r 6.00
8.00
Optimize heights
© Dario Bressanini
10.00
0.00
2.00
4.00
6.00
r (a.u.)
8.00
10.00
48
Y for
4He
2

Y literature (Rick & Doll) E = -0.00046 cm-1

Optimize 
Unbound

Optimize E (numerically)
E = -0.00075 cm-1

Y with Exp()
E = -0.00084 cm-1

Y using splines
E = -0.00081 cm-1

QMC
E = -0.00089(1) cm-1

Numerical
E = -0.00091 cm-1
© Dario Bressanini
49
Y for
4He
3

Y literature (Rick & Doll) E = -0.0798 cm-1

Optimize Energy
E = -0.0829(4) cm-1

Y with Exp()
E = -0.0851(2) cm-1

Y using splines
E = -0.0868(2) cm-1

QMC exact
E = -0.08784(7) cm-1
three-body terms are not important in Y for the trimer
© Dario Bressanini
50
Work in Progress

Various impurities embedded in a Helium cluster

Y for bigger clusters using splines

Optimize the energy instead of the variance of the
local energy.

Geometric structure of trimers and tetramers
© Dario Bressanini
51
Conclusions

The substitution of a 4He with a 3He leads to an
energetic destabilization.
 3He
weakly perturbes the 4He atoms distribution.
 3He
moves on the surface of the cluster.
 4He23He
bound, 4He3He2 unbound.
 4He33He
and 4He23He2 bound.

QMC gives accurate energies and structural information
© Dario Bressanini
52
Acknowledgments
Gabriele Morosi
Massimo Mella
Mose’ Casalegno
Giordano Fabbri
Matteo Zavaglia
© Dario Bressanini
53
A reflection...
A new method for calculating properties in nuclei, atoms, molecules, or
solids automatically provokes three sorts of negative reactions:
 A new method is initially not as well formulated or understood as
existing methods
 It can seldom offer results of a comparable quality before a
considerable amount of development has taken place
 Only rarely do new methods differ in major ways from previous
approaches
Nonetheless, new methods need to be developed to handle
problems that are vexing to or beyond the scope of the
current approaches
(Slightly modified from Steven R. White, John W. Wilkins and Kenneth G. Wilson)
© Dario Bressanini
54