Zero temperature qmc outline of lecture
1 / 57

Zero Temperature QMC Outline of lecture - PowerPoint PPT Presentation

  • Uploaded on

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. Notation. Individual coordinate of a particle r i All 3N coordinates R= (r 1 ,r 2 , …. r N )

I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
Download Presentation

PowerPoint Slideshow about 'Zero Temperature QMC Outline of lecture' - bachyen

An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
Zero temperature qmc outline of lecture
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


  • Individual coordinate of a particle ri

  • All 3N coordinates R= (r1,r2, …. rN)

  • Total potential energy = V(R)

  • Kinetic energy :

  • Hamiltonian :

Variational monte carlo vmc

Variational Principle. Given an appropriate trial function:


Proper symmetry


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)

Variational Monte Carlo (VMC)

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.

First Major QMC Calculation

  • Zero temperature (single state) method

  • Can be generalized to finite temperature by using “trial” density matrix instead of “trial” wavefunction.

The electron gas d m ceperley phys rev b 18 3126 1978
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?

Fermions antisymmetric trial function
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.

Jastrow factor for the e gas
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/dr|0=-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 t~1/k2

Wavefunctions beyond jastrow


Wavefunctions beyond Jastrow

  • 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.

Dependence of energy on wavefunction3d Electron fluid at a density rs=10Kwon, Ceperley, Martin, Phys. Rev. B58,6800, 1998

  • Wavefunctions

    • Slater-Jastrow (SJ)

    • three-body (3)

    • backflow (BF)

    • fixed-node (FN)

  • Energy <f |H| f> converges to ground state

  • Variance <f [H-E]2f> to zero.

  • Using 3B-BF gains a factor of 4.

  • Using DMC gains a factor of 4.



Twist averaged boundary conditions


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 r+L

  • If one integrates over all phases the momentum distribution changes from a lattice of k-vectors to a fermi sea.

  • Smaller finite size effects



Polarization of 3deg



Polarization of 3DEG

  • Twist averaging makes calculation possible--much smaller size effects.

  • Jastrow wavefunctions favor the ferromagnetic phase.

  • Backflow 3-body wavefunctions more paramagnetic

  • We see second order partially polarized transition at rs=52

  • Is the Stoner model (replace interaction with a contact potential) appropriate? Screening kills long range interaction.

  • Wigner Crystal at rs=105

Phase diagram
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.

Summary of variational vmc
Summary of Variational (VMC)

Simple trial function

Better trial function


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

Summary and problems with variational methods

What goes into the trial wave function comes out! “GIGO”

We need a more automatic method! Projector Monte Carlo

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:

Projector monte carlo variants green s function mc diffusion mc reptation mc
Projector Monte Carlo(variants: Green’s 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.

Diffusion monte carlo
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?

Monte carlo process
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 Trotter’s 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.

Basic dmc algorithm
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?



floor function

Importance sampling kalos 1970 ceperley 1979
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 ratio>99% . Determines time step.

  • We have three terms in the evolution equation. Trotter’s theorem still applies.

Commute  through H

Schematic of dmc
Schematic of DMC

Ensemble evolves according to

  • Diffusion

  • Drift

  • branching



  • 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.

  • What’s wrong with node release:

    • Because walks don’t 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

Transient estimate approach
Transient Estimate Approach

  • () converges to the exact ground state

  • E is an upper bound converging to the exact answer monotonically

Model fermion problem particle in a box

Positive walkers


Negative walkers

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.)

Scaling in released node
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 : t=projection time

    EF and EBare 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.

Exact fermion calculations
Exact fermion calculations

  • Possible for the electron gas for up to 60 electrons.

  • 2DEG at rs=1 N=26

  • Transient estimate calculation with SJ and BF-3B trial functions.

General statement of the fermion problem
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 ande?

  • If you can map the quantum system onto an equivalent problem in classical statistical mechanics then:

With 0 <a < 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.

Solved problems
“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: u<0 Hubbard model, nuclear Gaussian model.

  • Harmonic oscillators or systems with many symmetries.

  • Any problem with <i|H|j> 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.

The sign problem
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.

Fixed node method
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.

Nodal properties
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.)

Fixed phase method ortiz martin dmc 1993
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, …

Summary of t 0 methods variational vmc fixed node fn released node rn

Simple trial function



Better trial function



Summary of T=0 methods:Variational(VMC), Fixed-node(FN), Released-node(RN)

Reptation monte carlo rqmc
Reptation Monte Carlo (RQMC)

  • Similar technique to Diffusion MC:

    • Instead of imaginary time=computer 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.

Reptation monte carlo
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.

Reptation moves
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!





Tests on bcc hydrogen, N=16

Good results in a few slices p~20.

Coupled electron ion mc



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?)

Basics of the classical random walk method markov chain with rejections
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:


    E(S)=energy of ionic arrangement “S”,

    Z=partition function.

    Only the difference in energy enters in the acceptance probability.

Problem: QMC energy difference will be noisy

  • 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!!

Average energy of Lennard-Jones liquid

The penalty method dmc dewing j chem phys 110 9812 1998
The Penalty methodDMC & Dewing,J. Chem. Phys. 110, 9812(1998).

  • Assumeestimated energy difference e is normally distributed* with variance 2and the correct mean.

    < e > = E

    < [e- E]2 > = 2

    * OK because of central limit theorem for <

  • a(e; ) is acceptance ratio.

  • average acceptance A(E) = < a(e) >

  • 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)

  • Anexact 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.

Why hydrogen
Why Hydrogen?

  • Most abundant element

  • Theoretically clean:

    • 1 electron and 1 proton per atom

    • No pseudopotential required

  • A rich variety of properties, including:

    • Metal-insulator transition in fluid and solid: 70 year quest for metallic hydrogen

    • Possible liquid-liquid phase transition

    • Possibility of superconducting and superfluid phases

  • Equation of state not yet fully described but crucial in understanding planetary formation

  • Fusion applications

Qmc methods for dense hydrogen
QMC methods for dense hydrogen

Path Integral MC for

T > EF/10

Diffusion MC T=0

Coupled-electron Ion MC

Path Integral MC with an effective potential

Trial functions
Trial functions.

  • What do we choose for the trial function?

  • Standard choice: Slater-Jastrow function:

    with the orbital from a rescaled LDA calculation.

  • Reoptimization of trial functions during the CEIMC run is a major difficulty in time and reliability.

  • Requires an LDA calculation after each proton move.

  • It would be better to have a trial function which depends analytically on proton coordinates.

  • backflow + three body trial function are very successful for homogeneous systems

  • We generalize them to many-body hydrogen.

Energy difference methods

Nodal surfaces




“Reptile” space

Energy difference methods

  • We need an efficient way of computing difference: [E(S)-E(S*)]

  • Naïve (direct) method is to do separate (uncorrelated) samples of S and S*. Noise increases by 2.

  • Correlated methods map S walks into S* walks.

  • Simplest is “VMC re-weighting” (1-sided)

  • With fixed-node fermions, we need to worry about changes in the nodal surfaces. 1-sided methods can give the wrong answer because the distributions are not defined in the same regions of path space.

Optimal importance sampling
Optimal Importance Sampling

  • What distribution has the lowest variance for the energy difference? (ignoring autocorrelation effects)

  • Sum of squares is almost as good, and eliminates barriers which might be hard to cross.

  • Symmetric in the two distributions.

  • Generalizable to reptation MC, by including the action.

  • We are using to do other energy differences: energy of a molecular fragment.

Experimentally known high pressure phase diagram of H


H2 Bond-ordered phase

  • Wigner-Huntington (1935) predicted that eventually hydrogen should be a metal.

  • Experiments have not reached that pressure.

Two Possible Phase Diagrams

liquid H

liquid H

Solid H

Solid H

  • Ashcroft suggested a low temperature liquid metallic ground state.

  • Does the liquid go to T=0K?

  • How low in temperature is needed to see quantum protonic transitions?

  • How about electronic superconductivity?

Compare with DFT simulations

  • Temperature dependence in LDA is off by 100%.

  • This effect also seen in Natoli et al. calculation of various metallic hydrogen crystal structures and for liquid H2 structures.

  • In LDA (and some other functionals) energy landscape is too flat!

Ashcroft suggested a low temperature liquid metallic ground state.

We find that the solid is stable below 500K.

Results for state.

H2 liquid

  • Experiment from shock waves

  • Orbital was a bond-centered Gaussian

  • Comparison of g(r) between CEIMC and Car-Parinello MD.

  • Reasonable agreement between methods.

Trial wavefunction for h h 2 transition
Trial Wavefunction for H-H state.2 transition

  • Need a more flexible wavefunction to describe metal-insulator transition

  • We have to find within ~1 sec for 54 atom system.

  • Solve for {fi} in plane-wave basis (PWs)

  • Many basis functions, few {fi}  iterative diagonalization

  • Sampling e-p cusp difficult with PWs: coefficients decay slowly with wavevector

    • Remove e-p cusp from {fi} using

      approximate orbital on the real-space grid.

    • Cusp is fixed analytically in Jastrow term. No basis error.

    • High-k less important. Truncate basis-set after removal. Calculation of Slater determinant much faster.

Plasma phase transition
Plasma Phase Transition state.

  • Study nature of transition from molecular to non-molecular fluid using CEIMC

  • Simulations at T=2000K with P=50-200GPa

Vmc simulations
VMC Simulations state.

Energy from VMC (Jastrow-Slater)

  • 32 atoms

  • 216 twist angles

  • Look at proton-proton correlation function

  • Clear bonding peak.

  • Circles: simulations started from molecular fluid

  • Crosses: from non-molecular fluid

    Clear hysteresis in H-H2 transition.

  • Reptation qmc simulations
    Reptation QMC Simulations state.

    • VMC: Hysteresis; probably 1st order.

    • RQMC: No hysteresis; continuous transition

    An advantage of CEIMC: state.

    Averages are almost free.


    • Path Integrals for ions (particularly for protons or light ions)

      (M time slices to average over.)

    • k-point sampling (integrate over Brillouin zone of supercell). Twist averaged boundary conditions converge much faster than periodic boundary conditions for metals.

      (M k-points)

  • In explicit methods such as LDA, these extra variables will increase the CPU time by M.

  • With QMC there will be little increase in time if imaginary time and/or k are simply new variables to average over. Except for startup time, it just increases dimensionality of integral.

  • Summary
    Summary state.

    • Noisy Monte Carlo can be much more efficient than simply reducing noise to a small value. Add a penalty to the acceptance formula:

      a = min(1,exp[ - E(S)+E(S*) - v/2])

    • CEIMC simulations require accurate wavefunctions.

    • Reptation MC can compute energy differences efficiently.

    • CEIMC does not support liquid hydrogen ground state.

    • DFT-LDA temperature dependence seems seriously in error in atomic phase.

    • CEIMC shows no phase transition from H2 liquid to H liquid.

    • k-point sampling, ionic path integrals, disorder can be added cheaply, making a fast code for hydrogen, as fast as CPMD! Also converges much faster to thermodynamic limit.

    • Can we generalize the methods to heavier elements? Requires fast generation of orbitals-work in progress.