{'fraction_non_alphanumeric': 0.033134166214014125, 'fraction_numerical': 0.014122759369907659, 'mean_word_length': 5.028063610851263, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Chemical accuracy from quantum Monte Carlo for the Benzene Dimer

We report an accurate study of interactions between Benzene molecules using variational quantum Monte Carlo (VMC) and diffusion quantum Monte Carlo (DMC) methods. We compare these results with density functional theory (DFT) using different van der Waals (vdW) functionals. In our QMC calculations, we use accurate correlated trial wave functions including three-body Jastrow factors, and backflow transformations. We consider two benzene molecules in the parallel displaced (PD) geometry, and find that by highly optimizing the wave function and introducing more dynamical correlation into the wave function, we compute the weak chemical binding energy between aromatic rings accurately. We find optimal VMC and DMC binding energies of -2.3(4) and -2.7(3) kcal/mol, respectively. The best estimate of the CCSD(T)/CBS limit is -2.65(2) kcal/mol [E. Miliordos et al, J. Phys. Chem. A 118, 7568 (2014)]. Our results indicate that QMC methods give chemical accuracy for weakly bound van der Waals molecular interactions, comparable to results from the best quantum chemistry methods.

of those interactions involve a larger negative contribution from the electrostatic forces as in hydrogen bonded systems. If the electrostatic contribution is small, like in interactions of rare gas atoms with molecules, the minima depths are often below 1 kcal/mol. Therefore, from this point of view, a system such as the benzene dimer is also a van der Waals complex.
Finally, by including nonlocal terms in DFT correlation energy functional, vdW-DFs 25,26 include the long range nonlocal correlation energy obtained by the plasmon pole approximation. These functionals were originally applied on different van der Waals systems to obtain potential energy curve (PEC) 26 . DFT-vdW functionals result in significant improvements in equilibrium spacings between noncovalently bound complexes, as well as in binding energy of weak interacting systems. The efficiency and accuracy of different DFT-vdW functionals on bulk systems were analyzed [27][28][29] . The functionals are a clear improvement over semi-local functionals, although tests on a wider range of systems are desirable.
Numerous theoretical works have compared different approximate quantum chemistry based methods for noncovalent weakly bound systems [30][31][32][33] . Particularly the coupled-cluster theory through perturbative triplets CCSD(T) which is often considered as the gold standard for chemical accuracy 34,35 . However, due to its substantial computational cost, scaling as N 7 where N is the number of electrons, more efficient methods for vdW systems are highly desirable. These methods also cannot be applied to condensed matter. Thus we consider Quantum Monte Carlo (QMC), and compare with non-local density functionals.
Quantum Monte Carlo, which solves the electronic Schrödinger equation stochastically [36][37][38][39] , is an alternative approach to quantum mechanical methods. Diffusion quantum Monte Carlo (DMC) provides accurate energies for vdW systems [40][41][42][43][44][45][46] . DMC is also able to produce an accurate description of systems where many-body interactions play a key role 47,48 . In general QMC based methods are faster than the most accurate post-Hartree-Fock schemes for large number of particles N. The computational cost of QMC methods scales usually as N 3 -N 4 depending on the method.
The benzene dimer has become a benchmark system for electronic structure methods for systems where van der Waals interactions are important. Despite its simplicity the problem of identifying the global minimum structure is particularly challenging as there are only subtle differences in the binding energies of the different configurations. According to quantum chemistry results, two critical factors for the binding energy of the benzene dimer are basis set and electron correlation 49 . Our previous comprehensive study of benzene molecules 50 illustrates the importance of basis set in QMC energy calculations. Once the Jastrow factor is optimized by keeping fixed the Slater determinant, we obtained a good description of the atomization energy of the benzene molecule only when the basis of atomic orbitals is large enough and close to the CBS limit. In this work, we demonstrate that by using better trial wave functions and converged basis sets, we obtain a chemically accurate description of binding energy between aromatic rings.

II. COMPUTATIONAL DETAILS
We used the CASINO QMC code 51 with a trial function of Slater-Jastrow (SJ) form, where R is a 3N-dimensional vector that defines the positions of all N electrons, r ↑ i is the position of the i'th spin-up electron, r ↓ j is the position of the j'th spin-down electron, exp[J(R)] is the Jastrow factor, and det[ψ n (r ↑ i )] and det[ψ n (r ↓ j )] are Slater determinants of spin-up and spin-down one-electron orbitals. These orbitals were obtained from DFT calculations using the plane-wave-based Quantum Espresso code 52 . We used the local density approximation (LDA) to generate the orbitals in the Slater determinant for the trial wave function. We chose a very large basis-set cut-off of 200 Ry to guarantee convergence to the complete basis set limit 50 . The plane-wave orbitals were transformed into a blip polynomial basis 53,54 . The quality of the blip expansion , meaning the fineness of the blip grid, can be improved by increasing the grid multiplicity parameter and consequently results in a greater number of blip coefficients. The value of this parameter in our work is 2.0. The local density approximation (LDA) pseudopotentials are generated using the OPIUM pseudopotential generation program 55 . We also checked that the Kleinman-Bylander 56 transformation did not generate ghost states. In our DMC calculations, the pseudopotential energy was evaluated using a variational technique 57 . We used DMC time steps of 0.01 a.u. and 0.04 a.u. and extrapolated the results linearly to zero time step.
The Jastrow factor is a positive, symmetric, explicit function of interparticle distances.
We used a Jastrow factor consisting of polynomial one-body electron-nucleus (1B), two-body electron-electron (2B), and isotropic three-body electron-electron-nucleus (3B) terms. The main approximation in fermionic QMC is the fixed node approximation. To reduce this error, we used backflow transformation (BF) in our trial wave functions 58 . In the backflow transformation, the orbitals in the Slater determinant are evaluated not at the actual electron positions, but a quasi-electron positions that are functions of all the particle coordinates.
The backflow function, which describes the offset of the quasi-electron coordinates relative to the actual coordinates, contains free parameters to be determined by an optimization method. It allows the nodal surfaces to move within variational optimization, so with BF the QMC is no longer strictly fixed node. However, the subsequent DMC computations use the nodal surface that were determined during the VMC step.
In QMC calculations, correlated wave function can be obtained by replacing the single determinant by a sum over configuration state functions (CSFs), using BF transformations of the electronic coordinates or by using pairing wave functions [59][60][61][62][63] . Our BF transformation includes both electron-electron and electron-proton terms, given by where α ij (r ij ) and β iI (r iI ) are polynomial functions of electron-electron and electron-proton distance, respectively, containing optimizable parameters.
We use two methods for wave function optimization: variance minimization and energy minimization 64,65 . The parameters of Jastrow and backflow are first optimized by variance minimization at the variational Monte Carlo (VMC) level [66][67][68][69][70] . Since trial wave functions generally cannot exactly represent an eigenstate, the energy and variance minima do not co- incide. Therefore energy minimization should produce lower VMC energies. We have found that the lower VMC energies lead to lower DMC energies if we use backflow transformations.
This is due to improved many-body nodes as well as reduction in the errors induced by using non-local pseudopotentials.  tively. Free-energy landscape calculations using Car-Parrinello molecular meta-dynamics methods using the BLYP density functional with dispersion corrections predict that Tshape geometry is more stable at all temperatures 85 . However, the PD configuration with C 2h symmetry has been determined using optical absorption spectroscopy, whereas a polar V-shape configuration with C 2v symmetry has been suggested by multiphoton ionization mass spectroscopy 86,87 .    The VMC energies obtained with BF are lower than DMC energies without BF. This indicates that VMC with BF could be a useful level of theory to describe nonlocal long term interactions. In general, VMC calculations are significantly less expensive than DMC ones. Also, VMC has advantages for calculating expectation values of quantities more than the energy. Our QMC results clearly demonstrate that increasing the complexity of the wavefunction by including the BF correlations and 3B-Jastrow terms, plays key role in accurate describing of vdW interactions. the optimal values of binding energies at VMC level obtained by 1B+2B+3B-Jastrow and 1B+2B+3B-Jastrow plus BF correlations are -0.9(2), and -2.3(4) kcal/mol, respectively.
Employing 3B-Jastrow factor significantly improves the VMC binding energy, whereas it doesn't lower the DMC binding energy considerably. At the variational level, the inclusion of a 3B-Jastrow term provides additional dynamical correlation into the wavefunction and it is essentially useful for studying nonlocal vdW interactions.
Although a 3B-Jastrow factor improves the binding energy of benzene dimer in VMC, substantial enhancement is obtained by employing BF correlations. The optimal values indicate that VMC and DMC energies are improved by -1.4 and -0.9 kcal/mol, respectively.
It suggests that BF is effective at improving the nodal surface of benzene dimer in PD geometry. Considering that the LDA wave function often has too many nodal pockets, it is conceivable that BF coordinate transformations could modify the number of nodal pockets of a wave function. Table II shows the binding energy of benzene dimer in PD geometry obtained by different methods. Since the benzene dimer is a standard test for high-level quantum chemistry methods for proper characterization of vdW interactions there are many more results; a comprehensive comparison between different high-level ab initio approaches is reported recently 11 . They provide databases for noncovalent interactions. 49 bimolecular complexes in 345 geometry configurations are partitioned into subsets based on bonding motif. Benzene dimer is in the dispersion-dominated subset. Our DMC result is close to CCSD(T) (Table   II), which is considered the gold standard for chemical accuracy. The complete basis set (CBS) limit can now be estimated more precisely in the CCSD(T) framework. Our VMC results is comparable with those ones obtained by LRDMC method using Jastrow-AGP wave The T-shape configuration was studied before using QMC methods 46 . Using 1B, 2B, and 3B Jastrow factors, they found that the binding energy of T-shape configuration obtained by FN-DMC is -2.88(16) kcal/mol. An accurate study of noncovalent systems illustrates the  intermolecular interactions based on vdW dispersive forces. We find that adding 3B-Jastrow terms and BF transformations leads to significant improvement in the accuracy of the weak vdW interaction between aromatic rings. BF-VMC energies are significantly lower than SJ-VMC and therefore BF-VMC could be useful alternative for a SJ-DMC calculations , which are more expensive. The accuracy of our VMC results compared with DMC is evidence of the high accuracy of our trial wave function. BF correlations give substantial enhancement in trial wave function of aromatic rings. By improving the nodal surface of wave function, it leads to a significant reduction in binding energy between two benzene molecules. Improved trial wave functions will be useful in VMC calculations of quantities other than the energy, which are usually more difficult to obtain accurately than the energy.
We used single determinant wave functions in this work , but BF can also be combined with other types of wave functions such as multideterminant or pairing wave functions.

V. APPENDIX: ENERGY CURVE FITTING FUNCTIONS
In this appendix we report fitting parameters. We used piecewise polynomial fitting functions. Tables III and IV list fitting parameters including smoothing parameter α, the sum of squares due to error (SSE), R-square, Root Mean Squared Error (RMSE), and optimal parameters (opt) for each method.