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. Notation • Individual coordinate of a particle ri • All 3N coordinates R= (r1,r2, …. rN) • Total potential energy = V(R) • Kinetic energy : • Hamiltonian :

3. 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) Variational Monte Carlo (VMC)

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

5. 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?

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

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

8. smoothing 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.

9. 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. FN -SJ FN-BF

10. kx 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 PBC TABC

11. Polarization transition 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

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

13. Summary of Variational (VMC) Simple trial function Better trial function applications

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

15. 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:

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

17. 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?

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

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

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

21. Schematic of DMC Ensemble evolves according to • Diffusion • Drift • branching ensemble

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

23. Transient Estimate Approach • () converges to the exact ground state • E is an upper bound converging to the exact answer monotonically

24. Positive walkers Node 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.)

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

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

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

28. “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.

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

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

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

32. 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, …

33. Simple trial function VMC FN Better trial function applications RN Summary of T=0 methods:Variational(VMC), Fixed-node(FN), Released-node(RN)

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

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

36. 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!

37. VMC 300K DMC timestep Tests on bcc hydrogen, N=16 Good results in a few slices p~20.

38. electrons ions 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?)

39. 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”, Z=partition function. Only the difference in energy enters in the acceptance probability.

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

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

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

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

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

45. Nodal surfaces S* overlap s “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.

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

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

48. 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?

49. 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!

50. Ashcroft suggested a low temperature liquid metallic ground state. We find that the solid is stable below 500K.