Fractional charge by fixed-node diffusion Monte Carlo

TL;DR

This paper investigates how fixed-node diffusion Monte Carlo (FNDMC) can accurately reproduce the linear energy behavior with fractional electron numbers, demonstrating its effectiveness in electronic structure calculations involving fractional charges.

Contribution

It shows that FNDMC, when combined with mean-field trial wave functions, can restore the piecewise linearity of energy versus fractional electron number, improving charge localization accuracy.

Findings

FNDMC restores piecewise linearity in $E(N)$ for fractional charges.

FNDMC effectively localizes charge in systems with fractional electrons.

FNDMC outperforms other methods in charged noncovalent systems.

Abstract

Fixed-node diffusion Monte Carlo (FNDMC) is a stochastic quantum many-body method that has a great potential in electronic structure theory. We examine how FNDMC satisfies exact constraints, linearity and derivative discontinuity of total energy $E(N)$ vs. fractional electron number $N$, if combined with mean-field trial wave functions that miss such features. H and Cl atoms with fractional charge reveal that FNDMC is well able to restore the piecewise linearity of $E(N)$. The method uses ensemble and projector ingredients to achieve the correct charge localization. Water-solvated Cl$^-$ complex illustrates superior performance of FNDMC for charged noncovalent systems.

Equations

$$\Psi =\underset{\tau \to \infty }{\lim }\exp [-\tau (\hat{H}-E_T)]{\Psi }_T,$$

$$\hat{H}=\sum _{i=1}^8{\hat{H}}^i,$$

$$\hat{H}{\Psi }_k(M)=E(M){\Psi }_k(M),$$

$${\Psi }_k(M)=\hat{A}[\prod _{i=1}^8{\Psi }_{{\alpha }_i}^i(N_i)]$$

$${\hat{H}}^i{\Psi }_{\alpha }^i(N_i)=E^i(N_i){\Psi }^i{}_{{\alpha }_i}(N_i),$$

$$E(M)=\sum _{i=1}^8{E}^i(N_i)$$

$$\Psi (\mathbf{R},\tau )=\sum _{\vec{\lambda },\vec{\sigma }}{c}_{\vec{\lambda },\vec{\sigma }}|\vec{\lambda },\vec{\sigma }⟩,$$

$$\Psi (\mathbf{R},\tau )=\sum _{j=1}^K\delta [\mathbf{R}-{\mathbf{R}}_j(\tau )],$$

$$\Psi (\mathbf{R},\tau )=\sum _{\vec{\lambda },\vec{\sigma }}\sum _{j=1}^K|\vec{\lambda },\vec{\sigma }⟩⟨\vec{\lambda },\vec{\sigma }|{\delta }_{\mathbf{R},{\mathbf{R}}_j(\tau )}⟩.$$

$$E^i(N)\le \frac{1}{2}[E^i(N+1)+E^i(N-1)],$$

$$E(M)=(8-b)E^1(a)+b{E}^1(a+1),$$

Full text

Motivation: fixed-node diffusion Monte Carlo (FNDMC) is a many-body stochastic quantum Monte Carlo (QMC) projector methodAnderson1975 ; Reynolds1982 ; Umrigar1993 ; Mitas1991 that has been popular in electronic structure theory for its accuracy, scalability and versatilityKolorenc2011rev ; Luchow2011rev ; Austin2012rev ; Morales2014rev ; Dubecky2016rev . FNDMC acts in a continuum position space; for a Hamiltonian $\hat{H}$ it projects out an exact ground state (GS) $\Psi$ that has non-zero overlap with the antisymmetric trial state $\Psi_{T}$, in imaginary time $\tau$:

[TABLE]

where $E_{T}$ is an offset energy that keeps the norm of $\Psi$ asymptotically constant. The total FNDMC energy is an upper bound to the exact energyMoskowitz1982 and the related FN bias goes quadratically with the nodal displacement errorMitas1991 . $\Psi_{T}$ can be well sophisticated so that the related FN bias becomes negligibleMorales2012 ; Clay2015 , nevertheless, practical FNDMC feasible for large systems requires simple and efficient ansätze like the popular Slater-Jastrow wave function Ceperley1977 , $\Psi_{T}=\Psi_{S}J$, where $\Psi_{\mathrm{S}}$ is a single Slater determinant and $J$ is a positive-definite explicit correlation Jastrow termJastrow1955 .

Mean-field theories used to produce $\Psi_{S}$, including DFTPer2012 , as well as the states $\Psi_{S}$ and $\Psi_{T}$ (see below), however miss fundamental constraints on exact electronic structure theory: the total energy $E(N)$ as a continuous function of a particle number $N$ must show piecewise linear relationship with possible derivative discontinuities at integer $N$Perdew1982 ; Yang2000 ; Mori2009 . Lack thereof has been referred to as (de)localization error and it has severe consequences; it leads to an artificial charge adjustment and related spurious energy minimization which causes poor predictivity of charge-transfer/reaction barriers, band gaps, or noncovalent interactions (e.g., solvation of ions)Cohen2008 ; Mori2014rev . In addition, $E(N)$ must also satisfy constancy condition of fractional spin Mori2009 important for correct description of strong correlation.

Here we examine how FNDMC satisfies linearity of $E(N)$ for fractional $N$, and, derivative discontinuity at integer $N$, if combined with $\Psi_{T}$ based on popular spin-restricted mean-field $\Psi_{S}$ that show unphysical convex $E(N)$ dependencyMori2014rev . Such an analysis is important, for instance, to understand if the method is able to provide right energetics of charged noncovalent systemsCohen2008 for right reasons. We consider atoms with fractional charge, H and Cl, to show that FNDMC is well able to restore the piecewise linearity of $E(N)$ from the states that do not possess such a property. Walker population analyzes indicate that ensemble and projector features are both operative in this achievement. Insights gained from water-solvated Cl- complex suggest that accuracy and robustness of FNDMC in charge-involving noncovalent interactions relate to the accurate charge localization.

Theory: the effective fractional charge per atom was achieved for both, H and Cl, by modeling cube (system) composed of 8 atoms (subsystems) of the same type separated enough (infinite-separation limit) so that the interactions between the subsystems can be neglected. The Hamiltonian of such a system with $M=\sum_{i=1}^{8}N_{i}$ electrons,

[TABLE]

satisfies the Schrödinger equation,

[TABLE]

where

[TABLE]

is a $k$-th $g(M)$-times degenerate antisymmetric GS wave function $[g(M)=\prod_{i=1}^{3}g_{i}(N_{i})]$ assembled as an antisymmetrized product of individual subsystem $g_{i}(N_{i})$-times degenerate wave functions $\Psi^{i}_{\alpha_{i}}(N_{i})$ that each satisfies

[TABLE]

$\hat{A}$ is an $M$-electron antisymmetrization operator, and,

[TABLE]

is a corresponding degenerate total energy of the system. Addition of an integer charge $q$ to such a system adds effective charge $q/8$ per subsystem thus enabling fractional charge as well.

QMC simulations of systems in large separation limit require additional considerations for the lack of VMC thermalization within usual simulation times (see Methods), and, non-ergodicity (electron localization on specific subsystems) due to the branching term in FNDMC that eliminates configurations with small weights, including those with electrons possibly attempting to diffuse between the distant subsystems. We argue that the total energy expectation values $\langle E\rangle$ from such FNDMC simulations are energetically correct.

One needs to realize that within a properly initialized non-ergodic equilibrium FNDMC simulation (starting from $\Psi_{T}^{2}$), any walker $\mathbf{R}_{j}=(\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{r}_{3},\dots,\mathbf{r}_{N})_{j}$, irrespective of its specific localization $\{N_{i}\}_{j}$, i.e. a specific set of occupation numbers of the subsystems labeled by $i$, randomly samples one of the degenerate symmetry-broken GS wave functions $\Psi_{k}$ and contributes to $\langle E\rangle$ according to the so-called mixed distribution $\Psi_{k}\Psi_{T}$. For a totally symmetric $\Psi_{T}$, like, e.g., our $\Psi_{T}$ based on a spin-restricted Slater determinant, that satisfies $\langle\Psi_{T}|\Psi_{k}\rangle\neq 0$ for all $k$, it is not hard to accept that the population-based $\langle E\rangle$, that is of our sole interest here, is not affected by the symmetry breaking.

Hydrogen cube: first, we present the results obtained for H8 cubeMori2014rev with an edge of $a=$1000 Å, 0-16 electrons and singlet spin multiplicity ($N_{\alpha}=N_{\beta}$), allowing us to study essentially isolated H atom with an effective fractional charge in range $0\leq N\leq 2$, i.e., along H${}^{+}\rightarrow$H$\rightarrow$H*-* pathway and 0.25$e$ increments (Fig. 1). In agreement with the previous fractional-charge studies of HMori2009 ; Mori2014 ; Mori2014rev , we find that the considered mean-field methods, HF and DFT with PBEPerdew1996 functional, miss linearity and derivative discontinuity of $E(N)$ and show unphysical convex behavior instead (Fig. 1, red, blue).

Interestingly, FNDMC with related $\Psi_{T}$ (cf. Fig. 3) recovered piecewise linear $E(N)$ dependency within the statistical resolution (Fig. 1, black). The results were indistinguishable for HF and PBE (more $\Psi_{T}$ types were thus not considered). FNDMC with $\Psi_{T}$ based on PBE orbitals produced electron affinity EA(H)=0.713(1) eV to be compared vs. experimental value of 0.754 eVLykke1991 . Slight discrepancy of $\sim$0.04 eV is attributed to the use of effective core potential (ECP) and residual FN bias due to 1-determinantGasperich2017 . Note that if FNDMC method satisfies linearity condition, the quality of the nodal surface determines only the slope of the actual $E(N)$ segment.

Chlorine cube: next we consider Cl8 cube ($a=$1000 Å) in order to examine how FNDMC describes fractional charge in many-electron system. We consider only singlet states ($N_{\alpha}=N_{\beta}$). Effective fractional charge per Cl atom ranges between $16\leq N\leq 18$ (1Cl${}^{+}\rightarrow^{2}$Cl$\rightarrow^{1}$Cl*-*) with 0.25$e$ increments (Fig. 2). We observe that although the mean-field methods (HF, PBE, HSE06Heyd2006 ) miss linearity, and derivative discontinuityCohen2008 at point $N=17$ per atom (neutral Cl), FNDMC well recovers both features from the related $\Psi_{T}$ and produces accurate EA(Cl)=3.64(2) eV consistent with experiment (3.613 eVBerzinsh1995 ). It appears that FNDMC is well able to produce correct energetics in fractional-charge systems.

Wave function analysis: let us now focus on how FNDMC achieves linearity of $E(N)$ in fractional charge models. As mentioned above, the electrons within each walker localize on specific centers, leading to the correspondence $\mathbf{R}_{j}\leftrightarrow\{N_{i}\}_{j}$ between the explicit (alive) walker position (varying along simulation) and subsystem occupation vector in occupation representation (fixed during the non-ergodic simulation unless the walker dies out). Since each walker randomly samples one of the degenerate GS wave functions $\Psi_{k}$, and the subsystem occupations randomly vary between the walkers, we find it convenient to introduce irreducible subsystem occupation vector $\vec{\lambda}_{j}$ obtained as a descending-ordered set of subsystem occupations $\{N_{i}\}_{j}$, $\vec{\lambda}_{j}=(\mathrm{max}\{N_{i}\}_{j},\dots,\mathrm{min}\{N_{i}\}_{j})$. In addition, we introduce a corresponding vector of total subsystem spin momenta $\vec{\sigma}_{j}$ ordered in a descending way for each subset of identical $\lambda_{k}$. The pair $(\lambda_{k},\sigma_{k})_{j}$, labels the state of a subsystem $k$ (after reordering) with $\lambda_{k}$ electrons and a total spin momentum $\sigma_{k}$, as sampled by the walker $j$. We arrive at the more convenient correspondence, $\mathbf{R}_{j}\leftrightarrow|\vec{\lambda},\vec{\sigma}\rangle_{j}$, that enables equal-footing comparisons between the walkers irrespective of the specific localization (and local position), and, unambiguous analysis of the system state within FNDMC simulation. We are interested in coefficients $\{c_{\vec{\lambda},\vec{\sigma}}\}$ of an expansion of irreducible GS

[TABLE]

subject to $\sum_{\vec{\lambda},\vec{\sigma}}|c_{\vec{\lambda},\vec{\sigma}}|^{2}=1$, that could be obtained from a stochastic realization of $\Psi(\mathbf{R},\tau)$ sampled by $K$ walkers within equilibrium QMC simulation,

[TABLE]

by virtue of the projectors,

[TABLE]

For $M$ electrons in the system with a specific partition $\vec{\lambda}_{j}$ between the subsystems, the GS consistent with the subsystem total energy convexity conditionPerdew1982 ,

[TABLE]

requires that each distant subsystem contains $a=\mathrm{int}(M/8)$ electrons, and, for nonzero $b=M-8a$, $b$ subsystems contain 1 additional electron each. $\vec{\sigma}$ is such that it minimizes the total energy. For $M=10$, any 2 of the 8 subsystems must contain 2 electrons and the remaining subsystems 1 electron each, $\vec{\lambda}=(2,2,1,1,1,1,1,1)$, assuming a spin GS, e.g., $\vec{\sigma}=(0,0,0.5,0.5,0.5,-0.5,-0.5,-0.5)$ for a singlet H8 cube. Other partitions $\vec{\lambda}$ (or non-GS $\vec{\sigma}$) represent excited-state configurations, e.g., $\vec{\lambda}=(2,2,2,1,1,1,1,0)$, or, $\vec{\lambda}=(3,1,1,1,1,1,1,1)$. Therefore, for the special case of distant subsystems, an exact GS may be conveniently identified by 1-term eq. 7 expansions.

In Fig. 3, we compare the total energies $E(N)$ obtained for H cube by VMC sampling of $|\Psi_{\mathrm{S}}^{2}|$ and $|\Psi_{T}^{2}|$, and FNDMC using $\Psi_{T}$. Clearly, contrary to the piecewise linear FNDMC, the VMC total energies show convex behavior, although in case of SJ, it seems that the derivative discontinuity starts to develop due to the use of Jastrow termGoetz2017 . For $\Psi_{\mathrm{S}}$, the $E(N)$ curve is indeed consistent with HF in Fig. 1 within the statistical error.

Analysis of the related walker ensembles in terms of $|\vec{\lambda},\vec{\sigma}\rangle$ configurations revealed multi-configuration states containing large fractions of excited-states for VMC, and, only GS configurations for FNDMC equilibrium ensembles, for all $N$, orbital types, and Cl cube as well (cf. Fig. 5, Appendix I). For instance, $\Psi_{S}$/$\Psi_{T}$ of PBE orbitals of neutral H cube with $N=1$ electrons contained only 6%/11% of GS configurations and 6/5 types of configurations, whereas $\Psi$ sampled by FNDMC using $\Psi_{T}$ contained GS configurations only. For $N=1.25$, $\Psi_{S}$/$\Psi_{T}$ contained 22%/38% of GS configurations and 6/4 configuration types, whereas, again, FNDMC contained 100% of GS. Details of such analyzes for H and Cl cubes, DFT and HF orbitals, VMC and FNDMC, and all considered $N$, are reported in Appendix III.

The observation that FNDMC produces 100% of GS configurations in our fractional-charge models is very important as it may improve our understanding of the method robustness in chargedWineman-Fisher2019 and charge-transfer systems (e.g., reactionsSwann2017 ; Krongchon2017 , excited statesValsson2010 ; Dubecky2011 ). It has been shownYang2000 that the linearity of $E(N)$ for fractional $N$, and thus accurate charge localization, directly relates to the subsystem GS occupations. Namely, if the total energy of the system can be expressed, using $E^{i}(\lambda_{i})=E^{j}(\lambda_{i})$, as

[TABLE]

where $a$ and $b$ were defined above, then $E(N)=E(M)/8$ for $N=M/8$ is a piecewise linear function. The success of FNDMC in recovering the $E(N)$ piecewise linearity is thus attributed to the ensemble nature of the method causing that each electron resides on a specific subsystem at a time (fractional density is only possible via ensemble), and, projector property that apparently projects out proper (GS) configurations.

Cl-(H2O)6: Finally, we consider a practical application of the concepts discussed above for noncovalent interactions of charged systemsWineman-Fisher2019 . It was found that FNDMC provides correct $E(N)$ dependency for Cl cube contrary to mean-field (Fig. 2). In systems with a charged Cl atom, one would thus generally expect that FNDMC would be more robust than, say, DFT, and insensitive to the $\Psi_{T}$ orbitals. We have considered interaction energiesDubecky2016rev ($\Delta E$) of Cl-(H2O)6 complex by various DFT approximations and FNDMC using the same DFT functionals. Fig. 4 summarizes the results and illustrates superior robustness of FNDMC that reproduces the reference CCSD(T)/complete-basis-set results to within 1% for all $\Psi_{T}$. In contrast, $\Delta E$ from DFT and HF altogether varies by more than 20% due to the known (de)localization error.

This is now easy to understand. Randomly chosen DFT determines fractional charge $N$ from a certain interval according to its intrinsic delocalization error. In addition, it produces energy on a convex curve with its intrinsic curvature, that varies between functionalsCohen2006 . As a consequence, for a large set of DFT approximations, one obtains a large distribution of biases in the total energy of a complex, and, consequently in $\Delta E$, as observed. On the other hand, if FNDMC predicts a correct charge localization, as expected from the results reported above, the fractional charge $N$ in subsystems of a complex is accurate and well defined. Inaccuracy of $\Psi_{T}$ dictates only the slope of the linear dependency, which leaves only a limited interval of bias that can be realized with varying (reasonable) $\Psi_{T}$. The dispersion of the results is thus much smaller, or, in other words, the method is much more robust, as observed, and we attribute its success to the accurate charge localization.

Summary: Fractional-charge computations of H and Cl atoms revealed that FNDMC is well able to restore the correct piecewise linear $E(N)$ behavior from $\Psi_{T}$ that does not have such a property. Insights from the walker population analyses, in combination with formal discussions, indicate that FNDMC correctly describes $E(N)$ linearity thanks to its ensemble nature and projection property. Inaccuracy of $\Psi_{T}$ is expected to determine only the slope of $E(N)$ segments. This finding enables understanding of FNDMC robustness, as demonstrated in noncovalent Cl-(H2O)6 complex, from a new perspective. Study of fractional spin by FNDMC is underway.

Methods: Naive simulations of cube models with an edge of $a=1000$ Å would not be possible because the VMC equilibration times would vastly overcome the maximum simulation times available. In order to overcome this issue in an automated fashion (avoiding manual initialization of configurations that may be biased by user imagination), the walkers for large-cube FNDMC simulations were initialized by ergodic VMC thermalization and stepwise quasi-adiabatic extension of the system size starting from small $a$. We verified that the walker distributions produced in this way sample $\Psi_{T}^{2}$, and, in case of a Slater determinant with Hartree-Fock (HF) orbitals, VMC reproduced the HF total energy for both, H and Cl, and all considered $N$. This was not the case for walkers initialized directly in a large cube.

FNDMC computations used QMCPACKKim2018 code with imaginary time step of 0.005 a.u. and T-movesCasula2006 . Nuclei were represented by Burkatzki-Filippi-Dolg effective core potentialsBurkatzki2008 . Single-determinant Slater-Jastrow $\Psi_{T}$ were based on orbitals expanded in 1s-augmented VTZ one-particle basis sets without highest angular momentum channelsDubecky2017b obtained with a tight SCF convergence (Gaussian G09G09 ). Jastrows containing up to electron-electron-nucleus terms, with a cutoff radius of 10 Bohr, were optimized by the linear methodToulouse2007 , as usual, but cube simulations reused Jastrows optimized for H- and Cl, respectively. The structure of Cl-(H2O)6 (reported in Appendix II) was optimized at the MP2/aug-VTZ level (G09).

Discussions with A. Gendiar have been fruitful. Financial support by Czech Science Foundation (18-24321Y, 18-25128S), University of Ostrava (IRP201826), and ERDF (ITMS 26220220179), is gratefully acknowledged. The computations used IT4Innovations National Supercomputing Center (LM2015070).

Appendix

I Chlorine cube

II Structure of solvated chlorine-anion complex

The structure of the considered Cl-(H2O)6 complex optimized at the MP2/aug-VTZ level in G09 (xyz-format):

   19

Cl 0.000000 0.000000 0.000000 O 3.290138 0.000000 0.000000 O -3.290138 0.000000 0.000000 O 0.000000 3.264842 0.000000 O 0.000000 -3.264842 0.000000 O 0.000000 0.000000 3.353781 O 0.000000 0.000000 -3.353781 H 2.664747 0.734646 0.000000 H 2.664747 -0.734646 0.000000 H -2.664747 0.734646 0.000000 H -2.664747 -0.734646 0.000000 H 0.000000 2.636302 0.731619 H 0.000000 2.636302 -0.731619 H 0.000000 -2.636302 0.731619 H 0.000000 -2.636302 -0.731619 H 0.000000 0.732695 2.727551 H 0.000000 -0.732695 2.727551 H 0.000000 0.732695 -2.727551 H 0.000000 -0.732695 -2.727551

III Analysis of walkers

Walker ensemble analyzes for H and Cl cube models, various methods (VMC, FNDMC), and, trial states (Slater determinant $\Psi_{S}$, Slater-Jastrow $\Psi_{T}$). For each $N$ (in electrons per atom), reported are %, $\vec{\lambda}$ and $\vec{\sigma}$. GS configurations are indicated by the bold font.