Studies on the Transcorrelated Method

We investigate the possibility of using a transcorrelated (TC) Hamiltonian to describe electron correlation. A method to obtain TC wavefunctions was developed based on the mathematical framework of the bi-variational principle. This involves the construction of an effective TC Hamiltonian matrix, which can be solved in a self-consistent manner. This was optimized using a method we call second-order-moment minimization and demonstrate that it is possible to obtain highly accurate energies for some closed-shell atoms and helium-like ions. The effects of certain correlator terms on the description of electron–electron and electron–nuclear cusps were also examined graphically, and some TC wavefunctions were compared against near-exact Hylleraas wavefunctions.


■ INTRODUCTION
Capturing the effects of electron correlation is a central problem in electronic structure theory.A possible approach to tackle the problem involves the use of a similarity transformed Hamiltonian = H H e e , where the τ is a polynomial dependent on electronic positions and incorporates explicitly the correlation between various electron pairs.The use of such a Hamiltonian is known as the transcorrelated (TC) method.The inclusion of r 12 terms to describe electronic correlation can be dated back to Hylleraas 1 and later popularized by Kutzelnigg, 2 forming the basis of R12/F12 methodology 3,4 today.−11 This was done using a custom basis set and an optimized Jastrow factor.Hirschfelder, 12 Bartlett and Nooijen, 13 Klopper and co-workers, 14,15 and Ten-no 16−21 have also considered the use of such similaritytransformed Hamiltonians to eliminate the singularities associated with the r 1 ij term in the many-electron Hamiltonian.
The TC method joins in the family of explicitly correlated methods 22−24 and is closely related to canonical TC theory 25,26 and F12 methods.In canonical TC theory, the similarity transformation of the molecular Hamiltonian is not unitary.The formal difference leads to different approximations being made.F12 methods have now found widespread use in providing highly accurate electronic structure calculations.F12 methods employ geminal functions formed from the application of a rational generator (to impose cusp conditions) and a strong orthogonality projector on a two-electron determinant.The TC method attempts to recover dynamic correlation similarly by introducing cusp conditions and r 12 terms through its Jastrow factor.This represents a straightforward approach to considering the problem of explicitly correlated electrons.
There are two principal difficulties working with the TC Hamiltonian.First, the TC Hamiltonian will involve threeelectron operators, which can be expensive computationally.Second, the TC Hamiltonian is non-Hermitian.Unlike with Hermitian operators, the variational principle does not hold for non-Hermitian operators; hence, unphysical energies below that found from full configuration interaction (FCI) may be obtained.Furthermore, non-Hermitian matrices are difficult to work with due to the possibility for numerical instability in matrix computations.
However, with the introduction of variational Monte Carlo (VMC), the calculation became more computationally feasible and promising results were shown for a variety of atoms, molecules, 15,27−29 and periodic systems. 30,31To tackle the issue of non-Hermiticity, Luo proposed to replace the non-Hermitian ansatz with a Hermitian approximation so that a variational approach becomes viable. 32,33he TC Hamiltonian has also more recently been used with a variety of quantum chemistry methods with promising results.−37 Their numerical results show that the unboundedness of the non-Hermitian operator did not pose serious difficulties and that highly accurate results, even up to spectroscopic accuracies, 38 could be realized.More recently, they have used the TC Hamiltonian with coupled cluster, 39 and the energies found demonstrated better basis set convergence.On the other hand, Reiher and co-workers have developed a TC analogue of Density Matrix Renormalization Group and applied it to the Fermi-Hubbard model and some homo-nuclear diatomics. 40,41They have similarly found that the use of transcorrelation accelerates the convergence to the complete basis set limit.
Building on the recent successes of the TC method, this work shall attempt a deterministic approach and examine its viability as a computational tool to capture the effects of electron correlation.While Boys and Handy's formulation was deterministic, the computational resources of their time would have restricted the scope of their work.We therefore think that it would be useful to explore the possibilities of a deterministic approach with the computational resources today.Unlike Boys and Handy, however, we will solve the non-Hermitian Hamiltonian matrix self-consistently.Instead of parametrizing the correlator using the contraction equations, we propose to find it via what we shall call second-order-moment (SOM) minimization, an analogue of variance minimization that we have adapted for this work.
■ THEORETICAL BACKGROUND TC Hamiltonian.Given a Slater determinant Φ, Boys and Handy suggested 5 that the eigenfunction of the many-electron Hamiltonian can be approximately written as Φ TC = e τ Φ. τ is a function which accounts for electron correlation and is known as the correlator, while e τ is commonly referred to as the Jastrow factor.From the Schrodinger equation = H E e e (2) = H E (4) In the final step, we make the identification = H H e e .This is the TC Hamiltonian.The TC Hamiltonian can be expanded using the Baker−Campbell−Hausdorff expansion By using a correlator of the form τ = ∑ i<j u(r i , r j ), the third-and higher-order commutator terms vanish.The commutators can be further expanded to give where N e is the total number of electrons in the system studied.
Substituting τ = ∑ i<j u(r i , r j ), the TC Hamiltonian takes the form where ( , ) ( , ) The presence of the terms i j makes the TC Hamiltonian non-self-adjoint.This has been derived previously in several papers 6,14,15,35 but is recapitulated here for completeness.
One-Electron Effective Hamiltonian.The TC Hamiltonian is non-self-adjoint and will therefore have left-and right-eigenvectors.The left-and right-eigenvectors = ( ... ) are Slater determinants formed from molecular orbitals {ψ 1 ψ 2 ...ψ n } and {ϕ 1 ϕ 2 ...ϕ n }, respectively.The Slater determinants satisfy the following equations with E denoting the energy associated with the eigenvectors.A bi-orthonormal set of molecular orbitals, that is, ⟨ψ i |ϕ j ⟩ = δ ij can always be found via Loẅdin pairing 42 and hence we assume bi-orthonormality throughout this paper.
The TC energy can be identified as The denominator is unity due to the bi-orthonormality condition.The effective TC Hamiltonian can be found by taking the functional variation of the TC energy.The functional variation can be found by using the method of Lagrange multipliers.Forming the Lagrangian under the constraint of bi-orthonormal orbitals We seek the solution to to find a stationary point of the energy with respect to the constraint.We prove in the Appendix that using the condition = 0, we get the equation

Journal of Chemical Theory and Computation
We also introduce a notation = P ( 1) . S N is the symmetric group of degree N.For example 3 therefore gives all the possible permutations (with the correct parity) of the three-particle ket |ijk⟩.H eff is the effective TC Hamiltonian.It is a functional of the bi-orthogonal set of molecular orbitals {ψ i } and {ϕ i } and can thus be solved for iteratively through a self-consistent approach.
Jastrow Factor.The following form of the correlator was first introduced by Boys and Handy 8  (20)   Scaling of the inter-particle distances as r ̅ is known as the Padéform. 43Scaled distances are commonly used for Jastrow factors such that at large inter-particle distances, the terms in the Jastrow factors will approach a constant.There have been a number of scaling functions employed in literature. 44ollowing the work of Schmidt and Moskowitz, 45,46 we will use the Padéform with a = b = 1 due to the simplicity of implementation.
Optimizing Correlator Parameters.The correlators are a function of the set of parameters {c mno }.However, determination of these parameters is a non-trivial task.While the parameter c 001 in eq 18 has been determined previously to be 1/2 to satisfy the cusp condition, the other parameters c mno have yet to be determined.The unbounded nature of the nonself-adjoint TC Hamiltonian operator prevents the use of energy minimization for this.However, minimization of the local energy variance can be performed to find these parameters.Schmidt and Moskowitz applied VMC to calculate and minimize the variance.They performed this with correlators consisting of 7, 9, and 17 terms and found that with a 17 term correlator, 68−100% of the correlation energies for atoms helium through neon could be recovered using their variance-minimized parameters.
Handy also independently developed a variance minimization procedure to optimize the TC parameters. 47 The minimization of U TC was performed through the Davidson method and near-exact energy for the helium atom was calculated through this method, albeit with a slight modification of the Jastrow factor.However, helium is a twoelectron system; for any systems with more than two electrons, the three-electron operator in the TC Hamiltonian will in general give a non-zero term.As such, the calculation of the TC variance in eq 21 will require the evaluation of six-electron operators.The high computational cost and poor scaling has deterred research efforts along this line of inquiry.
SOM Minimization.While variance is well-defined for a self-adjoint operator, there is little literature for its non-selfadjoint counterpart.It is well known from linear algebra that a non-Hermitian matrix has left-and right-eigenvectors, which are not necessarily identical.The TC Hamiltonian is a nonself-adjoint operator and would similarly have left-and righteigenfunctions Ψ and Φ, respectively.We assume the use of a bi-orthogonal basis such that ⟨Ψ|Φ⟩ = 1.Instead of variance minimization, we propose the minimization of the SOM, a biorthogonal analogue of the variance for a non-self-adjoint Hamiltonian where . This is a bi-orthogonal extension to the usual definition of the variance (or second central moment in some papers 48,49 ).To the best of the authors' knowledge, the minimization of U SOM has not previously been performed.We shall first analyze some limiting cases to gain a better understanding of the quantity U SOM .
In the limit of = † H H (self-adjointness), the left-and righteigenfunctions become identical, Ψ = Φ.U SOM therefore reduces to the standard definition of the variance In the limit that Φ is an exact eigenfunction, that is, H̅ Φ = λΦ where where we have made use of the bi-orthogonality of the left-and right-eigenfunctions. Then A similar proof holds for the limit that Ψ is an exact eigenfunction.Since the exact eigenfunction has to satisfy the condition that U SOM = 0, the parameters should be varied such that the quantity U SOM becomes as close to zero as possible.The evaluation of U SOM would similarly require the evaluation

Journal of Chemical Theory and Computation
of six-electron terms (three from each H̅ ) and it is therefore as computationally challenging as Handy's TC variance.To sidestep this difficulty, the resolution of identity is employed.In the bi-orthogonal basis, the identity is given by where k runs through all of the possible Slater determinants for the given basis.
where the factor of a half was added to take into account double counting of ij and ab.The σ terms denote the various spins of electrons.The penultimate step is an approximation as we ignore the triple excitation terms when they are much smaller than the double excitation terms.In addition, the single excitation term has terms with: . By analogy to Brillouin's theorem for the Hartree−Fock method, single excitation determinants will not interact directly with the ground-state determinant, that is, . We can therefore ignore the single excitation terms and deduce that ab SOM (28)   The use of RI approximation in the context of transcorrelation is not new.The RI approximation had been used to evaluate the three-electron terms efficiently. 18In this work, the approximation was used to evaluate the expectation value of the square of the TC Hamiltonian.Without any approximations, the expectation value of a six-electron operator would have to be evaluated, which would have been computationally impractical.
Bi-Variational Principle.Having found the appropriate correlator parameters, we can construct the TC Hamiltonian and solve for its eigenfunctions.However, when the Hamiltonian is non-self-adjoint (as in the case for the TC Hamiltonian), the variational principle does not hold and the expectation value of the Hamiltonian is not bound from below.A naive minimization of the Hamiltonian's expectation value can therefore lead to values below the exact ground state energy, which are unphysical.However, one can formulate a different variational principle for a generic operator, which is not necessarily self-adjoint.While the mathematical exposition on the bi-variational principle has been previously undertaken by Loẅdin, 50−52 the essential parts of the proofs are reviewed here as it is crucial to the development of the TC method.
For a non-self-adjoint operator H̅ , we can define left and right eigenfunctions Ψ and Φ, respectively, such that We note that λ and μ are related by complex conjugation, that is, λ = μ* For a given pair of trial functions Ψ i and Φ i such that the expectation values λ i and μ i are given by The expectation values λ i and μ i have vanishing first-order variations (δλ i = 0), that is, they correspond to stationary points about the exact eigenvalues λ and μ, respectively.This is known as the bi-variational principle for a pair of adjoint operators. 50onversely, we can show that if δλ i = 0 for all δΦ and δΨ This implies that the trial function Φ i is an eigenfunction of H̅ with eigenvalue λ i and Ψ i is an eigenfunction of † H with eigenvalue λ i * = μ i .Equation 31 implies that if the trial functions Φ i and Ψ i are correct to the first order, the approximation of the eigenvalue λ i to the exact eigenvalue λ is correct to the second order.
Matrix Representation.The bi-variational equations can be recast in matrix form.In the following, the tensor notation of Head-Gordon et al. 53 shall be used.Given an atomic-orbital basis {χ 1 ...χ n }, we can expand any pair of trial functions Φ i and From eq 32, the bi-variation equations can then be expressed as The expressions in eq 35 were obtained through leftmultiplying by χ σ and integrating over all space.In the last step, we make the identification that Solving the TC Equation.We are now in a position to apply the bi-variational approach on the TC Hamiltonian.The effective TC Hamiltonian matrix has to be solved iteratively as the two-and three-electron terms are dependent on the trial functions Φ i and Ψ i .The following workflow was utilized:

Journal of Chemical Theory and Computation
1. Perform Hartree−Fock calculation and use the Hartree− Fock coefficients as a starting guess.2. Build the effective TC Hamiltonian matrix.3. Diagonalize the matrix to get new coefficients for the left-and right-eigenvectors.4. Repeat until convergence.
In doing so, we are simultaneously optimizing both the leftand right-eigenvectors.This is a different approach to that of Dobrautz, Luo, and Alavi 36 where only the right-eigenvector is optimized.While our approach requires the optimization of both left-and right-eigenvectors, which translates to a more expensive calculation, we gain the benefit of bounding the error of the calculation by the bi-variational principle.The bivariational approach proposed here shares similarity with that introduced by Ten-no. 18However, we treat the TC Hamiltonian in full, without approximations to the threeelectron term.The inclusion of three-electron contributions to the Fock matrix is also in contrast to Ten-no's approach 21 where the three-electron term is omitted.The error from our calculations comes from the error in using numerical grids and not from an approximation in our working equations.
Maximum Overlap Method.Convergence of the bivariational approach can be difficult in some cases.Taking inspiration from the work of Gilbert and co-workers, 54 we first assume that the Hartree−Fock coefficients are a good guess at our final coefficients.Therefore, at each iteration, the set of orbitals with the largest overlap to the occupied orbitals in the previous iterations will be picked.This process proceeds until convergence is reached.This is known as the maximum overlap method (MOM).Given the right coefficient matrix from the previous iteration C old , the left coefficient matrix from the current iteration D new , and the atomic orbital overlap matrix S, the maximum overlap matrix O MOM is given by The bi-orthogonal solutions from each iteration are determined only up to a phase factor, and hence the modulus is taken to ensure that the overlap remains positive.
Even with traditional implementations of MOM, it is found that it is possible for self-consistent field (SCF) iterations to converge onto unwanted solutions.This has led to the introduction of the initial MOM (IMOM), 55 where new orbitals in each iteration are picked based on their overlaps with the initial guess orbitals.This prevents the solutions from drifting away from the initial guess and has been shown to give better convergence to desired solutions.In this work, we adapt it for bi-orthogonal orbitals, such that the maximum overlap matrix O IMOM is given by where C initial is the initial left coefficient matrix.Both forms of the MOM were implemented for improved convergence of the iterative procedure.

■ COMPUTATIONAL DETAILS
The TC method is implemented in Python.The matrix elements relating to the correlator were found via numerical integration.These integrations were performed with grids found in the PySCF package.Throughout this work, we used Treutler−Ahlrichs grids with Becke partitioning.Q-Chem 5.3 was used for conventional non-orthogonal configurational interaction (NOCI) calculations and for finding Hartree−Fock solutions.Mathematica was used to plot Figures 2−8.
■ TC ENERGIES USING SCHMIDT−MOSKOWITZ PARAMETERS Schmidt and Moskowitz have previously found sets of 7, 9, and 17 correlator parameters for first-row atoms via variance minimization. 56Following Alavi and co-workers, we shall refer   to these sets as SM7, SM9, and SM17, respectively.The correlator has the form in eq 18.For ease of reference, the various terms incorporated in the correlator for SM7, SM9, and SM17 are tabulated in Table 1.
Using the correlator parameters found by Schmidt and Moskowitz, we solved the TC Hamiltonian for the first-row atoms with a series of correlation-consistent basis sets (cc-pVXZ, X = D, T, Q).This was done by using the corresponding Unrestricted Hartree−Fock orbitals as a starting guess and varying it until self-consistency.The data in Tables 2−4 show that for a fixed set of correlator parameters, the TC total energies of small atoms converge with basis set size.However, the TC energies do not necessarily decrease with an increasing number of parameters used.For example, in atoms from boron through neon, the TC energies increase going from 9 parameters to 17 parameters.This can be understood when we consider the origin of the parameters used.Schmidt and Moskowitz used Slater-type orbitals (STOs) in their optimization studies to obtain the parameters.On the other hand, this work employs Gaussian-type orbitals (GTOs).One major difference between the STOs and GTOs is that the electron-nuclear cusp condition is fulfilled while using STOs but not when using GTOs.Different corrections are therefore required for the Hartree−Fock solutions expressed with different orbital bases, leading to the need for different parameters.Hence, the parameters used in this study may not be optimal.
Alavi and co-workers have also used correlation-consistent bases in their work on the TC Hamiltonian.However, they are able to find energies in excellent agreement with experimental values (Tables 2−4).We believe that this is due to the effective multi-reference nature of the FCIQMC method such that any errors incurred from using these SM parameters are corrected for by adjusting the weight of each determinant.
■ SOM MINIMIZATION OF CORRELATOR PARAMETERS Singlet State Atoms.To improve upon the accuracy of our results, we allowed correlator parameters to vary alongside the orbitals.The correlator parameters were optimized by using SOM minimization.The parameters found from Schmidt and Moskowitz (the set of 7 parameters) were used as a starting guess for the optimization, with an additional mno = 001 term to correct for the electron-nuclear cusp conditions (SOM8 in Table 1).The starting guess for the mno = 001 term is zero.This set of 8 parameters will be referred to as SOM8.The BFGS algorithm was used for optimization.In practice, we have found it to be useful to optimize the parameters using a two-step SOM minimization procedure where we first keep the orbitals fixed through the optimization process and after the first round of optimization, we perform SOM minimization with orbital relaxation in each iteration of the second optimization cycle.In doing so, we are less likely to get caught in local minima after orbital relaxation.The TC energies found using the optimized parameters are tabulated in Table 5. Closed-shell singlet state atoms were used as examples as it side-steps the possible issue of spin-contamination since we would be able to use a Restricted Hartree−Fock determinant as the reference determinant.
The energies found did not appear to suffer from nonvariationality.For the 1 S states of Be, C, and O, very accurate energies could be found, which recover more than 98% of the correlation energy.While the results for helium and neon were not as encouraging, we note that a highly accurate energy of −2.9037E h for helium could be found by using a different starting guess (Table 6).A similarly good result of −2.9033E h (see Table 7) can also be found by using the SOM18 set of parameters.
Figure 4. Plotting the Hylleraas wavefunction against the TC wavefunction with c = −1.12.c was found by SOM minimization and was performed with a cc-pV5Z basis set.The Slater determinant used is found via Hartree−Fock calculation using a cc-pVQZ basis.While there is some discrepancy between the TC wavefunction and the Hylleraas wavefunction near the nucleus, the two wavefunctions are very similar further away from the nucleus.In practice, there were a number of local minima found during optimizations depending on the starting guess and hence a number of possible energies can in principle be found.A sample of possible solutions for helium are tabulated in Table 6.This suggests that the starting guess is very important to obtain the right correlator parameters and one should be cautious about the parameters found from such an optimization.
In the case of neon, the use of a larger set of parameters (SOM18) gave a non-variational energy of −129.0019Eh .This demonstrates the possibility of obtaining non-variational energies and highlights the potential pitfalls of using SOM minimization to obtain correlator parameters.
Helium-like Systems.Encouraged by the possibility of highly accurate energies using SOM minimization, we examined the approach on a series of helium-like systems (Table 7).To find the correlator parameters for this series, we used the set of 17 parameters from Schmidt and Moskowitz and with an additional mno = 001 term as a starting guess for the helium atom and performed SOM minimization to obtain the optimized parameters (SOM18).These optimized parameters were then used as starting guesses for each of these ions.
We found that it was important to use an augmented basis set (aug-cc-pVQZ) for the negatively charged hydride anion as more diffuse functions are required to describe the expanded orbitals.In contrast, a basis optimized for describing core−core correlations (cc-pCVQZ) was found to be useful to describe the contracted orbitals in cations.
For comparison, the same calculations were performed with a cc-pVQZ basis set and the results are tabulated in Table 8.The use of the cc-pVQZ basis increased the absolute error from SOM minimization and the TC energies found were   For example, 001 corresponds to the m = 0, n = 0, and o = 1 term, i.e., r ij term."SM" refers to Schmidt−Moskowitz parameters while "SOM" refers to parameters found via SOM minimization.The parameters were the same as those used by Schmidt and Moskowitz. 56FCIQMC (cc-pVQZ basis) data were found by Alavi and coworkers. 35Experimental values were found by Chakravorty and co-workers. 57he difference (in Hartrees) and the percentage of correlation energy found were similarly reported.*The exact energies of 1 S C and 1 S O were deduced from spectroscopic measurements. 58,59The exact energies of the other closed shell atoms were taken from experimental values found by found by Chakravorty and co-workers. 57The optimized parameters can be found in the Appendix (Table 10).All calculations were run with a cc-pVQZ basis set.
mostly non-variational.This shows that the choice of basis set is imperative to the accuracy of the TC method.
Using the appropriate basis sets for cations and anions, we found that the absolute error from SOM minimization is lower than that for HF for each ion (Table 7).From H − through N 5+ , SOM minimization recovers a large proportion of the correlation energy.From O 6+ through Ne 8+ , the percentage of correlation energy recovered drops considerably.This is likely due to the highly contracted nature of the 1s orbitals in these highly charged cations and a bigger basis with more contracted basis functions would be required to more accurately describe the electron correlation.However, it is possible that a better starting guess could similarly improve the correlation energy recovered.
It is also gratifying to note that most of the energies found from SOM minimization using appropriate basis sets do not exhibit non-variationality.Li + and Be 2+ were found to have non-variational energies .However, the error is small and within chemical accuracy (within ∼0.0016E h ).
Ionization Energies.A further investigation of the TC method could be made by determining the ionization energies within the SOM minimization framework.For consistency, we shall use the cc-pCVQZ basis set for the following ionization processes (Table 9): Be → Be 2+ , C → C 4+ , and O → O 6+ .
It is observed that the SOM18 energies for the atoms are higher than that found in Table 5.This is an indication that the near-exact energies found with 8 parameters previously may have been fortuitous.The comparatively poorer description of the atoms as compared to the ions consequently leads to poorer ionization energies.While the predicted ionization energies are an improvement from Hartree−Fock levels of theory, they exhibit appreciable (∼30 mE h ) errors when compared to the exact ionization energies.This shows that while SOM minimization could in principle provide near-exact energies (both absolute and relative), the optimization procedure may be difficult in practice.

■ GRAPHICAL ANALYSIS OF CORRELATION
Electron−Electron Cusp.To better understand the effects of the electron−electron and electron−nucleus terms in the correlator, we studied the effects of various Jastrow factors e τ on a Hartree−Fock solution Φ HF of a helium atom.We first attempted to study the effects of varying the angle θ between two electrons confined to the same electron−nucleus distance (Figure 1, left).Using the correlator = + c r r 1 12 12 and the Slater determinant found with a Hartree−Fock calculation with a cc-pVQZ basis, the TC wavefunction e τ Φ HF was plotted as a function of θ (Figure 2) for varying values of the parameter c.The set of parameters SOM8 was used, with starting guesses of 0 unless otherwise stated in the first column.Starting from SOM8 with c 001 = 0.5, c 100 = −2, and other parameters zero, a highly accurate energy of the helium atom could be found.The optimized parameters can be found in the Appendix (Table 11).All calculations were run with a cc-pVQZ basis set.The absolute error for both HF and SOM minimization methods increases with the magnitude of nuclear charge.The absolute error found from HF is consistently above that of those found from SOM minimization.The optimized correlator parameters found are tabulated in the Appendix (Tables 12 and 13).a The absolute error from SOM minimization is significantly higher than those found in Table 7.Most of the TC energies found are also non-variational.
The 6-term Hylleraas wavefunction, 60 which represents a good approximation to the exact wavefunction of helium, is also plotted for comparison.
For ease of reference, the Hylleraas wavefunction for He is given by 60 where The Jastrow factor's introduction of electron−electron cusps to the Hartree−Fock solution can be seen from Figure 2. The shape of the electron−electron cusp gets increasingly similar to that of the Hylleraas wavefunction as the coefficient increases from 0.1 to 0.5.This supports the use of = + r r to correct for the electron−electron cusps.The coefficient = c 1 2 is fixed in the TC calculations, which therefore necessitates the need for higher order terms in r ij to correct for the depth of the cusp.
Electron−Nuclear Cusp.To examine the effect of the electron−nuclear cusp, we use the special case whereby the nucleus and an electron are constrained to be 1.26 Bohr away and the other electron is free to move along the line defined by them (Figure 1, right).We use a different correlator where r is the variable electron−nuclear distance and r 12 = r − R is the electron−electron distance and vary the value of parameter c.The Hartree−Fock wavefunction most closely matches that of the Hylleraas wavefunction at the nucleus while the function with c = −1 is more similar further away from it (Figure 3).For comparison, the function with c = −2 was plotted as c = −2 and is what would be expected from a simple application of Kato's cusp conditions.This illustrates that the coefficient c need not be equal to the negative of the nuclear charge, −Z, as the electron−nuclear interaction term in the Jastrow factor affects the overall wavefunction and not only at the cusp.Performing SOM minimization with this we found a value of c = −1.12,and the TC wavefunction with c = −1.12 is plotted.The plot good agreement with the Hylleraas wavefunction (Figure 4).
A similar series of calculations were attempted for Li + .The Hylleraas wavefunction for Li + is given by From 5, can be seen that the Hartree−Fock wavefunction is very similar to that of the Hylleraas wavefunction near the nucleus.At regions further from the nucleus, the c = −1 wavefunction most closely resembles the Hylleraas wavefunction.SOM minimization gave c = −3.79,but in this case, we were unable to reproduce the Hylleraas wavefunction (Figure 6).There are several reasons why such a discrepancy can exist.
First, resolution of identity is an approximation, which may not be valid depending on the size of the basis set used.
Second, the Hartree−Fock wavefunction resembles the Hylleraas wavefunction well.While the addition of a Jastrow factor to it can provide a better description of the electron− nuclear cusp, it comes at the cost of affecting other parts of the wavefunction.More terms in the correlator may need to be added to more accurately describe the TC wavefunction at different points in space.
To address the latter point, the TC wavefunctions for He and Li + found using the SOM18 set of parameters (Table 7) were plotted against the respective Hylleraas wavefunctions (Figures 7 and 8).From Table 7, it can be observed that the TC energies are very close to the exact energies, and this suggests that the TC wavefunctions should look similar to that of the Hylleraas wavefunction.This is reflected graphically in the depiction of the electron−electron cusp (Figure 7) and electron-nuclear cusp (Figure 8), where each of the TC wavefunctions agree well with the corresponding Hylleraas wavefunction.The main difference between the TC wavefunction and Hylleraas wavefunctions appears to be the description of electron−electron interactions at larger electron−electron distances, which may hint at the use of a differently scaled form of the correlator to account for longer range effects.Overall, the plots show that SOM minimization can get highly accurate wavefunctions given sufficiently many correlator parameters, supporting the utility of SOM minimization when appropriate starting guesses are used.

■ CONCLUSIONS
A self-consistent method for solving the non-self-adjoint TC Hamiltonian has been implemented successfully.We have found that it is possible to obtain highly accurate energies of some first row atoms from this approach.The correlator parameters found in the literature are not optimized for the Gaussian orbital basis used in this current study and had to be re-optimized through a method we refer to as SOM minimization.This allowed us to find optimized parameters for any system, in principle.However, the optimization of multiple parameters is challenging and in practice, we have found it to be useful to optimize the parameters using a twostep SOM minimization procedure.SOM minimization has been found to give good energies for the first row atoms.However, the percentage of correlation energy recovered has been found to decrease with increased nuclear charge across a series of helium-like ions.We believe that this is due to the inability of the basis set to accurately describe highly charged cations, and a custom basis with more contracted basis The error from SOM minimization is smaller than that found by HF.The ionization energies are found by taking the energy difference between the closed shell singlet atoms and that of the two-electron ion.The optimized correlator parameters found are tabulated in the Appendix (Table 14).
functions would be helpful to describe the correlation in these systems.The analysis on ionization energies of select atomic systems suggests that there is still room to improve upon the SOM minimization method.We leave this for future work.A graphical analysis has also been done to illustrate the effects of some correlator terms on the overall wavefunction and demonstrated the importance of including higher-order correlator terms in the Jastrow factor to give a more accurate Developing methods to make SOM minimization more robust for calculating near-exact energies for a range of systems would also be desirable, and we hope that this work will inspire further efforts in this area.

■ APPENDIX A TC Hamiltonian
Expansion of Commutator Terms.We first separate the many-electron electronic Hamiltonian into a kinetic energy T ( ) and potential energy V ( ) term V is multiplicative and hence commutes with τ in the commutator [ H , , is The commutator terms in eq evaluated term by term as follows We can now make the identification Therefore Since [[ ] ] H , , is a multiplicative term, higher-order commutators of the form [[[ ] ] ] H , , ... vanish.Given that τ = ∑ i<j u(r i , r j ), we can further expand eq 6. symmetry u(r i , r j ) = u(r j , r i ).We therefore find where we relabeled by j and used the symmetry of u in the last line.

Journal of Chemical Theory and Computation
We similarly find Finally Substituting these commutator terms back into eq 5 recovers eq 6.
Lagrangian Approach to the TC Hamiltonian.We show here that the method of Lagrange multipliers can be used to derive the effective TC Hamiltonian.
In last line, we have renamed n-electron operators by the O .This is for brevity of notation and the understanding that the mathematics after is concerned only with the number of electrons the operators act upon.Taking an infinitesimal change of the Lagrangian with respect to variation of the spinorbitals We shall analyze each term of the above expression in turn.We will adopt the shorthand We assume that the bra always contains molecular orbitals from the set {ψ i } and that the ket always contains molecular orbitals from the set {ϕ i }.
One-Electron Term.

Journal Theory and Computation
where we used permutation symmetry of the integrals, e.g., .
Three-Electron Term.
where we have defined the following (3 3 We use the shorthand such that We seek for any arbitrary changes of ψ and ϕ independently.Hence, 0 and = independently.
Effective TC Hamiltonian.Consider the condition  These parameters are used to obtain the TC energies in Table 5.

Journal of Chemical Theory and Computation
Since the expression holds for any δψ i *, the terms in the square bracket must be zero.Hence for all i.Rearrangement of the equation recovers the equation given in eq 14.

■ B Correlator Parameters
Correlator Parameters for Closed-Shell Atoms (SOM8).The correlator parameters for closed-shell ions (Table 5) found from using SOM minimization on a set of 8 parameters are tabulated in Table 10.
Correlator Parameters for Helium from Different Initial Guesses (SOM8).The correlator parameters for helium (Table 6) found from using different initial guesses are tabulated in Table 11.
Correlator Parameters for Helium-like Ions (SOM18).The correlator parameters for helium-like ions found from using SOM minimization on a set of 18 parameters are tabulated in Tables 12 and 13.The parameters for helium were used as a starting guess for the helium-like systems in Table 7.

Figure 1 .
Figure1.(Left) two electrons constrained at a fixed electron−nuclear distance of R. The electron−electron distance is modulated by their angle of separation, θ. θ is measured in radians.(Right) one electron is fixed at distance R and the other is constrained to the (dotted) line defined by the nucleus and the fixed electron.R is found by the expectation value of the electron−nuclear separation in near-exact wavefunction given by Nakashima and Nakatsuji and co-workers.61

Figure 2 .
Figure 2. Graphical description of the electron−electron cusp of the TC wavefunction with c = 0.1 through c = 0.5 against Hartree−Fock and Hylleraas wavefunctions.The function with c = 0.5 most closely matches that of the Hylleraas wavefunction, but the cusp is shallow as compared to the Hylleraas wavefunction.

Figure 3 .
Figure 3. Graphical description of the electron−nuclear cusp of the TC wavefunction of He with c = −2, −1, and 1 against Hartree−Fock and Hylleraas wavefunctions.The position at which the other electron is fixed is shown by the dashed blue line.

Figure 5 .
Figure 5. Graphical description of the electron−nuclear cusp of the TC wavefunction of Li + with varying values of c against Hartree−Fock and Hylleraas wavefunctions.The Hartree−Fock wavefunction resembles the Hylleraas wavefunction near the nucleus.The TC wavefunction with c = −1 most closely matches that of the Hylleraas wavefunction further away from the nucleus.

Figure 6 .
Figure 6.Plotting the Hylleraas wavefunction against the TC wavefunction with c = −3.79 and c = −1.c was found by SOM minimization with a cc-pCVQZ basis set.The Slater determinant used is found via Hartree−Fock calculation using a cc-pVQZ basis.

Figure 7 .
Figure 7. Two plots describing the electron−electron cusp corresponding to Figure 1 (Left).(Left) plot of the TC wavefunction for He against the corresponding Hylleraas wavefunction.(Right) plot of the TC wavefunction for Li + against the corresponding Hylleraas wavefunction.Both plots show that the 18 parameter TC wavefunction reproduces the shape of the electron−electron cusp, but the cusp is shallower than the Hylleraas wavefunction.

Figure 8 .
Figure 8. Two plots describing the electron−nuclear cusp corresponding to Figure 1 (Right).(Left) plot of the TC wavefunction for He against the corresponding Hylleraas wavefunction.(Right) plot of the TC wavefunction for Li + against the corresponding Hylleraas wavefunction.Both plots show that the 18 parameter TC wavefunction reproduces the Hylleraas wavefunction well.
0.5, c 100 = +1 −2.8969 c 001 = 0.5, c 100 =−1 −2.8989 c 001 = 0.5, c 100 =−2 −2.9037 a wavefunction.Thus far, SOM minimization has been attempted for closedshell systems.Further work has to be done on open-shell systems where there is a possibility of spin-symmetry breaking, leading to an unphysical wavefunction.This could in principle be tackled with NOCI, 62−64 which forms the basis for future work.SOM minimization should also be attempted on larger systems to test if the method works more generally.However, as pointed out by Alavi and co-workers, 35 memory use is a bottleneck in TC calculations.This can be challenging especially the need to use large basis sets for SOM minimization as it relies on the resolution-of-the-identity approximation.A possible solution would be to use an auxiliary basis set instead, which would allow us to use a smaller basis set to represent the Slater determinant but a large auxiliary basis to satisfy the resolution-of-the-identity approximation.

Table 1 .
Summary of the Various Sets of Parameters Used, Where Each Number in the Second Column Has the Form mno a

Table 2 .
57mparison of the TC Total Energies (in Hartrees) Found with the Bi-Variational Approach Using 7 Parameters against Literature and Experimental Values aThe parameters were the same as that used by Schmidt and Moskowitz.56FCIQMC(cc-pVQZbasis)datawere found by Alavi and co-workers.35Experimentalvalueswere found by Chakravorty and co-workers.57 a

Table 3 .
Comparison of the TC Total Energies (in Hartrees) Found with the Bi-Variational Approach Using 9 Parameters against Literature and Experimental Values a 57he parameters were the same as those used by Schmidt and Moskowitz.56FCIQMC(cc-pVQZbasis) data have not been previously reported.Experimental values were found by Chakravorty and co-workers.57

Table 4 .
Comparison of the TC Total Energies (in Hartrees) Found with the Bi-Variational Approach Using 17 Parameters against Literature and Experimental Values a

Table 5 .
Comparison of Energies Found after SOM Minimization and the Exact Energies a

Table 6 .
TC Energies of the Helium Atom in Hartrees for Different Starting Guesses a

Table 7 .
Comparison of Absolute Error in the HF Energy against the Absolute Error in the Energy Found by SOM Minimization a

Table 8 .
Comparison of Absolute Error in the HF Energy against Those Found by SOM Minimization in the Cc-pVQZ Basis a

Table 9 .
Comparison of Ionization Energies Found in HF against Those Found by SOM Minimization in the Cc-pCVQZ Basis a

Table 11 .
Correlator Parameters for Helium Found from Using SOM Minimization with Different Initial Guesses a a These parameters are used to obtain the TC energies in Table6.

Table 10 .
Correlator Parameters for Closed-Shell Atoms Found from Using SOM Minimization a

Table 12 .
Correlator Parameters for Helium-like Ions (H − to B 3+ ) Found from Using SOM Minimization

Table 13 .
Correlator Parameters for Helium-like Ions (C 4+ to Ne 8+ ) Found from Using SOM Minimization

Table 14 .
Correlator Parameters for Closed-Shell Singlet States for Be, C, and O Found from Using SOM Minimization