Ab initio quantum Monte Carlo calculations of spin superexchange in cuprates: the benchmarking case of Ca$_2$CuO$_3$

In view of the continuous theoretical efforts aimed at an accurate microscopic description of the strongly correlated transition metal oxides and related materials, we show that with continuum quantum Monte Carlo (QMC) calculations it is possible to obtain the value of the spin superexchange coupling constant of a copper oxide in a quantitatively excellent agreement with experiment. The variational nature of the QMC total energy allows us to identify the best trial wave function out of the available pool of wave functions, which makes the approach essentially free from adjustable parameters and thus truly ab initio. The present results on magnetic interactions suggest that QMC is capable of accurately describing ground state properties of strongly correlated materials.

In view of the continuous theoretical efforts aimed at an accurate microscopic description of the strongly correlated transition metal oxides and related materials, we show that with continuum quantum Monte Carlo (QMC) calculations it is possible to obtain the value of the spin superexchange coupling constant of a copper oxide in a quantitatively excellent agreement with experiment. The variational nature of the QMC total energy allows us to identify the best trial wave function out of the available pool of wave functions, which makes the approach essentially free from adjustable parameters and thus truly ab initio. The present results on magnetic interactions suggest that QMC is capable of accurately describing ground state properties of strongly correlated materials. For decades, transition metal oxides have been amongst the most intriguing materials due to the complex correlated behavior of the 3d or 4d electrons of a transition metal ion. In particular, strong electronic correlations often give rise to non-trivial magnetism, such as quantum spin-liquid states in low-dimensional Mott insulating oxides. High-temperature superconductivity in copper oxides (cuprates) is also believed to originate from magnetic spin excitations that bind Cooper pairs [1][2][3]. Electronic correlations, however, also make this class of materials among the most difficult to describe theoretically, both from model as well as ab initio perspectives.
One of the practical challenges of ab initio electronic structure theory is to accurately predict the strength of magnetic coupling between localized spins of transition metal ions [4]. For solids, a natural method of choice is periodic density functional theory (DFT). DFT gives access to the system's ground state energy corresponding to different configurations of localized spins, which can be mapped onto the eigenstates of the Ising model to extract the magnetic couplings J. This approach relies on the accuracy of the description of the ground state. Unfortunately, the presently available approximations to the exchange-correlation functional in DFT either poorly account for the exchange and correlation effects [local density approximation (LDA)] or depend on empirical input parameters (LDA+U, hybrid functionals). Often, the only way to find an appropriate approximation for the system of interest is by comparing theoretical calculations with experiment, compromising the predictive nature of such calculations.
This problem is generic to a broad class of transition metal oxides. Let us exemplify the aforementioned limitations of DFT by considering the case of the Mott insulators Ca 2 CuO 3 and Sr 2 CuO 3 . These systems are one of the best realizations of the one-dimensional (1D) spin-1 2 antiferromagnetic Heisenberg chain model, demonstrating spin-liquid behavior and separation of spin and orbital degrees of freedom [10]. The crystal structure of Ca 2 CuO 3 and Sr 2 CuO 3 is similar to that of the superconducting two-dimensional cuprates, with the difference that in the CuO 2 plane the oxygen atoms along the crystallographic b direction are missing, so that the Cu chains run along the a direction [ Fig. 1 (a)]. The Cu-O-Cu bridge provides a favorable path for superexchange coupling between Cu spins, resulting in a particularly strong coupling constant J. The experimental estimate of J has been extracted from various probes, performed mostly on  ), but, generally speaking, their applicability to condensed phase systems is intrinsically limited. In view of this, developing a new, more universal and accurate, ab initio approach to computing magnetic interactions is critically important. Needless to say, a method capable of accurately describing magnetic interactions in transition metal oxides will also provide an improved ab initio description of many other ground state properties of these complex systems, and hence may yield new physical insights.
Here, we apply the diffusion Monte Carlo method [13] within the fixed-phase approximation (FP-DMC) to compute the value of the spin superexchange interaction constant in a transition metal oxide. To the best of our knowledge, this is among the first calculations of magnetic couplings in complex oxides that has been attempted with a method capable of chemical accuracy in full periodic boundary conditions. The fixed-phase error is controlled by scanning over a set of trial wave functions and using the variational nature of DMC energy, as explained below. This way, the variational principle determines the choice of the initial DFT functional to generate a trial wave function and thus eliminates empiricism from the calculations. We choose the 1D cuprate antiferromagnet Ca 2 CuO 3 as our test system because of the simplicity of its underlying spin model and relatively light constituent atoms, as compared to other cuprates, including Sr 2 CuO 3 , to minimize the potential role of relativistic effects. Our result for the nearest-neighbor Cu spin coupling, J = 0.159 ± 0.014 eV, is in excellent agreement with the value extracted from the temperature dependence of magnetic susceptibility (Refs. 6-8, Table I).
Historically, quantum Monte Carlo has played a fundamental role in advancing electronic structure theory. Released-node DMC calculations on the homogeneous electron gas by Ceperley and Alder [15] were used to construct the local density approximation to the exchangecorrelation functional which is at the core of modern DFT methods. In FP-DMC, the Schrödinger equation is rewritten in a form of an integral diffusion equation, which is stochastically solved via quantum Monte Carlo sampling by iteratively propagating the wave function in imaginary time. As a result, the ground state (GS) wave function is projected out. In order to handle the fermion sign problem, the complex phase of the target GS wave  [11], while the dark gold band represents the magnetic susceptibility measurements [7,8].
function is fixed to those of an input trial wave function. The fixed-phase approximation introduces a variational systematic error in the DMC energy which can be assessed and controlled by comparing results obtained with different trial wave functions. Exchange and correlation effects are fully accounted for within the fixed phase approximation. Even though the computational costs of FP-DMC grow as N 3 , where N is the number of particles, the costs of calculating quantities per particle grow only as N 2 , which is a great advantage in the present case since J is a characteristic of a Cu-Cu bond. Due to improvements in algorithms and available computer power, QMC has achieved a growing success in accurately predicting the properties of complex materials [16][17][18][19][20][21]. We used the DFT plane-wave code Quantum Espresso [22] in order to self-consistently generate trial wave functions in a single Slater determinant form. The LDA+U method was applied, with the values of the on-site Coulomb repulsion between Cu 3d electrons, U , varying as 2, 3.5, 5.5, and 7 eV. The energy cut-off was set to 500 Ry due to the use of very hard pseudopotentials and inclusion of semi-core electrons in the valence (described below). The corresponding DMC energies were subsequently compared to determine the best trial wave function. The spin superexchange coupling constant J is calculated by mapping the FP-DMC total energies of the ferromagnetic (FM) and antiferromagnetic (AFM) Cu spin configurations onto the eigenstates of the 1D chain Ising model, since in DMC, similarly to DFT, the Hamiltonian commutes with S z , not S 2 . For all considered U values, both FM and AFM solutions are insulating in DFT. FP-DMC calculations were performed with QMCPACK [23]. The DMC imaginary time step and the number of walkers have been converged to, respectively, 0.005 Ha −1 and 2000 per boundary twist (in the 2 × 1 × 1 supercell). To assess finite size errors we considered 2×1×1 and 2×2×1 supercells, defined with respect to the conventional Immm crystallographic unit cell of Ca 2 CuO 3 [ Fig. 1 (a)]. The 2×1×1 supercell contains four formula units, i.e., 28 atoms with 228 electrons. We also performed averaging over twisted boundary conditions on a 2×4×1 k-point grid for the 2×1×1 supercell (2×2×1 for the 2 × 2 × 1 supercell). The necessity of twist averaging indicates that coupled cluster calculations, such as those in Ref. [9], are under converged with respect to finite size effects. The ionic potentials were approximated by employing pseudopotentials (PPs). Through an extensive investigation, we have found the inclusion of semi-core electrons to be essential for high quality results. Core sizes used for the pseudopotentials are as follows: Hecore for oxygen atoms (6 electrons in valence), Ne-core for calcium atoms (10 electrons in valence), and Ne-core for copper atoms (19 electrons in valence). The quality of the PPs has been carefully tested within both DFT and DMC, as reported in detail in the Supplementary Information.
We first present the FP-DMC results obtained for the 2 × 1 × 1 supercell, which contains two Cu atoms along the chain direction a. Figure 1 (b) displays the FP-DMC total energies of the FM and AFM states as a function of the trial wave function, characterized by the LDA+U parameter U . In these calculations we stress that the U is simply a convenient optimization parameter for generating FP-DMC wave functions. Both the FM and AFM curves roughly follow a parabolic U dependence, which suggests a small fixed-phase error and high quality of the trial wave functions. The minima of both curves are around U =3.5 eV. From the difference between the AFM and FM FP-DMC total energies we compute the spin superexchange constant J, shown in Fig. 1 (c) together with the respective LDA+U J values for comparison. Also indicated are the ranges of experimentally determined J values of Sr 2 CuO 3 : the light gold band represents all reported experimental estimates [11], while the dark gold band represents susceptibility measurements [6][7][8], which is one of the most reliable probes. Since, unfortunately, no equivalent experimental data on Ca 2 CuO 3 are available, we can only compare our theoretical calculations with the experiments performed on Sr 2 CuO 3 . This is a valid approach as the spin exchange couplings of the two cuprates should differ by no more than a few percent [24]. From Fig. 1 (c), one readily sees that at U =3.5 eV, corresponding to the minimal FP-DMC energies in Fig. 1 (b), the FP-DMC result J = 0.16(3) eV is in an excellent agreement with the susceptibility data. Noteworthily, the J values obtained with LDA+U deviate considerably from experiment at all U .
We now assess the finite size error associated with these results by performing FP-DMC calculations on a 2×2×1 supercell, obtained by doubling the original 2 × 1 × 1 supercell in the direction perpendicular to the Cu chains. The U =3.5 eV LDA+U trial wave function is used here as the one to provide the best complex phase for FP-DMC, as has been established above. The resulting J = 0.159 (14) eV is to be compared with the 2 × 1 × 1 result of 0.16(3) eV. From this, we conclude that the finite size error must be within 0.03 eV, which is the statistical accuracy of the 2 × 1 × 1 calculations of Fig. 1 (c).
In conclusion, we have presented a theoretical determination of the value of the spin superexchange constant in a transition metal oxide with the FP-DMC method. Our results for the 1D antiferromagnetic cuprate Ca 2 CuO 3 are in excellent agreement with experiment. Moreover, this is a purely ab initio approach, where the fixed phase error is controlled via the variational principle, with no empirical adjustable parameters. In this sense, FP-DMC is superior to DFT where in order to improve the description of exchange and correlations one often resorts to LDA+U or hybrid functionals and chooses the "best" functional empirically. The success of FP-DMC in the present case implies that this method is capable of accurately describing the complicated spin superexchange processes between the correlated Cu 3d orbitals and oxygen 2p orbitals, involving on-site Coulomb correlations and p−d orbital hybridization. We hope that our present successful application of FP-DMC will stimulate future studies of magnetic and other properties of strongly correlated transition metal oxides with this highly competitive ab initio method.
Note: At the time of submitting this manuscript, we learned of similar DMC calculations, performed independently, published on arXiv.org [25].
The authors would like to thank L. Shulenburger for sharing expertise in pseudopotential construction and for providing access to pseudopotential datasets prior to publication. The work was supported by the Materials Sciences & Engineering Division of the Office of Basic Energy Sciences, U.S. Department of Energy. PRCK was supported by the Scientific User Facilities Division, Office of Basic Energy Sciences, U.S. Department of Energy. Computational time used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725.

SUPPLEMENTARY INFORMATION
For oxygen and calcium we used the pseudopotentials (PPs) optimized by Shulenburger and Mattsson[S1], who also demonstrated their good quality by performing numerous tests. Of a much greater concern in the present study was the proper performance of the Cu PP since the magnetic properties of Ca 2 CuO 3 are largely determined by the behavior of the Cu 3d electrons. Therefore, we subjected our candidate Cu PPs to a comprehensive selection process, as presented below. The candidate Cu PPs were generated using the Opium code[S2]. Our final choice is the hard Ne-core Cu PP [ Fig. S2 (a), (c)] that satisfies the most stringent accuracy criteria and thus ensures the validity of the bulk calculations.
A. Bulk DFT calculations and rejection of Ar-core Cu pseudopotentials Using our hard Ne-core Cu PP, we are able to accurately reproduce with Quantum Espresso the LDA results of the all-electron (AE) code WIEN2k[S4] for bulk Ca 2 CuO 3 in an antiferromagnetic (AFM) and a ferromagnetic (FM) states. Thus, in Fig. S1 the Ca 2 CuO 3 densities of states (DOS) obtained from PP and AE calculations are compared, for an AFM [ Fig. S1 (a)] and an FM [ Fig. S1 (b)] Cu spin configurations. For both configurations, the agreement with the AE code is very good. Interestingly, we were able to equally well reproduce the AE DOS when also using properly optimized Mg-core Cu PPs. This, however, did not hold for any of the Ar-core Cu PPs we tried: in this case, the bandwidth of the PP states as well as the conduction gap (AFM configuration) are systematically larger than in the AE calculations. This allowed us to discard Ar-core Cu PPs already at this stage of PP validation.
As for the spin superexchange coupling constant J, Quantum Espresso with the hard Ne-core Cu PP gives 0.64 eV in LDA, while WIEN2k gives 0.72 eV. In LDA+U, the WIEN2k J is rapidly decreasing as 1/U which is not observed in the PP DFT results of Fig. 1 (c). This may be a peculiarity of the implementation of the LDA+U scheme within the plane-wave basis method of Quantum Espresso.

B. FP-DMC atom ionization energies and rejection
of Mg-core Cu pseudopotentials As we have pointed out, in DFT calculations for bulk Ca 2 CuO 3 , the Ne-core and the Mg-core Cu PPs appear to be equally good, provided proper optimization has been carried out. However, this does not necessarily mean that they will perform well in diffusion Monte Carlo calculations. In order to test the latter, we calculated the One of the reasons of the poor performance of the Mg-core PP in this test, is that the 3s orbital, which is treated as core here, has a significant overlap in space with the 3p orbital, treated as valence [see Fig. S2 (c)]. This causes less trouble in DFT which is formulated in terms of Kohn-Sham orbitals so that such a division based on orbital character is natural. DMC, on the other hand, operates with a full manyelectron wave function where a removal of the 3s electrons negatively affects the representation of the motion of the nearby 3p electrons. This issue has also been discussed in the context of GW [S5, S6]. Although the hard Ne-core PP has been proven to be of a good quality for both DFT and DMC, it has a disadvantage in terms of computational load and memory demands as it requires a minimum of 500 Ry energy cutoff for the Quantum Espresso plane-wave basis. This results from the quite small cut-off radii of the s-, p-, and d-channels, as shown in Fig. S2 (a). In view of this, we constructed and tested an alternative Ne-core Cu PP with a soft core and low energy cut-off requirements of less that 200 Ry [ Fig. S2 (b)]. It demonstrated excellent characteristics in DFT tests but, unfortunately, gives worse results in DMC than the hard-core Cu PP. In particular, with the soft-core PP the equilibrium interatomic separation distance in a CuO molecule is overestimated by more than 3% in FP-DMC [holds for LDA, LDA+U (U =3.5, 6 eV), and B3LYP wave functions], whereas with the hard-core PP it is overestimated by only 0.6% (Fig. S3). Also the Cu atom ionization energy is slightly underestimated: 7.548(42) eV. Thus all calculations in DMC reported in the paper were obtained with the hard Ne-core pseudopotential. The surprising sensitivities to pseudopotential formulation that we find even for small core pseudopotentials indicates that careful testing is essential and results from large core pseudopotentials must be treated with caution.