Essentially exact ground-state calculations by superpositions of nonorthogonal Slater determinants.

An essentially exact ground-state calculation algorithm for few-electron systems based on superposition of nonorthogonal Slater determinants (SDs) is described, and its convergence properties to ground states are examined. A linear combination of SDs is adopted as many-electron wave functions, and all one-electron wave functions are updated by employing linearly independent multiple correction vectors on the basis of the variational principle. The improvement of the convergence performance to the ground state given by the multi-direction search is shown through comparisons with the conventional steepest descent method. The accuracy and applicability of the proposed scheme are also demonstrated by calculations of the potential energy curves of few-electron molecular systems, compared with the conventional quantum chemistry calculation techniques.

Background

In recent years, there have been many significant achievements regarding electronic structure calculations in the fields of computational physics and chemistry. However, theoretical and methodological approaches for providing practical descriptions and tractable calculation schemes of the electron–electron correlation energy with continuously controllable accuracy still remain as significant issues [1-15]. Although density functional theory (DFT) supplies a computationally economical and practical method, there are many unexplored problems raised by unreliable results obtained for some systems in which highly accurate electron–electron correlation energy calculations are required, since results by DFT depend significantly on the exchange-correlation energy functional used to perform the calculation [16-18].

The available quantitatively reliable methods require higher computational costs than the DFT method [18]. Although quantum Monte Carlo methods [19-23] can be applied to molecular and crystal systems and show good quantitative reliability where extremely high-accuracy calculations are required, difficulties in calculating forces for optimizing atomic configurations are a considerable disadvantage and inhibit this method from becoming a standard molecular dynamics calculation technique. Configuration interaction (CI), coupled cluster, and Møller-Plesset second-order perturbation methods, each of which use a linear combination of orthogonalized Slater determinants (SDs) as many-electron wave functions, are standard computational techniques in quantum chemistry by which highly accurate results are obtained [24], despite suffering from basis set superposition and basis set incompleteness errors. The full CI calculation can perform an exact electron–electron correlation energy calculation in a space given by an arbitrary basis set. However, it is only applicable for small molecules with modest basis sets since the required number of SDs grows explosively on the order of the factorial of the number of basis.

The required number of SDs in order to determine ground-state energies can be drastically decreased by employing nonorthogonal SDs as a basis set. The resonating Hartree-Fock method proposed by Fukutome utilizes nonorthogonal SDs, and many noteworthy results have been reported [25-30]. Also, Imada and co-workers [31-33] and Kojo and Hirose [34,35] employed nonorthogonal SDs in path integral renormalization group calculations. Goto and co-workers developed the direct energy minimization method using nonorthogonal SDs [36-39] based on the real-space finite-difference formalism [40,41]. In these previous studies, steepest descent directions and acceleration parameters are calculated to update one-electron wave functions on the basis of a variational principle [25-30,36-39]. Although the steepest descent direction guarantees a secure approach to the ground state, a more effective updating process might be performed in a multi-direction search.

In the present study, a calculation algorithm showing an arbitrary set of linearly independent correction vectors is employed to optimize one-electron wave functions with Gaussian basis sets. Since the dimension of the search space depends on the number of linearly independent correction vectors, a sufficient number of correction vectors ensure effective optimization, and the iterative updating of all the one-electron wave functions leads to smooth convergence to the ground states. The primary purpose of this article is to demonstrate the advantage of using multiple correction vectors in searching for the ground state over the conventional steepest descent search in which only one correction vector is used. As a demonstration of the accuracy and applicability of the proposed calculation algorithm, essentially exact potential energy curves of few-electron molecular systems with long interatomic distances are described for cases where the conventional calculation methods of quantum chemistry fail.

The organization of the article is as follows. In the ‘Optimization algorithm’ section, the proposed calculation algorithm for constructing a basis set of nonorthogonal SDs by updating one-electron wave functions with multiple correction vectors is described. The expression of the conventional steepest descent direction with a Gaussian basis set is also given for comparison. The convergence characteristics to the ground states of few-electron systems for calculations using single and multiple correction vectors are illustrated in the ‘Applications to few-electron molecular systems’ section. As demonstrations of the proposed calculation procedure, the convergence properties to the ground states of few-electron atomic and molecular systems are also shown. Finally, a summary of the present study is given in the ‘Conclusions’ section.

Optimization algorithm

The calculation procedures for constructing a basis set consisting of nonorthogonal SDs for N-electron systems using single and multiple correction vectors are described here. An N-electron wave function ψ(r1, σ1, r2, σ2,…, r, σ) is expressed by a linear combination of nonorthogonal SDs as follows:

Here, rand σdenote the position and spin index of the ith electron, respectively. L is the number of SDs, and Φ(r1, σ1, r2, σ2,…, r, σ) is the Ath SD, given by

with ϕ(r) and γ(σ) being nonorthogonal and unnormalized one-electron basis functions and spin orbital functions, respectively. The one-electron wave function ϕ(r) is constructed as a linear combination of Gaussian basis functions x(r) [24] as

Here, M and D are the number of basis functions and the sth expansion coefficient for the ith one-electron wave function ϕ(r), respectively.

The steepest direction is implemented in the expression of the total energy functional E of the target system on the basis of the variational principle, without the constraints of orthogonality and normalization on the one-electron wave functions. The updating procedure of the pth one-electron wave function belongs to the Ath SD which is represented as

where a is the acceleration parameter, which is determined by the variational principle with respect to the total energy E, i.e., [28]

The component of the steepest descent vector K is given by

where

and

Here, denotes the element of the jth row and ith column of the matrix . The overlap integral matrix Sis defined by

where the elements Sare the overlap integrals between the one-electron basis functions, i.e.,

and are the matrix elements of the Hamiltonians,

and

respectively. Here, V(r) stands for an external potential.

The proposed calculation procedure employs linearly independent multiple correction vectors for updating the one-electron wave function. The pth one-electron wave function in the Ath SD is updated by

where C(j = 1, 2,…, L + N) and N are the expansion coefficient and the number of correction vectors, respectively. The components of the correction vectors G determine N linearly independent correction functions ξ(r) which are defined as linear combinations of Gaussian basis functions as

Since the linearly independent correction vectors can be given arbitrarily, randomly chosen values are employed in the present study. A larger number of correction vectors N realize a larger volume search space; however, the number of the linearly independent vectors N is restricted to the dimension of the space defined by the basis set used.

Thus, we have a linear combination of L + N SDs as the new N-electron wave function

where

Figure 1 illustrates the flow of the present calculation procedure. Unrestricted Hartree-Fock (UHF) solutions for a target system are used for initial one-electron wave functions. The coefficients of Equation 17 are given by solving the generalized eigenvalue equations obtained by employing the variational principle applied to the total energy, and we can have a new N-electron wave function as a linear combination of L SDs as shown in Equation 17. Iteration of the above updating process for all the one-electron wave functions of all SDs increasing the number of the SDs’ L leads to an essentially exact numerical solution of the ground state.

Figure 1

Flow of the present algorithm.

Applications to few-electron molecular systems

Convergence performances for searching for the ground state of a C atom with the 6-31G** basis set are shown in Figure 2. The UHF solutions are adopted as initial states, and the number of employed SDs is 30. The steepest descent direction and acceleration parameter are adopted for the calculation using one correction vector (N =1), and seven randomly chosen linearly independent correction vectors are added to the steepest descent correction to create a calculation with eight correction vectors (N =8). An indispensable advantage of the multi-direction search over the single steepest descent direction search is clearly demonstrated. Although the steepest descent vector gives the direction with the largest gradient, it does not necessarily point toward the global energy minimum state. On the contrary, a linear combination of multiple correction vectors can be used to obtain the minimum energy state within the given space by adopting the variation principle.

Figure 2

Effectiveness of multi-direction search on total energy convergence. Effectiveness of multi-direction search on total energy convergence as a function of the number of iterations for a C atom with the 6-31G** basis set is shown.

Figure 3 illustrates the convergence performance of the proposed method for the electron–electron correlation energy of a HF molecule with the 6-31G** basis set as a function of the number of employed SDs. Calculated correlation energies are shown by ratios to exact ones obtained by full CI. The convergence performance to the exact ground state is improved by increasing the number of correction vectors, since the volume of the search space for a one-electron wave function increase with increasing N. The essentially exact ground-state energy is obtained using less than 100 nonorthogonal SDs with an error of 0.001%, compared with the exact value in which 99.5% of the electron–electron correlation energy is counted. The obtained convergence is so smooth that the accuracy of the total energy is controllable by adjusting the number of employed SDs. On the other hand, the full CI method requires over 108 orthogonal SDs, and thus the reduction in the numbers of SDs is a significant advantage of adopting nonorthogonal SDs. The ground-state energy obtained by the proposed method does not depend on the components of the correction vectors; however, the rate of convergence does depend on the number of employed correction vectors N.

Figure 3

Convergence performance of the proposed method for the correlation energy. Convergence performance of the proposed method for the correlation energy of a HF molecule with the 6-31G** basis set as a function of the number of employed SDs is shown.

The potential energy curve calculated when a single H atom is extracted from a CH4 molecule as shown in Figure 4. Calculations are performed using the 6-31G* basis set. Although the bond lengths are close to the equilibrium one, the errors in the energies obtained by coupled-cluster theory with singles and doubles (CCSD) plus perturbative triples (CCSD(T)) are a few milliHartree; at longer bond lengths, the accuracy of the results appears to deteriorate [42]. In contrast, the proposed calculation procedure ensures essentially exact ground states at all bond lengths, since no approximations are employed.

Figure 4

Potential energy curve of a CHmolecule obtained using the proposed algorithm with 6-31G* basis set.

Figure 5 illustrates the potential energy curve along the symmetric stretching coordinate of a H2O molecule in the 3-21G basis set. The angle between the O-H bonds is fixed at 107.6°. These results shown for the proposed calculation method, CCSD and CCSD(T) exhibit the same trends as for a CH4 molecule. The results for near the equilibrium bond length demonstrate comparable performance between the four methods, whereas results for long bond lengths indicate only that the proposed method has comparable performance with full CI not producing the same unphysical energy curves as CCSD and CCSD(T) around 2.3 Å [42].

Figure 5

Potential energy curve of a HO molecule obtained using the proposed algorithm with 3-21G basis set.

Conclusions

A reliable and tractable technique for constructing the ground-state wave function by the superposition of nonorthogonal SDs is described. Linear independent multiple correction vectors are employed in order to update one-electron wave functions, and a conventional steepest descent method is also performed as a comparison. The dependence of convergence performance on the number of adopted correction vectors is also illustrated. The electron–electron correlation energy converges rapidly and smoothly to the ground state through the multi-direction search, and an essentially exact ground-state energy is obtained with drastically fewer SDs (less than 100 SDs in the present study) compared with the number required in the full CI method. For the few-electron molecular systems considered in the present study, essentially exact electron–electron correlation energies can be calculated even at long bond lengths for which the standard single-reference CCSD and CCSD(T) show poor results, and the practicality and applicability of the proposed calculation procedure have been clearly demonstrated. In future studies, calculations employing periodic boundary conditions and effective core potentials (ECPs) [43] will be performed. A new procedure to reduce the iteration cost should be found in order to increase the applicability of the proposed algorithm for the calculation of essentially exact ground-state energies of many-electron systems.

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

HG conceived, planned this study, carried out the coding of the computation program, and drafted the manuscript. MK and KH participated in the discussions on the basic theory of the present method. AS performed tunings of the code and made all of calculations. All authors read and approved the final manuscript.