Neural network variational Monte Carlo for positronic chemistry
Quantum chemical calculations of the ground-state properties of
positron-molecule complexes are challenging. The main difficulty lies in
employing an appropriate basis set for representing the coalescence between
electrons and a positron. Here, we tackle this problem with the recently
developed Fermionic neural network (FermiNet) wavefunction, which does not
depend on a basis set. We find that FermiNet produces highly accurate, in some
cases state-of-the-art, ground-state energies across a range of atoms and small
molecules with a wide variety of qualitatively distinct positron binding
characteristics. We calculate the binding energy of the challenging non-polar
benzene molecule, finding good agreement with the experimental value, and
obtain annihilation rates which compare favourably with those obtained with
explicitly correlated Gaussian wavefunctions. Our results demonstrate a generic
advantage of neural network wavefunction-based methods and broaden their
applicability to systems beyond the standard molecular Hamiltonian.
G. Cassella, W. M. C. Foulkes, D. Pfau, J. S. Spencer
# Neural network variational Monte Carlo for positronic chemistry

###### Abstract

Gamma rays produced by positron annihilation are used as a sensitive probe of matter at atomic length scales. Technologies for manipulating positrons and studying their interactions with ordinary matter are rapidly progressing. This motivates the development of accurate _ab initio_ methods for modelling positronic interactions with molecular matter. Here, we apply the recently developed Fermionic neural network (FermiNet) wavefunction ansatz to the problem of finding the ground-state properties of mixed positron-electron systems. We find that FermiNet produces highly accurate, in some cases state-of-the-art, ground-state energies across a range of atoms and small molecules with a wide range of qualitatively distinct positron binding characteristics. We highlight the capabilities of our method by calculating the positron binding energy of the challenging non-polar Benzene molecule.

Since the existence of the positron - the positively charged anti-particle of the electron - was first postulated by Dirac [1], many predictions have been made concerning the formation of bound states between positrons and ordinary matter. The formation of bound electron-positron states results in greatly enhanced annihilation rates. This process has recently been successfully exploited in positron annihilation spectroscopy experiments [2; 3; 4; 5; 6; 7; 8; 9], where the energy-dependent annihilation rate enables sensitive measurements of the positron binding energy [10].

Enhanced positron annihilation _in materio_ has valuable applications in medical physics [11], materials science [12; 13], and astrophysics [14]. Experimental apparatus for trapping large numbers of positrons continues to grow in sophistication, offering a glimpse into a world of exotic antimatter chemistry [15]. These advances motivate the development of a strong _ab initio_ description of positronic bound states to accelerate the continued development of new antimatter-based technology.

This problem has been addressed using many standard computational chemistry tools. Positron binding to atoms and molecules has been studied using wavefunction expansions in explicitly correlated Gaussians [16; 17; 18; 19] (ECG), configuration interaction (CI) methods [20; 21; 22; 23], and quantum Monte Carlo (QMC) methods [24; 25; 26; 27; 28; 29; 30; 31]. The annihilation rate and lifetime of positrons in solids, particularly in the presence of defects, has been studied using density functional theory [32; 33] and quantum Monte Carlo methods [34; 35].

Despite the intense theoretical interest, describing the positronic wavefunction remains challenging for several reasons. Due to the repulsive potential between nuclei and positrons, expanding the positronic wavefunction in a basis often requires including diffuse basis functions and basis functions with large angular momenta. In many cases, it is not appropriate to consider the positronic wavefunction to be centred on the nuclei of a molecule, introducing additional difficulties in choosing appropriately centred basis functions. These difficulties have historically resulted in very slow convergence of CI calculations and limited accuracy for QMC methods which utilise wavefunction ansatze based on Hartree-Fock orbitals.

The most successful description of positron binding to molecules comes from a recent work that develops a many-body theory of positron binding to molecules and produces positron binding energies in close agreement with experimental measurements [36]. In their work, Hofierka et al. highlight the shortcomings of QMC methods applied to positron binding - particularly, the lack of calculations for large non-polar molecules, which constitute the majority of experimentally studied systems. Here, we address this shortcoming.

We propose a new approach to calculating the ground state properties of molecular positronic bound states, based on recently developed neural network wavefunction ansatze for QMC [37]. The Fermionic neural network (FermiNet) models the many-body wavefunction without referencing a set of basis functions. This conveniently sidesteps a number of aforementioned difficulties in describing positronic wavefunctions. We extend FermiNet to represent the positronic component of the wavefunction on an even footing with the electronic component. With a minimal alteration to the neural network architecture, we obtain a flexible and accurate ansatz for mixed electron-positron wavefunctions. We present results for a range of systems with qualitatively distinct mechanisms for positron binding and obtain state-of-the-art accuracy for the ground-state energy in these systems. Finally, we show that our method obtains the positron binding energy for Benzene in close agreement with the experimental value and the many-body theory of Hofierka et al [36].

## I Results

We benchmark the accuracy of FermiNet for a series of well-studied positronic systems. These are presented here in (approximate) order of increasing complexity. Unless otherwise specified, all results presented herein were obtained using the network architecture and training protocol detailed by Spencer [38]. We pre-train the electron orbitals to the Hartree-Fock solution of the bare molecule and do not pre-train the positron orbitals. Errors in energy expectation values are evaluated using a reblocking approach [39] to account for sequential correlations in the Metropolis-Hastings sampling.
_Positronium hydride_ - The positronium (Ps) atom (consisting of a bound electron and positron) and hydrogen form a stable molecule. Here we work within the Born-Oppenhimer approximation, neglecting the proton's motion. Near exact ECG calculations are available for this system [16], yielding a ground-state energy of \(E_{0}=-0.7891794\) Hartrees. We obtain a ground state energy of \(E_{0}=-0.789144(3)\) Hartrees.

_Sodium and magnesium atoms_ - The first ionisation energy of sodium, 0.1886 Hartrees, is smaller than the binding energy of the Ps atom, 0.25 Hartrees. The positronic sodium atom is then more accurately described as a bound complex of a positronium atom and a sodium cation. The binding energy of this complex is calculated as \(\epsilon=E\left([\mathrm{Na}^{+},\mathrm{Ps}]\right)-E\left(\mathrm{Na}^{+} \right)-0.25\). We fail to predict binding without variance matching, obtaining \(\epsilon=-0.37\) milliHartrees. Utilizing a variance matching procedure, we obtain \(\epsilon=0.32(15)\) milliHartrees, predicting binding in agreement with previous FCSVM (\(\epsilon_{\mathrm{FCSVM}}=0.4\) milliHartrees) and CI (\(\epsilon_{\mathrm{CI}}=0.21\) milliHartrees) calculations.

We obtain a positron affinity of 0.01618(9) Hartrees for the magnesium atom. This agrees with previous FCSVM (0.015612 Hartrees) and CI (0.01615 Hartrees) results.

_Lithium hydride and beryllium oxide_ - We have calculated the ground-state energy of LiH and its positronic complex for a range of interatomic separations (see Data Tables in Supplementary Material). Fitting these potential energy surfaces using Nesterov's algorithm, implemented in the MOLCAS package, we obtain equilibrium bond distances of 3.0196 Bohr, with ground-state energy -8.07050(1) Hartrees for LiH, and 3.371 Bohr with a ground state energy of -8.10774(1) Hartrees for [LiH, e\({}^{+}\)]. These deviate very slightly from the widely accepted literature values of 3.015 Bohr (LiH) and 3.348 Bohr ([LiH, e\({}^{+}\)]). We have calculated ground-state energies at the canonical separations for comparison with previous results, shown in Table 1. Our calculations are not only in excellent agreement with previous work but are seen to yield the most accurate variational result for [LiH, e\({}^{+}\)] reported to date.

We have calculated the ground-state wavefunction of BeO over a range of interatomic separations - from below the equilibrium separation to dissociation. We fit the potential energy surfaces using MOLCAS and obtain an equilibrium bond distance of 2.515 Bohr with a ground-state energy of -89.90572(4) Hartrees for BeO, and 2.530 Bohr with a ground-state energy of -89.93082(3) Hartrees for [BeO, e\({}^{+}\)]. We compare our results against previous calculations in Table 1 and find that FermiNet yields the most accurate variational result for [BeO, e\({}^{+}\)].

The electronic ground state of BeO transitions from the spin-singlet configuration at the equilibrium interatomic separation to the spin-triplet configuration at the dissociative limit. We enforce the appropriate ground-state spin configuration by choosing \(S_{z}=1\) ( as FermiNet is a spin-assigned wavefunction) at wide interatomic separations. The resulting potential energy surfaces are plotted in Fig. 1. At \(\sim 4\) Bohr, the electronic ground-state transitions between the spin-singlet and spin-triplet, causing a sharp change in the ground-state dipole moment and the resulting positron binding energy with the molecular ground-state, which vanishes almost completely. At \(\sim 6\) Bohr, we see a smooth transition between two qualitatively distinct binding modes between the molecule and the positron - binding to the molecular dipole field below this separation and binding exclusively to the 'lone' beryllium atom beyond this separation. This is readily seen by visualizing the ground-state positron density on either side of the maximum, as shown in Fig. 1 for the triplet state, and corroborated by the dipole moment falling below the critical threshold for binding (\(\sim 1.625\) Debye) at the transition.

_Dilithium_ - For the dilithium molecule, we obtain a vertical positron affinity of 66.75(2) milliHartrees at an equilibrium bond length of 5.015 Bohr. The molecular binding energy of Li\({}_{2}\) from our calculations is \(\sim 38\) milliHartrees, and the positron affinity of a lone lithium atom is known from the literature to be \(\sim 2\) milliHartrees, meaning this system is very stable against the \([\mathrm{Li}_{2},\mathrm{e}^{+}]\rightarrow[\mathrm{Li},\mathrm{e}^{+}]+ \mathrm{Li}\) dissociation channel. The one-particle density of the positronic ground-state of dilithium is shown in Fig. 2(a).

_Benzene_ - The experimental positron binding energy of Benzene is 5.51 milliHartrees [10]. Utilizing our variance-matching procedure (described in detail for Benzene in the Supplementary Material), we obtain a finite positron binding energy of 4.1(3) milliHartrees. This falls in very close agreement with the binding energy obtained by Hofierka et al. of \(\sim\)4.26 milliHartrees [36]. The one-particle density of the positronic ground-state of benzene is shown in Fig. 2(b).

## Discussion

The selection of systems presented here spans a broad range of positron binding phenomena: positronium formation, binding with an induced atomic dipole moment, binding with a static molecular dipole moment, and binding due to correlations with covalent bonding electrons in molecules. In all cases where benchmarks are available, FermiNet VMC produces excellent and, in some cases, state-of-the-art results for the positron affinity. From the density plots presented, it is clear that this performance is consistent between wavefunctions with very different qualitative characteristics. An identical ansatz is used for every system studied - we have not employed any fine-tuning or system-specific treatments. An important aspect of our approach is that, due to being a basis-set-free method, we do not provide any information about the nature of the positron binding in the input to the calculation (e.g. via the placement of basis set functions). Rather, the location and nature of the positron binding emerges naturally during optimization of the wavefunction. The high level of accuracy achieved across various systems shows that FermiNet offers a flexible and accurate ansatz for mixed electron-positron wavefunctions.

To date, quantum Monte Carlo calculations of positron binding have focused on small, polar molecules. As pointed out by Hofierka et al. [36], accurate quantum Monte Carlo results for large non-polar organics, which comprise the majority of experimentally relevant systems, are lacking. We believe that the present work amends this gap. Our results for positronium hydride, sodium and magnesium atoms,
Figure 1: **(a)** Ground-state positronic density, projected into the molecular plane, of a positron attached to a beryllium oxide molecule over a range of interatomic distances, accumulated via MCMC sampling. Image scale has been normalized by bondlength. The scale bar in the bottom left of each panel indicates the relative size of one Bohr radius, \(a_{0}\). **(b)** Energy (of the bare molecule and positronic complex), dipole moment (of the bare molecule), and positron affinity of beryllium oxide over a range of interatomic distances. Green (black) markers in the dipole moment and positron affinity plots indicate values accumulated for the electronic singlet (triplet) projected wavefunction. Solid lines and dots in the dipole moment plot indicate the corresponding values obtained by Buenker et al. [23]. Horizontal red lines in the dipole moment and positron affinity plots indicate the mean-field critical value for positron binding to the molecular dipole field and the positron affinity of a lone beryllium atom, respectively.
and small diatomic molecules demonstrate that our approach can achieve state-of-the-art accuracy compared with previous work. Further, our results for the non-polar dilithium and benzene molecules demonstrate that this accuracy is retained when describing modes of positron binding governed entirely by strong electron-positron correlation effects. The results in Fig. 2 offer an intuitive understanding of the binding mechanism between non-polar molecules and positrons: correlation-dominated binding is facilitated by the presence of a centre of increased electronic density away from the atomic nuclei of a molecule. In dilithium, this is the covalent bond, and in benzene, this is the increased electronic density in the centre of the molecule from the delocalisation of the \(\pi\)-bonds in the ring.

Our result for the positron binding energy of benzene fall within 'chemical accuracy' (\(\sim\) 1.6 milliHartrees) of the experimental value, but this level of accuracy is insufficient for other species of experimental interest. Many chemical species possess positron binding energies below chemical accuracy, and obtaining this level of accuracy represents a significant challenge to computational chemistry methodsis a signi. We believe that the necessary improvement in accuracy, and obviation of our need to utilise variance matching, can be achieved by adopting the recently introduced PsiFormer architecture [40], an improvement upon the FermiNet architecture utilising Transformer networks [41]. The PsiFormer obtains more

\begin{table}
\begin{tabular}{l c c c c c c}  & \multicolumn{3}{c}{LiH} & \multicolumn{3}{c}{BeO} \\ \cline{2-7} Method & Bare & Positronic & Binding energy & Bare & Positronic & Binding energy \\ \hline FermiNet-VMC & -8.07051(1) & **-8.10775(1)** & 0.03723(2) & **-89.90572(4)** & **-89.93802(3)** & 0.02510(7) \\ SJ-VMC & -8.06307(3) [25] & -8.08034(4) [25] & 0.01727(7) & -89.3173(25) [30] & -89.3365(13) [30] & 0.0192(38) \\ SJ-FN-DMC & -8.0704(1) [25] & -8.10718(11) [25] & 0.03678(2) & -89.7854(13) [30] & -89.8134(12) [30] & 0.02800(25) \\ CISD & -8.03830 [20] & -8.05530 [20] & 0.017 & - & - & - \\ MRD-CI & -8.06827 [21] & -8.09764 [21] & 0.02937 & -89.759352 [23] & -89.773133 [23] & 0.013781 \\ ECG & **-8.07054**[18] & -8.10747 [18] & 0.03693 & - & - & - \\ GW & - & - & 0.039(1) [36] & - & - & - \\ \end{tabular}
\end{table}
Table 1: Ground state energy of LiH and BeO, and their positronic complexes, obtained via various computational methods at equilibrium bond-length. Statistical errors are omitted where they are smaller than the reported precision or otherwise omitted in the referenced source. The lowest variational energy in each column is bolded.

Figure 2: Orthographic projections of the ground-state one-particle density for positronic **(a)** dlithium and **(b)** benzene molecules. Left/right columns show the electronic/positronic density, \(\rho_{e}\)-/\(\rho_{e^{+}}\). In the dilithium molecule, the positron is seen to be strongly localized to a torus wrapping the covalent bond. In the benzene molecule, the positronic orbital is extremely diffuse, resulting in a much noisier Monte Carlo estimate of the one-particle density. The positronic orbital resembles a \(p\)-orbital ‘sandwiching’ the aromatic ring. The scale bar in the bottom left of each subfigure indicates the relative size of one Bohr radius, \(a_{0}\).
accurate total energies for the ground-state energy of the bare benzene molecule than FermiNet. Investigating the accuracy of the PsiFormer applied to positron binding energy calculations is a promising avenue for future work.

As with previous work utilizing FermiNet, computational scaling remains an issue. Calculations involving \(\gtrsim 50\) particles are too expensive for presently available computational resources, prohibiting the application of the method to large molecules. One approach to calculating the ground-state wavefunction of large molecules is utilizing pseudopotentials to remove core electrons from the calculations. The positronic component of the wavefunction is often small near the atomic nuclei, so we would expect the correlation effect between the core electrons and the positron to be small. Pseudopotentials have been successfully employed in previous QMC calculations of positron lifetimes in solids [35]. As the architecture employed herein only extends the original FermiNet architecture to treat different particle types separately, we expect that future advances in the computational efficiency of FermiNet and related neural network wavefunctions will be able to be readily combined with our approach to model positron binding.

## Conclusion

We have shown that the FermiNet ansatz for VMC calculations can be extended to include positrons naturally, treating positrons on an equal footing to electrons. Because this ansatz does not depend on a basis set, our treatment sidesteps traditional methods' issues in selecting and converging an appropriate basis set for describing positronic wavefunctions. Our method produces highly accurate results for several molecules with various binding mechanisms for positrons without any system-specific tuning. We expect that the simplicity of this method will lend itself to many challenging applications beyond those presented here, e.g. calculations involving multiple positrons. With additional computational effort, this method can provide accurate predictions for positron annihilation experiments.

## Methods

We find the ground-state wavefunction and corresponding ground-state energy of the many-body Coulomb Hamiltonian in the Born-Oppenheimer clamped nuclei approximation,

\[\mathcal{H}=-\frac{1}{2}\sum_{i}\nabla_{i}^{2}+\sum_{i,j}\frac{q_{i}Z_{j}}{ |\mathbf{r}_{i}-\mathbf{R}_{j}|}+\sum_{i>j}\frac{q_{i}q_{j}}{|\mathbf{r}_{i}- \mathbf{r}_{j}|}+\sum_{i>j}\frac{Z_{i}Z_{j}}{|\mathbf{R}_{i}-\mathbf{R}_{j}|}, \tag{1}\]

where \(q_{i},\mathbf{r}_{i}\) are the particle charges and positions, and \(Z_{i},\mathbf{R}_{i}\) are the charges and positions of fixed nuclei. Here, and throughout, we utilise Hartree atomic units: \(\hbar=e=m_{e}=1\). For the mixed electron-positron systems considered here, \(q_{i}=\pm 1\). We solve for the many-body ground-state wavefunction using the variational Monte Carlo (VMC) method: [42] a many-body wavefunction \(\Psi_{\theta}\), parameterized by \(\theta\), is continuously updated via a gradient descent procedure to minimize the energy expectation value,

\[\left\langle E\right\rangle_{\theta}=\frac{\int\Psi_{\theta}^{*}(\mathbf{r}) \mathcal{H}\Psi_{\theta}(\mathbf{r})d\mathbf{r}}{\int\Psi_{\theta}^{*}(\mathbf{r})\Psi_{ \theta}(\mathbf{r})d\mathbf{r}}, \tag{2}\]

where \(\mathbf{r}=(\mathbf{r}_{1},\dots,\mathbf{r}_{N})\). This integral, and its gradient with respect to \(\theta\), are evaluated via Monte Carlo integration.

The FermiNet represents the many-body wavefunction as a sum of block-diagonal determinants containing many-particle orbitals which depend upon the coordinates of all particles in a permutation equivariant manner. This is written

\[\Psi(\mathbf{r})=\sum_{k}^{n_{\text{fold}}}\prod_{\chi}\text{det}\left[\psi_{i}^{ \chi\chi}(\mathbf{r}_{j}^{\chi};\{\mathbf{r}_{j}^{\chi}\};\{\mathbf{r}_{j}^{\chi}\})\right], \tag{3}\]

where the set \(\{\mathbf{r}_{/j}\}\) includes all particle coordinates except \(\mathbf{r}_{j}\), and \(\chi=(\sigma,q)\) labels species of particles which are distinguished by their spin \(\sigma\in(\uparrow,\downarrow)\) and charge \(q\in(+,-)\). Here we have made a slight abuse of notation for the sake of brevity: permutation invariance for the set \(\{\mathbf{r}_{j}^{\chi}\}\) is only maintained between particles of the same species. We emphasise that these are not the dense determinants discussed in recent works extending FermiNet [43], with the exception of the Benzene calculation for which dense determinants were used for the electronic component of the wavefunction. The many-particle orbitals \(\psi_{i}\) are represented by a deep neural network [37] (architecture described in the Supplementary Material). Multiplicative coefficients are omitted from the sum as they are trivially absorbed into the orbitals. Gradient descent is performed via the Kronecker-factored approximate curvature (KFAC) algorithm, an approximation of natural gradient descent [44] which scales well to large neural networks [45]. Natural gradient descent is closely related to the stochastic reconfiguration method, well-known in the quantum Monte Carlo literature [46]. The present work introduces two alterations to the original FermiNet architecture. Firstly, we have included positronic orbitals as additional species, i.e. additional blocks in the determinant. Secondly, we utilise distinct weights in the neural network layers for every species (unlike the original FermiNet, where spin up and down electronic orbitals shared weights).

A single determinant of the form in Eq. (3) can represent any fermionic many-body wavefunction [47]. In practice, the argument for the universality of FermiNet determinants depends upon the representation of discontinuous functions which cannot be constructed using realistic neural networks. Despite this, FermiNet-VMC calculations obtain state-of-the-art accuracy in ground state energy calculations for a range of molecules and solids [48; 49; 37; 43; 43; 50; 38].

We only consider calculations involving a single positron in the present work. We have discussed the treatment of the positronic spin coordinate only to demonstrate how our technique may be extended to calculations involving many positrons, as such systems have recently attracted theoretical interest.

Ground-state wavefunctions for bare and positronic molecules are not guaranteed to be similarly converged af
ter an equal number of gradient descent steps. This introduces uncontrolled error in calculating positron binding energies via VMC. Previous work has shown that Fermi-VMC calculations yield ground-state energies within chemical accuracy (\(\sim 1.5\) milliHartrees) of exact results for many small molecules [37; 38]. With this level of accuracy, the uncontrolled error will be negligible for molecules with a large positron binding energy. However, for molecules with very small binding energies, or large molecules for which the uncontrolled error may become large compared to the positron binding energy, there is no guarantee that an accurate estimate of the positron binding energy will be obtained by comparing ground-state calculations of different quality. In these cases, we employ the variance matching technique described by Entwistle _et al._[51], addressed in the Supplementary Material.

###### Acknowledgements.

 This work was undertaken with funding from the UK Engineering and Physical Sciences Research Council (EP/T51780X/1) (GC). Calculations were carried out with resources provided by the Baskerville Accelerated Compute Facility through a UK Research and Innovation Access to HPC grant, and we acknowledge PRACE for awarding us access to JUWELS at GCS@FZJ, Germany. Via his membership of the UK's HEC Materials Chemistry Consortium, which is funded by EPSRC (EP/R029431), Foulkes used the UK Materials and Molecular Modelling Hub for computational resources, MMM Hub, which is partially funded by EPSRC (EP/T022213).
