Optimizing single Slater determinant for Electronic Hamiltonian with Lagrange multipliers and Newton-Raphson methods as an alternative to ground state calculations via Hartree-Fock self consistent field
Considering the emblematic Hartree-Fock (HF) energy expression with single Slater determinant and the ortho-normal molecular orbits (MO) in it, expressed as a linear combination (LC) of atomic orbits (LCAO) basis set functions, the HF energy expression is in fact a 4 th order polynomial of the LCAO coefficients, which is relatively easy to handle. The energy optimization via the Variation Principle can be made with a Lagrange multiplier method to keep the ortho-normal property and the Newton-Raphson (NR) method to find the function minimum. It is an alternative to the widely applied HF self consistent field (HF-SCF) method which is based on unitary transformations and eigensolver during the SCF, and seems to have more convenient convergence property. This method is demonstrated for closed shell (even number of electrons and all MO are occupied with both,  and  spin electrons) and restricted (all MOs have single individual spatial orbital), but the extension of the method to open shell and/or unrestricted cases is straightforward.
INTRODUCTION<!>Energy optimization with Lagrange multipliers<!>APPENDIX<!>Appendix.2 The ratios of energy terms:<!>Appendix.4 Approximate Coulomb energies:
<p>The widely applied standard HF-SCF approximate solution [1-3] of the non-relativistic electronic Schrödinger equation for an N-electron molecular system (containing A=1,…,M atoms with ZA nuclear charges at positions RA [4][5][6]) approximates the total ground state electronic energy with single Slater determinant 0 S0), and the energy expression for closed shells with restricted MOs in S0 is E0,total HF-SCF = 2 i=1 N/2 (i|h|i) + ai=1 N/2 j=1 N/2 (2(ii|jj)-(ij|ji)) + Vnn (1) (see p.68, equation 3.127 on p.134 in) ref. [1]. For example, using i(1) i(r1) and j(1) j(r1) spatial functions for pairwise occupation: if N/2=1  S0= |1i(1), 2i(2)>= ()i(1)i(2), if N/2=2  S0= |1i(1), 2i(2), 3j(3), 4j(4)>, generally i(1) i(r1) the i th MO, with i=1,2,…N/2 as a function of position vector r1=(x1,y1,z1) of electron 1, as commonly abbreviated in computational quantum chemistry; unusual, (but useful in sums) is the function notation i(1) instead of f(r1). The h(r1)= -(1/2)1 2 -A=1…M ZA/RA1 with sum and RA1 |r1-RA| is the kinetic energy (T) plus Coulomb nuclear-electron attraction (Vne) operator, the ij term approximates the Coulomb electron-electron repulsion energy (Vee (N(N-1)/2)0 * 0r12 -1/2 n=1 N dsndrn, for which a particular closed shell restricted MOs case is written out), and Vnn= A=1…MB=A+1…M ZAZB/RAB is the Coulomb nuclear-nuclear repulsion energy. The parameter 'a' in Eq.1 is the "coupling strength parameter": a=1 has a physical meaning, while a=0 means that there is no electronelectron interaction useful for mathematical purposes [7][8]. Generally, (i|h|j) i * (1)hj(1)dr1 and (ij|kl) i * (1)j(1)r12 -1 k * (2)l(2) dr1dr2, which simplifies a bit for real functions (i(1) * =i(1)) what one has for MOs, as well as the integration by parts provides for u=x,y,z: -i(1)( 2 /u1 2 )j(1)dr1= (i(1)/u1)(j(1)/u1)dr1, because the MOs are well behaved (square integrable and tend to zero as |r1|). The ground state oneelectron density (, with sn=  or  spin state of electron n=1…N) and its closed shell approximation is (r1)N0 * 0 ds1n=2 N dsndrn and  HF-SCF (r1) i=1 N/2 i(1) 2 with i(1)= k=1…K cikGk(1).</p><p>(</p><p>For open shell case (N is odd or  not-closed MO), the energy expression is similar to Eq.1 and can be discussed analogously; (1) is an approximate for , e.g. (1):=  HF-SCF (r1) or else (e.g. Appendix.4).</p><p>The LCAO approximation for the i=1,2…N/2 MOs is given in Eq.2, where the set {cik} contains the LCAO coefficients to the basis set {Gk}k=1,…K, and K ≥ N/2 must be. With MO energy, increasing upward, picture the known arrangement of the rectangular matrix of LCAO coefficients as The number of elements in set {cik} is NK/2. Generally, the atom centered Gk(1) (x1-RAx) nx (y1-RAy) ny (z1-RAz) nz exp(-b|r1-RA| m ) Gaussian type orbits (GTO, m=2) for which analytical integration is available, or the more effective, but not easy to handle Slater type orbits (STO, m=1) are chosen. (Contracted GTO are also used as Gk, as more powerful basis set, etc. not detailed here, also the concepts like minimal basis, STO-3G (STO approximated with three GTO), etc..) For the demonstration we use all the atomic STO basis functions (p.88 and pp.92-94 in ref.</p><p>[9]) up to (n,l)= (3,+2) quantum numbers for all individual atoms up to ZA=18 (Ar), (non-relativistic cases), and numerical integration for all cross products (see below).</p><p>For initial LCAO values in the iteration for Eqs.4-5 below, one can use the sophisticated "Harris approximation" used in HF-SCF practice, or simply, mainly because stationary (equilibrium geometry) close to neutral molecules are in focus, our choice is as follow: e.g. for a quasi-neutral C atom (NA=6 from N= A=1 M NA) in molecule the cik= 1.0 for occupied 1s, 2s, the cik= 1.0 (or the finer 2/3) for partially occupied 2px, 2py, 2pz and cik= 0.0 for higher, unoccupied excited AOs (3s, 3px, etc.), based on the configuration of ground state atomic C (1s 2 2s 2 2px 1 2py 1 ). Based on some concepts of chemical bonds (core electrons, etc.), this means that in a proper matrix arrangement (Appendix.1) the cik= (i,k)= 1 if i=k and 0 if i≠k, the Kronecker delta (see the values 1 as initial values for the {cik} matrix above), and some cii=1 are zeroed out as explained in the case of a neutral C atom; a chemically plausible choice for initial {cik}. In case of a well chosen "minimal basis", the {Gk} can be partitioned to atoms it centered as UA=1 M {Gk(centered on A)}, wherein all subset contains ortho-normalized AOs (that is, solutions of the Schrödinger equation with M=1, that is, these initial {cik} satisfy all the diagonal and some off-diagonal equations in Eq.5 below), and as seen in the practice, in this hierarchy, the values in minimizing set {cik} are in or not far (Apendix.1) from the interval [-1,1], a convenient property for the iteration below.</p><!><p>The idea is based on the fact that inserting the LCAO approximation for i(1) in Eq.2 into Eq.1 yields a 4 th order multivariable polynomial of {cik} if a=1, 2 nd order only if a=0, and the Lagrangian (keeping the</p><p>(3) The ij are the Lagrange multipliers [10], only an upper diagonal matrix (i.e. j runs in [i, N/2] only), and (i,j) is the Kronecker delta. Picture the multipliers {ij} as an upper triangle, square matrix as</p><p>The number of elements in set {ij} is the number of diagonal plus upper diagonal elements in an (N/2)x(N/2) matrix: N(N+2)/8, the proof is elementary. Energy minimization is to find the minimum of L with respect to parameters {cik} and {ij}, where i,j=1…N/2 with j≥i runs for MOs and k=1…K runs for basis set functions. The polynomial order of parameters {LCAO coefficients}U{Lagrange multipliers}= {cik}U{ij} of four main terms in Eq.3 is 2 nd ,4 th , 3 rd and 0 th , the last Vnn term is only an additive constant.</p><p>The minimum of L is obtained by solving the nonlinear system {L/cik=0, L/ij=0 for all i,j,k}, e.g. with NR method (p.263 of ref.</p><p>[11]), for which the second derivatives of L, (the  2 L/cikcjk',  2 L/iji'j' and  2 L/ciki'j) are also necessary. The derivation is straightforward, but can be done simpler if the double sums are shifted inside to have terms like i=1 N/2 i(1) 2 j=1 N/2 j(2) 2 along with (i=1 N/2 i(1) 2 )/cik= 2i(1)Gk, etc.:</p><p>(5) Notice the hectic double indexing in ij in Eq.4: e.g. 21,22,23 for i=2 and N/2=3, it should be considered as 12,22,23, since j≥i; i.e. picking from upper diagonal. (Careless doubling the terms with ji beside ij is not allowed, because that would cause singular matrix in Eq.13 below. For the Jacobian, hij/i'j'=0, hij/cik=Gk(1)j(1)dr1 and hij/cjk=Gk(1)i(1)dr1 if j>i, hii/cik=2Gk(1)i(1)dr1, hij/ci'k=0 if i≠i'≠j, and similarly,  2 L(a=0)/cik 2 = gik(a=0)/cik= 2Gk[2h+ii]Gkdr1, etc..) The initial parameters for the set {ij} can be calculated from Eq.4 using the initial parameters chosen for the set {cik} above, because Eq.4 is a simple linear system for {ij}, as well as notice that, Eq.4 has more equations than necessary for this step: Pick N(N+2)/8 from NK/2. The grad(L)= (L/c11,..., L/N/2,N/2) via Eqs.4-5 shows the opposite direction in which L decreases most quickly, and |grad L| determines how fast the L changes in that direction in {cik}U{ij} space.</p><p>We minimize Eq.1 first with a=0 via the LCAO coefficients {cik}, that is solving Eqs.4-5 for zero 1 st derivatives with NR method (with the help of Eqs.7-13 below) for {cik}U{ij}, and use it as initial parameters for the wanted case when a=1 (but doing the same computation procedure). The reason for these two steps is that Eq.3 is only 2 nd order in {cik} and 3 rd order in {cik}U{ij} if a=0, and the zero 1 st derivatives (Eqs.4-5) can be found in stable ways, as well as the stationary set {cik} does not change too much [7][8] in a[0,1] (see Appendix.1). Strictly saying, the "perturbation theory" comprises mathematical methods for finding an approximate solution to a problem, by starting from the exact solution of a related, simpler problem, and in fact this pre-calculation trick with parameter 'a' can also be considered as a perturbation, see Appendices 1-2. A weaker approximation [1-2] than the fully minimized Eq.1 with a=1 can be calculated by LCAO coefficients from Eqs.4-13 at a=0, and using directly in Eq.1 along with a=1. We focus on ground state here, but the MOs at a=0 can be used as a basis for "configuration interactions (CI)" methods for ground-and excited states, wherein the off-diagonal elements in the CI matrix depend on operator r12 -1 only, for this, one must introduce at least one/two additional virtual (1+N/2) th electrons [7][8].</p><p>To solve the system in Eqs.4-5, we need the first derivatives for gik and hij, which is straightforward. The hij/i'j'= 0 since L is linear in ij, a lucky situation, but the two groups of functions (g and h) with double indexing (ik and ij) along with j≥i is hectic, so the transformation of indexing is useful, see below. We see from Eqs.3-5 that the powers build up as L= cikci'k' + cikci'k'i''j + acikci'k'ci''k''ci'''k''' (6) gik= cik +ciki'j + acikci'k'ci''k'' and hij= cikci'k' (7) where the symbolic  means proper LC with indices run, and the coefficients in these individual LCs can be mixed up from the elementary integrals to be calculated in Eqs.3-5, which are the same as in standard HF-SCF procedure: Two center integrals Gk(1)Gk'(1)dr1 from Eq.5 and the two and three center</p><p>(first and second term in the right, resp.) integrals from Eq.4 along with the four center integrals  GkGk' Gk''Gk'''r12 -1 dr1dr2, (9) where Gk's are primitive (or contracted) Gaussians, see Appendix.3-4. Eq.7 shows the degree of multivariable polynomials (g is 3 rd order, h is 2 nd order) and the missing powers (g has no squares like cikci'k', h has no linear terms like cik or ij etc.). In the first step calculation with a=0, the system in Eq.7 is second order (parabolic) only to solve for minimizing L in Eq.3. (In the 1 st sum in Eq.6 i=i' in cikci'k' by Eq.3, etc..)</p><p>The set {cik}U{ij} of variables in the multi-variable polynomials in Eqs.4-5 is not easy to handle because the two indices, again, i=1…N/2, j=i…N/2 for (i) and (j) MOs where the number of electrons N is even, as well as k=1…K for basis set functions {Gk}. It is useful to transfer them to one index variables: Let fn= gik for functions and xn= cik for LCAO coefficients, where n=(i-1)K + k (10) ran as n=1, 2, …NK/2. For the other functions hij and Lagrange multipliers 11,…, 1,N/2, 22…, 2,N/2, 33…, 3,N/2, …, N/2,N/2 the n in fn= hij and xn= ij continues as n= NK/2+1, NK/2+2, NK/2+3……P NK/2+N(N+2)/8= (N/2)(K+N/4+1/2), (11) respectively. (N(N+2)/8 is the cardinality of the set {ij}.) Eqs.10-11 also show that the number of variables (dimensionality) in the 4 th order polynomial L(c11,…,N/2,N/2)=L(x1,…,xP) in Eq.3 is P in Eq.11 with even N. The left hand side of the equation system in Eqs.4-5 is a zero column vector, but the i=j=1,…N/2 cases of Eq.5 rearrange it as 1= i(1) 2 dr1, providing non-singular system. Finally, the system in Eq.7 transfers (with the analogue meaning for symbolic ) to 0= fn= xn1 + xn1 xn2 + axn1 xn2 xn3 (12) where fn= fn(x1,….xP) with n, n1, n2, n3=1…P, with P in Eq.11. This transfer of indices for Eq.7, which can even start at Eq.6 as L= xn1 xn2 + xn1 xn2 xn3 + axn1 xn2 xn3 xn4 with 0= fn L/xn, avoids the caution j≥i for ij mentioned after Eq.4. In this way, the polynomial generation in Eq.12, for example, for the term anx1 2 x4 comes from b114x1x1x4+ b141x1x4x1+ b411x4x1x1 with an=b114+b141+b411. This transfer of indices in Eqs.10-11 can be conveniently done with e.g. FORTRAN "do-loops" in the programming.</p><p>Let [xn(0)] be the column vector of the initial values for xn, the function value at this point is the column vector [fn( 0 (13) It avoids the eigensolver and unitary transformation what e.g. HF-SCF uses to find the minimizing LCAO coefficients. Eq.13 should be applied for pre-optimization (a=0 in Eqs.4-5, 2 nd order polynomial) first, then a re-optimization for a=1 in Eqs.4-5 (4 th order). The choices for [M] in Eq.13 are: 1.: Unit matrix, called "iterative method", however, a strong condition in this case requires the initial [xn(0)] to be very close to the optimum, generally not feasible (may be in the re-optimization step), 2.: The NR as [M(m)]= -[W(m)] -1 , where [W] [fn1/xn2] is the Jacobian matrix of Eqs.4-5 (or Hessian of L in Eq.3), 3.: -1/[fn1/xn1]m diagonal and zero off-diagonal elements in [M(m)], called "diagonal NR", particularly,  2 L/ij 2 =0 from Eq.5 does not allow this here, 4.: The fast "gradient (steepest descent)" method converging along -grad</p><p>fn], (even matrix inversion is not necessary, only transpose, more, since W is Hessian  W is symmetric  W=W T ), 5.: The "diagonal gradient" neglects the off-diagonal part of W yielding b=(n=1 P (fn/xn) 2 fn 2 )/(n=1 P (fn/xn) 4 fn 2 ), but with slower convergence, however, the sum in the denominator avoids the problem of e.g.  2 L/ij 2 =0. Notice that if P=1, the "diagonal gradient" reduces to b=(f1/x1) -2 and 1x1 size [W]=f1/x1, so -b[W] T [f1]= -f1/(f1/x1), the 1 dimensional NR. In cases 1, 3 and 5, even the transfer of variables (Eqs.10-11) would not be necessary and the use of PxP matrix [M] in Eq.13 would reduce with the use of Px1 column vector [fn].</p><!><p>Appendix.1 The initial LCAO coefficients: For example, in case of LiH molecule, the HF-SCF expands the 2x6 matrix {cik} to 6x6 matrix (virtual orbits) and uses unitary transformation, etc., while the method with Eqs.3-5 uses the 2x6 matrix {cik} with 2x2 upper triangle matrix {ij}: Simple initial LCAO (atomic H config.: 1s 1  c12:=1.0, atomic Li config.: 1s 2 2s 1 c21:=1.0, N=4):</p><p>basis Li(1s) H(1s) Li(2s) Li(2px) Li(2py) Li(2pz) 2.MO:</p><p>1.0 0 0 0 0 0 1.MO: 0 1.0 0 0 0 0 Converged LCAO by HF-SCF/STO-3G/a=0 for Eq.1 (E0,total= -11.456970 hartree, hypothetic bound state): 2.MO:</p><p>1.00550 0.00435 -0.02504 0.0 0.0 -0.00159 1.MO:</p><p>-0.10689 0.66300 0.30999 0.0 0.0 -0.32185 Converged LCAO by HF-SCF/STO-3G/a=1 for Eq.1 (E0,total= -7.860313 hartree, real bound state): 2.MO:</p><p>-0.99129 -0.00349 -0.03290 0.0 0.0 -0.00603 1.MO: 0.16462 -0.55155 -0.45875 0.0 0.0 0.34460 Notice that, the two sets of converged {|cik|} are similar (apart from -1 phase factor).</p><!><p>The term E0 E0,total-Vnn= T + Vne + Vee is called electronic energy in a real case (a=1 in Eq.1), and the ratios in the right side is interesting in relation to the parameter 'a'. With the help of small H2 and larger C10H8 equilibrium geometry molecules we demonstrate that, the magnitude of V&T energy ratios in Eq.3 does not depend strongly on N. The calculation below represents the involvement of Vee interesting for convergence in view of manipulation with parameter 'a' in Eq.4; notice that, operator 'h' in Eq.3 counts for T+Vne, the |Vne| > T, Vee > 0 and Vne < 0 always in real systems, as well as the virial ratio -(Vne+Vee+Vnn)/T=2 has the strict value 2 (a=1 or a≠1), not the listed ones: 2)=K(K-1)/2) cross products for all possible combinations of basis function pairs (Gk, Gk'), but in case of Coulomb energy, the magnitude of the number of cross products, i.e. the four center integrals in Eq.9 is K 4 (more exactly ( K 2)( K 2)), which can be decreased via symmetry since r1 and r2 are equivalent, etc., but the K 4 magnitude is huge: In case of a small system, like water dimer (H2O)2, there are K ~100 contracted Gaussians (using 6-311++G(3d,2p) basis set), so about K 4 = 10 8 cases have to be evaluated. After calculating these cross products in Eqs.8-9 with basis set {Gk}, the stationer LCAO coefficients have to be found (by an HF-SCF procedure or the procedure described here) along with the multiplications and sums in Eqs.3-5 during the iteration, to obtain finally the desired energy (L) in Eq.3 to complete Eq.1. In view of computation, the origin of the difficulty with cross products can be symbolized as follows: Consider the expansion of (1 K ak)(1 K bk')= a1b1 + a1b2 + … aKbK, the left hand side has K-1 additions in both sums plus 1 multiplication, all together 2K-1 operations, while the right hand side has K 2 multiplications plus K 2 -1 additions, all together 2K 2 -1 operations, finally about 2K vs. 2K 2 , and similarly, for products of four sums: 4K vs. 4K 4 . The cardinality value K 2 of cross products in Eq.8 is feasible, but the large K 4 value counting for Eq.9 is a reason to look for good approximations to make a shortcut: Numerical integration in certain circumstances [6] or approximate formulas (Appendix.4) work.</p><!><p>We mention the important point in standard HF-SCF that, in case of GTO functions the terms in Eqs.8-9 can be evaluated analytically (using that for example, products of GTO is LC of GTO's which is not true for STO, etc.). In HF-SCF [1-3], the emblematic Vee J-K= (j(1,2)-k(1,2))r12 -1 dr1dr2 (14) approximation is used, see its particular case in Eq.1 and its derivative with cik in Eq.4. The J and K integrals contain the cross terms in Eq.9, which can be evaluated analytically and before the iteration. (Eq.14 is only an approximation, so it needs "correlation energy (1 %)" calculation, that is, correcting the not-enough single Slater determinant (S0) approximation, however, that is another question); the expression in Eq.14 has been a milestone equation in calculating Coulomb interactions on the quantum level. These famous J and K are called "Coulomb-J-integral" and "exchange-K-integral", resp.</p><p>Another (milestone) approximation (suffering also from the necessity of "exchange and correlation energy (1 %)" calculation) in density functional theory (DFT) [2-3] is Vee (1/2) (1)(2) r12 -1 dr1dr2 .</p><p>(15) Its derivatives (/cik)(1)(2)r12 -1 dr1dr2= 8Gki (2)r12 -1 dr1dr2 modifies Eq.4 accordingly.</p><p>Approximations in Eqs.14-15 also suffer from the above mentioned "K 4 operations" problem in their analytical integration, but these analytical evaluations in Eqs.8-9 are necessary only once before the iteration starts. In ref. [4] the approximation with approximate ground state one-electron density ()</p><p>Vee CJ[ 6/5 dr1] 5/3 + j=1,…nCxj[ [1+1/(3j)] dr1] j (16) is reviewed, wherein the 1 st term accounts for main value, and the rest is correction, even the correlation energy can be included. If Eq.16 is used for Eqs.1-5, the price is that: 1., the integrals can be evaluated numerically only (for both, GTO and STO) e.g. with the scheme used in ref. [6], 2., integrals must be evaluated in every iteration step with the improved {cik}U{ij}, along using (/cik)[ 6/5 dr1] 5/3 ~ [Gki 1/5 dr1] 2/3 continuing with numerical integration in the right, but another choice is a numerical derivation for /cik also, since many similar integrands are needed anyway. The benefit from approximation in Eq.16 is that it belongs to the numerical Vee f((1))dr1 forms, and importantly, in fact the above mentioned "K operations" vs. analytical "K 4 operations" () r12 -1 dr1dr2 in Eqs.14-15, so for larger molecules Eq.16 may be useful, but a less studied area than the widely used and tested Eqs.14-15.</p>