Zero Temperature QMC Outline of lecture

1
Zero Temperature QMC Outline of lecture

• Variational Monte Carlo
• Diffusion Monte Carlo
• Fermion systems
• Coupled Electron-Ion Monte Carlo
• Applications
• Electron gas
• High pressure hydrogen

2
Variational Monte Carlo

• Historically first quantum simulation method
• Slater-Jastrow trial function
• Examples liquid helium and electron gas.
• Wavefunctions for Quantum solids
• WaveFunctions beyond Slater-Jastrow back flow
  and 3-body
• Twist Averaged Boundary Conditions

3
The Variational Theorem

• Assume is a trial function
  where R are the quantum degrees of freedom
  (positions, spin) a are variational
  parameters.
• Conditions matrix elements exist, symmetries and
  boundary conditions are correct.

4

• Expand trial function in terms of the exact
  eigenfunctions
• Energy and variance are second order in
  (1-overlap).
• Other properties are first order.
• Temple lower bound

?
E0 E
5
Linear Basis approach

• Assume trial function is a linear combination of
  known functions a basis fn(R).

Problem that MC solves Unless we use 1-particle
basis, integrals are too slow to perform.
6
Properties of solution to GEP
E

• For a basis of size m, there exist m
  eigenvalues and orthonormal eigenfunctions
• McDonalds theorem the nth eigenvalue in a basis
  is an upper bound to the nth exact eigenvalue.
• We can always lower all the energies by
  augumenting the basis
• When basis is complete, we get exact answers!

7
Symmetry reduces complexity

• If an operator P commutes with H P,H0 we can
  reduce complexity by working in a basis with that
  symmetry.
• rotational symmetry use Ylm.
• Translation symmetry use plane waves.
• Inversion symmetry even/odd functions
• Matrix elements are non-zero only for states
  within the same sector.
• By reordering we can block diagonalize
• Reduces complexity from M3 to k(M/k)3M3/k2.
• McDonalds theorem applies to each sector
  individually.

8
Notation

• Individual coordinate of a particle ri
• All 3N coordinates R (r1,r2, . rN)
• Total potential energy V(R)
• Kinetic energy
• Hamiltonian

9
Liquid heliumthe prototypic quantum fluid

• Interatomic potential is known more accurately
  than any other atom because electronic
  excitations are so high.
• A helium atom is an elementary particle. A weakly
  interacting hard sphere.

• Two isotopes
• 3He (fermion antisymmetric trial function, spin
  1/2)
• 4He(boson symmetric trial function, spin zero)

10
Helium interaction

• Repulsion at short distances because of overlap
  of atomic cores.
• Attraction at long distance because of the
  dipole-induced-dipole force. Dispersion
  interaction is
• c6r-6 c8 r-8 .
• He-He interaction is the most accurate. Use all
  available low density data (virial coefficients,
  quantum chemistry calculations, transport
  coefficients, .) Good to better than 0.1K (work
  of Aziz over last 20 years).
• Three body interactions are small but not zero.

11
Helium phase diagram

• Because interaction is so weak helium does not
  crystallize at low temperatures. Quantum exchange
  effects are important
• Both isotopes are quantum fluids and become
  superfluids below a critical temperature.
• One of the goals of computer simulation is to
  understand these states, and see how they differ
  from classical liquids starting from
  non-relativistic Hamiltonian

12
Variational Monte Carlo (VMC)

• Variational Principle. Given an appropriate trial
  function
• Continuous
• Proper symmetry
• Normalizable
• Finite variance
• Quantum chemistry uses a product of single
  particle functions
• With MC we can use any computable function.
• Sample R from ?2 using MCMC.
• Take average of local energy
• Optimize ? to get the best upper bound
• Error in energy is 2nd order
• Better wavefunction, lower variance! Zero
  variance principle. (non-classical)

13
First Major QMC Calculation

• PhD thesis of W. McMillan (1964) University of
  Illinois.
• VMC calculation of ground state of liquid helium
  4.
• Applied MC techniques from classical liquid
  theory.
• Ceperley, Chester and Kalos (1976) generalized to
  fermions.

• Zero temperature (single state) method
• Can be generalized to finite temperature by using
  trial density matrix instead of trial
  wavefunction.

14
Trial function for helium 4Jastrow or pair
product

• We want finite variance of the local energy.
• Whenever 2 atoms get close together wavefunction
  should vanish.
• The pseudopotential u(r) is similar to classical
  potential
• Local energy has the form
• G is the pseudoforce
• If v(r) diverges as ?r-n how should u(r) diverge?
  Assume
• u(r)?r-m
• Keep N-1 atoms fixed and let 1 atom approach
  another and analyze the singular parts of the
  local energy.
• Gives a condition on u at small r.
• For Lennard-Jones 6-12 potential, Jastrow goes as
  m5

15
Periodic boundary conditions

• Minimum Image Conventiontake the closest
  distance
• rM min ( rnL)
• Potential is cutoff so that V(r)0 for rgtL/2
  since force needs to be continuous. Remember
  perturbation theory.
• Image potential
• VI ? v(ri-rjnL)
• For long range potential this leads to the Ewald
  image potential. You need a back ground and
  convergence method.

16
(No Transcript)
17
Optimization of trial function

• Try to optimize u(r) using reweighting
  (correlated sampling)
• Sample R using P(R)??2(R,a0)
• Now find minima of the analytic function Ev(a)
• Or minimize the variance (more stable but
  wavefunctions less accurate).
• Statistical accuracy declines away from a0.

18
Quantum Crystal Trial Function

• Jastrow trial function does not freeze at
  appropriate density.
• Solution is to break spatial symmetry by hand.
• Introduce a bcc lattice
• Zi
• bcc has the lowest Madelung energy, but others
  may have lower zero point energy.
• Introduce localized one-body terms (Wannier
  functions).
• Make a Slater determinant possibly with spin
  ordering.
• More complicated trial functions and methods are
  also possible.

• C is a variational parameter to be optimized.

19
Scalar Properties, Static Correlations and Order
Parameters

• What do we get out of a simulation? Energy by
  itself doesnt tell you very much.
• Other properties
• do NOT have an upper bound property
• Only first order in accuracy
• EXAMPLES
• Static properties pressure, specific heat etc.
• Density
• Pair correlation in real space and fourier space.
• Order parameters and broken symmetry How to tell
  a liquid from a solid
• Specifically quantum the momentum distribution

20
Other quantum properties

• Kinetic energy
• Potential energy
• Pair correlation function
• Structure function
• Pressure (virial relation)
• Momentum distribution
• Non-classical showing effects of bose or fermi
  statistics
• Fourier transform is the single particle
  off-diagonal density matrix
• Compute with McMillan Method.
• Condensate fraction 10

Like properties from classical simulations No
upper bound property Only first order in accuracy
21
Order parameters

• A system has certain symmetries e.g.translation
  invariance.
• At high temperatures one expects the system to
  have those same symmetries at the microscopic
  scale. (e.g. a gas)
• BUT as the system cools those symmetries can be
  broken. (a gas condenses or freezes).
• At a liquid gas-transition the density is no
  longer fills the volume droplets form. The
  density is the order parameter.
• At a liquid-solid transition, both rotational
  symmetry and translational symmetry are broken
• The best way to monitor the transition is to look
  how the order parameter changes during the
  simulation.

22
Microscopic Density

• ?(r) lt ?i ?(r-r i) gt
• Or you can take its Fourier Transform
• ? k lt ?i exp(ikri) gt
• (This is a good way to smooth the density.)
• A solid has broken symmetry (order parameter).
  The density is not constant.
• At a liquid-gas transition the density is also
  inhomogeneous.
• In periodic boundary conditions the k-vectors are
  on a grid k2?/L (nx,ny,nz) Long wave length
  modes are absent.
• In a solid Lindemanns ratio gives a rough idea
  of melting
• u2 lt(ri-zi)2gt/d2

23
Snapshots of densities

• Liquid or crystal or glass? Blue spots are
  defects

24
Pair Correlation Function, g(r)

• Primary quantity in a liquid is the probability
  distribution of pairs of particles. Given a
  particle at the origin what is the density of
  surrounding particles
• g(r) lt ?iltj ? (ri-rj-r)gt (2 ?/N2)
• Density-density correlation function
• From g(r) you can calculate all pair quantities
• (potential, pressure, )
• A function gives more information than a number!

25
g(r) in liquid and solid helium

• First peak is at inter-particle spacing. (shell
  around the particle)
• goes out to rltL/2 in periodic boundary
  conditions.

26
(The static) structure factor S(k)

• The Fourier transform of the pair correlation
  function is the structure factor
• problem with (2) is to extend g(r) to infinity
• S(k) is measured in neutron and X-Ray scattering
  experiments.
• Can provide a direct test of the assumed
  potential.
• Used to see the state of a system
• liquid, solid, glass, gas? (much better than g(r)
  )
• Order parameter in solid is ?G where G is a
  particular wavevector (reciprocal lattice vector).

27
Bragg peak

• In a perfect lattice S(k) will be non-zero only
  on the reciprocal lattice vectors G S(G) N
• At non-zero temperature (or for a quantum system)
  this structure factor is reduced by the
  Debye-Waller factor
• S(G) 1 (N-1)exp(-G2u2/3)
• To tell a liquid from a crystal we see how S(G)
  scales as the system is enlarged. In a solid,
  S(k) will have peaks that scale with the number
  of atoms.

• The compressibility is given by
• We can use this is detect the liquid-gas
  transition since the compressibility should
  diverge as k approaches 0. (order parameter is
  density)

28
Crystal liquid
29
Momentum Distribution

• Momentum distribution
• Classically momentum distribution is always a
  Gaussian
• Non-classical showing effects of bose or fermi
  statistics
• Fourier transform is the single particle
  off-diagonal density matrix
• Compute with McMillan Method.
• For fermions we need to use the determinant
  update formulas to find the effect of the
  movement of 1 electron.

30
Derivation of momentum formula

• Suppose we want the probability nk that a given
  atom has momentum hk.
• Find wavefunction in momentum space by FT wrt all
  the coordinates and integrating out all but one
  electron
• Expanding out the square and performing the
  integrals we get.
• Where
• (states occupied with the Boltzmann
  distribution.)
• For a homogeneous system, n(r,s)n(r-s)

31
Calculation of n(r)

• Naïve procedure
• Generate sample with VMC
• Take one particle at random and displace by a
  distance r.
• Find ratio of trial function for old and new
  position.
• O(1) work/1 distance

• McMillan procedure
• Generate sample with VMC
• Determine change of trial function if each
  particle is destroyed. ak
• Insert a new particle at a random position in the
  box and determine trial function for the
  insertion b
• For each k perform the average b/ak and add to
  n(r-rk).
• Repeat for O(N) insertions.
• O(N) work/O(N2) distances

32
The electron gasD. M. Ceperley, Phys. Rev. B 18,
3126 (1978)

• Standard model for electrons in metals
• Basis of DFT.
• Characterized by 2 dimensionless parameters
• Density
• Temperature
• What is energy?
• When does it freeze?
• What is spin polarization?
• What are properties?

33
Charged systems

• How can we handle charged systems?
• Just treat like short-ranged potential cutoff
  potential at rgtL/2. Problems
• Effect of discontinuity never disappears (1/r)
  (r2) gets bigger.
• Will violate Stillinger-Lovett conditions because
  Poisson equation is not satisfied
• Image potential solves this
• VI ? v(ri-rjnL)
• But summation diverges. We need to resum. This
  gives the ewald image potential.
• For one component system we have to add a
  background to make it neutral.
• Even the trial function is long ranged and needs
  to be resummed.

34
Ewald summation method

• Key idea is to split potential into k-space part
  and real-space part. We can do since FT is
  linear.
• Bare potential converges slowly at large r (in
  r-space) and at large k (in k-space)

35
Classic Ewald

• Split up using Gaussian charge distribution
• If we make it large enough we can use the minimum
  image potential in r-space.
• Extra term for insulators

36
How to do it
O(N3/2) O(N) O(N3/2) O(N3/2) O(N1/2) O(N3/2) O(
1)

• r-space part same as short-ranged potential
• k-space part
• Compute exp(ik0xi)(cos (ik0xi), sin (ik0xi)),
  k02?/L ? i
• Compute powers exp(i2k0xi) exp(ik0xi
  )exp(ik0xi) etc. This way we get all values of
  exp(ik . ri) with just multiplications.
• Sum over particles to get ?k all k.
• Sum over k to get the potentials.
• Forces can also be done by taking gradients.
• Constant terms to be added.
• Checks perfect lattice V-1.4186487/a (cubic
  lattice).

37
Optimized EwaldJ. Comput. Physics 117, 171
(1995).

• Division into Long-range and short-ranged
  function is convenient but is it optimal? No
• Trial functions are also long-ranged but not
  simply 1/r. We need a procedure for general
  functions.
• Natoli-Ceperley procedure. What division leads
  to the highest accuracy for a given radius in r
  and k?
• Leads to a least squares problem.
• FITPN code does this division.
• Input is fourier transform of vk
• on grid appropriate to the supercell
• Output is a spline of vsr(r)
• and table of long ranged function.

38
Problems with Image potential

• Introduces a lattice structure which may not be
  appropriate.
• Example a charge layer.
• We assume charge structure continues at large r.
• Actually nearby fluid will be anticorrelated.
• This means such structures will be penalized.
• One should always consider the effects of
  boundary conditions, particularly when
  electrostatic forces are around!
• You need to have a continuum model to understand
  the results of a microscopic simulation.

39
Fermions antisymmetric trial function

• At mean field level the wavefunction is a Slater
  determinant. Orbitals for homogenous systems are
  a filled set of plane waves.
• We can compute this energy analytically (HF).
• To include correlation we multiply by a
  pseudopotential. We need MC to evaluate
  properties.
• New feature how to compute the derivatives of a
  deteminant and sample the determinant. Use tricks
  from linear algebra.
• Reduces complexity to O(N2).

Slater-Jastrow trial function.
40
Jastrow factor for the e-gas

• Look at local energy either in r space or
  k-space
• r-space as 2 electrons get close gives cusp
  condition du/dr0-1
• K-space, charge-sloshing or plasmon modes.
• Can combine 2 exact properties in the Gaskell
  form. Write EV in terms structure factor making
  random phase approximation. (RPA).
• Optimization can hardly improve this form for
  the e-gas in either 2 or 3 dimensions. RPA works
  better for trial function than for the energy.
• NEED EWALD SUMS because potential trial function
  is long range, it also decays as 1/r, but it is
  not a simple power.

• Long range properties important
• Give rise to dielectric properties
• Energy is insensitive to uk at small k
• Those modes converge t1/k2

41
Derivation of the e-gas Jastrow

• For simplicity, consider boson trial function

42
Generalized Feynman-Kacs formula

• Lets calculate the average population resulting
  from DMC starting from a single point R0 after a
  time t.

43
Wavefunctions beyond Jastrow
smoothing

• Use method of residuals construct a sequence of
  increasingly better trial wave functions.
  Justify from the Importance sampled DMC.
• Zeroth order is Hartree-Fock wavefunction
• First order is Slater-Jastrow pair wavefunction
  (RPA for electrons gives an analytic formula)
• Second order is 3-body backflow wavefunction
• Three-body form is like a squared force. It is a
  bosonic term that does not change the nodes.

44
Backflow wave function

• Backflow means change the coordinates to quasi-
  coordinates.
• Leads to a much improved energy and to
  improvement in nodal surfaces. Couples nodal
  surfaces together.
• Kwon PRB 58, 6800 (1998).

3DEG
45
Dependence of energy on wavefunction 3d Electron
fluid at a density rs10 Kwon, Ceperley,
Martin, Phys. Rev. B58,6800, 1998

• Wavefunctions
• Slater-Jastrow (SJ)
• three-body (3)
• backflow (BF)
• fixed-node (FN)
• Energy ltf H fgt converges to ground state
• Variance ltf H-E2 fgt to zero.
• Using 3B-BF gains a factor of 4.
• Using DMC gains a factor of 4.

FN -SJ
FN-BF
46
Comparison of Trial functions

• What do we choose for the trial function in VMC
  and DMC?
• Slater-Jastrow (SJ) with plane wave orbitals
• For higher accuracy we need to go beyond this
  form.
• Need correlation effects in the nodes.
• Include backflow-three body.
• Example of incorrect physics within SJ

47
Analytic backflowHolzmann et al, Phys. Rev. E
68, 0467071-15(2003).

• Start with analytic Slater-Jastrow using Gaskell
  trial function
• Apply Bohm-Pines collective coordinate
  transformation and express Hamiltonian in new
  coordinates
• Diagonalize resulting Hamiltonian.
• Long-range part has Harmonic oscillator form.
• Expand about k0 to get backflow and 3-body
  forms.
• Significant long-range component to BF
• OPTIMIZED BF
  ANALYTIC BF
• 3-body term is non-symmetric

rs1,5,10,20
48
Results of Analytic tf

• Analytic form EVMC better for rslt20 but not for
  rs?20.
• Optimized variance is smaller than analytic.
• Analytic nodes always better! (as measured by
  EDMC)
• Form ideal for use at smaller rs since it will
  minimize optimization noise and lead to more
  systematic results vs N, rs and polarization.
• Saves human machine optimization time.
• Also valuable for multi-component system of
  metallic hydrogen.

49
Twist averaged boundary conditions

• In periodic boundary conditions (? point), the
  wavefunction is periodic?Large finite size
  effects for metals because of shell effects.
• Fermi liquid theory can be used to correct the
  properties.
• In twist averaged BC we use an arbitrary phase ?
  as r ?rL
• If one integrates over all phases the momentum
  distribution changes from a lattice of k-vectors
  to a fermi sea.
• Smaller finite size effects

PBC TABC
50
Twist averaged MC

• Make twist vector dynamical by changing during
  the random walk.
• Within GCE, change the number of electrons
• Within TA-VMC
• Initialize twist vector.
• Run usual VMC (with warmup)
• Resample twist angle within cube
• (iterate)
• Or do in parallel.

51
Grand Canonical Ensemble QMC

• GCE at T0K choose N such that E(N)-?N is
  minimized.
• According to Fermi liquid theory, interacting
  states are related to non-interacting states and
  described by k.
• Instead of N, we input the fermi wavevector(s)
  kF. Choose all states with k lt kF (assuming
  spherical symmetry)
• N will depend on the twist angle ?. number of
  points inside a randomly placed sphere.
• After we average over ? (TA) we get a sphere of
  filled states.
• Is there a problem with Ewald sums as the number
  of electrons varies? No! average density is
  exactly that of the background. We only work with
  averaged quantities.

52
Single particle size effects

• Exact single particle properties with TA within
  HF
• Implies momentum distribution is a continuous
  curve with a sharp feature at kF.
• With PBC only 5 points on curve
• rs4 N33 polarized
• No size effect within single particle theory!
• Kinetic energy will have much smaller size
  effects.

53
Potential energy

• Write potential as integral over structure
  function
• Error comes from 2 effects.
• Approximating integral by sum
• Finite size effects in S(k) at a given k.
• Within HF we get exact S(k) with TABC.
• Discretization errors come only from non-analytic
  points of S(k).
• the absence of the k0 term in sum. We can put
  it in by hand since we know the limit S(k) at
  small k (plasmon regime)
• Remaining size effects are smaller, coming from
  the non-analytic behavior of S(k) at 2kF.

54
3DEG at rs10
TABC
TABC1/N
PBC
GC-TABC1/N
We can do simulations with N42! Size effects now
go like We can cancel this term at special
values of N! N 15, 42, 92, 168, 279,
55
Brief History of Ferromagnetism in electron gas

• What is polarization state of fermi liquid at low
  density?
• Bloch 1929 got polarization from exchange
  interaction
• rs gt 5.4 3D
• rs gt 2.0 2D
• Stoner 1939 include electron screening contact
  interaction
• Herring 1960
• Ceperley-Alder 1980 rs gt20 is partially
  polarized
• Young-Fisk experiment on doped CaB6 1999 rs25.
• Ortiz-Balone 1999 ferromagnetism of e gas at
  rsgt20.
• Zong et al Redo QMC with backflow nodes and
  TABC.

56
Ceperley, Alder 80
T0 calculations with FN-DMC

• 3d electron gas
• rslt20 unpolarized
• 20ltrslt100 partial
• 100ltrs Wigner crystal

Energies are very close together at low density!
More recent calculations of Ortiz, Harris and
Balone PRL 82, 5317 (99) confirm this result but
get transition to crystal at rs65.
57
Polarization of 3DEG

• We see second order partially polarized
  transition at rs52
• Is the Stoner model (replace interaction with a
  contact potential) appropriate? Screening kills
  long range interaction.
• Wigner Crystal at rs105

• Twist averaging makes calculation possible--much
  smaller size effects.
• Jastrow wavefunctions favor the ferromagnetic
  phase.
• Backflow 3-body wavefunctions more paramagnetic

58
Phase Diagram

• Partially polarized phase at low density.
• But at lower energy and density than before.
• As accuracy gets higher, polarized phase shrinks
• Real systems have different units.

59
Recent calculations in 2D
Tanatar,Ceperley 89 Rapisarda, Senatore 95 Kwon
et al 97
T0 fixed-node calculation Also used high
quality backflow wavefunctions to compute energy
vs spin polarization. Energies of various phases
are nearly identical Attaccalite et al PRL 88,
256601 (2002)

• 2d electron gas
• rs lt25 unpolarized
• 25lt rs lt35 polarized
• rs gt35 Wigner crystal

60
Polarization of 2D electron gas

• Same general trend in 2D
• Partial polarization before freezing
• Results using phase averaging and BF-3B
  wavefunctions

rs20
rs10
rs30
61
Linear response for the egas

• Add a small periodic potential.
• Change trial function by replacing plane waves
  with solutions to the Schrodinger Eq. in an
  effective potential.
• Since we dont care about the strength of
  potential use trial function to find the
  potential for which the trial function is
  optimal.
• Observe change in energy since density has mixed
  estimator problems.

62
Fermi Liquid parameters

• Do by correlated sampling Do one long MC random
  walk with a guiding function (something
  overlapping with all states in question).
• Generate energies of each individual excited
  state by using a weight function
• Optimal Guiding function is
• Determine particle hole excitation energies by
  replacing columnsfewer finite size effects this
  way. Replace columns in slater matrix
• Case where states are orthogonal by symmetry is
  easier, but non-orthogonal case can also be
  treated.
• Back flow needed for some excited state since
  Slater Jastrow has no coupling between unlike
  spins.

63
Summary of Variational (VMC)
Simple trial function
Better trial function
applications
64
Summary and problems with variational methods

• Powerful method since you can use any trial
  function
• Scaling (computational effort vs. size) is almost
  classical
• Learn directly about what works in wavefunctions
• No sign problem

• Optimization is time consuming
• Energy is insensitive to order parameter
• Non-energetic properties are less accurate. O(1)
  vs. O(2) for energy.
• Difficult to find out how accurate results are.
• Favors simple states over more complicated
  states, e.g.
• Solid over liquid
• Polarized over unpolarized

What goes into the trial wave function comes out!
GIGO We need a more automatic method! Projector
Monte Carlo
65
Projector Monte Carlo

• Originally suggested by Fermi and implemented in
  1950 by Donsker and Kac for H atom.
• Practical methods and application developed by
  Kalos

66
Projector Monte Carlo(variants Greens function
MC, Diffusion MC, Reptation MC)

• Project single state using the Hamiltonian
• We show that this is a diffusion branching
  operator. Maybe we can interpret as a
  probability. But is this a probability?
• Yes! for bosons since ground state can be made
  real and non-negative.
• But all excited states must have sign changes.
  This is the sign problem.
• For efficiency we do importance sampling.
• Avoid sign problem with the fixed-node method.

67
Diffusion Monte Carlo

• How do we analyze this operator?
• Expand into exact eigenstates of H.
• Then the evolution is simple in this basis.
• Long time limit is lowest energy state that
  overlaps with the initial state, usually the
  ground state.
• How to carry out on the computer?

68
Notation

• Individual coordinate of a particle ri
• All 3N coordinates R (r1,r2, . rN)
• R can depend on imaginary time, time slice or
  Trotter index t or on iteration number Rt.
• Total potential energy V(R)
• Kinetic energy
• Hamiltonian

69
The Greens function

• Operator notation
• We define the coordinate greens function (or
  density matrix by
• Roughly the probability density of going from R0
  to R in time t. (but is it a probability?)
• Properties

70
Froebinius Theorem

• When can we consider the wavefunction as a
  probability? First how about the Greens
  function?
• But if we start with a non-negative function it
  will stay non-negative, and can be interpreted as
  a p.d.f.
• Not true for all Hamiltonians (require
  off-diagonal matrix elements to be non-positive.)
  (not pseudopotentials, not magnetic fields.)
• Only true for the bosonic ground state.

71
Monte Carlo process

• Now consider the variable t as a continuous
  time (it is really imaginary time).
• Take derivative with respect to time to get
  evolution.
• This is a diffusion branching process.
• Justify in terms of Trotters theorem.
• Requires interpretation of the wavefunction as a
  probability density.
• But is it? Only in the boson ground state.
  Otherwise there are nodes. Come back to later.

72
Trotters theorem

• How do we find the solution of
• The operator solution is
• Trotters theorem (1959)
• Assumes that A,B and AB are reasonable
  operators.
• This means we just have to figure out what each
  operator does independently and then alternate
  their effect. This is rigorous in the limit as
  n??.
• In the DMC case A is diffusion operator, B is a
  branching operator.
• Just like molecular dynamics At small time we
  evaluate each operator separately.

73
Evaluation of kinetic density matrix
74
Putting this together

• n is number of time slices.
• ? is the time-step
• V is diagonal
• Error at finite n comes from commutator is
  roughly
• Diffusion preserves normalization but potential
  does not!

75
Basic DMC algorithm

• Construct an ensemble (population P(0)) sampled
  from the trial wavefunction. R1,R2,,RP
• Go through ensemble and diffuse each one
  (timestep ?)
• number of copies
• Trial energy ET adjusted to keep population
  fixed.
• Problems
• Branching is uncontrolled
• Population unstable
• What do we do about fermi statistics?

ndrn uprn floor function
76
Sampling the normal distribution

• Inverse mapping is a little slow, also of
  infinite range.
• Trick generate 2 at a time r(x,y)
• Or sample angle using rejection technique
• Sample (x,y) in square
• Accept if x2y2 lt1
• Normalize to get the correct r.

77
Code to sample normal distribution

• Normal distribution ltxgtx0 and lt(x-x0)2gt?2

1 xsprng()-0.5 ysprng()-0.5 r2xxyy
if (r2gt0.25) go to 1 radius sigmasqrt
(-2ln(sprng())/r2) xnormalx0xradius
ynormaly0yradius

• No trig functions, 1 log, 1 sqrt, 1 divide
• Mixes up regularity of random numbers
• Efficiency of angle generation is 4/?.
• Gets 2 ndrns each time.

78
Population Bias

• Having the right trial energy guarantees that
  population will on the average be stable, but
  fluctuations will always cause the population to
  either grow too large or too small.
• Various ways to control the population
• Suppose P0 is the desired population and P(t) is
  the current population. How much do we have to
  adjust ET to make P(tT)P0?
• Feedback procedure

• There will be a (small) bias in the energy caused
  by a limited population.

79
Importance SamplingKalos 1970, Ceperley 1979

• Why should we sample the wavefunction? The
  physically correct pdf is ?2.
• Importance sample (multiply) by trial wave
  function.
• Evolution diffusion drift
  branching
• Use accept/reject step for more accurate
  evolution.
• make acceptance ratiogt99 . Determines time
  step.
• We have three terms in the evolution equation.
  Trotters theorem still applies.

Commute ? through H
80
Brownian Dynamics

• Consider a big molecule in a solvent. In the high
  viscosity limit the master equation
  (Smoluchowski or Fokker-Planck eq.) is

Diffusion Quantum Monte Carlo without branching
is the same as Brownian Dynamics. Use same
techniques.
81
Greens function for a gradient

• What is Greens function for the operator?
• This operator just causes probability
  distribution to drift in the direction of F.
• Smoluchowski equation for Brownian motion it was
  the effect of gravitational field on the motion
  of colloids.
• In practice, we limit the gradient so the walk is
  not pushed too far.

82

• To the pure diffusion algorithm we have added a
  drift step that pushes the random walk in
  directions of increasing trial function
• Branching is now controlled by the local energy
• Because of zero variance principle, fluctuations
  are controlled.
• Cusp condition can limit infinities coming from
  singular potentials.
• We still determine ET by keeping asymptotic
  population stable.
• Must have accurate time evolution. Adding
  accept/reject step is a major improvement.

83

• Importanced sampled Greens function
• Exact property of DMC Greens function
• We enforce detailed balance to decrease time step
  errors.
• VMC satisfies detailed balance.
• Typically we choose time step to have 99
  acceptance ratio.
• Method gives exact result if either time step is
  zero or trial function is exact.

84
Schematic of DMC

• Ensemble evolves according to
• Diffusion
• Drift
• branching
• ensemble

85
(No Transcript)
86
Local Markov process

• Let us make a course approximation to a molecular
  dynamics of a large molecule in the presence of a
  solvent that we want to eliminate.
• Particle will follow a continuous trajectory
  R(t).
• Let us assume that the movement is Markovian and
  local.
• Describe in terms of the Greens function
  ?(r,r,t)
• This is probability density of a particle at
  (r,t) starting at (r,0).
• The master equation describes its time evolution
• Evolution
• Initial condition

87
General form of evolution

• What could we have for O?
• Not hard to prove that O must be a sum of first
  and second derivative operators. The most
  general form is
• Clearly probability is conserved.
• Steady state solution is

88

• How do we prove dynamics is unique?
• If it is a continuous walk-it has to be made of
  little random steps.
• Central limit theorem says that the only things
  that survive adding together many little steps
  are
• Mean value this is the drift term
• The variance this is the diffusion term
• Hence the smart MC gaussian is the most general
  Markovian form, if we allow a general diffusion
  matrix.
• If we choose it to be isotropic and independent
  of position then D(R) is a constant.
• it could vary in some media both in position and
  direction
• also it could couple particles far apart,
  hydrodynamical effects.

89
Moment Expansion
For the variance
90
Mixed estimators

• Problem is that PMC samples the wrong
  distribution.
• OK for the energy
• Linear extrapolation helps correct this
  systematic error
• Other solutions
• Maximum overlap
• Forward walking
• Reptation/path integrals

91
Forward WalkingKalos et al. 1974.

• Lets calculate the average population resulting
  from DMC starting from a single point R0 after a
  time t.
• We can estimate the correction to the mixed
  estimator by weighting with the number of
  descendants of a given configuration.
• Problem the fluctuations in the weights
  eventually diverge.
• Dont make t too large.

92
Fusion sticking coefficientPhys. Rev. A 31, 1999
(1985).

• Consider the 3 body system (? d t)
• For the sticking coefficient, we need the exact
  wavefunction at the point where 2 nuclei are at
  the same position. (this is a singular point)

93
Other projector functions can be used

• Common effect on long-time (iteration) limit.
• 3rd choice generates a Krylov sequence. Only
  works for bounded spectra such as a lattice model.

G
E
94
Greens Function Monte CarloKalos, Levesque,
Verlet Phys. Rev. A9, 2178 (1974).

• It is possible to make a zero time-step-error
  method
• Works with the integral formulation of DMC
• Sample time-step from Poisson distribution
• Express operator in a series expansion and sample
  the terms stochastically.
• Recent Revival Continuous time Monte Carlo for
  lattice models.

95
Fermions?

• How can we do fermion simulations? The initial
  condition can be made real but not positive (for
  more than 1 electron in the same spin state)
• In transient estimate or released-node methods
  one carries along the sign as a weight and
  samples the modulus.
• Do not forbid crossing of the nodes, but carry
  along sign when walks cross.
• Whats wrong with node release
• Because walks dont die at the nodes, the
  computational effort increases (bosonic noise)
• The signal is in the cancellation which dominates
• Monte Carlo can add but not subtract

96
Transient Estimate Approach

• ?(?) converges to the exact ground state
• E is an upper bound converging to the exact
  answer monotonically

97
VMC
300K
DMC
timestep
bcc hydrogen, N16
98
Model fermion problem Particle in a box

• Symmetric potential V(r) V(-r)
• Antisymmetric state f(r)-f(-r)

Initial (trial) state
Final (exact) state
Sign of walkers fixed by initial position. They
are allowed to diffuse freely. f(r) number of
positive-negative walkers. Node is dynamically
established by diffusion process. (cancellation
of positive and negative walkers.)
99
Scaling in Released-Node
Initial distribution
Later distribution

• At any point, positive and negative walkers will
  tend to cancel so the signal is drown out by the
  fluctuations.
• Signal/noise ratio is tprojection time
• EF and EB are Fermion, Bose energy (proportional
  to N)
• Converges but at a slower rate. Higher accuracy,
  larger t.
• For general excited states
• Exponential complexity!
• Not a fermion problem but an excited state
  problem.
• Cancellation is difficult in high dimensions.

100
Exact fermion calculations

• Possible for the electron gas for up to 60
  electrons.
• 2DEG at rs1 N26
• Transient estimate calculation with SJ and BF-3B
  trial functions.

101
General statement of the fermion problem

• Given a system with N fermions and a known
  Hamiltonian and a property O. (usually the
  energy).
• How much time T will it take to estimate O to an
  accuracy e? How does T scale with N and e?
• If you can map the quantum system onto an
  equivalent problem in classical statistical
  mechanics then

With 0 lta lt 4

• This would be a solved quantum problem!
• All approximations must be controlled!
• Algebraic scaling in N!
• e.g. properties of Boltzmann or Bose systems in
  equilibrium.

102
Solved Problems

• 1-D problem. (simply forbid exchanges)
• Bosons and Boltzmanons at any temperature
• Some lattice models Heisenberg model, 1/2 filled
  Hubbard model on bipartite lattice (Hirsch)
• Spin symmetric systems with purely attractive
  interactions ult0 Hubbard model, nuclear Gaussian
  model.
• Harmonic oscillators or systems with many
  symmetries.
• Any problem with ltiHjgt ?0
• Fermions in special boxes
• Other lattice models
• Kalos and coworkers have invented a pairing
  method but it is not clear whether it is
  approximation free and stable.

103
The sign problem

• The fermion problem is intellectually and
  technologically very important.
• Progress is possible but danger-the problem maybe
  more subtle than you first might think. New ideas
  are needed.
• No fermion methods are perfect but QMC is
  competitive with other methods and more general.
• The fermion problem is one of a group of related
  problems in quantum mechanics (e.g dynamics).
• Feynman argues that general many-body quantum
  simulation is exponentially slow on a classical
  computer.
• Troyer Wiese show that some quantum sign
  problems are NP hard.
• Maybe we have to solve quantum problems using
  analog quantum computers programmable quantum
  computers that can emulate any quantum system.

104
Fixed-node method

• Initial distribution is a pdf.
• It comes from a VMC simulation.
• Drift term pushes walks away from the nodes.
• Impose the condition
• This is the fixed-node BC
• Will give an upper bound to the exact energy, the
  best upper bound consistent with the FNBC.

• f(R,t) has a discontinuous gradient at the nodal
  location.
• Accurate method because Bose correlations are
  done exactly.
• Scales well, like the VMC method, as N3.
  Classical complexity.
• Can be generalized from the continuum to lattice
  finite temperature, magnetic fields,
• One needs trial functions with accurate nodes.

105
Proof of fixed-node theorem

• Suppose we solve S.E. in a subvolume V determined
  by the nodes of an antisymetric trial function.

106
Nodal Properties

• If we know the sign of the exact wavefunction
  (the nodes), we can solve the fermion problem
  with the fixed-node method.
• If f(R) is real, nodes are f(R)0 where R is the
  3N dimensional vector.
• Nodes are a 3N-1 dimensional surface. (Do not
  confuse with single particle orbital nodes!)
• Coincidence points ri rj are 3N-3 dimensional
  hyper-planes
• In 1 spatial dimension these points exhaust the
  nodes fermion problem is easy to solve in 1D
  with the no crossing rule.
• Coincidence points (and other symmetries) only
  constrain nodes in higher dimensions, they do not
  determine them.
• The nodal surfaces define nodal volumes. How many
  nodal volumes are there? Conjecture there are
  typically only 2 different volumes ( and -)
  except in 1D. (but only demonstrated for free
  particles.)

107
Nodal Picture 2d slice thru 322d space

• Free electron
• Other electrons
• Nodes pass thru their positions
• Divides space into 2 regions
• Wavelength given by interparticle spacing

108
SPIN?

• How do we treat spin in QMC?
• For extended systems we use the Sz
  representation.
• We have a fixed number of up and down electrons
  and we antisymmetrize among electrons with the
  same spin.
• This leads to 2 Slater determinants.
• For a given trial function, its real part is also
  a trial function (but it may have different
  symmetries), for example momentum
• For the ground state, without magnetic fields,
  spin-orbit interaction we can always work with
  real functions.
• However, in some cases it may be better to work
  with complex functions.

109
Fixed-Phase methodOrtiz, Martin, DMC 1993

• Generalize the FN method to complex trial
  functions
• Since the Hamiltonian is Hermitian, the
  variational energy is real
• We see only one place where the energy depends on
  the phase of the wavefunction.
• If we fix the phase, then we add this term to the
  potential energy. In a magnetic field we get also
  the vector potential.
• We can now do VMC or DMC and get upper bounds as
  before.
• The imaginary part of the local energy will not
  be zero unless the right phase is used.
• Used for twisted boundary conditions, magnetic
  fields, vortices, phonons, spin states,

110
Problem with core electrons

• Bad scaling in both VMC and DMC
• In VMC, energy fluctuations from core dominate
  the calculation
• In DMC, time step will be controlled by core
  dynamics
• Solution is to eliminate core states by a
  pseudopotential
• Conventional solution semi-local form
• Ensures that valence electrons go into well
  defined valence states with the wavefunction and
  energy for each angular momentum state
  prescribed.
• PP is non-local OK for VMC. Leads to an extra MC
  integral. But DMC uses a locality approximation
  and good trial functions. Extra approximation.

111
Solid state / chemical applications

• Use LDA derived pseudopotentials
• Take orbitals from other methods
• Gaussian-xx give orbitals for molecules
• DFT-PW codes give orbitals for extended systems
• HF is slightly better because of self-interaction
  effects within DFT
• Multiply by a Jastrow function (electron gas or
  otherwise). Can include higher order e-e-n terms
• Must add a compensating e-n term in order to
  cancel out purely repulsive character of e-e
  correlation.
• Assuming LDA density is correct, this can be done
  by making sure VMC electron densityLDA electron
  density. (Fahy correction).

112
Summary of T0 methodsVariational(VMC),
Fixed-node(FN), Released-node(RN)
113
Problems with projector methods

• Fixed-node is a super-variational method
• DMC dynamics is determined by Hamiltonian
• Zero-variance principle allows very accurate
  calculation of ground state energy if trial
  function is good.
• Projector methods need a trial wavefunction for
  accuracy. They are essentially methods that
  perturb from the trial function to the exact
  function. (Note if you dont use a trial
  function, you are perturbing from the ideal gas)
• Difficulty calculating properties other than
  energy. We must use extrapolated estimators or
  forward walking.
• Bad for phase transitions and finite temperature,
  complex systems.
• Path Integral MC solves some of these problems.

114
Reptation Methods
115
Reptation Monte Carlo (RQMC)

• Similar technique to Diffusion MC
• Instead of imaginary timecomputer time, keep
  entire path in memory
• Update with a Metropolis based method instead of
  branching diffusing random walks
• Get exact properties no extrapolation or mixed
  estimators
• Good for energy differences.
• How to move the particles? Reptation means move
  like a snake. This is how polymers can move.

116
Reptation Monte Carlo

• ?(?) converges to the exact ground state as a
  function of imaginary time. More accurate than
  VMC
• E is an upper bound converging to the exact
  answer monotonically
• Do Trotter break-up into a path of p steps a la
  PIMC.
• Bosonic action for the links
• Trial function at the end points.
• For fixed-phase add a potential to avoid the
  sign problem. Exact answer if potential is
  correct.

117
Reptation moves

• Let d be the direction of the move
• -1 tail move
• 1 head move
• Standard method.
• Choose d at random
• Acceptance probability is
• Takes O(p2) steps to decorrelate.
• One way reptation gives the wrong answers.
• Bounce method
• add d to the state.
• Change d only on rejections.
• Use same acceptance formula!!
• Does not satisfy detailed balance but still gives
  correct answer since it is an eigenfunction of T.
• Moves are 1/(rejection rate) times more
  effective!

118
Link action

• PIMC method pair potential plus nodal action
• DMC method. Use importanced-sampled evolution to
  suggest action (good for accurate trial
  functions)
• Symmetrize with respect to R and R to get higher
  accuracy.

119
VMC
300K
DMC
timestep
Tests on bcc hydrogen, N16 Good results in a
few slices p20.
120
Coupled Ionic-Electronic Simulations Carlo
Pierleoni, Kris Delaney and DMCIllinois ,
LAquila
As computers get more powerful, we can use more
accurate intermolecular potentials
121
Quantum Monte Carlo

• We need to use simulation techniques to treat
  many-body quantum problems, just as in classical
  liquids, lattice spin models
• Various QMC methods for electronic structure
  problem
• Variational Monte Carlo
• Projector Monte Carlo e.g. Diffusion MC,
  Reptation MC
• Path Integral Monte Carlo
• Let us assume that QMC gives most accurate
  energies for general quantum many-body systems?
  we want to use it for the energy surface of an
  ionic system.
• We want to use the Born-Oppenheimer
  approximation a big separation in electron and
  ionic energies.
• Try out CEIMC method on a simple but difficult
  many-body system warm dense hydrogen.

122
Coupled electron ion MC

• To calculate thermodynamic properties, we sample
  the classical (or quantum) Boltzmann
  distribution.
• The electronic energy E(S) is obtained by solving
  the electronic Schroedinger equation for a given
  position, S, of the ions.
• MC is simpler and more rigorous than MD
• No ergodic problems.
• Flexibility in choosing transition moves.
• Possibility of canceling noise. (how do we do
  this in MD?)

123

• In Coupled Electron-Ionic Monte Carlo (CEIMC )
• Do classical MC or Path Integral MC for the ions
  at Tgt0.
• Let electrons be at zero temperature, a
  reasonable approximation for room temperature
  simulations.
• Do MC, not MD, for the ions.
• Assumptions are convenient, not necessary.
• Advantages to having T0 electrons
• Zero variance principle in ground state can
  reduce computation time for high quality wave
  functions.
• Quality of zero temperature wavefunctions is
  higher than for Tgt0 implies fewer time slices.
• Possibility of using existing DFT/SCF/QMC
  wavefunctions at T0.
• Energy difference methods can be efficient.

124
Plan of talk

• Describe some features of the CEIMC method
• The penalty method
• New trial functions
• Reptation MC
• Energy difference method
• Path Integrals for protons
• Twist averaged boundary conditions
• Hydrogen Phase diagram
• Crystal-liquid transition
• Molecular-atomic transition
• Advantages of Coupled-Electron Ion MC

125
A single CIEMC step

• Move ions from S to S
• Pre-reject moves based on empirical potential
• Determine new trial function ?T (RS)
• Sample electron coordinates from P(RS,S)
• Determine energy difference E(S)-E(S).
• Use penalty method to accept or reject

R
electrons
ions
S ? S
126
Basics of the classical random walk methodMarkov
chain with rejections

• The Metropolis, Rosenbluth, Teller (1953) method
• Move from S to S with probability
  T(S?S)T(S?S)
• Accept move with probability
• a(S?S) min 1 , exp ( - (E(S)-E(S))/kBT )
• Repeat many times
• Given ergodicity, the distribution of the state,
  S, will be
• exp(-E(S)/kBT)/Z
• E(S)energy of ionic arrangement S,
• Zpartition function.
• Only the difference in energy enters in the
  acceptance probability.

127
Problem QMC energy difference will be noisy
Average energy of Lennard-Jones liquid

• Ignoring noise gives a systematic increase in the
  energy because high energy moves are occasionally
  wrongly accepted.
• Acceptance formula is non-linear min 1 , exp (
  - ?E/kBT )
• But it is possible get the exact distribution,
  independent of noise level!!

128
Approaches to deal with this noise

• Ignore problem will smooth out distributions.
• Compute energy difference very accurately takes
  a lot of computer time!!
• Langevin Dynamics (Grossman and Mitas PRL, 2005)
• Solve MD with random component to forces.
• Noise must be small with respect to time step and
  temperature.
• Dynamics is artificial since noise doesnt
  correspond to a physical process.
• Linear Noise (Kennedy and Kuti, PRL 54, 2473,
  1985).
• Used in lattice gauge simulations where noise is
  small.
• Sign problem for large fluctuations
• Non-optimal it increases rejections
  unnecessarily, especially for large noise.
• Our approach modify the formula for accepting
  moves.
• Similar to fuzzy MC by Krajcí and Hafner, PRL
  74, 5100 (1995) but with higher efficiency and
  fewer approximations.
• Recent development (more general distributions)
• R. C. Ball et al. PRL 91, 030201 (2003).

129
The Penalty method DMC Dewing, J. Chem. Phys.
110, 9812(1998).

• Assume estimated energy difference ?e is normally
  distributed with variance ?2 and the correct
  mean.
• lt ?e gt ?E
• lt ?e- ?E2 gt ?2
• OK because of central limit theorem for ?lt?
• a(?e ?) is acceptance ratio.
• average acceptance A(?E) lt a(?e) gt
• Markov chain goes to the correct distribution if
  flux of transitions from S to S is symmetric
  detailed balance A(?E) exp (- ?E ) A(-?E)
• An exact solution is to use a modifed acceptance
  formula a(x,?) min 1, exp(-x- ?2/2)
• ?2/2 is penalty additional rejections caused
  by the noise.

130
Estimating the error

• Problem both the mean and variance are unknown
  and are spatially varying.
• Both energy and variance are estimated from the
  data and must be included in the acceptance rate
  formula a(?,v)
• Sample M independent estimates of ?E ek
  1?k?M.
• ?? ek/M ? ?E
• v? ek - ?2/(M(M-1)) ? ?2
• 2 parameter integral equation, normal in ? and
  ?2(vM) in v.
• Unique smooth solution is a Bessel function
  JM/2(?Ms). Expand
• a(?,s) exp - ? - v/2 - v2/4(M1)
  -v3/3(M1)(M3)
• Error of noise causes an additional penalty!
• Asymptotic series only good for vlt4M. But very
  accurate for vlt0.1M.

131
2 Parameter Integral Equation
The estimated variance, s, is distributed
according to a chi-squared distribution with
(N-1) degrees of freedom.
Detailed Balance Equation
to be satisfied for all values of y and ?2 .
132
Optimization of Efficiency w.r.t. Noise Level

• As the noise level increases, the acceptance
  ratio decreases because of the penalty.
• Assume it is usual MC rule to reduce the
  variance. That is
• ?2t/M where M data points ? CPU
  time ? (?/?)2 constant
• ? error of the energy per step.
• For double well, optimal noise turns out to be
  3kBT!
• Large noise (order k B T) is much more efficient
  than small noise.
• Trying many cheap steps, rejecting most of
  them, is more successful than spending much
  computer time crafting a step which will be
  accepted.
• Nonetheless there is n upper limit to the noise
  because you need to estimate the variance also.

Average time to achieve a given error in mean
energy
133
Other uses of the Penalty method

• Consider a problem where there are a finite but
  large number of terms in the potential.
• Example long range potential (order N terms
  after 1 particle is moved.) Ewald sums for
  coulomb potential scale as N1/2 .
• Fast multipole method O(1) is only appropriate
  if all particles are moved together as in
  molecular dynamics.
• We must compute the neighboring term
  directly--otherwise error will not be normally
  distributed. Split pair potential into
  short-range and long-range terms v(r)
  vlr(r)vsr(r)
• We can sample terms at position, r, of the
  long-range part. Since it does not diverge
  anywhere, its variance is finite.
• The optimal sampling probability is

134
Pre-rejection

• For efficiencys sake we want to make moves as
  far as possible.
• We use an empirical potential (Vm) to decide if a
  move is at all reasonable (first level
  rejection) choose Vm self-consistently.
• 2-stag