Science, Faculty of
Chemistry, Department of
DSpace
UBCV
Hsu, Bernie (Po Cheng)
Master of Science - MSc
University of British Columbia
In recent decades, Ab initio quantum chemistry methods have played a crucial role in elucidating the electronic structures of molecules and materials. Despite their widespread application, the constant emergence of new methods lacking rigorous mathematical foundations raises concerns about their reliability. This thesis seeks to bridge this gap by addressing the electronic structure problem through a mathematical lens, leveraging existing numerical analysis results to address the electronic structure problem. Focusing on wavefunction methods, specifically Variational Monte Carlo, our investigation delves into the distinctive features of the Schrodinger equation across different chemical systems, encompassing variations between molecules and materials. We conduct a numerical analysis to derive the convergence of the Ewald Summation of the simulated potential. Subsequently, we employ the Atomic Cluster Expansion to parameterize the trial wavefunction, demonstrating its adaptability to different chemical systems. Building upon existing numerical analyses of the Atomic Cluster Expansion, we argue that our methodology not only provides superior parameterization flexibility compared to existing methods but also ensures a systematic and mathematically sound process for deriving the parameterizations. Finally, we implement simulations to assess the effectiveness of our proposed wavefunction architectures in the 1D uniform electron gas, a system allowing us to probe electronic correlation effects. Through this comprehensive approach, our research aspires to contribute to the development of more robust and mathematically justified methods in the field of Ab initio quantum chemistry.
Wavefunction Parameterizations Using The Atomic Cluster ExpansionbyBernie (Po Cheng) HsuHBSc, The University of Toronto, 2021A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMaster of ScienceinTHE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES(Chemistry)The University of British Columbia(Vancouver)December 2023© Bernie (Po Cheng) Hsu, 2023The following individuals certify that they have read, and recommend to the Faculty of Graduate and Post-doctoral Studies for acceptance, the thesis entitled:Wavefunction Parameterizations Using The Atomic Cluster Expansionsubmitted by Bernie (Po Cheng) Hsu in partial fulfillment of the requirements for the degree of Master ofScience in Chemistry.Examining Committee:Christoph Ortner, Professor, Mathematics, UBCCo-SupervisorYan Alexander Wang, Associate Professor, Chemistry, UBCCo-SupervisorRoman Krems, Professor, Chemistry, UBCSupervisory Committee MemberAlireza Nojeh, Professor, Electrical and Computer Engineering, UBCSupervisory Committee MemberiiAbstractIn recent decades, Ab initio quantum chemistry methods have played a crucial role in elucidating theelectronic structures of molecules and materials. Despite their widespread application, the constant emer-gence of new methods lacking rigorous mathematical foundations raises concerns about their reliability.This thesis seeks to bridge this gap by addressing the electronic structure problem through a mathematicallens, leveraging existing numerical analysis results to address the electronic structure problem.Focusing on wavefunction methods, specifically Variational Monte Carlo, our investigation delves intothe distinctive features of the Schro¨dinger equation across different chemical systems, encompassing vari-ations between molecules and materials. We conduct a numerical analysis to derive the convergence of theEwald Summation of the simulated potential. Subsequently, we employ the Atomic Cluster Expansion toparameterize the trial wavefunction, demonstrating its adaptability to different chemical systems. Buildingupon existing numerical analyses of the Atomic Cluster Expansion, we argue that our methodology not onlyprovides superior parameterization flexibility compared to existing methods but also ensures a systematicand mathematically sound process for deriving the parameterizations. Finally, we implement simulations toassess the effectiveness of our proposed wavefunction architectures in the 1D uniform electron gas, a sys-tem allowing us to probe electronic correlation effects. Through this comprehensive approach, our researchaspires to contribute to the development of more robust and mathematically justified methods in the field ofAb initio quantum chemistry.iiiLay SummaryWavefunction methods continue to stand out as one of the most influential strategies for addressing theelectronic structure problem. Given that this problem is fundamentally rooted in the Schrodinger equation,we contend that any effort to inspire novel methodologies must be firmly grounded in sound mathematicalfoundations - a facet often absent in many recent proposals. This work strives to bridge this gap by establish-ing a robust mathematical framework, providing a solid footing for the development of new wavefunctionmethods. We will also draw comparisons with well-established methods to underscore the importance ofour approach.ivPrefaceThis thesis is an original, unpublished, independent work by the author, Bernie Hsu, under the supervi-sion of Dr. Christoph Ortner and Dr. Alexander Wang.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Electronic Structure Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.1 Variational Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.1.1 1D Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.2 Variational Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2.1 Outline and Disclaimer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.3 Miscellaneous Math . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 Periodic Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93.1 Periodization of the Potential . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93.2 Ewald Summation (3D) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123.3 Periodization of 1D UEG and Its Ewald Sum . . . . . . . . . . . . . . . . . . . . . . . . . 133.4 Proofs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213.4.1 (3.13) Satisfies Poisson Equation (3.9) . . . . . . . . . . . . . . . . . . . . . . . . . 213.4.2 Derivation of 1D Ewald Sums . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.4.3 Proof of Theorem 3.3.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 264 Wavefunction Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31vi4.1 Conventional Quantum Chemistry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 314.1.1 Hartree-Fock Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 334.1.2 Cusp Condition and The Helium Atom . . . . . . . . . . . . . . . . . . . . . . . . 344.2 Explicitly Correlated Wavefunctions - A Brief Survey . . . . . . . . . . . . . . . . . . . . . 374.3 Explicitly Correlated Wavefunctions - A Closer Look . . . . . . . . . . . . . . . . . . . . . 384.4 ACE and Its Adaptation to Correlated Wavefunctions . . . . . . . . . . . . . . . . . . . . . 404.4.1 ACE Parameterization of Jastrow Factor . . . . . . . . . . . . . . . . . . . . . . . . 404.4.2 Rotationally Invariant Jastrow Factor . . . . . . . . . . . . . . . . . . . . . . . . . 424.4.3 Explicitly Incorporating Kato Cusps . . . . . . . . . . . . . . . . . . . . . . . . . . 434.4.4 Backflow Transformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 444.4.5 A Comparison Between CASINO and ACE Wavefunction . . . . . . . . . . . . . . 474.5 Proofs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 494.5.1 Diagonalizability of Hartree-Fock UEG . . . . . . . . . . . . . . . . . . . . . . . . 494.5.2 Useful Tools for Cusp Constraint Calculation . . . . . . . . . . . . . . . . . . . . . 504.5.3 Cusp Constraints Calculation for (4.45) . . . . . . . . . . . . . . . . . . . . . . . . 505 Numerical Implementation of VMC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 565.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 565.2 Recap and Basic Set-ups . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 575.3 Wavefunction Architectures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 585.3.1 ACE-parameterized Backflow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 595.3.2 ACE-parameterized Backflow with Coordinate Difference . . . . . . . . . . . . . . 595.3.3 Wigner Ansatz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 605.4 Hartree-Fock Theory of 1D UEG - Benchmarking ACE . . . . . . . . . . . . . . . . . . . . 615.5 Ground State Computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 625.5.1 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 635.6 Proofs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 705.6.1 Minimal HF Ground State of 1D UEG . . . . . . . . . . . . . . . . . . . . . . . . . 706 Conclusion and Outlooks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74viiList of FiguresFigure 3.1 Comparison between Vb as naively periodized in (3.27) and Ewald summation . . . . . . 15Figure 3.2 Vb as defined in (3.27). We treat each summand as O(1) operation, absorbing also thequadrature. Theoretical bound of error is computed via the formula 12L(12b2L2+ 1)1(M−1)2in Theorem 3.3.2. Note that the formula predicts that, for a fixed L, 12b2L2+ 1 ∼ 1 forany small enough b’s. The fact that all three lines on the right look essentially identicalverify this prediction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16Figure 3.3 Real space part Vb − Vb ⋆L gb. Each summand treated as O(1)-cost. Theoretical boundof error is computed via b2(2 + 4√π)1L31(M−1)2 from Theorem 3.3.2 . . . . . . . . . . 17Figure 3.4 Reciprocal space part Vb ⋆L gb. Each summand treated as O(1)-cost. Theoretical boundof error is computed via 12πb[Lb2πK]3exp(− [ b2πKL ]2) from Theorem 3.3.2. We havealso plotted 12πb[Lb2π]3exp(− [ b2πKL ]2) as the dashed-dotted-line, which demonstratesthat the polynomial dependence K−3 has negligible influence on the theoretical error. . . 18Figure 3.5 Convergence of Total Evaluation of Ewald Summation . . . . . . . . . . . . . . . . . . 20Figure 4.1 The error in the electronic energy (in Hartrees) of the ground-state helium atom. Theerror is plotted as a function of the number of terms in the expansions. This image isadapted from Helgaker et al. [1] Figure 7.6. . . . . . . . . . . . . . . . . . . . . . . . . 35Figure 5.1 VMC Convergence for rs = 0.5, b = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . 65Figure 5.2 VMC Convergence for rs = 1, b = 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . 66Figure 5.3 VMC Convergence for rs = 1, b = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 67Figure 5.4 VMC Convergence for rs = 5, b = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 68Figure 5.5 VMC Convergence for rs = 5, b = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 69viiiList of AbbreviationsACE Atomic Cluster ExpansionCC Coupled ClusterCI Configuration InteractionDFT Density Functional Theoryee electron-electronen electron-nucleuseen electron-electron-nucleusFCI Full Configuration InteractionHF Hartree-FockLHS Left Hand SideNN Neural-NetworkPBC Periodic Boundary ConditionRHS Right Hand SideSTO Slater Type OrbitalVMC Variational Monte CarloixAcknowledgementsI would like to thank Dr. Christoph Ortner for giving me such a fun project. It is truly a privilege that Ihave been able to explore as much as I selfishly want, while getting helpful guidance whenever I am stuckin the research.I would like to thank Dr. Alex Wang for patiently teaching me chemistry. It must have been difficult tocommunicate with me when I don’t even know what sp-hybridization is.I would like to thank Nancy Qiu for being such a great friend ever since highschool, and for stepping inat that crucial moment in my academic career.I would like to thank Jerry Ho and DeXuan Zhou for all the coding help. This project literally wouldhave been impossible if I had to write the scripts myself.Special thanks to Dr. Ariel Zhitnitsky for introducing me to the world of quantum field theory. I aspireto one day fit your favorite topic into my research.Finally, I thank Dr. Roman Krems and Dr. Alireza Nojeh for agreeing to be my examining committeemembers. The suggestions provided for my thesis have been enjoyably thought-provoking.xChapter 1IntroductionThe accurate prediction of electronic properties is essential across various scientific domains, encom-passing chemistry, physics, materials science, and biochemistry. Ab initio quantum chemistry methods,rooted in solving the Schro¨dinger equation without resorting to empirical parameters, have proven indis-pensable in unraveling the electronic structures of molecules and materials. Over recent decades, substantialstrides have been made in advancing ab initio methods [2, 3], resulting in noteworthy enhancements incomputational accuracy and efficiency.Computational methods in this realm can be categorized into two main schemes: the wavefunction andthe density functional approach. The former directly addresses the Schro¨dinger equation and its solution -the wavefunction. The earliest and simplest yet foundational wavefunction method is the Hartree-Fock (HF)model [2, 4], approximating the many-electron wavefunction as a single determinant (see Section 4.1.1).HF was groundbreaking in establishing a systematic framework for comprehending molecular structuresand properties. Nevertheless, HF falls short of achieving the gold standard of chemical accuracy (1.6milli-Hartree) [2] and is inadequate in describing systems with strong electron correlation, bond break-ing, and dispersion interactions, constraining its applicability to specific chemical phenomena. The secondscheme, known as Density Functional Theory (DFT), focuses on electron density [3, 5] rather than thewavefunction. Represented as a map of three spatial variables (fixing spins), as opposed to the wavefunc-tion’s 3 × (number of electrons), DFT often offers sufficient accuracy at a computational cost lower thanwavefunction methods. However, similar to HF, DFT faces challenges in accuracy for specific systems andstruggles in practice to capture strong correlation phenomena.Hartree-Fock and Density Functional Theory serve as foundational pillars in investigating electronicproperties, and the pursuit of enhanced accuracy within each scheme constitutes the core of contemporaryab initio research. In response to HF’s limitations, novel post-Hartree-Fock methods, surpassing the single-determinant architecture, have emerged. Simultaneously, to compute electron density with greater accuracy,there is a continuous stream of functionals or methodologies governing their formulation. Nevertheless,the introduction of a new wavefunction architecture or density functional often lacks a satisfactory mathe-matical analysis. Proposals are typically guided by physical intuition and performance statistics, and thissemi-empirical approach raises concerns regarding reliability and transferability across different chemicalsystems.1This thesis narrows its focus to one wavefunction method, specifically Variational Monte Carlo (VMC),and undertakes a mathematical approach to the electronic structure problem. It involves an analysis ofdifferences and similarities arising from the Schrodinger equation in distinct chemical systems, includingvariations between molecules and materials, and between 3-dimensional and 1-dimensional systems. Sub-sequently, we will motivate new wavefunction architectures leveraging established mathematical findings.In particular, we will parameterize the trial wavefunction using the Atomic Cluster Expansion (ACE) [6, 7],and show how it adapts to the technicalities arising from the aforementioned differences and similarities. Tosubstantiate this motivation, we draw comprehensive comparisons between our proposed parameterizationand existing wavefunction architectures such as Coupled Cluster [1, 4], CASINO code [8], and Neural-Network (NN) wavefunctions [9, 10]. Utilizing existing numerical analyses that have been done on ACE[7, 11], we argue that it not only offers the parameterization flexibility on par with Neural-Network wave-functions, surpassing the CASINO code, but also that the procedures leading to the parameterizations issystematic and mathematically sound, like conventional quantum chemistry methods (e.g, FCI, CoupledCluster, etc). In the final chapter, we will explore various wavefunction architectures based on the AtomicCluster Expansion, subjecting them to empirical tests to assess their accuracy. Through the amalgamationof physical insights and mathematical rigor, we aim to conceive wavefunction architectures from whichrigorous accuracy results can be stated and proven.2Chapter 2Electronic Structure Problem2.1 Variational PrincipleThe problem under consideration is the computation of electronic ground state wavefunction and itsenergy. That is, letting the Hamiltonian H be defined byH = −12N∑i=1∆i + V (r1, . . . , rN ),V (r1, . . . , rN ) =12N∑i=1N∑j=1j ̸=i1|ri − rj | −N∑i=1M∑A=1ZA|ri −RA| +12M∑A=1M∑B=1B ̸=AZAZB|RA −RB| , (2.1)we want to find a function Ψ0 of electrons coordinates and the lowest number E0 such thatHΨ0 = E0Ψ0, (2.2)and Ψ0(xp1 , . . . ,xpN ) = (−1)pΨ0(x1, . . . ,xN ), ∀p ∈ SN .Here, N denotes the number of electrons, M number of nuclei, and ∆i =∇2i . The electron coordinates arewritten in the form xi = (ri, σi) ∈ R3×{↑, ↓}. ri andRA are spatial positions of respectively the electronsand nuclei. The Hartree units are assumed throughout the entire thesis. The second condition is PauliExclusion principle. SN is the collection of all permutations p of the tuple (1, . . . , N), and (−1)p denotesthe sign of the permutation. Under the Born-Oppenheimer approximation[1, 4], the positions R1, . . . ,RMare fixed parameters that we need to initialize prior to every numerical simulation. This work also considersthe uniform electron gas (UEG) model, where the discrete ionic contribution {(RA, ZA)|A = 1, . . . ,M} isreplaced by a uniform distribution {(R, Z)|R ∈ Ω} in a box Ω that contains all the nuclei. The Hamiltonian3is then given byH = −12N∑i=1∆i + VUEG(r1, . . . , rN ),VUEG(r1, . . . , rN ) =12N∑i=1N∑j=1j ̸=i1|ri − rj | −N∑i=1∫ΩdR|Ω|Z|ri −R| +12∫ΩdR|Ω|∫ΩdR′|Ω|Z2|R−R′| , (2.3)where |Ω| denotes the volume of the box. To avoid confusion, we note that the nomenclature “UniformElectron Gas” refers strictly to the Hamiltonian (2.3), and is extracted from standard condensed mattertextbooks [12, 13]. We are not making physical claims about the electrons, and the electron density neednot be uniform.Let us remark that there is no explicit removal of self-interactions in the ion-ion contribution. To contrast,in the general case (2.1), the sum requires B ̸= A. Indeed, the functionr 7→∫ΩdR1|r −R| (2.4)is continuous as remarked in [14, VI.2], hence no need for such removal. This non-removal of self-interaction will later play a role in our analysis upon periodic boundary condition (see (3.19)) in the UEGcase. For the rest of this thesis, Ω = [−L/2, L/2]3 or [−L/2, L/2] depending on the context, and thecontext will be clear.In practice, it is difficult to solve the eigenvalue problem directly, and practically impossible as thenumber of electrons increases. We therefore need to resort to some form of approximation. One way to doso is to construct a wave-function ansatz and appeal to the following theorem[1, 4].Variational Principle. Let Ψ be any trial wave-function and E0 be the ground state energy of H . ThenE0 ≤ E(Ψ) := ⟨Ψ|H|Ψ⟩⟨Ψ|Ψ⟩ =∫dx1 . . . dxNΨ∗(x1, . . . ,xN ) (HΨ) (x1, . . . ,xN )∫dx1 . . . dxN |Ψ(x1, . . . ,xN )|2 . (2.5)Moreover,E0 = minΨE(Ψ).In this thesis, we only consider the cases where spatial dimension is either be 3 or 1 (and we replaceri → xi in the later case). The integral thus means∫dx =∑σ=↑,↓∫dr. We leave the case of 2D modelfor future research.2.1.1 1D ModelIn this subsection we will introduce the 1D model implemented in the existing VMC simulations [15,16], and later in Chapter 5.By a 1D model, we mean particle motions being confined on one axis, each hence possessing only 14spatial coordinate (r → x). They interact through an effective potential (instead of the bare Coulomb 1/r)motivated by the following physical intuition: A harmonic potential V⊥ = (y2 + z2)/8b4, where b controlsthe “width” of the confinement, is imposed along the directions transverse to the x axis. If the electrons aresufficiently confined in the transverse plane, the transverse component of electrons can then be assumed tobe in the ground state,θb(y, z) =1√2πb2exp(−y2 + z24b2),of this potential, leading to a wavefunction of the form Ψ(x, y, z) = Ψ(x)θb(y, z). The effective 1D inter-action vb is obtained by integrating out the transverse part of the wavefunction,vb(x) =∫dr⊥dr′⊥|θb(r⊥)|2|θb(r′⊥)|2[x2 + |r⊥ − r′⊥|2]1/2=√π2bexp(x24b2)erfc( |x|2b)(2.6)This yields the effective 1D interaction between electrons. One can further show thatvb(x)≈ 1/|x| when x is large= √π/2b x = 0. (2.7)That is, the interaction is Coulombic when a pair of electrons are far apart, and finite when they are close. Itis also worthwhile to note that the original derivation [17] in fact first worked out the Fourier transform ofthis model,vˆb(k) = E1([bk]2) exp([bk]2). (2.8)Note that, due to its 1/|x| behavior at large x, the Fourier transform is defined in the L2 sense.This summarizes all the background we need from this model. For a much more involved “derivation”of this electron-electron interaction and discussion of the physics, we refer the readers to [17, 18].2.2 Variational Monte CarloProvided that a wave-function Ψθ is parameterized, the subsequent step is then to minimize the energyfunctional E(Ψθ) with respect to parameters θ. Plenty of methods have been proposed to numerically carryout such minimization. In this thesis, we consider the Variational Monte-Carlo (VMC) method.The VMC method first reformulates the variational problem by writingE(Ψθ) =∫dx1 . . . dxN(HΨθ) (x1, . . . ,xN )Ψθ(x1, . . . ,xN )|Ψθ(x1, . . . ,xN )|2⟨Ψθ|Ψθ⟩=∫dx1 . . . dxNEL(x; Ψθ)Π(x), (2.9)5where the local energy is defined byEL(x; Ψ) :=(HΨ) (x)Ψ(x). (2.10)Equation (2.9) can then be seen as the average of local energy EL over the probability distribution Π :=|Ψθ|2/ ⟨Ψθ|Ψθ⟩. Monte-Carlo integration [19] then allows us to approximate the integral by a sum oversamples. That is, ∫dx1 . . . dxNEL(x; Ψθ)Π(x) ∼ 1|S|∑x∈SEL(x; Ψθ), (2.11)and the error would fall as O(1/|S|1/2). In the above, S is a finite collection of points randomly sampledfrom Π. The larger the number |S| of samples, the better the accuracy of the approximation. Typically,such sample S of arbitrary distribution is generated by Metropolis-Hastings Algorithm [19, 20]. StochasticGradient Descent [21] or a variant thereof is then applied to the approximate expression (2.11). That is, theoptimization is facilitated by the approximate gradient∇θE(Ψθ) ∼ 1|S|∑x∈S∇θEL(x; Ψθ), (2.12)where∇θEL(x; Ψθ) is an expression that one can write down using ordinary rules of calculus.2.2.1 Outline and DisclaimerSection 2.2 summarizes all the necessary knowledge on numerical sampling and optimization that thisthesis concerns. The rest of the thesis is dedicated to “everything else” in VMC. We will discuss in Chapter3 additional details, such as periodic boundary conditions and Ewald summation, not present in molecularsimulation but indispensable to treat extended systems. In Chapter 4, an overview of wavefunction ansatz isgiven. Going from the point of view of conventional quantum chemistry, to existing VMC code in the liter-ature, we will systematically motivate new parameterizations of Ψθ. Finally, Chapter 5 exhibits numericalimplementation and appraises the performances of the proposed trial wavefunctions.Much of the emphasis of this thesis is focused on theory. By an nD code, we mean the numerical imple-mentation of an n-dimensional model, and as of the time of writing, the 3D code is still under developmentand the author plays no role in this front of the endeavor. As such, although the theoretical portion is generalwith respect to spatial dimensions, all the numerical implementation is done with 1D code.The rest of the thesis also contains several theorems to substantiate our theory. Proofs are written tojustify them, unless the results can be readily found in literatures and such are gathered in the subsequentSection 2.3. However, these proofs are long enough to disturb the flow of presentation. We therefore gatherthem at the end of each chapter. On another note, all the “proofs” in this thesis are only rigorous up to andexcluding distribution theory (see for instance [22, ch.9]) in order to avoid excessive technicalities. We leavethe refinement of rigor for future work.62.3 Miscellaneous MathThroughout this thesis we use the following definitions/facts/conventions. The special functions we useextensively areerf(x) =2√π∫ x0dte−t2, (2.13)erfc(x) = 1− erf(x) = 2√π∫ ∞xdte−t2, (2.14)E1(x) =∫ ∞xdte−ttfor x > 0. (2.15)Given a function f , L-periodic in each axis, the Fourier coefficient is defined viaf˜(k) =∫Ωdre−ik·rf(r) for k ∈ (2π/L)n Zn, (2.16)Ω := [−L/2, L/2]n, (2.17)where n is the spatial dimension of interest. In this thesis, n = 3 or 1. The Fourier series is given byf(x) =1Ln∑k∈(2π/L)nZnf˜(k)eik·r. (2.18)Given any function f , the Fourier transform is defined byfˆ(k) :=∫Rndre−ik·rf(r) for k ∈ Rn, (2.19)and the inverse Fourier transform byf(x) =1(2π)n∫Rndkfˆ(k)eik·r. (2.20)We define the convolution between f and g over the box Ω byf ⋆L g(r) :=∫Ωdr′f(r′)g(r − r′). (2.21)Poisson Summation. Suppose a continuous f , not necessarily periodic, satisfies |f(r)| ≤ C(|r|+ 1)−n−ϵand |fˆ(k)| ≤ C(|k|+ 1)−n−ϵ for some C, ϵ > 0. Then∑m∈Znf(r −mL) = 1Ln∑k∈(2π/L)nZnfˆ(k)eik·r. (2.22)Fourier Transform of Convolution. Given an L-periodic function f , and a function g that is not neces-7sarily L-periodic, their convolution evaluated at k ∈ (2π/L)n Zn is given by(f ⋆L g) (k) = f˜(k)gˆ(k) (2.23)Given a function f that has domain in 3D, let us write the position in spherical coordinates r = rrˆ. Wedefine its spherical average at radius r byf¯(r) := −∫S2drˆf(rrˆ) :=14π∫ π0sin(θ)dθ∫ 2π0dϕf(rrˆ) (2.24)Dominated Convergence Theorem (DCT). Consider f : X × [a, b] → R, (−∞ < a < b < ∞) andsuppose that f(·, t) : X → R is integrable for each t ∈ [a, b](a) Suppose also that there exists an integrable g such that |f(x, t)| ≤ g(x) for all x, t, and if limt→t0 f(x, t) =f(x, t0) for every x, thenlimt→t0∫Xdxf(x, t) =∫Xdx limt→t0f (x, t) (2.25)(b) If ∂f/∂t exists and there exists an integrable g such that | ∂f/∂t | ≤ g(x) for every x and t. Thenddt∫Xdxf(x, t) =∫Xdx∂f∂t(x, t) (2.26)As X that shows up in this thesis is usually a closed and bounded set, the g that shows up in this thesiscan most of the time be taken as constantWe say a function f is C1 around r if all its partial derivatives exist and are continuous in a smallneighborhood containing r.8Chapter 3Periodic Boundary Conditions3.1 Periodization of the PotentialAny realizable extended system, be it gas, liquid, or solid state, involves an Avogadros’s number ofparticles. Numbers of electrons and nuclei, N and M in equation (2.1), are thus O(1023). In practice, it isimpossible to simulate that many particles in ab initio calculations. Instead, we define a simulation supercellΩ and impose the Periodic Boundary Condition (PBC) on the potential term V . Afterwards, we run severalnumerical simulations with increasing number of particles (instantiated in Ω), and estimate the energy ofO(1023) particles by an extrapolation N,M → ∞. In principle, such extrapolation can also be computedwith open or vacuum boundary, the canonical choice for molecular computation. The advantage of periodicboundary over open or vacuum boundary on V is that it alleviates the spurious boundary effects and theextrapolation error, incurred due to the discrepancy between real and simulation particle numbers, whereN,M = O(10) to O(102). To detail this periodization procedure, we will follow the treatment in [23, II.B]and fill in some details along the way. Note that this scheme is general and not restricted to VMC.The first step is to rewrite the potential V , as given by equation (2.1), into a more convenient form. Thisis exhibited byV (r1, . . . ,R1, . . .) =12(N∑i=1(−1)ϕe(ri) +M∑A=1ZAϕA(RA)), (3.1)ϕe(r) = limr′→r(ϕ(r′)− (−1)|r′ − r|), (3.2)ϕA(R) = limr′→R(ϕ(r′)− ZA|r′ −R|), (3.3)ϕ(r′) =N∑j=1(−1)|r′ − rj | +M∑B=1ZB|r′ −RB| . (3.4)ϕ as expressed in (3.4) is physically interpreted as the electrostatic potential afforded by the total chargedistribution. We have made explicit the number that represents charge. Namely, −1 and ZA. By removingself-interaction, (3.2) and (3.3) each computes the interaction of the labelled point charge with its surround-9ings, constituted of other charges. (3.1) then sums each point electrostatic contribution. The particularconvenience in ϕ is founded in the fact that it is also the solution to the electrostatic (Poisson) equation ofpoint sources  −∆ϕ = 4π(∑Nj=1(−1)δrj +∑MA=1 ZAδRA)ϕ(r)→ 0 as r →∞,(3.5)where δrj denotes the point source density (a.k.a Dirac delta measure) localized at the particle coordinate rj .This solution can furthermore be viewed, through principle of superposition, as constructed by the Green’sfunction G of the (signed) Laplace operator. That is, G is defined by −∆G = 4πδ0G(r)→ 0 as r →∞, (3.6)and it is well-known that the unique solution to this equation is given by the Coulomb potential,G(r) =1|r| . (3.7)One natural way to introduce the PBCs then arises. We replace the open boundary condition in equation(3.6) by PBC, thus redefining the electrostatic point source. This shifts our interest to the Green’s functionGPBC defined by  −∆GPBC = 4πδ0 in ΩGPBC(r′ + nL) = GPBC(r′) ∀n ∈ Z3. (3.8)However, applying Fourier transform to this equation would give k2G˜PBC(k) = 4π for any k ∈ (2π/L)3Z3.This is a contradiction at k = 0 and hence the solution does not exist. The remedy [14, IV.2][24, AppendixA] is to screen the point charge by a cancelling uniform background equal in charge magnitude. Mathemat-ically, this means adding a constant term to the Poisson source term,−∆GPBC = 4π(δ0 − 1/L3)in ΩGPBC(r′ + nL) = GPBC(r′) ∀n ∈ Z3.(3.9)The solution now exists and is unique up to a constant. This can again be seen through applying the Fourier10Transform.G˜PBC(k) = 4π/k2 k ∈ (2π/L)3 Z3\{0}undetermined k = 0,(3.10)=⇒ GPBC(r′) = 1L3G˜PBC(0) +∑k ̸=04πk2eik·r′ . (3.11)It is mentioned in Lieb and Simons [14, IV.2] that different choices of such (Madelung) constant termsentail different physical interpretation. Rather than delving into the physics debate, we will simply take thecanonical choice (as did Lieb and Simons [14])G˜PBC(0) = 0. (3.12)With this choice, it can be demonstrated that the Green’s function (3.11) can be re-expressed in a morephysically appealing form,GPBC(r′) =∑n∈Z3(1|r′ − nL| −∫Ωdr′′L31|r′ − r′′ − nL|). (3.13)Two arguments are given in 3.4. This form will later be the starting point of the periodization of the 1Dmodel. Here we particularly note that (3.13) is simply the potential of single point charge (3.6) (and itsperiodic images), screened by that of the uniform background (and its periodic images.) Moreover, we cannow periodize ϕ in (3.4) through replacing G by its periodized analog,ϕPBC(r′) =N∑j=1(−1)GPBC(r′ − rj) +M∑B=1ZBGPBC(r′ −RB), (3.14)as done in [23]. Substituting back into V , we then obtainV (r1, . . . ,R1, . . .) =12N∑i=1N∑j=1j ̸=i(−1)2GPBC(ri − rj) +N∑i=1M∑A=1(−1)ZAGPBC(ri −RA) (3.15)+12M∑A=1M∑B=1B ̸=AZAZBGPBC(RA −RB) + N(−1)2VMad2+VMad2M∑A=1Z2A,VMad = limr→0(GPBC(r)− 1r). (3.16)Let us note that the net charge∑Ni (−1) +∑MA=1 ZA in the simulation box is assumed to be 0 due to thecancelling uniform background introduced to periodize V .11The treatment of the UEG case follows the same strategy. We rewrite potential (2.3) viaVUEG(r1, . . . , rN ) =12(N∑i=1(−1)ϕe(ri) +∫ΩdRL3Zϕ(R))(3.17)ϕe(r) = limr′→r(ϕ(r′)− −1|r′ − r|)(3.18)ϕ(r′) =N∑j=1−1|r′ − rj | +∫ΩdRL3Z|r′ −R| . (3.19)Notice here we do not stipulate a removal of Z-self-interaction as in (3.3) precisely due to the remarkfollowing (2.3). Carrying out the same replacement procedure as in (3.14), we arrive atϕPBC(r′) =N∑j=1(−1)GPBC(r′ − rj) +∫Ωd3RL3ZGPBC(r′ −R), (3.20)where GPBC is as given in (3.13). A simplification is now in play. Due to our choice of constant term in(3.11), ∫ΩdRGPBC(r′ −R) = G˜PBC(0) = 0.This givesϕPBC(r′) =N∑j=1(−1)GPBC(r′ − rj), (3.21)VUEG(r1, . . . ,R1, . . .) =12N∑i=1N∑j=1j ̸=i(−1)2GPBC(ri − rj) + N(−1)2VMad2. (3.22)3.2 Ewald Summation (3D)Now that we have studied the theoretical groundwork of the periodized potential, the next step is toconsider the practical aspect in its numerical evaluation. Since the technicalities are rather involved in the3D case, we will only sketch the steps leading to the results, and give a more careful account of the 1Dcounterpart later. The 3D results here are mainly drawn from [25, 7.5.3].The evaluation of the periodized potential V pivots on the evaluation of GPBC. This can be done eitherin real space via (3.13), or in reciprocal space via (3.11). Both expressions are infinite sums that need tobe truncated for numerical evaluations, thus committing errors. As mentioned in [14], the truncated (3.13)which we may write as∑|n|<N only ensures a convergence of some inverse power ofN . This is consideredslow, and on top of summation truncation, the integral of 1/| · | needs to be evaluated through numericalquadrature. This places further computational burden. The Fourier transformed version in (3.11) wouldnot help either because the O(1/k2) integrand also amounts to an error decay of O(1/k). The goal is toreach some evaluation scheme with exponential efficiency. This calls for the alternative evaluation scheme -12Ewald summation [26][25, 7.5.3]. The idea is to decompose the potential byGPBC = GPBC ⋆L gκ + (GPBC −GPBC ⋆L gκ) (3.23)gκ(r) :=4κ3√π∑n∈Z3exp(−κ2|r − nL|2), (3.24)where “⋆L” denotes the convolution over Ω (see (2.21)), and κ is a positive parameter to be tuned. Noticethat gκ is a sum of normalized Gaussians and the summation is over multiples of box dimension. This makesgκ Ω-periodic.We can think that we have further separated charge distribution δ0 − 1/L3 into gκ − 1/L3 (Gaussianscreened by uniform background) plus δ0−gκ (point source screened by Gaussian). Indeed,G = GPBC⋆Lgκis a solution to PBC Poisson equation with source gκ− 1/L3, while G = GPBC−GPBC ⋆L gκ is a solutionto that with source δ0 − gκ.Through some techniques from Fourier theory, one can expand GPBC ⋆ gκ in reciprocal space andGPBC −GPBC ⋆L gκ in real direct lattice space to writeGPBC ⋆L gκ(r) =1L3∑k∈(2π/L)3Z34πk2exp(−k2/(4κ2) + ik · r) (3.25)GPBC(r)−GPBC ⋆L gκ(r) = − πκ2L3+∑n∈Z3erfc (κ|r − nL|)|r − nL| , (3.26)where erfc is a special function defined in (2.14). Notice now all the terms involved in (3.25) and (3.26)no longer require explicit construction of numerical quadrature to evaluate1. Moreover, it can be shownthat the summand in GPBC ⋆L gκ behaves asexp(−|k|2/(4κ2))√π|k|2 , while those of GPBC −GPBC ⋆L gκ behavesasexp(−κ2|nL|2)√π|nL|2 . This says that the sum over n converges fast when κ is large, whereas the sum over kconverges better when κ is small. The parameter κ is to then be chosen to give good compromise betweenthe convergence between the two sums. The Ewald summation thus outperforms the original evaluationscheme upon finding an appropriate such compromise.3.3 Periodization of 1D UEG and Its Ewald SumIn this section, we will periodize the 1D pairwise interaction introduced in (2.1.1) and present the Hamil-tonian of 1D UEG in (3.28). We will perform a numerical analysis of the peiodized potential of 1D UEG.1D Ewald sums mirroring (3.25) and (3.26) are derived in (3.30) and (3.29).To set foot in the computation of 1D UEG model, like the 3D setting, we need a periodized pairwiseinteraction. To that end, we start by interpreting the bare particle interaction vb(x) (2.6) as given rise by1or at least they can be efficiently evaluated by packages designed for special functions, e.g [27]13point charge. Mirroring how G in (3.5) is periodized into GPBC as given by (3.13), we naively stipulateVb(x) =∑m∈Z(vb(x−mL)−∫ L/2−L/2dyLvb(x− y −mL))(3.27)and the Hamiltonian, mirroring (3.22), is given byH = −12N∑i=1∂2∂x2i+ VUEG(x1, . . . , xN ),VUEG(x1, . . . , xN ) =12N∑i=1N∑j=1j ̸=iVb(xi − xj) + NVMad2,VMad = limx→0(Vb(x)− vb(x)) = limx→0Vb(x)−√π2b. (3.28)This is precisely the Hamiltonian the two existing Monte-Carlo treatments [15, 16] implement. Again as inthe 3D case, we decompose Vb into Ewald sums.Theorem 3.3.1. Let gb be the periodized 1D Gaussian defined in (3.41), E1 be the exponential integralgiven in (2.15). By writing Vb = Vb ⋆L gb + (Vb − Vb ⋆L gb), we can evaluate each term byVb ⋆L gb =2L∞∑k=1E1([2πbkL]2)cos(2πkxL), and (3.29)Vb − Vb ⋆L gb =∑m∈Z(vb(x−mL)− 1|x−mL| erf( |x−mL|2b)). (3.30)Proof. Section 3.4.2.Equations (3.29) and (3.30) together comprise precisely the evaluation carried out in [15, 16]. As asanity check, we verify that Ewald summation does agree with the naive periodization (3.27) by directlyplotting in Figure 3.1, where (3.27) is plotted as dots to compromise the costly quadrature. Notice thatthe confinement strength b now additionally takes on the role of κ in (3.25) and (3.26). However, since italso possesses physical meaning, we cannot exactly tune it at will and extract the most optimal compromisebetween the reciprocal and real space sums. Nonetheless, the Ewald summations do still avoid the need touse numerical quadratures, and provide efficient evaluation scheme.Recall what is missing in Section 3.2 is not only the derivation of the Ewald sums, but also a clearstatement on the rate of decay of error of each infinite sum with respect to truncation, as well as the gainsby going from the original to Ewald sums. Such considerations were avoided due to the extra technicaldifficulties that 3D grid space can cause (see for instance [28].) Instead, we compensate with the 1D version,as it is what we will implement in this work anyway. With this in mind, we hereby state the following.Theorem 3.3.2. Let K,M be any integers larger than 1. Then the following inequalities hold for every14Figure 3.1: Comparison between Vb as naively periodized in (3.27) and Ewald summationx ∈ [−L/2, L/2].∣∣∣∣Vb(x)− ∑|m|<M(vb(x−mL)−∫ L/2−L/2dyLvb(x− y −mL))∣∣∣∣ ≤ 12L(12b2L2+ 1)1(M − 1)2 , (3.31)∣∣∣∣Vb ⋆L gb(x)− 2LK∑k=1E1([2πbkL]2)cos(2πkxL)∣∣∣∣ ≤ 12πb[Lb2πK]3exp(−[b2πKL]2), (3.32)∣∣∣∣(Vb(x)− Vb ⋆L gb(x))− ∑|m|<M(vb(x−mL)− 1|x−mL|erf( |x−mL|2b)) ∣∣∣∣≤ 1Lb2L2(2 +4√π)1(M − 1)2 . (3.33)Proof. Section 3.4.3To reach a satisfactory accuracy in the evaluation of the infinite sum, a sufficiently large truncation isnecessary. However, as each summand amounts to an O(1) evaluation, this truncation must not be arbitrarilylarge and should be a judicious choice based on the available computational resource. Theorem 3.3.2 servesas that guideline, informing the error bounds when the stated truncation schemes are chosen.Numerical tests are performed against these analytical error bounds to yield empirical confirmation toTheorem 3.3.2 for several sets of b and L. Figure 3.2 shows the numerical test for Vb as periodized in(3.27). The integral is computed through numerical quadrature, and collectively with other terms in thesummand, treated as O(1)-cost operation. The markered-line are errors directly computed via (3.27) andtruncated accordingly. The dashed-lines are the plots of analytical values 12L(12b2L2+ 1)1(M−1)2 given inTheorem 3.3.2. The exact error requires the evaluation of the limit point∑m∈Z, which is not available. Asa consequence, we resort to plotting the approximated error by∑|m|≤4M −∑|m|≤M . In other words, we15Figure 3.2: Vb as defined in (3.27). We treat each summand as O(1) operation, absorbing also thequadrature. Theoretical bound of error is computed via the formula 12L(12b2L2+ 1)1(M−1)2 inTheorem 3.3.2. Note that the formula predicts that, for a fixed L, 12b2L2+ 1 ∼ 1 for any smallenough b’s. The fact that all three lines on the right look essentially identical verify this prediction.expect 4M -term to approximate the limiting point as M →∞. This is theoretically sound due to Cauchy’scriterion of sequence convergence, and founded because∣∣∣∣ ∑|m|≤4M−∑|m|≤M∣∣∣∣ ≤ ∣∣∣∣ ∑|m|>4M∣∣∣∣+ ∣∣∣∣ ∑|m|>M∣∣∣∣ ≤ CL,b 1716M2 ∼ CL,b 1M2where CL,b is the prefactor as given by the theorem, and we tacitly estimated M − 1 ≈ M . Identicaltechniques are employed for other two infinite sums. Figure 3.3 and 3.4 display numerical tests against,respectively, Vb − Vb ⋆L gb and Vb ⋆L gb. For Vb − Vb ⋆L gb, the error improves as b decreases (with Lfixed), contrasting the case of Vb as displayed in Figure 3.2. For all three plots, the dashed-line representingtheoretical error bound lies above the numerically evaluated error (markered-line) at each data point. Thisverifies that formulae in Theorem 3.3.2 are indeed upper bounds.As there are two infinite sums involved in Ewald summations, some error balancing shall be carried outin order to help us visualize how much gain we actually yield in converting (3.27) into its Ewald sums. This16Figure 3.3: Real space part Vb − Vb ⋆L gb. Each summand treated as O(1)-cost. Theoretical bound oferror is computed via b2(2 + 4√π)1L31(M−1)2 from Theorem 3.3.2will demonstrate how one such aforementioned judicious choice of truncation can be obtained. Indeed, thefollowing technique of error balancing is by no means unique or optimal, but it is robust enough to yielddesirable gain. First, we start by choosing an accuracy tolerance ϵ. Then, let us setϵ =12L(12b2L2+ 1)1M2naive(3.34)for the naive evaluation in (3.27) and (3.31), truncated at Mnaive, and set for Ewald sumsϵ/2 =12πb[Lb2π]3exp(−[b2πKL]2)(3.35)ϵ/2 = b2(2 +4√π)1L31M2, (3.36)where the former represents the reciprocal space error bound (3.29) and (3.33) truncated at K and the latter17Figure 3.4: Reciprocal space part Vb ⋆L gb. Each summand treated as O(1)-cost. Theoretical bound oferror is computed via 12πb[Lb2πK]3exp(− [ b2πKL ]2) from Theorem 3.3.2. We have also plotted12πb[Lb2π]3exp(− [ b2πKL ]2) as the dashed-dotted-line, which demonstrates that the polynomialdependence K−3 has negligible influence on the theoretical error.18represents real space error bound (3.30) and (3.32) truncated at M . The 1/2 factor comes from demandingϵ = errorEwald(K,M) := errorrecip(K) + errorreal(M) (3.37)Certainly, the net error does not have to be evenly allocated to the two summations. Notice also that wehave made the estimations M − 1 ∼ M and exp(−K2)/K3 ∼ exp(−K2), the later of which is justifiedthrough Figure 3.4. These estimations allow us to invert the equations easily:M =[b2(2 +4√π)2L3ϵ]1/2(3.38)K =L2πb[log(L34ϵ(πb)4)]1/2(3.39)Mnaive =[12L(12b2L2+ 1)1ϵ]1/2(3.40)Each truncation solving its equation represents the cost, and assures the corresponding summation are nu-merically evaluated within the error tolerance. By varying ϵ, we can extract the resulting M,K and Mnaive’sand evaluate the empirical errors using the full summation formulae: (3.27), (3.29), and (3.30). The resultis Figure 3.5. For Ewald summation, the net cost to evaluate within ϵ is P ≡ M + K, reflecting the factthat there are two infinite sums to truncate. We again provide both the plot of analytical formula (3.37) and(3.34) - dashed line, and plot of empirical evaluations - markered-line. The agreement is clear. For a givencost, we can see that Ewald summation yields better accuracy. However, this advantage diminishes as b isincreased. If b were like κ, not possessing physical meaning, we could tune it so that we are always workingwithin the top-left of Figure 3.5.This concludes the technical overview for potential evaluation. For every sample x ∈ S within oneoptimization step during a VMC simulation, the potential V (or Vb) is evaluated in calculating local energy(2.10) and its derivative. Typically, |S| = O(103) so an efficient evaluation scheme for PBC is necessary.Lee and Drummond [15, 29] in fact surrogates the evaluation with splines to speed up simulations for severalsets of L and b’s. We do not follow suit for our implementation in Section 5.19Figure 3.5: Convergence of Total Evaluation of Ewald Summation203.4 Proofs3.4.1 (3.13) Satisfies Poisson Equation (3.9)Proof By Direct SubstituitionHere we will demonstrate that GPBC as given by (3.13) is the unique solution to Poisson equation (3.9)with PBC and G˜(0) = 0. The argument relies on formal algebraic manipulation. For a mathematicaldiscussion beyond the scope of this thesis, we refer the readers to Lieb and Simon[14].As it is evident that GPBC is periodic over the box size by definition, it suffices to directly compute theLaplacian and the 0-th Fourier coefficient, ignoring convergence.(−∆GPBC)(r) =∑n∈Z3(−∆r 1|r − nL| −1L3∫Ω(−∆r) dr′|r − r′ − nL|)=∑n∈Z3(δ(r − nL)− 1L3∫dr′1Ω(r′)δ(r − r′ − n))(1Ω = indicator)=∑n∈Z3(δ(r − nL)− 1L31Ω(r − nL))= δ(r)− 1L3in Ω.This shows the desired result for Laplacian. As for Fourier coefficient,G˜PBC(0) =∫Ωdr′GPBC(r′)=∫Ωdr′∑n1|r′ − nL| −1L3∫Ωdr′∫Ωdr∑n1|r − r′ − nL|=∫Ωdr′∑n1|r′ − nL| −1L3∫Ωdr′∫Ω−r′du∑n1|u− nL|=∫Ωdr′∑n1|r′ − nL| −1L3∫Ωdr′∫Ωdu∑n1|u− nL|=∫Ωdr′∑n1|r′ − nL| −1L3(L3)∫Ωdu∑n1|u− nL| = 0.In the third equality, a change of variable u = r − r′ is used. In the 4th equality, the periodicity of theintegrand is invoked. This shows the desired result.21Proof By Fourier TransformWe will make use of the formula (̂1r)(k) =4πk2.Let k ∈ (2π/L)3Z3. Then G˜PBC(k) =∫Ωdre−ikxGPBC(x) =∑k∈(2π/L)3Z3∫Ωdre−ik·r(1|r −mL| −1L3∫Ωdr′1|r − r′ −mL|)=∑k∈(2π/L)3Z3∫Ω−mLdr”e−ik·(r”+mL)(1|r”| −1L3∫Ωdr′1|r”− r′|)=∑k∈(2π/L)3Z3∫Ω−mLdr”e−ik·r”(1|r”| −1L3∫Ωdr′1|r”− r′|)where in the second last equality a change of variable r” = r −mL is applied. This gives e−ik·r =e−ik·(r”+mL) = e−ik·r” · 1 in the last equality. Now continuing, we get=∫R3dr”e−ik·r”(1|r”| −1L3∫Ωdr′1|r”− r′|)=4πk2− 1L3∫Ωdr′∫R3dr”e−ik·r”1|r”− r′|=4πk2− 1L3∫Ωdr′e−ik·r′∫R3dre−ik·r1|r| (change variable r = r”− r′)=4πk2− δk04πk2,and we see the Fourier coefficients satisfy Poisson equation in reciprocal space (3.10). Further elucidationis perhaps required for when k = 0. We remark that this proof is identical to 3.42.3.4.2 Derivation of 1D Ewald SumsThis subsection provides a derivation of Ewald summation of the 1D model, Vb = Vb ⋆L gb + (Vb −Vb ⋆L gb), where the first term is expanded in reciprocal space (3.29), and the second in direct lattice space(3.30). We remark that Lee and Drummond [15] has also provided a derivation, but their approach deviatesfrom the 3D framework (in terms of convolution ⋆) earlier introduced in Section 3.2. We will stay close tothe framework previously introduced so that we can better comprehend the corresponding 3D analog.Other than those from Appendix 2.3, we will put in line several other tools before embarking on thederivations. Firstly, the aforementioned periodic 1D Gaussian is given bygb(x) =∑m∈Z12b√πexp(− 14b2(x−mL)2), (3.41)which is to be compared to the 3D counter part (3.24). This is again viewed as charge density that extends22over the whole space while being L-periodic. Secondly, Casula et al. [16] without justification relates theFourier coefficient V˜b to the Fourier transform vˆb, given in (2.8). We hereby state the relation and providethe proof.Proposition 3.4.1.V˜b(k) =vˆb(k) k ∈ (2π/L)Z\{0}0 k = 0 . (3.42)Proof. Fix k ∈ (2π/L)Z. Then V˜b(k) =∫ L/2−L/2dxe−ikxVb(x) =∑m∈Z∫ L/2−L/2dxe−ikx(vb(x−mL)− 1L∫ L/2−L/2dyvb(x− y −mL))=∑m∈Z∫ L/2−mL−L/2−mLdue−ik(u+mL)(vb(u)− 1L∫ L/2−L/2dyvb(u− y))=∑m∈Z∫ L/2−mL−L/2−mLdue−iku(vb(u)− 1L∫ L/2−L/2dyvb(u− y)),where in the second last equality a change of variable u = x − mL is applied. This gives e−ikx =e−ik(u+mL) = e−iku · 1 in the last equality. Now continuing, we get=∫ ∞−∞due−iku(vb(u)− 1L∫ L/2−L/2dyvb(u− y))= vˆb(k)− 1L∫ L/2−L/2dy∫ ∞−∞due−ikuvb(u− y)= vˆb(k)− 1L∫ L/2−L/2dye−iky∫ ∞−∞dxe−ikxvb(x) (change variable x = u− y)= vˆb(k)− δk0 vˆb(k).Note that vˆb(0) =∞ so the last line needs further interpretation just like the 3D case (3.4.1).We are in position to derive Ewald sum. We shall begin with the reciprocal space sum Vb ⋆ gb. Let us23write hb(x) =exp(− x24b2)2b√πfor convenience. ThenVb ⋆L gb(x) =∫ L/2−L/2dyVb(y)gb(x− y)dy=∑m∈Z∫ L/2−L/2dyVb(y)hb(x−mL− y)=∑m∈Z(Vb ⋆L hb)(x−mL)=1L∑k∈(2π/L)Ẑ(Vb ⋆L hb)(k)eikx (Poisson summation (2.22))=1L∑k∈(2π/L)ZV˜b(k)hˆb(k)eikx (by (2.23))=1L∑k∈(2π/L)Z\{0}vˆb(k)hˆb(k)eikx (by (3.42))=2L∞∑k=1E1([b2πkL])cos(2πkxL),where in the last equality, we have used the fact that both vˆb and hˆb are even. This completes the proof forreciprocal space sum.24For the real space sum, we will formally perform arithmetic on the infinity vˆb(0).(Vb ⋆L gb)(x) =∑m∈Z∫ L/2−L/2dyVb(y)hb(x−mL− y)=∑m∈Z∫ L/2+mL−L/2+mLduVb(u−mL)hb(x− u) (change variable u = y +mL)=∑m∈Z∫ L/2+mL−L/2+mLduVb(u)hb(x− u) (L-periodicity of Vb)=∫ ∞−∞duVb(u)hb(x− u)=∑m∈Z∫ ∞−∞duvb(u−mL)hb(x− u)−∑m∈Z∫ ∞−∞du1L∫ L/2−L/2dyvb(u−mL− y)hb(x− u)=∑m∈Z1|x−mL|erf( |x−mL|2b)(integrate directly, as did in [15])−∑m∈Z∫ ∞−∞du1L∫ L/2+mL−L/2+mLdwvb(u− w)hb(x− u) (change variable w = y +mL)=∑m∈Z1|x−mL|erf( |x−mL|2b)−∫ ∞−∞du(1L∫ ∞−∞dwvb(u− w))hb(x− u)=∑m∈Z1|x−mL|erf( |x−mL|2b)− 1Lvˆb(0)∫ ∞−∞duhb(x− u)=∑m∈Z1|x−mL|erf( |x−mL|2b)− 1Lvˆb(0) (Gaussians hb are normalized)25Thus,Vb(x)− Vb ⋆L gb(x) =∑m∈Z(vb(x−mL)− 1L∫ L/2−L/2dyvb(x− y −mL)− 1|x−mL|erf( |x−mL|2b))+1Lvˆb(0)=∑m∈Z(vb(x−mL)− 1|x−mL|erf( |x−mL|2b))− 1L∫ ∞−∞dyvb(x− y) + 1Lvˆb(0)and the last two terms cancel, hence the desired result.3.4.3 Proof of Theorem 3.3.2Let x be in the simulation box Ω = [−L/2, L/2].Convergence Analysis of Naive Summation (3.27)Using the definition as given by (2.14), we obtain through integration by parts thaterfc(x) =2√π∫ ∞x(−1)(2t)−1 ddt(e−t2)dt =1√πe−x2x− 1√π∫ ∞xt−2e−t2dt (3.43)=⇒ vb(x) = 1|x| −12bexp(x24b)∫ ∞|x|/2bt−2e−t2dt (3.44)As we are truncating the series, it suffices to consider just the tail (i.e terms labelled beyond truncationM > 1) of it, and show that the summands there are small enough. In this spirit, assume that m > M > 0.The first goal is to witness the cancellation of 1/|x| term in the m-th summand. This is crucial becauseabsolute convergence is what we are after. To begin, recall from (3.27) that the m-th summand readsvb(x±mL)− 1L∫ L/2−L/2vb(x− y ±mL)dy.Both ±m terms here are seen as the “m-th summand.” With accord to earlier result via integration by parts,define R byvb(x) =1|x| +R(x). (3.45)26We can then directly expand and computevb(x±mL)− 1L∫ L/2−L/2vb(x− y ±mL) = 1mL± x +R(mL± x)−1L∫ L/2∓mL−L/2∓mLdyR(x− y)− 1L[LmL± x +R+(L2(mL± x))−R−(L2(mL± x))]= R(mL± x)− 1L∫ L/2∓mL−L/2∓mLdyR(x− y)− 1L[R+(L2(mL± x))−R−(L2(mL± x))],(3.46)where R+(x) and R−(x) are (Lagrange type) Taylor remainders of respectively ln(1 + x) and ln(1− x)centred at 0 as given byR+(x) = − 12(1 + ξx)2x2; R−(x) = − 12(1− ηx)2x2This gives the desired cancellation.We now proceed to derive the bound on the RHS of (3.31). To do so, we collect the following results.Lemma 1.|R(x)| ≤ 2b2 1|x|3 .Lemma 2. ∣∣∣∣R+( L2(mL± x))−R−( L2(mL± x))∣∣∣∣ ≤ 4m− 52(2m− 3)2(L2(mL± x))2.Lemma 3. ∣∣∣∣ ∫ L/2∓mL−L/2∓mL dyR(x− y)∣∣∣∣ ≲ 4b2 L(mL± x)3 .Let us first assume these results. Since x ∈ Ω, we can then use triangle inequality to derive∣∣∣∣vb(x±mL)− 1L∫ L/2−L/2vb(x− y ±mL)∣∣∣∣ ≤ ∣∣∣∣R(mL± x)∣∣∣∣+ 1L∣∣∣∣ ∫ L/2∓mL−L/2∓mL dyR(x− y)∣∣∣∣+1L∣∣∣∣R+( L2(mL± x))−R−( L2(mL± x))∣∣∣∣≤ 6b2 1(mL± x)3 +(2m− 1)L4(mL± x)2≤ 6b2L31(m− 1)3 +1L12(m− 1)3 .In the second inequality we use m − 1 < m − 1 + (m − 2) = 2m − 3 =⇒ 4m−52(2m−3)2 ≤ 4m−52(m−1)2 ≤4m−42(m−1)2 =2(m−1) . In the last inequality we use the assumption that x is in the first cell.This shows that the m-th summand is ∼ O(1/m3), so the effective potential V as originally stipulated27is indeed absolutely convergent. To estimate the rate, we will simply integrate over the tail∣∣∣∣ ∑m∈Z(vb(x−mL)− 1L∫ L/2−L/2vb(x− y −mL))−M∑m=−M(vb(x−mL)− 1L∫ L/2−L/2vb(x− y −mL))∣∣∣∣≤∑m>M∣∣∣∣vb(x−mL)− 1L∫ L/2−L/2vb(x− y −mL)∣∣∣∣+ ∑−M>−m∣∣∣∣vb(x+mL)− 1L∫ L/2−L/2vb(x− y +mL)∣∣∣∣≤ 2∑m>M(6b2L31(m− 1)3 +1L12(m− 1)3)=1L(12b2L2+ 1) ∑m>M(1(m− 1)3)≤ 1L(12b2L2+ 1)∫ ∞Mdm(m− 1)3 =12L(12b2L2+ 1)1(M − 1)2 .This demonstrates the convergence (3.31) for the naively periodized Vb stated in Theorem 3.3.2.We last but not least present the proofs to previous lemmas.Proof of lemma 1.|R(x)| = 12bexp(x24b2)∫ ∞|x|/2bt−2e−t2dt≤ 12bexp(x24b2)( |x|2b)−3 ∫ ∞|x|/2bte−t2dt=12bexp(x24b2)( |x|2b)−3 12exp(− x24b2).Proof of lemma 2. The proof for both ± will be identical, so let us assume that we are dealing with theremainders evaluated at δ± = L/2(mL± x). For the Lagrange remainder, we have the relation0 < ξδ± , ηδ± ≤12(m− 1)where the last inequality follows from the assumption that x is in the first cell [−L/2, L/2] so that δ± ≤L/[2(m− 1)L]. Hence,1− ξδ± ≥ 1− δ± > 1−12(m− 1) =2m− 32m− 2 .28Thus,R+(δ±) =12(1 + ξδ±)2δ2± ≤12δ2±R−(δ±) =12(1− ηδ±)2δ2± ≤12(2m− 22m− 3)2δ2±=⇒ |R+(δ±)−R−(δ±)| ≤ (2m− 2)2 − (2m− 3)22(2m− 3)2 δ2± ≤4m− 52(2m− 3)2 δ2±.Proof of lemma 3. Let δ∓ be as defined in the proof to Lemma 2. By lemma 1, we have∣∣∣∣ ∫ L/2∓mL−L/2∓mL dyR(x− y)∣∣∣∣ ≤ 2b2 ∫ L/2∓mL−L/2∓mL dy 1|y − x|3= 2b2[1(−L/2 +mL∓ x)2 −1(L/2 +mL∓ x)2]= b2[1(mL∓ x)2] [1(1− δ∓)2 −1(1 + δ∓)2]≤ b2[1(mL∓ x)2] [4δ∓(1− δ2∓)2]∼ b2[1(mL∓ x)2]4δ∓.Convergence Analysis of Vb ⋆ gb (3.33)Via the definition of exponential integral (2.15), we have for y > 0 that∣∣E1(y)∣∣ ≤ 1y∫ ∞ye−t ≤ 1ye−yThus the tail of the sum can be bounded easily by a similar procedure.∣∣∣∣ 2L ∑k>KE1([b2πkL]2)cos(2πkxL)∣∣∣∣ ≤ 2L ∑k>KE1([b2πkL]2)≤ 12πb[Lb2πK]3exp(−[b2πKL]2),where the last inequality follows from |E1(y2)| ≤ e−y2y3. This demonstrates that reciprocal space sumconverges at the rate as given in Theorem 3.3.2.29Convergence Analysis of Vb − Vb ⋆ gb (3.32)We immediately have by 2.14 the dersired 1/|x| cancellationvb(x)− 1|x|erf( |x|2b)=1|x| +R(x)−(1|x| −1|x|erfc( |x|2b)).Passing to absolute values, we derive∣∣∣∣vb(x)− 1|x|erf( |x|2b) ∣∣∣∣ ≤ |R(x)|+ ∣∣∣∣ 1|x|erfc( |x|2b) ∣∣∣∣≤ 2b2 1|x|3 +1|x|2b√π|x| exp(−x24b2)≤ 2b2 1|x|3 +(2b)2√π|x|3 = b2(2 +4√π)1|x|3 .The last inequality invokes exp(−x2/4b) < 2b/|x|. With finally an integration on the tail of the real spacepart we see ∣∣∣∣ ∑|m|>M(vb(x−mL)− 1|x−mL|erf( |x−mL|2b)) ∣∣∣∣≤ 2b2(2 +4√π)∫ ∞Mdm(mL− x)3 = b2(2 +4√π)1L1(ML− x)2≤ b2(2 +4√π)1L31(M − 1)2 ,where the last inequality relies on the fact that x ∈ [−L/2, L/2]. This demonstrates the convergence forVb − Vb ⋆ gb as stated in Theorem 3.3.2.30Chapter 4Wavefunction Methods4.1 Conventional Quantum ChemistryIn this section we review the wavefunction ansatz more commonly used in the field of quantum chem-istry. Traditional ansatz pivots on 1) the single electron basis set, and 2) the truncation scheme for determi-nant expansion. Given a basis set1 {χi}Ki=1, the full determinantal expansion of these basis functions, alsoconventionally known as the Full Configuration Interaction (FCI), can be seen as an ansatz of the minimiza-tion problem (2.9),ΨFCI =∑ICIΨI , (4.1)where CI are the variational parameters, and ΨI is a Slater determinant. For I = (i1, . . . , iN ), we haveΨI =1√N !∣∣∣∣∣∣∣∣∣∣χi1(x1) χi2(x1) . . . χiN (x1)χi1(x2) χi2(x2) . . . χiN (x2)....... . ....χi1(xN ) χi2(xN ) . . . χiN (xN )∣∣∣∣∣∣∣∣∣∣, (4.2)and the summation in (4.1) runs over all ordered combinations of the indices i. Evaluation of an FCIwavefunction (4.1) quickly gets prohibitive as the number of determinants increase, yet there are situationswhere high accuracy is not of central importance. As such, hierarchical structures are designed for thedeterminantal expansion ansatz so that we can, but systematically, further truncate the summation. Herewe will give two well-known examples, the single-reference Configuration Interaction (CI) and CoupledCluster (CC) methods.1by basis set, we mean a finite subset of a complete basis, which does not have to be finite.31Single Reference Configuration Interaction (CI)Viewing at face value the expression in (4.1), the most straight forward hierarchy would be to decreewhat determinants to be included in the linear parameterization. CI is such hierarchical structure, definedwith respect to a reference. Let us assume the basis {χi} is endowed with some order, in such a way thatthe first N of them χ1, . . . , χN constitute the reference determinant ΨI0 . Then the hierarchical strcuture inCI is exhibited byΨCI =1 + ∑k>NN∑i=1Cki Eki +∑k2>k1i2>i1Ck1k2i1i2 Ek1k2i1i2+∑k3>k2>k1i3>i2>i1Ck1k2k3i1i2i3 Ek1k2k3i1i2i3+ . . .ΨI0 , (4.3)where the C’s are variational parameters, and Ek1k2...i1i2... is a linear operator defined by replacing χi1 in thedeterminant it acts on by χk1 , χi2 by χk2 , and so on (but return 0 if the χ to be replaced is not contained inthe determinant)2. Formula (4.3) is extracted directly from [1] and is termed intermediate normalized formin [4], but sometimes the reference determinant is also multiplied by a parameter, say C0, instead of 1. Thelatter convention is also discussed in the conventional textbooks [1, 4]. The collection of all Eki |I0⟩’s arecalled singles (with respect to I0), Ek1k2i1i2 |I0⟩ called doubles, etc. Thus, one can declare via the acronym,CISD, CISDT, CISDTQ, etc, to specify the truncation at the double, triple, quadruple, etc level. Albeitstraight forward, CI fails to possess size-extensivity [1] that the exact ground state wavefunction satisfies.Coupled Cluster (CC)Motivated by the lack of size-extensivity of CI, we have the CC ansatzΨCC =[1 +∑k>NN∑i=1Cki Eki]1 + ∑k2>k1i2>i1Ck1k2i1i2 Ek1k2i1i21 + ∑k3>k2>k1i3>i2>i1Ck1k2k3i1i2i3 Ek1k2k3i1i2i3 . . .ΨI0 (4.4)where the notations are analogous to the case of CI. Notice the bracketed terms commute, so there is noambiguity in the product definition.As opposed to the linear truncation in CI, the hierarchical structure of CC manifests in the length ofproduct. Thus, one can declare via the acronym, CCSD, CCSDT, CCSDTQ, etc, to specify the producttruncation at the double, triple, quadruple, etc level. However, schemes above CCSDT get impracticallycostly. Note also that the parameterization is nonlinear. By expanding the product, nevertheless, one canrewrite the expression as summation of determinants, it is in this sense that one retrieves FCI.2This operator is more formally defined via second quantization by Ek1k2...i1i2... := a†k1a†k2 . . . ai1ai2 . . . . However, recalling thatlinear maps are completely specfied by their actions upon basis elements - determinants in this case, the loose definition in the maintext suffice.324.1.1 Hartree-Fock TheoryThroughout the preceding discussion, we have been obscure about the single determinant referenceΨI0 upon which CI and CC rely. In principle, ΨI0 can be any determinant constructed through reasonablechemical or physical intuition.3 There is, however, a canonical choice in Hartree-Fock (HF) model [4, 30],the cornerstone of wavefunction based electronic structure methods. As FCI, HF is also a wavefunctionansatz, but predicates on a single Slater determinantΨHF, given by (4.2). It is used a lot but is known to oftenhave insufficient accuracy. Nevertheless, not only does the single determinant HF energy ⟨ΨHF|H |ΨHF⟩impose a reasonable variational upper bound for benchmarking, it provides a natural way to construct thepreviously unspecified basis {χi} that also exhibits qualitatively correct behavior. By assuming the ansatzΨ = ΨHF, we can then analytically minimize the energy functional with 1-electron functions subjected toan orthonormality constraint, ∫χi(x)χj(x)dx = δij .Such calculation, detailed in [4], will lead to the Hartree-Fock equations. Assuming the notations in Szabo[4], where the shorthands f(1) = f(r1), r12 = r1 − r2 are introduced, the HF equations are given byF(1)χi(1) = ϵiχi(1), (4.5)F(1) := h(1) +N∑j=1j ̸=iJj(1)−N∑j=1b̸=aKj(1), (4.6)h(1) = −12∆1 −M∑A=1ven(r1A), (4.7)Jj(1) =∫dx2χj(2)vee(r12)χj(2), (4.8)Kj(1)χ(1) =∫dx2χj(2)vee(r12)χ(2)χj(1), (4.9)where we have been deliberately ambiguous about the precise form of the two-particle interaction vee andven, as the referenced derivation based on variational calculus does not necessitate the use of Coulombinteraction. That is, the usual case is given byvee(rij) =1|rij | =1|ri − rj | and ven(riA) =ZA|riA| =ZA|ri −RA| .The solutions to the HF equations, i.e orbitals, constitute a basis {χi} ordered by their eigenvalues. Withinthe HF single determinant regime, the N orbitals giving the lowest ⟨ΨHF|H|ΨHF⟩ form the best approx-imate wavefunction, usually the ones corresponding to the N lowest eigenvalues. Due to the simplicityafforded by the single determinant ansatz, for any N orbitals χa (not necessarily corresponding to the N3There exists frameworks to include multiple reference [1], though we will not pursue such scheme in this thesis.33lowest eigenvalues), the HF energy can be computed by⟨ΨHF|H |ΨHF⟩ =N∑a[χa|h|χa] + 12N∑aN∑b[χaχa|χbχb]− [χaχb|χbχa], (4.10)[χa|h|χb] =∫dx1χ∗a(1)h(1)χa(1), (4.11)[χaχb|χcχd] =∫dx1dx2χ∗a(1)χb(1)vee(r12)χ∗c(2)χd(2). (4.12)For almost all practical scenarios, the HF equations themselves have to be solved numerically through theSelf-Consistent Field method [4]. This turns out to be unnecessary for UEG, and we conclude this sectionwith the following.Theorem 4.1.1. The HF equations of the UEG model, as introduced in (3.22), can be diagonalized byplane-waves. More precisely, leth(1) = −12∆1,vee(r) = GPBC(r),index i = (k, σ) ∈ (2π/L)3Z3 × {↑, ↓},χi(x) = χk,σ(r, s) = eik·rδσs /L3/2,∫dx =∑σ=↑,↓∫Ωd3r.ThenF(1)χk,σ(1) = ϵ(k)χk,σ(1),ϵ(k) =12k2 − 1L3∑q ̸=kG˜PBC(k − q).Proof. Section 4.5.14.1.2 Cusp Condition and The Helium AtomIn practice, the orbital usually takes on the form χiσ(r, s) = fi(r)δsσ with its spatial part fi in turnbeing specified as linear combination of atomic orbitals. For typical molecular calculations, each atomicorbital is centered on some nucleus. This implies the ansatz for conventional quantum chemistry does notdepend explicitly on electron-electron relative displacement rij = ri − rj . This is disadvantageous since itis known when two electrons get near, cusp behavior manifests in the exact wavefuinction. In fact, we havethe following precise mathematical statement due to Kato [31].Kato Cusp Condition. Let rij = ri−rj and riI = ri−RI , where ri and rj are electron coordinates, andRI nucleus coordinate. Then the solution Ψ to Schrodinger equation is everywhere continuous, possesses34Figure 4.1: The error in the electronic energy (in Hartrees) of the ground-state helium atom. The erroris plotted as a function of the number of terms in the expansions. This image is adapted fromHelgaker et al. [1] Figure 7.6.uniformly bounded first derivatives, and satisfies the cusp conditions,limrij↓0∂Ψ¯∂rij=12Ψ∣∣∣∣rij=0, andlimriI↓0∂Ψ¯∂riI= −ZΨ∣∣∣∣riI=0, (4.13)provided that there is no third particle or more coalescing. Ψ¯ denotes the spherical average (see (2.24)) ofΨ over the displacement (rij or riI ) coordinate system. Z is the atomic number of nucleus I .The partial derivative in Kato’s theorem is taken along the positive radial component before sending theradial distance between two particles ↓ 0. Intuitively, provided that the wavefunction is non-zero uponcoalescence of a pair of particles, Kato’s condition entails that the exact wavefunction then exhibits a non-zero radial slope, cusp-like behavior.4 The deficiency in using traditional quantum chemistry methods toconverge to electron-electron cusp can be witnessed through a practical example - the Helium atom. By4For a pair of electrons of parallel spin, the RHS of the condition vanishes, implying that we do not have first order cusp. Packand Brown [32, 33] showed that we may still have the next order condition,limrij↓0(1Ψ¯∂Ψ¯∂rij)=14. (4.14)35choosing Slater Type Orbital (STO) as basis, Helgaker et al. [1, Ch7.3] showed that the spatial part of theFCI ansatz for approximating 1S-ground state of Helium is equivalently given by the parameterization5ΨFCI = exp(ζ(r1 + r2))∑ijk∈Z≥0Cijk(ri1rj2 + rj1ri2)r2k12 , (4.15)where r1 and r2 are distances of electron 1 and 2 from the single nucleus, and r12 is the distance between thetwo electrons. Note that r12 is made explicit here only after some argument involving change of coordinates.The STOs, prior to and change of variables, depend only on electron-nucleus displacement. The poorconvergence of this parameterization is illustrated Figure 4.1, specifically by single-ζ as well as numericalCI. The central problem lies in the absence of the linear r12 term. No matter how the truncation is carriedout, differentiation with respect to r12 always yields zero. That is,∂ΨFCI∂r12∣∣∣∣r12=0= 0,hence clearly a violation to the cusp condition. Remarkably, the following simple modified ansatz withlinear r12 term included,ΨFCI−r12 =(1 +12r12)ΨFCI, (4.16)leads to the correct cusp behavior6∂ΨFCI−r12∂r12∣∣∣∣r12=0=12ΨFCI(r12 = 0) =12ΨFCI−r12(r12 = 0), (4.17)as well as the improved convergence as demonstrated by single-ζ CI-R12 Figure 4.1.The best performing scheme, labelled as Hylleraas on Figure 4.1, is simply the generalization of (4.15)to include all powers of r12,ΨHylleraas = exp(ζ(r1 + r2))∑ijk∈Z≥0Cijk(ri1rj2 + rj1ri2)rk12. (4.18)It demonstrates that one can indeed improve accuracy and variationally learn the cusp condition by employ-ing parameterizations beyond the determinantal expansion. However, such increase in parameters inevitablyalso increases the computational cost. In contrast, ΨFCI−r12 has just as many parameters as ΨFCI, yetachieves chemical accuracy (1mEh) with significantly fewer terms than FCI. As such, to produce the opti-mal accuracy-cost ratio, a prudent ansatz that directly incorporates cusp conditions while being more flexiblethan determinantal expansion would indubitably pave a better route.5In order to satisfy the spin and spatial angular momentum symmetries decreed by the atomic term symbol 1S, namely singletfor both, the FCI here is expanded in precomputed linear combinations of Slater determinants. Under this construction, the ansatz1) turns out to be rotationally invariant, and 2) factors into a symmetric spatial part and an antisymmetric spin part. The later ofwhich serves as a constant, hence ignored in (4.15).6Note that 1S demands antiparallel spins, hence the (Kato) 1/2 prefactor.364.2 Explicitly Correlated Wavefunctions - A Brief SurveyThe discussion from the last section foreshadows the notion of explicitly correlated methods, meaningan ansatz that explicitly depends on rij . From Fermi Hypernetted Chain [34], to Coupled Cluster R12 [35],to Variational Monte Carlo that this thesis is concerned with, the methodology has indeed already had somehistory. Common to all these mentioned is the starting point, usually known as the Slater-Jastrow ansatzΨSJ = eJΨ0, (4.19)where Ψ0 is the antisymmetric (Slater) part and J is the symmetric (Jastrow) part. The former is noneother than the Slater determinant (4.2). On the other hand, the Jastrow part depends explicitly on rij and ismultiplied to the antisymmetric part so that the ansatz satisfies Kato’s condition. This form of ansatz wasfirst proposed in [36]. An immediate generalization of this ansatz is a linear combination of terms eachassuming the same form as ΨSJ. However, as we saw in the Helium example that the main purpose is abetterment with as few determinants involved in a summation as possible, we will restrict our attention to asingle ΨSJ.Of all the mentioned existing methods, VMC has the most flexibility in the design one can employbecause the optimization procedure associated to it is virtually blind (at least in theory) to the specific formof the ansatz (See equation (2.11)). Expectation values computed within CC-R12, for instance, pivot onthe evaluations of two-centered Gaussian integrals, highlighting the method’s reliance on specific basis. Toexploit the flexibility of VMC further, the antisymmetric Slater part of the ansatz had been further gener-alized. Conventional determinantal wavefunctions rely on the parameterization of single particle orbitalsχ(i) = χ(xi). VMC, however, can accomodate generalized orbitals of the formχ(i) = χ(xi, x̸=i). (4.20)Here, χ is understood as a function of all N particles, but with the additional property that it be symmetricwith respect to all but particle-i’s coordinate. Indeed, one can show the expressionΨBF(x1, . . . ,xN ) =∣∣∣∣∣∣∣∣∣∣χ1(x1,x ̸=1) χ2(x1,x ̸=1) . . . χN (x1, x̸=1)χ1(x2,x ̸=2) χ2(x2, x̸=2) . . . χN (x2, x̸=2)....... . ....χ1(xN ,x ̸=N ) χ2(xN , x̸=N ) . . . χN (xN , x̸=N )∣∣∣∣∣∣∣∣∣∣(4.21)remains antisymmetric. This next order generalization on the orbital is now known as the Backflow Trans-formation, whoes origin is usually credited to Feynman and Cohen [37]. Nowadays, it has become commonfor VMC software [8][38] to include Slater-Jastrow, Backflow, and Slater-Jastrow-Backflow as choices dueto their effectiveness.With the development of Stochastic Gradient Descent and NN, the utmost flexibility of VMC by far iswitnessed just recently [9, 10]. With no Jastrow factor, FermiNet [9] uses NN to parameterize the symmetricportion of backflow transformed orbitals χj and includes several determinants. On the other hand, also with37multiple determinants, PauliNet [10] parameterizes both J and χj . However, the former model gives riseto ∼ 106 parameters and the later ∼ 105 parameters., both of which drastically exceed the typical amountrequired by FCI, which is ≲ 103 parameters, making physical insight less apparent in the model. Whatremains nevertheless notable is thatJPauliNetθ = Kato + JNNθ , andχPauliNetθ (i) = χHF(xi)χNNθ (x1, . . . ,xN ), (4.22)where θ denotes explicit dependence upon optimizable parameters. In particular, PauliNet uses explicit cuspfunction “Kato” and Hartree-Fock orbitals χHF, both free of optimizable parameters. This seems to suggestthat the reason the model in PauliNet results in an order of magnitude less parameters than FermiNet isbecause more conventional physical wisdom is incorporated.4.3 Explicitly Correlated Wavefunctions - A Closer LookWe have mentioned that the main purpose of Jastrow J factor in (4.19) is so that the ansatz satisfiesthe Kato cusp condition (4.13). However, we have never explained how this is achieved. We will start fillingthe missing gaps following calculations to that presented in Drummond et al. [33]. We first present the mainresult.Cusp Condition 1 on Jastrow. Suppose J is everywhere continuous, C1 away from coalescences, and itsfirst derivative is uniformly bounded and possesses a limit upon coalescence. On the other hand, suppose theSlater partΨ0 satisfies all mentioned regularities J satisfies, but is additionallyC1 around the coalescences.Then the conditions (∂J¯∂rij)rij↓0=12(∂J¯∂riI)riI↓0= −Z (4.23)guarantee the desired cusp condition (4.13), provided there is no third particle or more coalescing.It is now evident what it means to say that “the Jastrow factor is used to satisfy the cusp condition”:we simply equate the limit of its derivative to the proportion factor from Kato’s condition. In the above, theregularity assumptions imposed on J are simply those taken from what Kato (4.13) deduced for the eigen-function of the Schrodinger equation. The additional smoothness allowed for Ψ0 really is for convenience,on which we will comment more. We can nevertheless derive the condition without the additional regularity.Cusp Condition 2 on Jastrow. Suppose J and Ψ0 each is everywhere continuous, C1 away from coales-38cences, and has first derivative that is uniformly bounded and possesses a limit upon coalescence. ThenΨ(rij = 0) · 12= Ψ(rij = 0)(∂J¯∂rij)rij↓0+ eJ(rij=0)(∂Ψ¯0∂rij)rij↓0Ψ(riA = 0) · (−Z) = Ψ(riI = 0)(∂J¯∂riI)riA↓0+ eJ(riI=0)(∂Ψ¯0∂riI)riI↓0(4.24)provided there is no third particle or more coalescing.We shall first consider Cusp Condition 2 (4.24), as Cusp Condition 1 (4.23) will follow as a corollary.Concretely, let us consider electron 1 and electron 2 which define the coordinatesr12 = r1 − r2,r¯12 =r1 + r22, (4.25)where r¯12 denotes the center of mass. The first property to notice is that on the LHS of (4.13), the orderof limr↓0 ∂/∂r and spherical averaging can be swapped (proven in Section 4.5). Thus, by a change ofcoordinates,(r1, r2, . . . , rN ) = (12(2r¯12 + r12),12(2r¯12 − r12), r3, . . . , rN ) (4.26)= (12(2r¯12 + r12rˆ12),12(2r¯12 − r12rˆ12), r3, . . . , rN ), (4.27)we may write Ψ = Ψ(r12, . . .) and deduce from Kato’s condition that−∫S2drˆ12 limr12↓0∂Ψ∂r12=12Ψ(r12 = 0). (4.28)Inserting Ψ = eJΨ0, the product rule gives−∫S2drˆ12 limr12↓0∂Ψ∂r12= −∫S2drˆ12(Ψ∂J∂r12+ eJ∂Ψ0∂r12)r12↓0= Ψ(r12 = 0)−∫S2drˆ12(∂J∂r12)r12↓0+ eJ(r12=0)−∫S2drˆ12(∂Ψ0∂r12)r12↓0,(4.29)that is,Ψ(r12 = 0)12= Ψ(r12 = 0)(∂J¯∂r12)r12↓0+ eJ(r12=0)(∂Ψ¯0∂r12)r12↓0, (4.30)which is precisely Cusp Condition 2 (4.24). Simplification towards Condition 1 (4.23) can be obtained if fur-ther regularity is allowed on Ψ0. Let us now assume that Ψ0 is C1 around r12 = 0, then limr12↓0 ∂Ψ¯0/∂r1239vanishes by (4.57) and we obtain the conditionΨ(r12 = 0)12= Ψ(r12 = 0)(∂J¯∂r12)r12↓0, (4.31)which is satisfied by simply equating as done in (4.23). The above calculation can be adapted to other two-particle coalescences with some modifications7. This completes the justification of both Cusp Conditions(4.23) and (4.24).To satisfy the original cusp condition (4.13), one is allowed to use Condition 1 (4.23) when designinga Jastrow factor, but that also means Ψ0 is to be made C1 around coalescences. Upon extensive literaturereviews, however, it seems to be a common practice to take (4.23) as the Kato cusp condition withoutverifying whether Ψ0 is C1 or not. For the most basic form where Ψ0 is the usual Slater determinantconstructed with orbitals χ(i) = χ(xi), analysis can plausibly be done on the conventional basis (e.g STOs)that expand the orbitals. For when backflow transformation is adopted as done in traditional VMC software(e.g CASINO [8]) or NN-based architectures [9, 10], χ(i) = χ(xi,x ̸=i) could still be analyzed based on thefunctions of which it is composed, albeit the inevitable increase in complexity. Regardless, such empiricismhas already had some history and does yield successes, so we shall exercise the same negligence whenCondition 1 (4.23) is invoked.4.4 ACE and Its Adaptation to Correlated Wavefunctions4.4.1 ACE Parameterization of Jastrow FactorThe Jastrow factor is believed to admit a many-body like expansion [20, 39],J = Jee + Jen + Jeen + . . . , (4.32)and the three terms as exhibited are empirically justified as the three most important contributions [33,40, 41]. Namely, electron-electron (ee), electron-nucleus (en), electron-electron-nucleus (een). The taskis now to devise an approximation scheme, substituting each body-ordered term with functions containingoptimizable parameters. A possible form motivated by the CASINO code [8, 33] is given viaJee =N∑i ̸=ju(rij)Jen =N∑iM∑IηI(riI)Jeen =N∑i ̸=jM∑IfI(riI , rjI), (4.33)7For nucleus, the mass is ∞. For electrons with anti-parallel spins, we refer the readers to Drummond’s work [33]40where u, η, f are free of electron labels to favor the fact that electrons are indistinguishable. We could makethe parameterization spin dependent by augmenting, for instance, the ee term to uσiσj but we will not addressit in this thesis for the sake of clarity. Note that we have chosen RI as the center of the coordinatization forthe een term. A straight forward expansion yields the parameterizationu(r) =∑νCνϕν(r)ηI(r) =∑νCIνϕν(r)fI(r1, r2) =∑ν1ν2CIν1ν2ϕν1(r1)ϕν2(r2), (4.34)with the C’s be optimizable parameters and {ϕν} some basis to be chosen.The cost associated to this evaluation scheme can be estimated. Let us denote the truncation implicitin the expansion (4.34) over ν by K. Then for the two-body terms Jee and Jen, we have O(N2K) andO(NMK) evaluations respectively (treating evaluation of basis ϕ(r) as O(1)). For Jeen, O(N2MK2). ForJeee (if included), written in the same form as Jeen but with electrons as centers, we have O(N3K2), and soon. The key trend to note is that for every increment of 1 body order, the exponent increases/worsens by 1.Certainly, it would be desirable to ameliorate such scaling as there are situations where contributions fromhigher order Jastrow terms are significant [40].Recently, Drautz introduced an efficient scheme, dubbed Atomic Cluster Expansion (ACE) [6], toparameterize interactomic potentials. The pivotal “trick” Drautz invoked is applicable to more generalcontexts, including the present parameterization scheme for Jastrow factors. The idea is to add and subtractself-interaction for terms with a common center. In the case for terms with nuclear centers,Jen + Jeen =N∑iM∑I∑νCIνϕν(riI)−N∑iM∑I∑ν1ν2CIν1ν2ϕν1(riI)ϕν2(riI)︸ ︷︷ ︸self interaction+N∑iM∑I∑ν1ν2CIν1ν2ϕν1(riI)ϕν2(riI)︸ ︷︷ ︸self interaction+N∑i ̸=jM∑I∑ν1ν2CIν1ν2ϕν1(riI)ϕν2(rjI)=N∑iM∑I∑νC¯Iνϕν(riI) +N∑ijM∑I∑ν1ν2CIν1ν2ϕν1(riI)ϕν2(rjI),where the coefficient C¯Iν is obtained by realizing that the self interaction term (with a minus sign) is effec-tively a two body term itself, hence can be linearly spanned by ϕν(riI). Notice that the sum of the higherbody terms is now an unrestricted sum in electron indices, which allows us to writeJeen =M∑I∑ν1ν2CIν1ν2(N∑iϕν1(riI))· N∑jϕν2(rjI) . (4.35)41The cost of evaluation of this new three-body Jeen is O(NMK) + O(MK2), which is now linear in thenumber of electrons. By including higher body-order terms and applying the same trick, the ACE parame-terization can be easily generalized.J(single n) = Jen + Jeen + Jeeen + . . .=M∑I(∑νCIνAIν +∑ν1ν2CIν1ν2AIν1AIν2 +∑ν1ν2ν3CIν1ν2ν3AIν1AIν2AIν3 . . .), (4.36)whereAIν =N∑iϕν(riI). (4.37)With slight modification, we can also derive for all-electron Jastrow terms,J(all e) = Jee + Jeee + Jeeee + . . .=N∑i(∑νCνAiν +∑ν1ν2Cν1ν2Aiν1Aiν2 +∑ν1ν2ν3Cν1ν2ν3Aiν1Aiν2Aiν3 . . .), (4.38)whereAiν =N∑j ̸=iϕν(rij). (4.39)Since the parameterization is linear, we can write the Jastrow factor as dot product of two (perhaps long)vectors,J = J(all e) + J(single n) = C ·A, (4.40)where C is constituted by all the parameters, and A by all the features∏ν Aν .Absent in the foregoing discussion is a Jastrow factor with interactions involving more than 1 nucleus,and one may wonder whether they need to be included. Recall that the electronic Schrodinger equationunder Born-Oppenheimer approximation has clamped nuclei. Therefore, close range interaction betweenmore than 1 nucleus is impossible, provided we have suitably chosen the length scale. Hence, we arejustified to neglect them.For a rigorous analysis on the completeness and complexity of ACE parameterization of general many-body like expansion, we refer the readers to [7, 11].4.4.2 Rotationally Invariant Jastrow FactorOne notable trait of the Jastrow factor implemented in the CASINO code [8, 33] is its O(3)-invariance,i.e. rotation and reflection invariance. As a particle (e.g. electron or nucleus) is usually thought of as42a spherically symmetric entity, such trait is natural. The Jastrow factor as defined by (4.36), (4.38), and(4.40) is permutation invariant, but one can also build O(3)-invariance into ACE parameterisations, and theprocedure of which has been discussed in detail in [6, 7]. Since this is not the focus on the thesis we onlygive the briefest outline.First, we choose ν = (n, l,m) and the one-particle basisϕnlm(r) = Rn(r)Ylm(rˆ). (4.41)Although seemingly related to atomic orbitals, our only reason for choosing this basis is that it gives arepresentation of O(3). That is, the one-particle basis satisfiesϕnlm(Qr) =∑m′Dlm′m(Q)Rn(r)Ylm′(rˆ) ∀Q ∈ O(3), (4.42)where Dlm′m denotes the Wigner-D matrices. I was shown in [6, 7] that from this property alone one canconstruct a sparse coupling matrix C (containing generalized Clebsch-Gordan coefficients) and a resultingtransformed feature defined byB = C ·A,which is O(3)-invaraint. We therefore change the Jastrow factor from (4.40) toJ = C ·B, (4.43)to also gurantee O(3)-invariance of J .4.4.3 Explicitly Incorporating Kato CuspsIn this section, we explicitly demonstrate how cusp conditions can be included in a wavefunctionarichtecture. In principle, there are many ways to carry out the task. We will nevertheless specialize to themethod implemented by Drummond [33], which leads to CASINO [8].Drummond’s method offers a way to build exact cusp conditions into the ansatz while maintaining thelearnability of the Jastrow factors. The main idea is to introduce cusp conditions via linear constraints onthe parameters in the Jastrow factor. Such a formulation requires us to choose a basis. Taking Section 4.4.2as a motivation, we henceforth chooseν = (n, l,m),ϕν(r) = ϕnlm(r) = Rn(r)Ylm(rˆ). (4.44)We can now stateConstraints on ACE-parameterized Jastrow Factor. Set J = Jee + Jen + Jeen, and let each of the threeterms be paramterized by ACE as described in Section 4.4.1. The cusp constraint equations on J are then43given by12=∑nCn00R′n(0+)1√π−ZI =∑nCIn00R′n(0+)1√4π+∑n1n2lmCIn1lmn2lm12πR′n1(0+)Rn2(0+) ∀I0 =∑n1(CIn100n2lm+ CIn2lmn100)R′n1(0+) ∀n2lm. (4.45)Proof. Section 4.5.3The calculation is done for J expanded only up to body-order Jeen because Drummond also assumessuch J in [33]. Generalization to higher body-order is in principle doable, as CASINO demonstrates in [40].However, we will not explore it in this thesis because 1) the calculation is more involved, and 2) the maingoal is just to demonstrate one method (of many) to incorporate cusp conditions.In deriving (4.45), each body ordered term is expanded in the chosen basis, and Jeen is particularlymodified by ACE as written in (4.35). To include these constraints in the implementation, first notice theycan be collectively written asV =MC,where v is LHS of constraint equation written as column vector, C as in (4.40) is the parameter vector, andM is the matrix constituted by coefficients in (4.45). From basic linear algebra, it is known that any solutionto such linear system can be written generally in the formC = c0 +NC˜,where c0 is the particular solution satisfying V = Mc0 and will stay fixed during the optimization. N isa matrix formed by placing the basis vectors of null(M) as columns. C˜ is a new column of independentparameters, and presumably, is smaller than C in size. We can now setJ = C ·A = c0 ·A+ (NC˜) ·A = c0 ·A+ C˜ · (NTA), (4.46)and we simply need to optimize with respect to C˜.This concludes the practical aspect of the cusp implementation. The rest of this subsection will bedevoted to the derivation of the constraints. Several calculations to intermediate results, i.e. lemma, arerelegated to Section 4.5.4.4.4 Backflow TransformationIn the discussion of the Jastrow factor, we have seen that ACE arises as a cost-reducing parameter-ization of functions admitting many-body structure. By assuming a one particle basis (e.g. that given in(4.44)), we can construct a complete representation of J . Using a similar technique (add self-interaction),we will see in this section that ACE can also be adapted to a parameterization of backflow transformed44orbitals χ(i) = χ(xi,x ̸=i) that are widely used in VMC calculations.For the sake of simplicity, let us assume that the “highlighted” particle be i = 1, and that N = 3.Recall that the backflow transformation requires the orbital χ(1) = χ(x1, x̸=1) be symmetric with respectto the non-1-coordinates x ̸=1 = (x2,x3). Let us take a one-particle basis {gµ}, then by fixing x1, we canexpand in terms of direct symmetrization of the chosen basis,χ(x1,x2,x3) =∑µ2≤µ3Cµ2µ3(x1)Sym[gµ2gµ3 ](x2,x3)=∑µ2≤µ3Cµ2µ3(x1) (gµ2(x2)gµ3(x3) + gµ3(x2)gµ2(x3)) .Assume additionally that g0 = 1. That is, the constant 1 is part of the basis. Then we can rearrangeχ(x1,x2,x3) = 2C00(x1)+∑µ2<0Cµ20(x1) (gµ2(x2) + gµ2(x3))+∑0<µ3C0µ3(x1) (gµ3(x3) + gµ3(x2))+µ2,µ3 ̸=0∑µ2≤µ3Cµ2µ3(x1) (gµ2(x2)gµ3(x3) + gµ3(x2)gµ2(x3))= C ′00(x1) · · · redefinition of C00+∑µ̸=0C ′µ(x1) (gµ(x2) + gµ(x3))+µ2,µ3 ̸=0∑µ2≤µ3Cµ2µ3(x1) (gµ2(x2)gµ3(x3) + gµ3(x2)gµ2(x3)) ,where we have definedC ′00 := 2C00C ′µ :=Cµ0, µ < 0C0µ, µ > 0 .This now looks like a many-body expansion. We now invoke the self-interaction trick, and absorb the45self-interaction term into the lower body-ordered expansion.χ(x1,x2,x3) = C′00(x1)+∑µ̸=0C ′µ(x1) (gµ(x2) + gµ(x3))− self+µ2,µ3 ̸=0∑µ2≤µ3Cµ2µ3(x1) (gµ2(x2)gµ3(x2) + gµ2(x3)gµ3(x3))︸ ︷︷ ︸self+µ2,µ3 ̸=0∑µ2≤µ3Cµ2µ3(x1) (gµ2(x2)gµ3(x3) + gµ3(x2)gµ2(x3))= C ′00(x1)+∑µ̸=0C ′′µ(x1) (gµ(x2) + gµ(x3))+µ2,µ3 ̸=0∑µ2≤µ3Cµ2µ3(x1) (gµ2(x2) + gµ2(x3)) · (gµ3(x2) + gµ3(x3)) ,This is now the ACE representation for a symmetric function [11] (in variables x2 and x3), analogous toJeen(4.35). To finally arrive at the fully expanded orbital, we simply “unfix” x1 and expand the coefficient,χ(x1,x2,x3) =∑νCνϕν(x1)+∑νµ ̸=0Cνµϕν(x1) (gµ(x2) + gµ(x3))+µ2,µ3 ̸=0∑νµ2≤µ3Cνµ2µ3ϕν(x1) (gµ2(x2) + gµ2(x3)) · (gµ3(x2) + gµ3(x3)) . (4.47)Notice that we do not assume {ϕν} and {gµ} to be the same bases. A generalization to the complete N -particle representation is now evident and can be exhibited byχ(xi, x̸=i) =N∑B=1µ2,...,µB ̸=0∑νµ2≤···≤µBCνµ2...µBϕν(x1)Aiµ2 · · ·AiµB ,Aiµ :=N∑j ̸=igµ(xj). (4.48)In [21, 42],∑NB=1 in (4.48) is replaced by∑BB=1, where B is a tunable hyperparameter. This gives amodel that systematically resolve correlation order, a hierarchical feature that NN-based wavefunctionsdoes not possess. Most remarkably, by setting B = 1 (and J = 0), we observe that χ(xi, x̸=i) =χ(xi) =∑ν Cνϕν(x1), and the ACE-parameterized Slater-Jastrow-Backflow wavefunction is reduced to46the Hartree-Fock model. This suggests that ACE-parameterized wavefunction is a computationally efficientmany-body generalization of the Hartree-Fock model.Remark 1. From the above derivation, it is critical to assume 1 be included in the basis {gν}. Such assump-tion is adopted in Bachmayr et al. [11, Assumption3.1] in proving completeness of ACE-parameterizationof symmetric function. Without the 1, the self-term trick would not work as there is no lower order terms toabsorb −self , unless the function originally admits many-body structure. On the other hand, [42] proposesthat {ϕν} and {gν} each could be chosen to be quantum chemistry basis sets. However, such basis does notcontain a 1 for molecular systems and consist of linear combinations of atomic orbitals, a family of localizedfunctions (either Gaussian or exponential decay.) We therefore suspect there is a way to rigorously showthat 1 is not needed in the basis for certain circumstances, such as molecular calculation, or for when 1 isincluded in the span of the basis.Remark 2. As remarked in [11], for any symmetric function f , the representation by symmetrized productbasisf(x1, . . . ,xN ) =∑µ1...µNCµ1...µNSym[g1 . . . gN ](x1, . . . ,xN )=∑µ1...µNCµ1...µN∑p∈SNg1(xp(1)) . . . gN (xp(N))is in theory viable, and reduces the number of parameters by a factor of almost N ! compared to a naivetensor product. However, the cost of direct symmetrization is also O(N !), cancelling the earlier factorialsaving during evaluation. Thus, ACE offers again a complete representation of a function at significantlylower evaluation cost.4.4.5 A Comparison Between CASINO and ACE WavefunctionAs the form of the ACE wavefunction ansatz draws inspiration from the CASINO code, it is only fairto give a comparison. We compare the two main components of the explicitly correlated ansatz.Jastrow and cuspThe Jastrow factor proposed by Drummond et al.[33] for the CASINO code [8] reads,Jee =N∑i ̸=juCASINO(rij),Jen =N∑iM∑IηCASINOI (riI),Jeen =N∑i ̸=jM∑IfCASINOI (riI , rjI , rij). (4.49)47Thus, one can see that our proposed ACE Jastrow factor generalizes the CASINO Jastrow as follows.uACE(rij)set= uCASINO(|rij |),ηACEI (riI)set= ηCASINOI (|riI |),fACEI (r1I , r2I)set= fCASINOI (|r1I |, |r2I |, |r1I − r2I |). (4.50)The basis in terms of which uCASINO, ηCASINOI , fCASINOI are expanded is chosen to be simple pow-ers rnij , rniI . Recall from our discussion on cusp implementation (e.g Section 4.4.3) that uACE, ηACEI onlycontribute their Cn00 components. Thus, up to second body-order, our formalism agrees with CASINO’s,provided the radial part of the basis is also chosen to be simple powers. Our formulation differs slightlyfrom CASINO’s in the parameterization of fI . Since fCASINOI explicitly contains the linear term rij (andin fact all its odd powers), ∂fCASINOI/∂rij contributes to the constraint 1/2 = (∂J¯/∂rij )rij↓0, unlessexplicitly set to 0. Indeed, CASINO dose impose the latter. Our formulation does not have this issue as(∂f¯ACEI/∂rij )rij↓0 = 0 exactly (see (4)).Backflow TransformationCASINO’s backflow transformation [8, 43] starts with the pure Slater-Jastrow ansatz. A coordinatetransformation,yi = ri + ξeei + ξeni + ξeeni , (4.51)ξeei =N∑j ̸=iu(rij)rij ,ξeni =M∑IηI(riI)riI ,ξeeni =N∑j ̸=iM∑I(fI(riI , rjI , rij)rij + θI(riI , rjI , rij)riI),is then applied to the electron coordinates ri in the Slater part, upon expanding u, ηI , fI , θI in simplepowers of radial coordinates. The reason for employing such formulation is two-fold. On the one hand,(4.51) is an obvious analogy to their Jastrow parameterization. On the other, backflow transformation asorginally introduced in Feynman and Cohen’s variational approach for Hellium-4 [37] is also such a coor-dinate transformation. Substituting yi into the HF orbital of the form χ(ri, σi) = ϕ(ri)δσis , we see thatχ(yi, σi) = ϕ(yi)δσis is symmetric with respect to x̸=i. Not only permutation SN−1-symmetry, if onedecides to expand and parameterize χ in terms of symmetry adapted basis elements of an irreducible rep-resentation of the molecular point group, χ(yi, σi) then satisfies the same symmetry as χ(ri, σi). LaterNN-approaches [9, 10] takes away only the permutation symmetry of backflow transformed orbitals andparameterize by symmetric neural network. Our ACE-parametrized backflow transformation as investigatedin Section 4.4.4 shares the same spirit. However, ACE has been used to parameterize vector quantity such as48magnetic momenta [44, 45]. This means, contrary to the case of NN-architecture, one can in principle alsoformulate an ACE-parameterized backflow coordinate transformation. Such exploration will corroboratethe robustness of ACE, on top of being a cost-efficient and complete parameterization. We will, however,leave that for future work.4.5 Proofs4.5.1 Diagonalizability of Hartree-Fock UEGAssuming the notation from subsection 4.1.1. We calculate every term, particularly the Coulomb andexchange integral, in HF equations directly. Let us first fix (k, σ). Thenh(1)χk,σ(1) =12k2χk,σ(1), (4.52)∑b ̸=(k,σ)Jb(1)χk,σ(1) =∑b̸=(k,σ)[∫dx2|χb(2)|2GPBC(r1 − r2)]χk,σ(1) (4.53)=[2(N − 1) 1L3∫d3yGPBC(y)]χk,σ(1) = 0 (by (3.12)),∑(q,α)̸=(k,σ)Kq,α(1)χk,σ(1) =∑(q,α)̸=(k,σ)[∫dx2χ∗q,α(2)χk,σ(2)GPBC(r1 − r2)]eiq·r1δσ1α /L3/2︸ ︷︷ ︸χq,α(1)=∑(q,α)̸=(k,σ)∫Ωdr2L3ei(k−q)·r2∑σ2=↑,↓δασ2δσ2σ GPBC(r1 − r2) eiq·r1δσ1α /L3/2=∑q ̸=k[∫Ω+r1dyL3e−i(k−q)·yGPBC(y)]ei(k−q)·r1︸ ︷︷ ︸change variable y=r1−r2eiq·r1δσ1σ /L3/2=∑q ̸=k∫Ω︸︷︷︸periodicitydyL3e−i(k−q)·yGPBC(y) eik·r1δσ1σ /L3/2=∑q ̸=k1L3G˜PBC(k − q)eik·r1δσ1σ /L3/2=∑q ̸=k1L3G˜PBC(k − q)χk,σ(1). (4.54)In the foregoing derivations for both J and K, a variable change y = r1 − r2 is applied. We also used thatGPBC is periodic over L. As a reminder, G˜PBC denotes the Fourier transform of GPBC. This shows the49direct diagonalization of F(1) with eigenvalues(1/2)k2 − 1L3∑q ̸=kG˜PBC(k − q). (4.55)4.5.2 Useful Tools for Cusp Constraint CalculationSwaps Of Limits. Assume f : R3 → R is C1 except possibly at 0. Suppose that its derivative is uniformlybounded, and limr↓0 df¯(r)/dr exists. Thenlimr↓0df¯dr= limr↓0−∫S2drˆ∂f∂r= −∫S2drˆ limr↓0∂f∂r. (4.56)Proof. Since f is C1 on R3\{0}, so is its spherical average f¯ on r > 0. Since the derivative is uniformlybounded, we have by DCT(2.3) that exchange of differentiation and integration is allowed,df¯dr= −∫S2drˆ∂f∂r.Taking lim↓0 of the above yields the first eqality. The second follows from another application of DCT.Regularity Implies No Cusp. Assume f is C1 everywhere with uniformly bounded derivative, thenlimr↓0df¯dr= 0. (4.57)Proof. As f is C1 everywhere, all its directional derivatives exist everywhere and we have in sphericalcoordinates that∂f∂r(rrˆ) = (∇f(rrˆ)) · rˆ=⇒ limr↓0∂f∂r(rrˆ) = (∇f(0)) · rˆ.In particular, the RHS is a dot product between vector constant and angular dependence. Thus,limr↓0df¯dr= −∫S2drˆ limr↓0∂f∂r· · · Swap of Limits above= ∇f(0) · −∫S2drˆrˆ = 0.4.5.3 Cusp Constraints Calculation for (4.45)The key assumption to all the constraints is the preclusion of coalescences beyond 2 particles. Kato’scondition (4.13) tells us nothing when the exact wavefunction assumes such coalescences. Thus, we will50simply assume that when, say, electron 1 and electron 2 are approaching each other (r12 ↓ 0) in such a waythat their CM stays fixed, all other particles (including nuclei) are fixed at positions distinct from r¯12 andfrom each other. With this assumption in mind, we can derive the following intermediate results.Lemma 4. (∂J¯en∂r12)r12↓0=(∂J¯een∂r12)r12↓0=(∂J¯ee∂r1I)r1I= 0. (4.58)Lemma 5. (∂J¯ee∂r12)r12↓0=∑nR′n(0+)Cn00√π. (4.59)Lemma 6. (∂J¯en∂r1I)r1I↓0=∑nCIn00R′n(0+)1√4π. (4.60)Lemma 7.(∂J¯een∂r1I)r1I↓0=∑n1n2lmCIn1lmn2lm12πR′n1(0+)Rn2(0+)+∑n1n2lm(CIn100n2lm+ CIn2lmn100)R′n1(0+)1√4π∑j ̸=1Rn2(rjI)Ylm(rˆjI). (4.61)Let us postpone the proofs of these lemmas, and first derive the final result (4.45) using them. The firstequation in (4.45) is given by12=(∂J¯∂r12)r12↓0=(∂J¯ee∂r12)r12↓0+(∂J¯en∂r12)r12↓0+(∂J¯een∂r12)r12↓0︸ ︷︷ ︸=0=∑nR′n(0+)Cn00√π.For second and third, cusp condition (4.23) tells us−ZI =(∂J¯ee∂r1I)r1I↓0︸ ︷︷ ︸=0+(∂J¯en∂r1I)r1I↓0+(∂J¯een∂r1I)r1I↓0=∑nCIn00R′n(0+)1√4π+∑n1n2lmCIn1lmn2lm12πR′n1(0+)Rn2(0+)+∑n1n2lm(CIn100n2lm+ CIn2lmn100)R′n1(0+)1√4π∑j ̸=1Rn2(rjI)Ylm(rˆjI)51By linear independence of{∑j ̸=1Rn2(rjI)Ylm(rˆjI)|n2lm}(which looks like ACE basis), constants aregrouped together−ZI =∑nCIn00R′n(0+)1√4π+∑n1n2lmCIn1lmn2lm12πR′n1(0+)Rn2(0+),and coefficients of non-constants are zeros0 =∑n1(CIn100n2lm+ CIn2lmn100)R′n1(0+).The cusp constraint (4.45) thus follows. We now prove the lemmas.Proof of Lemma 4. Consider Jen. It is a linear combination of ϕnlm(riI). In the limit r12 ↓ 0, riI is eitherunmoved or goes to r1I , depending on what i is. Because we assume no more than two particles overlapand riI involve a nucleus, the limit either way is non-zero, around which ϕnlm is C∞. We can thus apply(4.56) and (4.57) to deduce−∫S2drˆ12(∂Jen∂r12)r12↓0=∑iInlmCInlm−∫S2drˆ12(∂∂r12ϕnlm(riI))r12↓0= 0.In addition, ∂J¯een/∂r12 = ∂J¯ee/∂r1I = 0 by the same argument.Proof of Lemma 5. (∂J¯ee∂r12)r12↓0=∑i ̸=j−∫S2drˆ12(∂∂r12u(rij))r12↓0Now notice, by (4.56) and (4.57), for (i, j) ̸= (1, 2) or (2, 1) we have−∫S2drˆ12(∂∂r12u(rij))r12↓0= 0.Then (∂J¯ee∂r12)r12↓0= −∫S2drˆ12(∂∂r12u(r12))r12↓0+−∫S2drˆ12(∂∂r12u(r21))r12↓0.52Consider the two terms in the sum separately. We utilize the orthonormality of Ylm to get−∫S2drˆ12(∂∂r12u(r12))r12↓0=∑nlmCnlmR′n(0+)−∫S2drˆ12Ylm(rˆ12)=∑nlmCnlmR′n(0+)δ0l δ0m1√4π(orthogonality)=∑nCn00R′n(0+)1√4π.Similarly,−∫S2drˆ12(∂∂r12u(r21))r12↓0=∑nlmCnlmR′n(0+)−∫S2drˆ12Ylm(rˆ21)=∑nlmCnlmR′n(0+)−∫S2drˆ12 (−1)lYlm(rˆ12)︸ ︷︷ ︸parity symmetry of Ylm=∑nCn00R′n(0+)1√4π.Adding the two yield the desired.Proof of Lemma 6. Again using the fact that no more than 2 particles overlap, (4.56) and (4.57), and or-thonormality, (∂J¯en∂r1I)r1I↓0=∑iB−∫S2drˆiB(∂∂r1IηB(riB))r1I↓0=(∂ηI(r1I)∂r1I)r1I↓0=∑nlmCInlmR′n(0+)−∫S2drˆ1IYlm(rˆ1I)=∑nCIn00R′n(0+)1√4π.Proof of Lemma 7. We first separate the ACE parameterized een contribution into self-term and (i ̸= j)-53term. (∂J¯een∂r1I)r1I↓0=∑ij∑n1l1m1n2l2m2CIn1l1m1n2l2m2(∂∂r1Iϕn1l1m1(riI)ϕn2l2m2(rjI))r1I↓0=∑i∑n1l1m1n2l2m2CIn1l1m1n2l2m2(∂∂r1Iϕn1l1m1(riI)ϕn2l2m2(rjI))r1I↓0+∑i ̸=j∑n1l1m1n2l2m2CIn1l1m1n2l2m2(∂∂r1Iϕn1l1m1(riI)ϕn2l2m2(rjI))r1I↓0.Note that terms involving nuclei other than I vanish due to strict two-particle coalescence as before. Workingseparately, we deriveN∑i∑n1l1m1n2l2m2CIn1l1m1n2l2m2(∂∂r1Iϕn1l1m1(riI)ϕn2l2m2(riI))r1I↓0=∑n1l1m1n2l2m2CIn1l1m1n2l2m2(R′n1(0+)Yl1m1(rˆ1I)Rn2(0+)Yl2m2(rˆ1I) +Rn1(0+)Yl1m1(rˆ1I)R′n2(0+)Yl2m2(rˆ1I))=∑n1l1m1n2l2m22CIn1l1m1n2l2m2R′n1(0+)Rn2(0+)−∫S2drˆ1IYl1m1(rˆ1I)Yl2m2(rˆ1I); (relabel)=∑n1l1m1n2l2m22CIn1l1m1n2l2m2R′n1(0+)Rn2(0+)14πδl1l2δm1m2 ,54which is the first term in RHS of Lemma 7. On the other hand,∑i ̸=j∑n1l1m1n2l2m2CIn1l1m1n2l2m2(∂∂r1Iϕn1l1m1(riI)ϕn2l2m2(rjI))r1I↓0=∑i ̸=j∑n1l1m1n2l2m2CIn1l1m1n2l2m2Yl1m1(rˆiI)Yl2m2(rˆjI)(δ1iR′n1(0+)Rn2(rjI) +Rn1(riI)δ1jR′n2(0+))=∑n1l1m1n2l2m2CIn1l1m1n2l2m2N∑iδ1iR′n1(0+)Yl1m1(rˆiI)N∑jj ̸=iRn2(rjI)Yl2m2(rˆjI)+∑n1l1m1n2l2m2CIn1l1m1n2l2m2N∑jδ1jR′n2(0+)Yl2m2(rˆjI)N∑ii ̸=jRn1(riI)Yl1m1(rˆiI)=∑n1l1m1n2l2m2CIn1l1m1n2l2m2R′n1(0+)Yl1m1(rˆ1I)N∑jj ̸=1Rn2(rjI)Yl2m2(rˆjI)+∑n1l1m1n2l2m2CIn1l1m1n2l2m2R′n2(0+)Yl2m2(rˆ1I)N∑ii ̸=1Rn1(riI)Yl1m1(rˆiI)=∑n1l1m1n2l2m2CIn1l1m1n2l2m2R′n1(0+)(−∫S2drˆ1IYl1m1(rˆ1I)) N∑jj ̸=1Rn2(rjI)Yl2m2(rˆjI)+∑n1l1m1n2l2m2CIn1l1m1n2l2m2R′n2(0+)(−∫S2drˆ1IYl2m2(rˆ1I)) N∑ii ̸=1Rn1(riI)Yl1m1(rˆiI)=∑n1n2l2m2CIn100n2l2m2R′n1(0+)1√4πN∑j ̸=1Rn2(rjI)Yl2m2(rˆjI)+∑n1l1m1n2CIn1l1m1n200R′n2(0+)1√4πN∑i ̸=1Rn1(riI)Yl1m1(rˆiI)=∑n1n2lm(CIn100n2lm+ CIn1lmn200)R′n1(0+)1√4πN∑j ̸=1Rn2(rjI)Ylm(rˆjI),which is the second term in RHS of Lemma 7.55Chapter 5Numerical Implementation of VMC5.1 MotivationThe system we will use in our simulations is the 1D UEG introduced in Sections 2.1.1 and 3.3. Wehave chosen 1D UEG for 3 reasons.1. This will be the first implementation of ACE-parameterized wavefunctions on extended systems, andwe want to explore what extra considerations, not present in molecular simulations, there are. Thisembodies Section 3.2. ACE-parameterized wavefunctions have been applied to 1D molecules in the work by Zhou et al.[21]. Notably, their numerical demonstration showcased the sufficiency of a single determinant withgeneralized orbitals to converge to the 1D molecular ground state. We want to see if that remains thecase for an extended system.3. Parameters in the UEG model can be tuned to reach low and high correlation regime. In particular, itis known in the literature that when the number density ρ0 = N/V of UEG increases, the kinetic termin the Hamiltonian becomes dominant, and the electrons tend to the low correlation (HF) regime. Onthe other hand, when ρ0 gets low, the electrons tend to localize, giving rise to the so-called WignerCrystal [46–48]. In this regime, the electrons are highly correlated. Altogether, UEG offers a stringentassessment of the efficacy of ACE-parameterized wavefunction across tunable correlation regime.The primary objective of the first point is to acquaint us with Ewald Summation and the formalities of poten-tial periodization. Consequently, our simulation will predominantly focus on the second aspect. To commentmore on the exclusion of Jastrow factor (see (4.23)), we note that there is no Kato cusp condition provenfor 1D system. Furthermore, as demonstrated in [21, 49], it turns out that a single determinant with gener-alized orbitals constitute a complete representation of 1D wavefunctions, provided that the wavefunction issufficiently smooth.565.2 Recap and Basic Set-upsLet us recall from (3.28) that the Hamiltonian is given byH = −12N∑i=1∂2∂x2i+ VUEG(x1, . . . , xN ),VUEG(x1, . . . , xN ) =12N∑i=1N∑j=1j ̸=iVb(xi − xj) + NVMad2,VMad = limx→0(Vb(x)− vb(x)) = limx→0Vb(x)−√π2b,where the periodized interaction is defined byVb(x) :=∑m∈Z(vb(x−mL)−∫ L/2−L/2dyLvb(x− y −mL)),vb(x) :=√π2bexp(x24b2)erfc( |x|2b),but in actuality evaluated by Ewald SummationVb = Vb ⋆L gb + (Vb − Vb ⋆L gb),Vb ⋆L gb =2L∞∑k=1E1([2πbkL]2)cos(2πkxL),Vb − Vb ⋆L gb =∑m∈Z(vb(x−mL)− 1|x−mL| erf( |x−mL|2b)).The truncation K and M for, respectively, Vb ⋆L gb and Vb−Vb ⋆L gb are chosen according to Theorem 3.3.2so that the error incurred in evaluation of Vb is O(10−8), as implemented in [15].To specify regimes for which we simulate, first recall from (2.6) that b is interpreted as the width ofthe wire. Abiding by the convention in existing works [15, 16, 48], let us also introduce the Wigner-Seitzradius (for 1D),rs := L/(2N), (5.1)where L = |Ω| is the size of the simulation cell and N is the number of electrons. We now declare that werun simulations with N = 30, on four ambience regimes specified by(rs, b) = (0.5, 1.0), (1.0, 2.0), (1.0, 1.0), (5.0, 1.0). (5.2)We will also only deal with paramagnetic states. That is, the number of ↑-electrons is identical to the numberof ↓-electrons.57According to [48], high rs (low density) and low b (thin wire) would correspond to strong correlationeffect. On the other hand, low rs and high b correspond to weak correlation. Intuitively then, we expect theHartree-Fock model to be more accurate for the first two regimes in (5.2) than the last two regimes. We willexplore this intuition, and see if increasing the correlation order B of the parameterization (4.48) resolve theexact ground state better. The performance of our VMC simulation will be compared against [48], whichsimulates the same systems with the CASINO code [8].5.3 Wavefunction ArchitecturesCommon to all the trial wavefunctions to be considered in this chapter is the form of generalizedorbitals. Writing xi = (xi, σi), we will assumeχj,s(xi,x ̸=i) = φj(xi, x̸=i)δσis =φj(xi, x̸=i)δσi↑ j ≤ 15φj(xi, x̸=i)δσi↓ j > 15. (5.3)This is essentially Unrestricted-Hartree-Fock theory for when B = 1 (in (4.48)). By convention, VMC isimplemented with spin-assigned wavefunction. This means, we can assumeσi =↑ i ≤ 15↓ i > 15 . (5.4)As a result, the trial wavefunction factors. That is,ΨACE =∣∣∣∣∣∣∣∣∣∣φ1(x1,x ̸=1) φ2(x1, x̸=1) . . . φ15(x1,x ̸=1)φ1(x2,x ̸=2) φ2(x2, x̸=2) . . . φ15(x2,x ̸=2)....... . ....φ1(x15,x ̸=15) φ2(x15, x̸=15) . . . φ15(x15,x ̸=15)∣∣∣∣∣∣∣∣∣∣︸ ︷︷ ︸highlighted electrons with ↑∣∣∣∣∣∣∣∣∣∣φ16(x16, x̸=16) φ17(x16, x̸=16) . . . φ30(x16,x ̸=16)φ16(x17, x̸=17) φ17(x17, x̸=17) . . . φ30(x17, x̸=17)....... . ....φ16(x30, x̸=30) φ17(x30,x ̸=30) . . . φ30(x30, x̸=30)∣∣∣∣∣∣∣∣∣∣︸ ︷︷ ︸highlighted electrons with ↓(5.5)As suggested by the expression (4.48), we will need to choose a basis. Since UEG is periodic and planewaves turn out to be the exact solution of the Hartree-Fock model (discussed in Section 5.4), the basis ischosen to be{ϕν} = {gµ} = {sin(kx), cos(kx) | k ∈ (2π/L)Z≥0},where g and ϕ are those in (4.48).We now declare three wavefunctions with which we will numerically explore.585.3.1 ACE-parameterized BackflowModulo the choice of basis, this is the wavefunction used throughout the molecular examples in Zhouet al. [21]. It turns out for them that a total body-order B = 2 suffices to converge in almost all simulations.Since it also helps with clarity, let us assume body-order B = 2 for now. (4.48) thus readsφj(xi,x ̸=i) =∑νµ∑β=↑,↓Cjνµβgν(xi)Aiµβ=∑νµ(Cjνµ↑gν(xi)Aiµ↑ + Cjνµ↓gν(xi)Aiµ↓),Aiµs :=30∑j ̸=igµ(xj)δσjs , (5.6)so we see Aiµs is a sum of 15 terms, rather than 30, depending on the spin s. Note that every spatial degreeµ (or ν) gives rise to two basis labels (µ, ↑) and (µ, ↓). Also, the highlighted electron xi and its surrounding{xj}j ̸=i are embedded into the same basis functions “gµ· spin”.The case for general B is exhibited byφj(xi,x ̸=i) =∑νµ∑β=spinsCjνµβgν(xi)Aiµβ, (5.7)Aiµβ := Aiµ1s1 · · ·AiµBsB ,where µ,β are now multi-indices of length B.5.3.2 ACE-parameterized Backflow with Coordinate DifferenceWe replace xj in Aiµs by xj − xi. This givesφj(xi, x̸=i) =∑νµ(Cjνµ↑gν(xi)Aiµ↑ + Cjνµ↓gν(xi)Aiµ↓), (5.8)Aiµs :=30∑j ̸=igµ(xj − xi)δσjs .Note that periodicity is still preserved. That is, Aiµs(. . . , xj+L, . . .) = Aiµs(. . . , xj , . . .), ∀j ̸= i. Moreover,if one redoes the derivation as given in Section 4.4.4 and begins withχ(x1, x2, x3) =∑µ2≤µ3Cµ2µ3(x1)Sym[gµ2gµ3 ](x2 − x1, x3 − x1),one will conclude that this variant (5.8) of ACE parameterization is also a complete representation usingverbatim arguments. The case for general B is given by the same notation as in (5.7).595.3.3 Wigner AnsatzShortly after Fock proposed the single-determinant ansatz [30], Wigner [50] proposed a modificationthat turns out to be a special case of the preceding ansatz with coordinate differences, where features aredropped. That is,φj(xi, x̸=i) =∑νµCjνµ↓gν(xi)Aiµ↓ j ≤ 15∑νµCjνµ↑gν(xi)Aiµ↑ j > 15, (5.9)Aiµs :=30∑j ̸=igµ(xj − xi)δσjs . (5.10)In words, if the highlighted electron i is ↑, we correlate it with those of ↓. Similarly, for a highlightedelectron of ↓, we correlate it with those of ↓. The motivation in Wigner’s proposal lies in the fact that theHF ansatz factors into a product of two determinants,ΨHF =∣∣∣∣∣∣∣∣∣∣φ1(x1) φ2(x1) . . . φ15(x1)φ1(x2) φ2(x2) . . . φ15(x2)....... . ....φ1(x15) φ2(x15) . . . φ15(x15)∣∣∣∣∣∣∣∣∣∣︸ ︷︷ ︸highlighted electrons with ↑∣∣∣∣∣∣∣∣∣∣φ16(x16) φ17(x16) . . . φ30(x16)φ16(x17) φ17(x17) . . . φ30(x17)....... . ....φ16(x30) φ17(x30) . . . φ30(x30)∣∣∣∣∣∣∣∣∣∣︸ ︷︷ ︸highlighted electrons with ↓.This suggests that ↑-electrons and ↓-electrons are statistically independent in HF theory. Wigner’s modifi-cation alleviates this lack of correlation description.Coming back to the ACE-formalism, a possible generalization of Wigner’s ansatz to general correlationorder B isφj(xi,x ̸=i) =∑νµCjνµβgν(xi)Aiµ↓ j ≤ 15∑νµCjνµβgν(xi)Aiµ↑ j > 15, (5.11)Aiµ↑ := Aiµ1↑ · · ·AiµB↑; Aiµ↓ := Aiµ1↓ · · ·AiµB↓.However, it is less apparent whether this representation is theoretically complete or not.Remark 3. To our knowledge, there has not been an implementation of Wigner ansatz since Wigner’s firstproposal.605.4 Hartree-Fock Theory of 1D UEG - Benchmarking ACETheorem 5.4.1. The HF equations of the 1D UEG model, as reviewed in (5.2), can be diagonalized byplanewaves. More precisely, leth(1) = −12∆1,vee(x) = Vb(x),orbital index = (k, σ) ∈ (2π/L)Z× {↑, ↓},χk,s(x, σ) = eikxδσs /L1/2,∫dx =∑σ=↑,↓∫Ωdx.Then the diagonalization is given byF(1)χk,σ(1) = ϵ(k)χk,σ(1),ϵ(k) =12k2 − 1L∑q ̸=k|q|<kFV˜b(k − q),and the HF ground state energy is given by⟨ΨHF|H|ΨHF⟩⟨ΨHF|ΨHF⟩ =∑|k|<kFk2 − 1L∑|k|,|q|<kFvˆb(k − q), (5.12)where kF := 2π(N/4)/L is the Fermi wave number, and vˆb is given by (2.8).Proof. For diagonalizability, arguments used in 3D case (Theorem 4.1.1) can be reused verbatim so we willnot repeat it. For ground energy computation, see Section 5.6.1.This theorem tells us that planewaves constitute the optimal HF ground state when the basis sizeis minimal (i.e when we truncate at a spatial basis set of size N/2). For our 30-electron system, thereare thus 15 occupied planewaves (one for every 2 spins), hence 15 wave-numbers k being summed incalculating ⟨ΨHF|H|ΨHF⟩ / ⟨ΨHF|ΨHF⟩, namely k = 0,±1, . . . ,±7. We will denote these plane waves byw1, . . . , w15, and thus|ΨHF⟩ =∣∣∣∣∣∣∣∣∣∣w1(x1) w2(x1) . . . w15(x1)w1(x2) w2(x2) . . . w15(x2)....... . ....w1(x15) w2(x15) . . . w15(x15)∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣w1(x16) w2(x16) . . . w15(x16)w1(x17) w2(x17) . . . w15(x17)....... . ....w1(x30) w2(x30) . . . w15(x30)∣∣∣∣∣∣∣∣∣∣(5.13)A good initial guess can reduce the run-time of VMC simulation by starting the optimization at the“right” range of energy values. For our ACE model, since we retrieve HF model by setting body-order B to611 (as discussed after (4.48)), we can naturally use HF theory to set up the suitable initial trial wave function.We outline how this is done for our 1D UEG.Setting B = 1 (so no Aiµs) and truncating at ν ≤ 15, the ↑-determinant in (5.5) is given by∣∣∣∣∣∣∣∣∣∣φ1(x1) φ2(x1) . . . φ15(x1)φ1(x2) φ2(x2) . . . φ15(x2)....... . ....φ1(x15) φ2(x15) . . . φ15(x15)∣∣∣∣∣∣∣∣∣∣=∣∣∣∣∣∣∣∣∣∣g1(x1) g2(x1) . . . g15(x1)g1(x2) g2(x2) . . . g15(x2)....... . ....g1(x15) g2(x15) . . . g15(x15)∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣C11 C21 . . . C151C12 C22 . . . C152....... . ....C115 C152 . . . C1515∣∣∣∣∣∣∣∣∣∣.(5.14)Note that the matrix of φ factors into a product of g and coefficients C. Similarly for ↓-determinant. Nowsince the basis set {gµ}15ν=1 consisting of cos and sin can be linearly transformed to planewaves {wl}15l=1 viaa unitary transformation (which can be constructed by Euler’s formula: eiθ = cos θ + i sin θ), we get forsome α ∈ R ∣∣∣∣∣∣∣∣∣∣g1(x1) g2(x1) . . . g15(x1)g1(x2) g2(x2) . . . g15(x2)....... . ....g1(x15) g2(x15) . . . g15(x15)∣∣∣∣∣∣∣∣∣∣= eiα∣∣∣∣∣∣∣∣∣∣w1(x1) w2(x1) . . . w15(x1)w1(x2) w2(x2) . . . w15(x2)....... . ....w1(x15) w2(x15) . . . w15(x15)∣∣∣∣∣∣∣∣∣∣, (5.15)where the phase eiα comes from determinant of the unitary transformation. Similarly for ↓-determinant.Thus, for any coefficients C, we see that⟨ΨACE|H|ΨACE⟩⟨ΨACE|ΨACE⟩ =⟨ΨHF|H|ΨHF⟩ |C|2|eiα|2⟨ΨHF|ΨHF⟩ |C|2|eiα|2 =⟨ΨHF|H|ΨHF⟩⟨ΨHF|ΨHF⟩ . (5.16)The fact that the coefficients C can be arbitrary means that we can simply set those coefficients corre-sponding to basis functions not included in {gµ}15ν=1 to be 0, and expect the first iteration (expensive) of theoptimization to evaluate to the same energy as that given by (5.12) (cheap).5.5 Ground State ComputationFor each of the four regimes (5.2), we set all our trial wavefunctions to B = 2 with the maximalwave number at body-order one to be 7 (so there are 15 such basis functions), and that at body-order twoto be 2 (so 5 such basis functions). This amounts to 4650 parameters for ACE-parameterized backflowwavefunctions, and 2550 for parameters for Wigner ansatz. With accord to the paragraph succeeding (5.16),we initialize the trial wavefunction so that it computes the minimal HF determinant at the first epoch. TheVMC optimization subsequently proceeds through 200 iterations, and the convergence performance arepresented in Figures 5.1, 5.2, 5.3, and 5.4. Statistical error (standard deviation) per electron is plotted as aband.The hyperparameters, “burnin”, “N chain”, “∆t”, and “lr” are selected with accord to Zhou et al.[21] to strike a balance between accuracy and efficiency. Briefly,“burnin” is the number of steps (with-62out optimization) we let the VMC algorithm take before the first epoch. This ensures the distributionΠ(x) = |ΨACE(x)|2/ ⟨ΨACE|ΨACE⟩ (see (2.11)) is accurately represented in the initial epoch. “N chain”corresponds to the number S of samples (see (2.11)). ∆t is the stepwidth in Metropolis Hastings Algorithm[19, 20]. Finally, “lr” is the learning rate at which the gradient descent updates the (learnable) parameters θin (2.11).The HF energy with minimal basis, serving as a benchmark upperbound, is computed with the formula(5.12). The CASINO ground state energy, serving as a target value, is taken from the Appendix of [48].5.5.1 DiscussionTo reiterate, the results are given in Figures 5.1, 5.2, 5.3, and 5.4.Based on the variational principle (2.5), we see that the ACE-parameterized backflow with coordinatedifference prevails as the best performing trial wavefunction in all four regimes. However, it is also clearfrom the plots that as the electron correlation is increased, the performance deteriorate in comparison tothe CASINO code. We have tried increasing the degree and body-order for the (rs, b) = (1, 1) simulation,and it does not lead to noticeable improvement. We have not tried resolving the first two regimes (rs, b) =(0.5, 1), (1, 2) further with higher polynomial degree and body-order. As our computations are already closeto the target value and each of the two regimes exhibits a smaller gap between the HF energy and CASINO’sground state than the case (rs, b) = (1, 1), the improvement is expected to be minute when the degree andbody-order do increase.The last regime (rs, b) = (5, 1) exhibits a behavior distinct from the other three (see Figure 5.4). Onone hand, it does not seem like 200 iterations is enough for the algorithm to converge. On the other hand,the convergence of the Wigner ansatz resembles that of the ACE parameterization with coordinate differ-ence, whereas in the case of the other three where the correlation is lower (judging from the gap betweenHF and CASINO), the convergence resembles that of ACE-parameterized backflow without coordinate dif-ference. Since the electrons tend to localize more in (rs, b) = (5, 1), the wavefunction exhibits a sharperprofile/shape, and higher polynomial degree is thus required to resolve the ground state. We tried increas-ing the number of iterations for degree [15, 5], and afterwards increased the degree to [21, 9]. The result isplotted in Figure 5.5. We see that the Wigner ansatz and the ACE-backflow eventually converge to the sameenergy, though each is still unable to reach the target value. Increasing the degree further does improve theenergy in both cases, but the rate at which the energy lowers is significantly slower, so we have decided notto include that portion of the computation. On the other hand, increasing the body-order B does not improve,and in fact worsens the energy value.The computation performed within (rs, b) = (5, 1) seems to suggest that the electrons of a certain spin(say ↑) tend to correlate with the electrons of the opposite spin (↓) more than those surrounding electronsof the same spin in the Wigner crystal regime. On the contrary, the electrons tend to correlate with itssurrounding electrons equally rather than biasing one spin when the correlation effect is low, as the Wigneransatz exhibits suboptimal performances in Figures 5.1 and 5.2.Recall that Zhou et al. [21] utilizes the ACE-parameterized backflow (without coordinate difference)wavefunction As it exhibits suboptimal convergence in all regimes, it is evident that the wavefunction of631D UEG demands more considerations than the 1D molecular model, with which accurate ground states areobtained without coordinate difference [21]. We suspect that such deviation is rooted in the difference ofsmoothness between the electron-electron interaction Zhou et al. uses and the one we simulate. While our Vbexhibits cusp-like behavior as shown in Figure 3.1, Zhou et al. simulates 1D molecule with soft-Coulomb,vee(x) =1√1 + x2, (5.17)which is a smooth function. It is thus natural to expect the wavefunction for Vb manifests some singularity,and the higher the electron correlation, the more singular. Since all our single determinant wavefunctionarchitectures are smooth, such singularity can only be inefficiently converged. With the introduction ofJastrow factor, as included in the wavefunction architecture of CASINO code [15, 29, 48], the accuracyshould improve. (See Section 6 - Conclusion and Outlooks.)64Figure 5.1: VMC Convergence for rs = 0.5, b = 1.65Figure 5.2: VMC Convergence for rs = 1, b = 2.66Figure 5.3: VMC Convergence for rs = 1, b = 1.67Figure 5.4: VMC Convergence for rs = 5, b = 1.68Figure 5.5: VMC Convergence for rs = 5, b = 1.695.6 Proofs5.6.1 Minimal HF Ground State of 1D UEGOne just needs to substitute in spin-assigned plane-waves, and invoke the periodicity of Vb to do thealgebra. We will nevertheless compute the ground state energy, which is not computed for the 3D case, butthe steps used here can be mirrored. Recalh(1) = −12∆1vee(x) = Vb(x)index i = (k, σ) ∈ (2π/L)Z× {↑, ↓}χi(x) = χk,σ(x, s) = eik·xδσs /L1/2∫dx =∑σ=↑,↓∫ L/2−L/2dx.We now substitute plane-wave basis into (4.10) to obtain[χk,β| − 12∆1|χk,β] = 12L∑σ∫ L/2−L/2dxe−ikxδβσk2eikxδβσ =k22, (5.18)[χk,βχk,β|χq,αχq,α] = 1L2∫dx1dx2e−ikx1δβσ1eikx1δβσ1Vb(x12)e−iqx2δασ2eiqx2δασ2=1L2∫ L/2−L/2dx1∫ L/2−L/2dx2Vb(x1 − x2)=1L∫ L/2−L/2dyVb(y) =V˜b(0)L= 0 (by (3.42)), (5.19)[χk,βχq,α|χq,αχk,β] =δαβL2∫ L/2−L/2dx1e−i(k−q)x1∫ L/2−L/2dx2Vb(x1 − x2)ei(k−q)x2 (5.20)=δαβL2∫ L/2−L/2dx1∫ L/2+x1−L/2+x1dyVb(y)e−i(k−q)y (5.21)=δαβL2∫ L/2−L/2dx1∫ L/2−L/2dyVb(y)e−i(k−q)y (5.22)=δαβ V˜b(k − q)L. (5.23)70Adding these up, the ground state energy at the HF level is exhibited by⟨ΨHF|H|ΨHF⟩ = 12∑β∑kk2 +12∑βα∑kq(−δαβ V˜b(k − q))(5.24)=∑kk2 − 1L∑kqV˜b(k − q) (5.25)=∑kk2 − 1L∑kqvˆb(k − q) (by (3.42)) (5.26)This shows the desired formula.We have not specified what wave-number k and q to include in the computation of the ground stateenergy here, as well as the orbital energy (HF eigenvalues) for the 3D case earlier. In theory, they shouldcorrespond to the N/2 plane-waves yielding the lowest ⟨ΨHF|H|ΨHF⟩. However, there does not seem tobe a methodological approach to determine just what these N/2 plane-waves are. In fact, treatise of UEGpresented in classic condensed matter textbooks such as Ashcroft and Mermin [12] or Mardar [13] simplyassume the wave-numbers are those satisfying |k| < kF 1, where kF denotes the radius of Fermi sphere.Note that we can still, based on Variation Principle, say that any N/2 plane-waves constitute an HF upperbound to the exact ground state energy. We will thus abide by conventional practice and include in the sumthose wave-numbers satisfying |k| < kF to compute the variational upper bound later used to benchmarklater VMC calculation.1These wave-numbers correspond to exact HF ground state for non-interacting electrons\fermions.71Chapter 6Conclusion and OutlooksIn conclusion, this thesis has delved into the intricacies of wavefunction methods, with a primary fo-cus on Variational Monte Carlo (VMC). Chapter 3 provided a comprehensive exploration of the technicalaspects essential for simulating solid state systems, emphasizing the implementation of Ewald summation tonumerically evaluate potential energy with efficiency and controlled error. In Chapter 4, we scrutinized thelimitations of conventional quantum chemistry and advocated for the adoption of explicitly correlated wave-function methods. The integration of Atomic Cluster Expansion (ACE) into these methods was discussed,enhancing the flexibility, computational efficiency, and mathematical rigor of the wavefunction ansatz. No-tably, we parameterized the Jastrow factor and the backflow transformation, laying the groundwork forfurther developments. The application of VMC code to the 1D Uniform Electron Gas (UEG) in Chapter 5yielded suboptimal results, yet provided valuable insights into the interpretation of the physics underlyingthe various wavefunction architectures.Looking ahead, the theoretical insights gained from Chapters 3 and 4 will be refined and extended tothe ongoing implementation for 3D systems, whether molecular or solid state. These insights will serveas a mathematical foundation to substantiate the creative enhancements incorporated into our wavefunctionarchitectures. Chapter 4 also serves as a foundation that will help us identify connections between ACE-parameterized wavefunctions with existing wavefunction architectures. To give one such example, uponlearning via Chapter 5 that parameterization with body order beyond B = 2 never enhanced the accuracy ofthe model, one appropriate direction seems to be an investigation into Antisymmetric Geminal Power [51],where the “building blocks” are the pairs of orbitals instead of the individuals.Taking a more open-ended perspective, Chapter 5 underscored the ubiquitous use of the term “corre-lation” without a comprehensive understanding. Our future endeavors will involve a rigorous mathematicalexploration of this concept, seeking to quantify correlation and illuminate its role in electronic systems.Initially, the total body order B was considered a measure of the “strength of correlation,” but our findingsindicate a more nuanced approach is necessary, as parameterization beyond B = 2 did not improve ourmodel accuracy. This realization prompts the need for a more intricate scheme of quantification.In this work, the 1D UEG model merely served as a simple tool for elucidating the numerical aspectsof PBC and the efficacy of ACE-parameterization. Nevertheless, it would be particularly beneficial if thereverse can be achieved. That is, utilizing ACE-VMC to glean physical insights from the system of interest.72Throughout our examination of the 1D model Vb, we have come to appreciate that 1D physics lends itself tomore analytical techniques, such as bosonization [52] and the Bethe ansatz [53, 54], compared to 2D or 3Dphysics. These methodologies, grounded in physical insights like field theory, yield correlation functionsthat unveil additional nuances of the system. As ACE is also more amenable to analytical means than NNor CASINO, it therefore seems worthwhile to investigate its connections with the analytical frameworkoffered by 1D physics, in order to establish ACE as a practical method and thereby enhancing its broaderapplicability.73Bibliography[1] T. Helgaker, P. Jorgensen, and J. Olsen, Molecular electronic-structure theory. John Wiley & Sons,2013. → pages viii, 2, 3, 4, 32, 33, 35, 36[2] J. A. Pople, “Nobel lecture: Quantum chemical models,” Review of Modern Physics, vol. 71, pp.1267–1274, Oct 1999. [Online]. Available: https://link.aps.org/doi/10.1103/RevModPhys.71.1267 →page 1[3] W. Kohn, “Nobel lecture: Electronic structure of matter—wave functions and densityfunctionals,” Review of Modern Physics, vol. 71, pp. 1253–1266, Oct 1999. [Online]. Available:https://link.aps.org/doi/10.1103/RevModPhys.71.1253 → page 1[4] A. Szabo and N. S. Ostlund, Modern quantum chemistry: Introduction to advanced electronic structuretheory. Courier Corporation, 2012. → pages 1, 2, 3, 4, 32, 33, 34[5] P. Hohenberg and W. Kohn, “Inhomogeneous electron gas,” Physical Review, vol. 136, pp.B864–B871, Nov 1964. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRev.136.B864 →page 1[6] R. Drautz, “Atomic cluster expansion for accurate and transferable interatomic potentials,” PhysicalReview B, vol. 99, p. 014104, Jan 2019. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevB.99.014104 → pages 2, 41, 43[7] G. Dusson, M. Bachmayr, G. Csa´nyi, R. Drautz, S. Etter, C. van der Oord, and C. Ortner, “Atomiccluster expansion: Completeness, efficiency and stability,” Journal of Computational Physics,vol. 454, p. 110946, 2022. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0021999122000080 → pages 2, 42, 43[8] R. J. Needs, M. D. Towler, N. D. Drummond, P. Lo´pez Rı´os, and J. R. Trail, “Variational and diffusionquantum Monte Carlo calculations with the CASINO code,” The Journal of Chemical Physics, vol.152, no. 15, p. 154106, 04 2020. [Online]. Available: https://doi.org/10.1063/1.5144288 → pages2, 37, 40, 42, 43, 47, 48, 58[9] D. Pfau, J. S. Spencer, A. G. D. G. Matthews, and W. M. C. Foulkes, “Ab initio solution of themany-electron schro¨dinger equation with deep neural networks,” Physical Review Reseach, vol. 2, p.74033429, Sep 2020. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevResearch.2.033429→ pages 2, 37, 40, 48[10] J. Hermann, Z. Scha¨tzle, and F. Noe´, “Deep-neural-network solution of the electronic schro¨dingerequation,” Nature Chemistry, vol. 12, no. 10, pp. 891–897, 2020. → pages 2, 37, 38, 40, 48[11] M. Bachmayr, G. Dusson, C. Ortner, and J. Thomas, “Polynomial approximation of symmetric func-tions,” arXiv preprint arXiv:2109.14771, 2023. → pages 2, 42, 46, 47[12] N. W. Ashcroft and N. D. Mermin, Solid State Physics. Saunders College Publishing, 1976. → pages4, 71[13] M. P. Marder, Condensed matter physics. John Wiley & Sons, 2010. → pages 4, 71[14] E. H. Lieb and B. Simon, “The thomas-fermi theory of atoms, molecules and solids,” Advances inmathematics, vol. 23, no. 1, pp. 22–116, 1977. → pages 4, 10, 11, 12, 21[15] R. Lee and N. Drummond, “Ground-state properties of the one-dimensional electron liquid,” PhysicalReview B, vol. 83, no. 24, p. 245114, 2011. → pages 4, 14, 19, 22, 25, 57, 64[16] M. Casula, S. Sorella, and G. Senatore, “Ground state properties of the one-dimensional coulomb gasusing the lattice regularized diffusion monte carlo method,” Physical Review B, vol. 74, no. 24, p.245427, 2006. → pages 4, 14, 23, 57[17] W. Friesen and B. Bergersen, “Dielectric response of a one-dimensional electron gas,” Journal ofPhysics C: Solid State Physics, vol. 13, no. 36, p. 6627, 1980. → page 5[18] S. Bednarek, B. Szafran, T. Chwiej, and J. Adamowski, “Effective interaction for charge carriers con-fined in quasi-one-dimensional nanostructures,” Physical Review B, vol. 68, p. 045328, Jul 2003. →page 5[19] M. J. Evans and J. S. Rosenthal, Probability and statistics: The science of uncertainty. Macmillan,2004. → pages 6, 63[20] R. M. Martin, L. Reining, and D. M. Ceperley, Interacting electrons. Cambridge University Press,2016. → pages 6, 40, 63[21] D. Zhou, H. Chen, C. H. Ho, and C. Ortner, “A multilevel method for many-electron schr\”{o} dingerequations based on the atomic cluster expansion,” arXiv preprint arXiv:2304.04260, 2023. → pages6, 46, 56, 59, 62, 63, 64[22] G. B. Folland, Real analysis: Modern techniques and their applications. John Wiley & Sons, 1999,vol. 40. → page 6[23] L. M. Fraser, W. Foulkes, G. Rajagopal, R. Needs, S. Kenny, and A. Williamson, “Finite-size effectsand coulomb interactions in quantum monte carlo calculations for homogeneous systems with periodicboundary conditions,” Physical Review B, vol. 53, no. 4, p. 1814, 1996. → pages 9, 1175[24] M. P. Tosi, “Cohesion of ionic solids in the born model,” in Solid state physics. Elsevier, 1964, vol. 16,pp. 1–120. → page 10[25] X. Blanc, “A mathematical insight into ab initio simulations of the solid phase,” Mathematical modelsand methods for ab initio quantum chemistry, pp. 133–158, 2000. → pages 12, 13[26] P. P. Ewald, “Die berechnung optischer und elektrostatischer gitterpotentiale,” Annalen der physik, vol.369, no. 3, pp. 253–287, 1921. → page 13[27] S. Johnson and et al, “Specialfunctions.jl.” [Online]. Available: https://specialfunctions.juliamath.org/stable/ → page 13[28] S. W. de Leeuw, J. W. Perram, and E. R. Smith, “Simulation of electrostatic systems in periodicboundary conditions. i. lattice sums and dielectric constants,” Proceedings of the Royal Society ofLondon. A. Mathematical and Physical Sciences, vol. 373, no. 1752, pp. 27–56, 1980. → page 14[29] R. Lee, Application of quantum Monte Carlo methods to excitonic and electronic systems, vol. PhDThesis. → pages 19, 64[30] V. Fock, “Na¨herungsmethode zur lo¨sung des quantenmechanischen mehrko¨rperproblems,” Zeitschriftfu¨r Physik, vol. 61, pp. 126–148, 1930. → pages 33, 60[31] T. Kato, “On the eigenfunctions of many-particle systems in quantum mechanics,” Communicationson Pure and Applied Mathematics, vol. 10, no. 2, pp. 151–177, 1957. → page 34[32] R. T. Pack and W. B. Brown, “Cusp conditions for molecular wavefunctions,” The Journal of ChemicalPhysics, vol. 45, no. 2, pp. 556–559, 1966. → page 35[33] N. D. Drummond, M. D. Towler, and R. J. Needs, “Jastrow correlation factor for atoms,molecules, and solids,” Physical Review B, vol. 70, p. 235119, Dec 2004. [Online]. Available:https://link.aps.org/doi/10.1103/PhysRevB.70.235119 → pages 35, 38, 40, 42, 43, 44, 47[34] G. Ripka, “A practical guide to the hypernetted chain method applied to fermi and bose systems,”Physics Reports, vol. 56, no. 1, pp. 1–64, 1979. → page 37[35] L. Kong, F. A. Bischoff, and E. F. Valeev, “Explicitly correlated r12/f12 methods for electronicstructure,” Chemical Reviews, vol. 112, no. 1, pp. 75–107, 2012, pMID: 22176553. [Online].Available: https://doi.org/10.1021/cr200204r → page 37[36] R. Jastrow, “Many-body problem with strong forces,” Physical Review, vol. 98, pp. 1479–1484, Jun1955. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRev.98.1479 → page 37[37] R. P. Feynman and M. Cohen, “Energy spectrum of the excitations in liquid helium,” Physical Review,vol. 102, pp. 1189–1204, Jun 1956. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRev.102.1189 → pages 37, 4876[38] J. Kim, A. D. Baczewski, T. D. Beaudet, A. Benali, M. C. Bennett, M. A. Berrill, N. S. Blunt,E. J. L. Borda, M. Casula, D. M. Ceperley, S. Chiesa, B. K. Clark, R. C. Clay, K. T. Delaney,M. Dewing, K. P. Esler, H. Hao, O. Heinonen, P. R. C. Kent, J. T. Krogel, I. Kyla¨npa¨a¨, Y. W. Li,M. G. Lopez, Y. Luo, F. D. Malone, R. M. Martin, A. Mathuriya, J. McMinis, C. A. Melton, L. Mitas,M. A. Morales, E. Neuscamman, W. D. Parker, S. D. P. Flores, N. A. Romero, B. M. Rubenstein,J. A. R. Shea, H. Shin, L. Shulenburger, A. F. Tillack, J. P. Townsend, N. M. Tubman, B. V. D.Goetz, J. E. Vincent, D. C. Yang, Y. Yang, S. Zhang, and L. Zhao, “Qmcpack: An open sourceab initio quantum monte carlo package for the electronic structure of atoms, molecules and solids,”Journal of Physics: Condensed Matter, vol. 30, no. 19, p. 195901, apr 2018. [Online]. Available:https://dx.doi.org/10.1088/1361-648X/aab9c3 → page 37[39] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, “Quantum monte carlosimulations of solids,” Rev. Mod. Phys., vol. 73, pp. 33–83, Jan 2001. [Online]. Available:https://link.aps.org/doi/10.1103/RevModPhys.73.33 → page 40[40] P. Lo´pez Rı´os, P. Seth, N. D. Drummond, and R. J. Needs, “Framework for constructing genericjastrow correlation factors,” Physical Review E, vol. 86, p. 036703, Sep 2012. [Online]. Available:https://link.aps.org/doi/10.1103/PhysRevE.86.036703 → pages 40, 41, 44[41] C.-J. Huang, C. J. Umrigar, and M. P. Nightingale, “Accuracy of electronic wave functions in quantumMonte Carlo: The effect of high-order correlations,” The Journal of Chemical Physics, vol. 107, no. 8,pp. 3007–3013, 08 1997. [Online]. Available: https://doi.org/10.1063/1.474658 → page 40[42] R. Drautz and C. Ortner, “Atomic cluster expansion and wave function representations,” arXiv preprintarXiv:2206.11375, 2022. → pages 46, 47[43] P. Lo´pez Rı´os, A. Ma, N. D. Drummond, M. D. Towler, and R. J. Needs, “Inhomogeneous backflowtransformations in quantum monte carlo calculations,” Physical Review E, vol. 74, p. 066701, Dec2006. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevE.74.066701 → page 48[44] R. Drautz, “Atomic cluster expansion of scalar, vectorial, and tensorial properties includingmagnetism and charge transfer,” Physical Review B, vol. 102, p. 024104, Jul 2020. [Online].Available: https://link.aps.org/doi/10.1103/PhysRevB.102.024104 → page 49[45] M. Rinaldi, M. Mrovec, A. Bochkarev, Y. Lysogorskiy, and R. Drautz, “Non-collinear magnetic atomiccluster expansion for iron,” arXiv preprint arXiv:2305.15137, 2023. → page 49[46] E. Wigner, “Effects of the electron interaction on the energy levels of electrons in metals,” Transactionsof the Faraday Society, vol. 34, pp. 678–685, 1938. → page 56[47] E. Mostaani, M. Szyniszewski, C. H. Price, R. Maezono, M. Danovich, R. J. Hunt, N. D. Drummond,and V. I. Fal’ko, “Diffusion quantum monte carlo study of excitonic complexes in two-dimensionaltransition-metal dichalcogenides,” Physical Review B, vol. 96, p. 075431, Aug 2017. [Online].Available: https://link.aps.org/doi/10.1103/PhysRevB.96.07543177[48] A. Girdhar, V. Ashokan, R. O. Sharma, N. D. Drummond, and K. N. Pathak, “Wire-width andelectron-density dependence of the crossover in the peak of the static structure factor from 2kf→ 4kf in one-dimensional paramagnetic electron gases,” Physical Review B, vol. 107, p. 115414,Mar 2023. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevB.107.115414 → pages56, 57, 58, 63, 64[49] M. Hutter, “On representing (anti) symmetric functions,” arXiv preprint arXiv:2007.15298, 2020. →page 56[50] E. Wigner, “On the interaction of electrons in metals,” Physical Review, vol. 46, pp. 1002–1011, Dec1934. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRev.46.1002 → page 60[51] J. Linderberg, “The antisymmetrized geminal power, a simple correlated wave function for chemicalbonding,” Israel Journal of Chemistry, vol. 19, no. 1-4, pp. 93–98, 1980. → page 72[52] H. J. Schulz, “Wigner crystal in one dimension,” Physical Review Letter, vol. 71, pp. 1864–1867, Sep1993. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevLett.71.1864 → page 73[53] H. Bethe and I. Zur Theorie der Metalle, “Eigenwerte und eigenfunktionen der linearen atomkette,”Physik, vol. 71, pp. 205–231, 1931. → page 73[54] H. Schulz, “Fermi liquids and non–fermi liquids,” arXiv preprint, 1995. → page 7378
Thesis/Dissertation
eng
Chemistry
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Ortner, Christoph
Wang, Yan Alexander
Wavefunction parameterizations using the atomic cluster expansion
Text