Coupled Electron-Ion Monte Carlo Study of Hydrogen Fluid

1. Coupled Electron-Ion Monte Carlo Study of Hydrogen Fluid Kris Delaney, David Ceperley University of Illinois at Urbana-Champaign Carlo Pierleoni Universita del l’Aquila, l’Aquila, Italy Supported by NSF DMR 04-04853 and NSF DMR 03-25939 ITR

2. Overview • Features of hydrogen fluid • Simulation Method: • The CEIMC approach • Choosing and sampling the ensemble • Energy differences with noise • The trial wavefunction • Simulating phase transition in finite systems • Results • Plasma phase transition • Finite-size error assessment • Conclusions

3. Why Hydrogen? • Important in planetary science and high-pressure physics • 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 • Possible liquid-liquid phase transition • Possibility of superconducting and superfluid phases • Equation of state not yet fully described

4. Features of Hydrogen Fluid • Current phase diagram: • Open problems: • Liquid-liquid phase transition • Shape of melt curve • Possibility of quantum fluid

5. The CEIMC Approach • Assume Born-Oppenheimer (BO) valid: • Separates electronic and ionic degrees of freedom • Electronic subsystem adiabatically remains in ground-state • May split system into two coupled Monte Carlo simulations: • Classical or quantum (PI) simulation for ions at finite T • T=0K QMC for electrons

6. The CEIMC Approach • Assume Born-Oppenheimer (BO) valid: • Separates electronic and ionic degrees of freedom • Electronic subsystem adiabatically remains in ground-state • May split system into two coupled Monte Carlo simulations: • Classical or quantum (PI) simulation for ions at finite T • T=0K QMC for electrons

7. The CEIMC Approach • Assume Born-Oppenheimer (BO) valid: • Separates electronic and ionic degrees of freedom • Electronic subsystem adiabatically remains in ground-state • May split system into two coupled Monte Carlo simulations: • Classical or quantum (PI) simulation for ions at finite T • T=0K QMC for electrons • BO good for T<<TF(~150,000K)

8. Choice of Ensemble • For thermodynamic observables, first generate ionic configurations to a phase-space probability distribution of some ensemble. • All ensemble averages equivalent in thermodynamic limit • Far from simulation regime • Careful choice required • Use constant T: energy fluctuates in exchange with environmental heat bath. • Use constant N: theoretical simplicity, especially in wavefunction generation. • Choices: NPT or NVT.

9. The NVT Ensemble • We choose NVT for these studies • Pros: • Algorithmically simple • Fewer computations per MC step than NPT • Cons: • May bias transitions involving crystalline phases if simulation cell doesn’t conform OK for fluid studies • Doesn’t capture density fluctuations beyond length-scale of simulation cell: NPT or mVT (GCE) would. Must monitor error due to finite cell size

10. NVT Partition Function • The NVT partition function for classical ions is: • Monte Carlo: • The momentum part is analytically soluble: • Sample only configuration space, not phase-space. • Generate samples according to: • Metropolis: • Propose uniform move in 3D box. • Accept according to: Spatial configuration of ions BO energy difference of ion configurations S,S’

11. The CEIMC Approach • Assume Born-Oppenheimer (BO) valid: • Separates electronic and ionic degrees of freedom • Electronic subsystem adiabatically remains in ground-state • May split system into two coupled Monte Carlo simulations: • Classical or quantum (PI) simulation for ions at finite T • T=0K QMC for electrons • BO good for T<<TF(~150,000K)

12. QMC Energy Differences • Compute DU(S,S’) with either VMC or RQMC • RQMC: projector-based method with no mixed-estimator • More accurate than VMC • Unbiased for energy differences. • Finite-size errors on energy difference are reduced by using twist-averaged boundary conditions • Improves electron kinetic energy part of DU • Does not improve potential energy error • Remaining issues: • Estimate of DU will have a statistical uncertainty • A fast and accurate trial function is required

13. Dealing with Noise • 3 Approaches: • Run QMC for longer. DU noise lower • Reduce noise through correlated sampling methods • Tolerate noise with modified acceptance criterion • Correlated Sampling: • Direct energy difference; noise often large compared with DU • Importance sampling; QMC walker (R) samples combination of S1,S2

14. Penalty Method • Metropolis acceptance ratio is biased if DU has noise • Adjust using “Penalty Method”: • Satisfy detailed balance on average • If DU is Gaussian distributed (CLT) acceptance ratio becomes • Tolerates noisy estimates of energy differences • Extra rejections through noise. • In practice, s estimated → further terms • J. Chem. Phys. 110, 9812 (1999) • Correlated sampling methods still important  noise lowered  fewer noise rejections

15. Trial Wavefunction • BO energy is defined as: • Approximation for Y. Choose Slater-Jastrow: • The {fi}s obey a single-particle equation • Veff can be: • Some parameterized potential. Parameters optimized with EVMC. Fast, but free parameters • Full self-consistent Kohn-Sham potential with approximate Vxc (eg, LDA, GGA) No free parameters, but slow • Choices: • RPA • Optimizable: • 1-body • 2-body • Choices: • Analytic backflow • Gaussian orbitals • (molecular configurations) • DFT-LDA • Optimizable OEP bands where

16. Trial Wavefunction: PW basis • We employ a plane-wave basis to represent the single-body Hamiltonian: • Fourier transform of potential is often possible analytically. For example, OEP potential in PBCs: Cell structure factor. The only S-dependence

17. Trial Wavefunction: Cusp • Solve for {fi} in plane-wave basis (PWs) • Many basis functions, few {fi}  band-by-band iterative diagonalization • Sampling eln-ion cusp difficult with PWs: coefficients decay slowly with wavevector magnitude • Remove eln-ion cusp from {fi} using RPA on the incomplete FFT grid. • Cusp is returned analytically in Jastrow term. No basis error. • High-K less important. Truncate basis-set after removal. Each Slater determinant much faster.

18. CEIMC Summary • Propose move S→S’ • Compute new Y(S’) • Expand h in plane-wave basis • Diagonalize h • Iterate for new h if DFT-LDA • Cusp removal + basis truncation • Compute DU(S,S’) • Correlated sampling • VMC or RQMC (fixed node) • TABC – typically 216 twists, 15,000 total electron moves • Accept/reject • Penalty method • Repeat to sample Boltzmann distribution • Typically 8,000--20,000 moves needed

19. Optimizing OEP • Choose an electron-nuclear potential including effective screening with free parameters: • Yukawa • Gaussian • Yukawa + Gaussian • …? • Variationally optimize free potential parameters (eg. Yukawa screening) • Either on each CEIMC nuclear move (slow!) • A variety of static nuclear configurations, including: • Crystals • Disordered • Molecular + non-molecular states • Range of densities • Try to discover parameter trends

20. Potential Optimization • Parameters optimized for molecular & non-molecular crystals, and disordered systems at densities ranging from rs=1-4 • Bare Coulomb interaction was always optimal or near-optimal

21. Wavefunction Test • Test wavefunctions on static crystals: • Veff = Bare Coulomb or DFT-LDA • Pure Gaussian orbitals • Backflow • Variational principle: lowest energy  best wavefunction

22. Simulating Phase Transitions • 1st-order phase-transitions hard to locate precisely: • Finite-size cell: periodic correlations suppress interface between pure phases close to transition; forces system into a metastable state. • Finite simulation time: probability of changing from one phase to another is related to DF (Helmholtz F for NVT) and height of barrier. • Result: hysteresis; signature of first-order transition. • Higher order transitions: usually no hysteresis . • Close to a critical point: • System behaviour is singular • Fluctuations in density become unbounded. • Finite-cell simulations more unreliable, especially NVT.

23. The Plasma Phase Transition • Study nature of transition from molecular to non-molecular fluid using CEIMC • Simulations at T=2000K with P=50-200GPa

24. Results (DU from VMC) • Simulation details: • 32 atoms, NVT ensemble • T = 2000K • P = 50 – 200 GPa • VMC for DU • 216 twist angles • Circles: simulations started from molecular fluid • Crosses: from non-molecular fluid Clear hysteresis

25. Results (DU from RQMC) • Simulation details: • 32 atoms, NVT ensemble • T = 2000K • P = 50 – 200 GPa • RQMC for DU • 216 twist angles • Circles: simulations started from molecular fluid • Crosses: from non-molecular fluid No hysteresis

26. Molecular Order Parameter • Molecular order parameter, l, defined as: • VMC: Hysteresis; probably 1st order. • RQMC: No hysteresis; continuous transition.

27. Finite-Size Error Assessment VMC RQMC • Pronounced finite-size error • Consistent with first order behaviour • No significant finite-size error • Consistent with continuous transition Circles: 32-atom cell Crosses: 54-atom cell

28. Conclusions • More accurate solution of Hamiltonian changes the nature of the transition from molecular to non-molecular fluid • VMC results indicate first-order at T=2000K • Compatible with CPMD but at different pressure • RQMC results indicate continuous transition at T=2000K • FSE not large for PPT, even when 1st order. No long-range nucleation. • Future work: • Investigate metallization of fluid • Effects of quantum nuclei