Fermion Quantum Monte Carlo based on the idea of sampling

1
Fermion Quantum Monte Carlobased on the idea of
sampling graphsAli AlaviUniversity of
Cambridge
Alex Thom James Spencer EPSRC
2
Overview Introduction and motivation Paths
integrals and the Fermion sign problem FSP as a
problem in path counting A useful
combinatorial formula From path-sums to
graph-sums Applications to molecular
systems Towards application to periodic systems
3
Essence of idea Express the many-electron path
integral in a finite Slater Determinant basis
Resum the path integral
over exponentially large numbers of paths to
convert path-sums gt graph-sums The
graphs are much more stable entities which can
then be sampled.
k
l
i
j
4
A graphical, or diagrammatic, expansion of the
partition function
.





2-vertex
3-vertex
Each vertex is a Slater determinant Each graph
represents the sum over all paths of length P
which visit all verticies of the
graph Non-pertubative expansion
5
Path Integrals
Consider the (thermal) density matrix
In terms of the eigenstates of the Hamiltonian
The energy can be calculated from
6
The density matrix can be represented in real
space
For a single electron located at x
PE terms
KE termsharmonic spring
7
One can simulate an electron as a ring-polymer,
moving in the external field (which itself can be
dynamic). Polarons Parrinello Rahman
PE along the path
KE along path
8
For N electrons
Describes closed paths which can exchange
identical particle coordinates
Coulomb interaction
Odd permutations subtract from the Partition
function Fermion sign problem
As N or b increases, there is an exponential
cancellation of contributions arising from even
and odd paths.
9
Slater Determinant space
Let Di be a Slater determinant composed out of
N orthonormal spin-orbitals e.g. Hartree-Fock
orbitals, Kohn-Sham, etc chosen out of a set of
2M
The Di form a orthonormal set of antisymmetric
functions. They are solutions to a
non-interacting, or uncorrelated
(mean-field) Hamiltonian H0
Full problem
Exact solutions are linear superpositions of
uncorrelated determinants
10
Paths in Slater determinant space
A closed path in S.D. space
A given S.D. can occur multiple times along a path
11
(No Transcript)
12
Hamiltonian matrix elements (Slater-Condon rules)
Since H contains at most 2-body interactions
i j
a b
Hamiltonian connects only single and double
excitations
Maximum connectivity
Spin selection rule
Other symmetries may also exist Hubbard model
translational invariance Moleculespoint group
symmetry
Cost of calculation
13
Search for a power series in rii
Two-hop
3-hop
14
Rearranging
15
Define the nested sum
which appears in the h-hop term
16
The hop series
17
Using induction, one can show
AJW Thom and A Alavi, J Chem Phys, 123, 204106,
(2005)
18
Residue Theorem gives
19
Some useful properties of Z-sums
? Replace upper limit of sums over h to ?
Symmetry
20
From hop-expansion to vertex expansion
Consider the 4-hop terms
4-vertex
3-vertex
chain
Star
2-vertex
21
Analytic summation over alternating series
Cycle function
22
Eg. A 2-vertex graph
j
(For simplicity)
i
23
Next compute S2
24
(No Transcript)
25
Star graphs
j
j
k
k
i
i
l
l
For a star-graph with g-spokes, G1,G2,Gg
attached to i
26
Chains graphs
G3
G3
G1
G2
27
General 3-vertex graph
Unfolded representation Each spoke represents
an Independent circuit on the graph
Folded representation
j
k
i
Denominator is cubic polynomial in z i.e. there
are 3 residues
28
Unfolded 4-vertex graph
29
Denominator is a quartic polynomial in z
30
A graphical, or diagrammatic, expansion of the
partition function
.





2-vertex
3-vertex
Each vertex is a Slater determinant Each graph
represents the sum over all paths of length P
which visit all verticies of the graph
31
2,3, and 4-vertex graphs




32
Monte Carlo sampling of graphs The energy can be
obtained from
If graphs can be sampled with an un-normalised
probability given by w(n) G, then the energy
estimator is
33
i.e. the number of positive sampled graphs should
exceed the number of negative sampled graphs in
such a way that this difference is finite
and does not vanish.
Monitor fraction of sampled graphs which are
trees, positive cyclic and negative cyclic
graphs.
34
For graphs that contain the HF determinant
Hartree-Fock energy
Ground state energy
ApproximationTruncate series at 2-vertex,
3-vertex or higher-vertex graphs. 2 vertex
Double-excitations Number of graphs
N2M2 3 vertex Quadruple excitations
N4M4 4 vertex
Hexatuple excitations
N6M6
35
N2 molecule
36
N2 molecule in VDZ basis
Types of sampled graphs (4-vertex level)
37
(No Transcript)
38
N2 sampled energies (4-vertex level)
39
N2 binding curve sampling graphs which contain
the HF determinant
40
(No Transcript)
41
Applications to periodic systems

• Taking a plane-wave PP code (CPMD) which can
  solve for
• KS orbitals and potential
• KS virtuals
• -gt Use these as the basis for the vertex series
• KS Hamiltonian becomes the reference
  (single-excitations now
• contribute)
• Need 2-index and 4-index integrals, which are
  computed on-the-fly
• using FFTs (time consuming part)
• Advantage (i) Treatment of periodic systems
• (ii) No BSSE
• (iii) Can be used as a post-DFT
  method

42
Graphite (4 atom) primitive cell. 16el, BHS PP
(Ec90Ry)
2-vertex
43
(No Transcript)
44
Conclusions and outlook
Development of QMC methods based on graphs gives
a method to combat the Fermion sign
problem Proof of concept for small molecular
systems Major effort is now being expended on
developing a periodic code..
..perhaps to return to surface problems in due
course!
45
Advantage of graph-sampling algorithm
O(N2) scaling!

• The observed stability at the 4-vertex level is
  extremely encouraging.
• Current work
• Extension to higher order graphs
• Improved Monte Carlo sampling
• Applications to large systems

46
Graphs
A graph a set of n distinct elements (in no
particular order) with a given connectivity
n distinct determinants
Connectivity of graph is determined by rij
47
Compactly expressed
A set of n connected determinants
Sum over all paths which visit all the
determinants in G
Each graph represents a sum over exponentially
large numbers of paths its weight can be expected
to be much better behaved than that of
individual paths.
48
A graph, G, is an object on which we can
represent the paths which visit all the vertices
in G
The weight of a given graph is the sum over all
paths of length P which visit all the vertices
of the graph
The prime indicates that the summation indicies
must be chosen in Such a way that each vertex in
G is visited at least once.
49
This condition ensures that the weights of two
different graphs Ga and Gb (I.e. two graphs that
differ in at least one vertex) do not
double-count paths which visit only Ga ? Gb
will not double-count paths which visit
50
Quantum Chemical applications
Dissociation of diatomic molecules Multiple-bond
dissociation, e.g. the N2 molecule, is a major
challenge to any ab initio method. Use HF
orbitals generated from MOLPRO Gaussian basis
set cc-pVDZ or VTZ Two-electron primitive
integrals read in from MOLPRO output and ? matrix
constructed on the fly.
51
Cost of the calculations
2-vertex lt1 s 3-vertex
150 secs 4-vertex
1 week over 109 4-vertex graphs to sum On
a pentium 4 processor 2003 vintage
How to make 4-vertex (and eventually higher
vertex) calculations practical?
52
So if on step t of an MC simulation consisting of
K steps we are at graph Gt
In order to perform a Metropolis MC simulation,
one needs to ensure that microscopic
reversibility is satisfied. In the present
implementation, we generate fresh graphs at each
step according to an algorithm to be shortly
described. In addition one needs to compute the
generation probability of a graph using this
algorithm, in order to unbias the Metropolis MC
acceptance ratio.
53
Tree graphs are graphs that do not contain cycles
j
l
k
i
The weight of trees is positive definite at all b
54
Exactly diagonalised by Krogh, Olsen CPL 344,
578, (2001), and by Chan, Kallay and Gauss, JCP,
121, 610 (2004)
15820024220 determinants
55
To summarize
ApproximationTruncate series at 2, 3 or higher
vertex terms. 2 vertex Double-excitations
N2M2 3 vertex Quadruple excitations
N4M4 4 vertex Hexatuple excitations N6M6
By Comparison CCSD N2(M-N)4 Nit
CCSDT N3(M-N)5 Nit
CCSDTQ N4(M-N)6 Nit
CCSDTQ56
N6(M-N)8 Nit
56
An iIlustration of the Monte Carlo ?8x?8 Hubbard
lattice with 6 e-
Momentum- space basis
8 site system at or near half-filling is strongly
open-shell
(1,1)
(0,0)
(1,0)
(2,0)
57
(No Transcript)
58
(No Transcript)
59
Finite T
We wish to compute the energy at a finite b-1kT
as
Where the trace is taken over all Ndet
determinants. Problem is that these sums are not
Monte Carlo-able.
60
Sampling Slater determinants
61
Where the expectation value is taken over an
ensemble of determinants sampled with probability
pi.
62
Perform Metropolis sampling of Di chosen
according to wi
63
Problem the weight itself is a path-integral!
High-temperature DM
Define
Discrete path integral wildly oscillatory
integrand. Cant use Monte Carlo! Hopeless to
calculate by brute-force!
64
Generation of graphs with a computable generation
probability
We adopt a Markov chain algorithm is which
successive determinants are added to a list until
the desired size of graph is reached.
However, since the connectivity of each
determinant is not uniform, such an algorithm
can produce a non-uniform generation probability.
Start at i, and selected a connected determinant,
j, with probability pij. This results in a
2-vertex graph, Gi,j. Next, select k,
connected to j, with probability pjk. If k is
distinct, then add k to the list Gi,j,k.
Otherwise, select a new determinant from the
current Position (i.e. the last visited
determinant). Continue this process until n
distinct verticies have been visited.
65
The generation probability can be calculated by
examining all possible ways of generating G
according to this algorithm. For example, for a
3-vertex graph, Gi,j,k
This procedure can be generalised for a n-vertex
graph (ngt3).
66
The general case is most compactly expressed in
matrix notation. Let us call our n verticies
Gi1,i2,,in, all distinct, with
i1i. Consider the generation probability of G
in the given order (i1,i2,,in). According to
this algorithm, we visit i2 for the first time
from i1, i3 for the first time from either i1 or
i2, etc. In general we visit ik for the
first time from any of the previous visited k-1
verticies. The algorithm terminates when we
first visit the n-th vertex. This is a
first-passage problem in Markov chain theory.
67
We will construct a series of transition
probability matrices in which vertex ik is an
absorbing state
Note that
represents the probability of arriving at ik in
exactly n steps given we Started from ik-1,
passing through some or all of (i1,i2,., ik-1).
68
So therefore the total probability of arriving at
ik, is simply the geometric series
Therefore the probability to generate the graph G
in the sequence
69
The probability to generate G in any order is
given by the sum over all n! permutations
In current implementation we choose
j3
j1
i
j2
Where Ni is the number of determinants connected
to i. In other words we do not introduce an
energetic bias in the selection of determinants.
70
Conclusions

• A new approach for Fermion Monte Carlo is being
  developed, based on
• sampling Slater determinant space with weights
  computed according to
• a novel path counting scheme.
• The mathematics of path-counting needs further
  investigation.
• The scheme has been applied to the Hubbard model
  and the N2 problem with encouraging results.

71
Topology of graphs
j
j
k
k
i
i
Chain (tree)
Star (tree)
Cyclic
i,j, and k all must be single or
double- excitations of each other.
j and k all must be s- or d-excitations of i, but
not necessarily of each other
j must be a s- or d- excitation of i, and k is a
triple or quadruple of i.
72
Future work
Technical Sampling graphs to counter the scaling
problem Calculation of electron
density Parallelisation of code Systems
Hubbard models e.g. stability of striped phases
Dispersion interactions (e.g. graphitic
systems)
73
(No Transcript)
74
(No Transcript)
75
(No Transcript)
76
Contribution to the weights
77
Tentative conclusion is for the Hubbard model
with U4, the 3-vertex approximation is not
perfect, although it is nevertheless an
improvement over UHF Captures about 20-50 of
the correlation energy. Can we estimate the
contribution of the higher order graphs through
a MC sampling? gt work in this direction is in
progress
78
Distribution of terms among the 2,3 and 4 vertex
graphs
UHF basis
79
(No Transcript)
80
(No Transcript)
81
Convergence of Êi with the vertex approximation
82
(No Transcript)
83
?10??10 Hubbard Model
184756 determinants at N10. Exactly
diagonalisable with effort on P4
Half-filled system is closed-shell
4
1
-1
-4
84
Two important questions How good is the 3-vertex
approximation? What is the best one-particle
basis to use?
85
j
A 3-determinant star
i
k
Contour integral solution
86
Again in exact agreement with the diagonalisation
result
87
j
Fully connected 3-vertex graph
k
i
Via diagonalisation
88
Via the contour integral
multiply top and bottom by (z-1)3
Factorise
Cancel factors
Evaluate two residues at z1-b/a, and z12b/a
89
(No Transcript)
90
(No Transcript)
91

• 8-site Hubbard model with N6 electrons
• 3-vertex weights against exact weights

92
Some simple examples. (a) A two-determinant
system
Exact solution via diagonalisation
Solution via Contour integral formula First
define A(z)
93
Next compute S2
In exact agreement with diagonalisation result
94
Residue theorem
95
Exact ground-state energy, UHF and lowest Êi vs
particle number
96
The electron correlation problem How to account
for the fact that electrons move around in a
correlated fashion? Quantum chemistry approach
is Start from Hartree-Fock and try to improve
systematically Hartree-Fock mean-field theory,
N3 N4 ?HFD0 Coupled Cluster CCSD(T), N7
?eT ?HF perturbation
theory Full-configuration interaction eN
? ?HF ?j cjDj
Expansions in Antisymmetric functions
Essential feature of HF theory maintains an
orbital (one-particle) picture of electronic
structure
97
Quantum Monte Carlo

• QMC refers to stochastic methods to solve the
  Schrödinger Equation
• (or sample path integrals) based on interpreting
  the S.E. as a diffusion equation
• in imaginary time
• ? is interpreted as a probability distribution.
  Long-time stochastic propagation
• diffusionlife/death processes leads to
  sampling the nodeless eigenstateof H.
• Application to Fermion systems is severely
  hampered by sign problems.
• Unconstrained sampling of the configuration space
  of Fermions leads to Boson
• Catastrophe.
• QMC can be stabilised by the introduction of
  constraints
• Fixed node approximation in diffusion MC J
  Anderson
• -Restricted path integral MC (fixing nodes of
  density matrix) Ceperley
• -Positive projection and constrained path MC for
  auxillary field QMC
• Fahy and Hamman, Zhang

98
Why not use antisymmetric spaces?

• We would like to explore the possibility of using
  an antisymmetrized space
• as the basis for quantum monte carlo.
• Can we sample a set of Slater determinants in
  such way that we can
• extract meaningful physical quantities (eg
  energy) at the end of
• the simulation?
• This strategy avoids the Boson catastrophe from
  the outset without
• imposition of fixed-node type approximations.
• We will show that
• Such a method is indeed numerically stable
• (2) The MC weights are obtained by summing over
  many paths of
• fluctuating sign.
• (3) Method depends on combinatorial ideas for
  path counting
• -gt So far it is not exact
• (4) Applications to (a) Hubbard model and (b)
  Dissociating molecules
• A major conceptual advantage is that it allows to
  build directly on the

99
Whats the problem?

100
(No Transcript)
101
4su(2pz)
su(2pz)
su(2pz)
2pg
pg
pg
1pu
pu
pu
3sg (2pz)
sg (2pz)
sg (2pz)
2su
su
su
1sg
sg
sg
102
Formally speaking, in the eigenvalue basis of H
103
Hubbard Model
Model of itinerant magnetism for narrow-band
systems. Intensively studied since the mid-80s
in the context of High Tc.
104
(No Transcript)
105
Partition function
Note that the sign of w(P) is a very poorly
behaved quantity Depends on the product of P
matrix elements. Therefore small variations in
the path can lead to wild fluctuations in the
sign of the path.
106
Exact Diagonalisation
Exact solutions are expressed as linear
superpositions of uncorrelated Determinants.
107
Conjecture the n-vertex graph gives rise to a
polynomial of degree n in the denominator of the
contour integral
Contour integrals reduce to a sum over n residues
108
(No Transcript)
109
One can simulate an electron as a ring-polymer,
moving in the external field (which itself can be
dynamic). Polarons Parrinello Rahman
Harmonic springs hold together a ring polymer
110
Motivation
The development of stable fermion QMC algorithms
which do not require fixed-node approximations,
but which maintain a N2 or N3 scaling. Does
working in antisymmetric spaces (eg Slater
determinants) help? Intuitively, Slater
determinant spaces are the right spaces to be
dealing with fermions i.e. one should build in
the fermion antisymmetry in the outset in any
N-particle representation.
111
Computational cost of electronic structure methods
Hartree Fock N3
MP2 - MP4 N4-N7
Coupled Cluster CCSD-(T) N6-N7
FCI eN
DFT N3
QMC N2-N3
112
Electron correlation is ubiquitous in chemistry
e.g. in molecular dissociation
Incorrect dissociation
113
Configuration Interaction
Consider the doubly excited determinant
Correlated wavefunction
Correct dissociation
114
What is the problem with configuration
interaction?

• Slowly convergent with respect to short and
  intermediate
• range correlation.
• Must include many determinants the problem grows
  exponentially
• with number of electrons and the number of
  virtuals
• -(linear) Truncated CI lacks size consistency
• Coupled cluster methods are nowadays preferred.