An approach to first principles electronic structure calculation by symbolic-numeric computation

This article is an introduction to a new approach to first principles electronic structure calculation. The starting point is the Hartree-Fock-Roothaan equation, in which molecular integrals are approximated by polynomials by way of Taylor expansion with respect to atomic coordinates and other variables. It leads to a set of polynomial equations whose solutions are eigenstate, which is designated as algebraic molecular orbital equation. Symbolic computation, especially, Gr\"obner bases theory, enables us to rewrite the polynomial equations into more trimmed and tractable forms with identical roots, from which we can unravel the relationship between physical parameters (wave function, atomic coordinates, and others) and numerically evaluate them one by one in order. Furthermore, this method is a unified way to solve the electronic structure calculation, the optimization of physical parameters, and the inverse problem as a forward problem.


INTRODUCTION
This article is intended for an introduction of a new approach to first principles electronic structure calculation by way of symbolic-numeric computation [1].
There is a wide variety of electronic structure calculation cooperating with symbolic computation. The latter is mainly purposed to play the auxiliary role (but not without importance) to the former. In the field of quantum physics [1][2][3][4][5][6][7][8][9], researchers sometimes have to handle with complicated mathematical expressions, whose derivation seems to be almost beyond human power. Thus one resorts to intensive use of computer, namely, symbolic computation [10][11][12][13][14][15][16]. The example concerning this, as is given in the reference 16, ranges various topics: atomic energy levels, molecular dynamics, molecular energy and spectra, collision and scattering, lattice spin models and so on. How to obtain molecular integrals analytically or how to manipulate complex formulas in many-body interaction, is one of such problems. In the former, when one uses special atomic basis for a specific purpose, to express the integrals by the combination of already-known analytic functions may sometimes be very difficult. In the latter, one must rearrange a number of creation and annihilation operators in a suitable order and calculate the analytical expectation value. In usual, a quantitative and massive computation is to follow from a symbolic one; for the convenience of the numerical computation, it is necessary to reduce a complicated analytic expression into a tractable and computable form. This is the main motive for the introduction of the symbolic computation as a forerunner of the numerical one and their collaboration has won considerable successes up to now. The present work should be classified as one of such trials. Meanwhile, the use of symbolic computation in the present work is not limited to indirect and auxiliary part to the numerical computation. The present work can be applicable to a direct and quantitative estimation of the electronic structure, skipping conventional computational methods.
The basic equation of the first principle electronic structure calculations is the Hartree-Fock or the Kohn-Sham equation, derived from the minimum condition of the energy functional in the electron-nuclei system [1][2][3], which is expressed as follows. ).
The second term in the parse of the left side is the potential from nuclei with charge Za.
The third term is the Coulomb potential generated by the charge distribution ρ. The forth term V is the quantum dynamical interaction operating in the many-electron system, called the "exchange and correlation". It is rewritten into the matrix eigenvalue problem, by adopting wavefunctions expanded by a certain localized basis set. This is analytical base, such as STO (Slater type orbital) or GTO (Gaussian type orbital), to construct one-electron and two-electron molecular integrals, whose concrete expression can be derived from symbolic manipulation [4][5][6][7][8][9] by means of computer algebra systems [10,11].
The use of analytic basis in the HFR equation is effective in the achievement in the precision of the numerical computation, but causes some difficulty in the mathematical operations to the equation itself, because the analytic expression is, in general, very complicated. We, however, can rewrite and approximate the HFR equation by polynomials in order to obtain much simpler expressions. The concept of the polynomial approximation to the HFR equation is nourished by Yasui [6][7][8][9]. The equation becomes the set of algebraic polynomial equations expressed by atomic coordinates, orbital exponents, and quantum numbers. Based on this, we will able to unravel the relationship among parameters and clarify their dependence. This idea should be outlined here. We express a molecular orbital by the linear combination of atomic orbitals (LCAO), C χ R , n , l , m , ζ , r, θ, .

(I.3)
The variables R , n , l , m , ζ are the atomic position, quantum numbers, and the orbital exponents respectively. The variables r, θ, are of the atom centered coordinates.
The key components to the molecular orbital calculations are molecular integrals, which are the matrix elements of the each part of the Hamiltonian operator obtained by the use of LCAO, such as, the kinetic energy, the nuclear and electronic potentials, and the overlapping integrals. The approximation to molecular integrals is obtained by Taylor expansion; For example, the two-centered overlapping integral is defined as, S AB χ R A , n , l , m , ζ , r A , θ A , A χ R B , n , l , m , ζ , r B , θ B , B dr . (I.5) The integration generates an analytic function of two orbital exponents and inter-atomic distance R. The polynomial approximation is given as, The other molecular integrals can be expressed in the similar way. Once the molecular integrals are approximated as polynomials, the HFR equation and the energy functional take polynomial expressions. The orbital exponents and atomic coordinates can equivalently be regarded as parameters in the calculus of variations, as well as LCAO coefficients. It will be adequate to call this multi-variable polynomial expression "algebraic molecular orbital equation," or "algebraic molecular orbital theory." It is not necessary to regard those equations as pure numerical eigenvalue problems. Those equations are a set of polynomials, to which both symbolic manipulations and numerical solving are applicable.
One should note some inconveniency in the conventional methods, which may be surmounted by "algebraic molecular orbital equation." The standard electronic structure calculation is a "forward problem". We suppose the material structure, execute the electronic structure calculations, and optimize the structure so that the energy functional will be minimized. The foundation for this treatment is so-called "the adiabatic approximation", which enables us to separate the dynamics of the nuclei and wavefunctions into two independent models, ruled by the classical and quantum dynamics respectively. The conventional method iterates two alternative computational phases, one of which are the optimization for the wavefunctions and the other for the positions of the nuclei. It is believed that this way is numerically stable. But, in view of effectiveness, this may be a lengthy and roundabout one. This also results in some inefficiency. Owing to the separation of the degrees of freedom of wavefunctions and nuclei, it is difficult for the conventional method to coop with cases where the dynamics of nuclei and wavefunctions are strongly coupled with each other. Meanwhile, the "inverse problem" will be to search the material structure which shows the desirable electronic properties. To do this, the conventional method must go with trial and error.
At first we suppose the material structure to evaluate the electronic properties, and, by adjusting the structure, we search the direction in which the desired properties will be obtained. We have to solve forward problems repeatedly to obtain the solution of the inverse problems. The reason to this is as follows. In the conventional methods, the computation has the fixed order of numerical procedures, consisted from the eigenvalue solution, the self-consistent-field calculation and the relaxation of atomic structure, which is implemented as nested loops of independent phases of the optimizations. The unknown parameters are computed from inner loops to outer ones in order. The conventional method is obliged to determine unknown variables in a fixed order in any cases. As we will see later, the concept "algebraic molecular orbital equation" suggests a solution strategy to this circumstance.

METHOD
With the view of these circumstances, we propose the following method, named "Symbolic-numeric ab-initio molecular dynamics and molecular orbital method" [1].
It is summarized as follows. "At first, HFR equation is approximated as a set of multi-variable polynomial equations (Algebraic Molecular Orbital Equation), and by the symbolic computation, it is rewritten into a certain form more convenient for the numerical treatment. The eigenstates are evaluated by the root finding by means of symbolic-numeric procedure. " The question in the present work is how to obtain the numerical solutions of the equations after the polynomial approximation and derive the significant information.
For the purpose of rewriting and solving a set of polynomial equations, the several types of hybrid techniques, so-called "symbolic-numeric solving", are proposed. In them, the symbolic manipulation is applied as a preconditioning toward the set of equations to be solved. The equations are transformed into the other ones which have the same roots, to which the numerical computation will be easy and stable. From the form of the transformed equations, the character of the solution, such as, the existence and the geometrical structure, can be determined. Then the solving process is passed over to the numerical one. For the mathematical background, see ref. [12][13][14][15]. A review of application of symbolic computations in the field of the computational chemistry is given in ref . 16.
In the present work, we make use of the symbolic-numeric solving and rewriting HFR equation, approximated as a set of polynomial equations. As a strategy, the algorithm of the "decomposition of polynomial equations into triangular sets" is applied [12,13]. In this algorithm, the following transformations are applied.
The starting equations , … , ,..., , ( valuables to ones with more. However, the total number of polynomials in the Gröbner bases may grow more than that of the starting polynomials. Though we can search the root at this stage, we furthermore apply the decomposition to the Gröbner bases and obtain several "triangular" systems of equations { , … , }, the first entry of which has one unknown x , the second has two unknowns x , x , the third three unknowns, etc, until, the last n-th has n-unknowns x , x , … , x by turns. In order to obtain all roots of the starting equations , … , , we may need to construct several sets of triangular set of equations, all of which can be generated by the algorithm in ref .12 or 13. Single triangular system has fewer numbers of entries than that in the Gröbner bases before the transformation. This makes the numerical procedure easier. Once we can obtain the triangular systems of the equations, we can evaluate each unknown one by one. In the numerical solution, only a Quasi-Newton-like method, or its kindred for one variable, is necessary. We can execute another type of the symbolic numerical solving. The foundation to this is the theorem of Stickelberger. As above, we regard the HFR equations as the set of polynomial equations expressed by unknowns X1,…,Xm. The set of polynomial equations constructs a zero-dimensional ideal I in the polynomial ring R=k[X1,…,Xm], whose zeros corresponds to a residue ring A=R/I. Here k is the coefficient field, which, in our cases, is the rational number field or the real number field. The ring A is a finite dimensional vector space over k, whose bases are expressed by monomials of X , X , … , X . Thus, the multiplication by X , X , … , X on each base results in the linear combination of the bases in A. The product operations by X , X , … , X are expressed as linear transformation matrices m h X , X , … , X ). The bases in the residue ring A=R/I and the transformation matrices are obtained by means of Gröbner bases technique.
The theorem of Stickelberger asserts that there is a one-to-one correspondence between an eigenvector of the matrix m and the zero point ξ ξ , … , ξ of the ideal I. The correspondence is given by The one of the merits in this treatment is as follows. In the conventional method, the input data is the atomic structure and the output if the electronic structure. By contrast, in the present method, the possible input data are not limited to the atomic structures.
We can select arbitrary parameters in the HFR equation and set them as the input. If the problem to be solved is properly established, we can compute other unknown variables properly. As to the properness of the problem, i.e., the existence of the roots of the set of polynomial equations, it can be judged from the ideal theory in mathematics on the condition whether its Gröbner bases have zero points set or not.

COMPUTATIONAL FLOW
The task flow is listed as follows.
It is the two electron repulsion between 1s orbitals and classified as "Coulombic type", denoted as 1s A 1s A |1s B 1s B . Here the notation "1s (A)" and "1s (B)" mean the atomic orbitals centered on atoms A and B. The Slater orbital takes a form r; z / π / e . There are other repulsion integrals, classified as the "exchange type"

[1s(A)1s(B)|1s(A)1s(B)] and the "hybrid type" [1s(A)1s(A)|1s(A)1s(B)
]. In general, the molecular integrals have more complicated expressions than (R.1). If STO is used, these integrals take more lengthy expressions, including transcendental functions, such as exponentials or exponential integrations, and infinite series summations [5]. In spite of this complicacy, we can treat them easily after the symbolic processing, by approximating them as finite degree polynomials by means of the Taylor expansion. It is noted here that the STO base can describe the physical property of the localized atomic wavefunction more precisely than by GTO (Gaussian Type Orbital) base, both in the neighborhood of the nucleus and in the remote region from it. Thus the STO base becomes more advantageous for the purpose of expressing the molecular equations as the polynomials of the atomic coordinates. This is the reason why STO is adopted here.
However, the following recipes are also applicable to the GTO calculations, and possibly, to semi-empirical calculations, such as AM1 (Austine Model 1) [17] and PM3 (Parameterized Model number 3) [18], or tight-binding model, where the matrix elements are given as analytic formulas.
As an example of a forward problem in the first principles molecular dynamics, the optimization of the structure (the distance between two hydrogen atoms) and the UHF electronic structure calculation are simultaneously executed.
We execute the UHF (Unrestricted Hartree-Fock) calculations in which the trial wavefunctions for up-and down-spins are defined as in (R.2) and (R.3). The corresponding eigenvalues are denoted as ev and ew. It is noted here: these trial functions with the bases of the orbital exponent 1 are too primitive to assure good agreements with experiments: for accuracy, we must optimize the orbital exponent to a suitable value. It is only to reduce the computational cost in the symbolic computation why we adopt such primitive trial functions. Here, the positions of hydrogen A and B are denoted as R A and R B . The inter-atomic distance is r R A R B . We use the notation, such as x A |x R A | and x B |x R B |.
The energy functional is transformed into the polynomial form by way of the fourth order Taylor expansion of the inter-atomic distance r, centered at the position of R0=7/5 atomic unit. The functional is generated in a standard way of molecular orbital theory, which is the total energy of the electron-nuclei system (given in the atomic units) with the constraint condition of the ortho-normality of the wavefunctions. The Lagrange multipliers are eigenvalues. The coefficients of real numbers are truncated to third decimal places and approximated as rational numbers, as is shown in (R.4). (a,b) are the LCAO coefficients for the up-spin electron,(c,d) are those for the down-spin electron, ev is the eigenvalues of the up-spin electron, ew is that of the down-spin electron, r is the inter-atomic distance. It is rather a rough approximation to use numerical coefficients truncated to third decimal places, which causes the computational error (the order of a few percents) by the present method, compared to the conventional way with the double precision calculation. This approximation is only intended to reduce the computational cost in the symbolic processing. In this sense, examples in this section are mock-ups for realistic calculations, the aim of which is to illustrate the application of the present method, aside from the accuracy.
Then the HFR equation is given by the set of equations in (R.8), where (t,s) is the LCAO coefficient for up-spin, (u,v) is that for down spin, ev is the eigenvalue for up spin, ew is that for down spin, and r is the inter-atomic distance.
The entry as is shown by   The triangular decomposition to (R.10) is shown in (R.11), which involves five decomposed sets of equations. One decomposed set includes seven entries, into each of which, the seven variables is added one by one, with the order of r, ew, ev, v, u, t, s. There are numerical coefficients that are very lengthy ones. This is due to a problem in the algorithm in the Gröbner bases generation [19][20][21]. The computational procedure applies the Buchberger's algorithm, in which the addition, subtraction, multiplication and division are iterated to the polynomial system. In the intermediate expression through the computation, some polynomials with huge degrees may arise, some of whose coefficients have extreme difference in the numerical scale compared to others.
This difference in the scale of coefficients will remain in the final result. To assure the numerical accuracy, we must resort to the computations with the arbitrary precision.
Though the solutions include complex ones, the physically admissible real solutions are shown in Table. 1. We obtain four combinations, where the two electrons of up or down spins are located the symmetric or asymmetric wavefunctions. This means we obtain both of the ground and the excited states.
Atomic and electronic structural optimization executed simultaneously The inter-atomic distance r and the wavefunctions can be optimized at the same time, as is done in Car-Parrinello method. In this case the inter-atomic distance r is an unknown to be determined, not being a fixed constant.
In order to obtain the ground state alone, we add eq. (R.12) into the equations in (R.8).
(This is the trick applicable to this example only. In general, the ground state is given as a solution where the sum of the total occupied eigenvalue becomes minimum one. To specify the ground state, it is enough to compute eigenvalues alone. For this purpose, in making the triangular decomposition of the equations, we can prepare the equations including only eigenvalues as unknowns. We have only to solve them.) With this treatment, we can replace the equation to be solved with a simpler one. The part of the equations including r is shown in eq. (R.13), whose real solutions are shown in Table 2.
The discrepancy between the solution and the experimental value ( r ～1.4 ) is due to the numerical error caused by the roughness of the fourth order Taylor expansion and the rationalization of the numerical coefficients, being truncated. In addition, it is also due to the not-optimized orbital exponent in the trial wavefunctions.

A recipe for an inverse problem
It is demonstrated here how to solve a kind of inverse problem. Suppose a problem, where the energy difference between the occupied and the unoccupied states has a certain value; we should evaluate the inter-atomic distance r at which the energy difference shows this value. This example is a miniature of the inverse problem to find the lattice constants at which the band gap shows the desired width. In this case, we execute RHF (Restricted-Hartree-Fock) calculations. The wavefunctions of the occupied and the unoccupied are given in (R.14) and (R.15). The eigenvalues are denoted as eocc， The required equations are presented in (R.16), whose details are omitted here. The set of the equations for the occupied state can be obtained by the same way as was done as the example of UHF calculation. The orthogonality condition to the occupied and the unoccupied states is added to it. Then, in this case, is the structure is stable? If we evaluate the inter-atomic forces, it can be easily judged. On the other hand, with the view of symbolic-numeric solving, we can make use of the following judgment. To do this, in the set of the equations in (R. 16) we insert the condition of eq. (R.17), the minimization condition of the energy functional with respect to the inter-atomic distance r.  The transformation matrix corresponding to the multiplication by the variable t is  Table 2, given as the previous example.

DISCUSSION
The advantages by the present method are recapitulated here. The algebraic molecular orbital equation, expressed as polynomial sets, informs us the relationship among unknown variables. Thus, in order to evaluate those unknowns, we can divide suitable parts of them into the prepared inputs and the expected outputs, respectively. The calculation is not confined to the conventional framework, such as, whose the input is the structure and the output is the electronic the states. The distinction between the forward and the inverse problems are cleared away, and we can treat all of them as the forward problems in a unified way. In order to cope with the inverse problem, we should check whether it is well-posed or not. The present method affords us a key to this. After the transformation from the fundamental equation to the In the computations in previous section, we fixed the orbital exponent z at unity and converted molecular integrals into polynomials of inter-atomic distance R alone. The algebraic molecular equations can be extended to a more general case. We can make multi-variable Taylor expansions for the atomic coordinates and orbital exponents in order to prepare polynomials of those two kinds of variables. If this is done, the equation includes both parameters and we could optimize atomic coordinates and orbital exponents simultaneously.
It should be admitted that the present work lies in a primitive stage. At present, this study does not necessarily afford us sufficiently precise calculation. One reason to this is the constraint by the ability of the hardware. The cost in the symbolic computation sometimes tends to be so massive (especially on memory usage) that we are obliged to use simplified molecular integrals, reduce the degrees of the approximating polynomial and rationalize numerical coefficients in lower accuracies. As for the overall computation cost in the explanatory calculations in the previous section, the generation of the molecular integrals, especially that in two-electron repulsion integrals are the most demanding part. In the desktop pc (2.0GHz dual core CPU, 2.0GB memory), in the proposed now [19][20][21]. For example, Brickenstein proposes "slimgb" algorithm in order to keep the intermediate expressions as slim as possible, by regularly replacing a swelled polynomial with a shorter, equivalent one [19]. Lichtblau, based on an empirical study, discusses on important points which should be handled carefully in using approximate arithmetic for coefficients [20]. Arnold presents modular algorithms for the purpose of limiting the enormous growth in rational-numbered coefficients [21]. Those strategies, as well as other ones with the same intention, will be of importance in more complex calculations.
Other difficulties will arise in the application to more complex systems.
Firstly: the generation of the molecular integrals must be burdensome. The shown example, hydrogen molecule is very simple one. We need only one-or two-centered integrals. The molecular integrals are computed from the 1s orbital alone, so that the integrals take the simplest expressions. For a more complex molecule, we must use more general, more complicated atomic basis, which give much more complex molecular integrals. Besides, we have to evaluate three-or four-centered two-electron repulsion integrals. There is a problem: as for the analytic expression of multi-centered repulsion integrals, the derivation by the STO is a formidable one. From this reason, the use of STO has been a minority in quantum chemistry. But, since there are several merits in STO as stated before, the studies to obtain analytic expressions for multi-centered repulsion integrals by means of STO are still pursued [22]. In the present stage, it is more pragmatic to employ GTO, the de facto standard, in which the analytic technique for the derivation of multi-centered integrals has been well established.
Secondly: in the explanatory numerical computations given above, we can evaluate both ground and excited states. Looking from a different angle, this means that we must confront all possible electronic configurations and must single out the necessary one from the solutions. It contrasts to the conventional computations, where it is a trivial work to get a ground state; one has only to compute the necessary number of lower occupied eigenvalues in ascending order from the bottom and put an electron in each of them; in the iterative minimization of the total energy functional, one can gather occupied states only, omitting the calculation of unoccupied states. On the other hand, in the polynomial equations, if one solves them without care, the whole eigenstates (including both of necessary and unnecessary ones, and as many as the Hamiltonian dimension) will appear in the set of solutions indiscriminately. The number of the unnecessary configurations shall grow enormously in more complex systems, where more electrons and a larger basis set are included. In order to sift out the ground states, we can pick up the configuration in which the total sum of the eigenvalues in the occupied states becomes minimal, by firstly solving the equations for energy spectrum obtained by polynomial triangulations. Maybe this tactics will be of use for sorting and indexing excited states. Also, to impose certain point group symmetry on wave function will be effective in the reduction of the computational cost. The RHF structural optimization for hydrogen molecule in the previous section, in which the ground state is extracted from the symmetrized wavefunction, is a clumsy example of this. The switching between symmetrized and asymmentrized wavefunctions leads to the selection of ground or excited states.
Thirdly: the polynomial approximation to the exact functional causes nonsensical solutions, which should be checked with care. In more complicated cases, this situation will be more troublesome. The admissible solutions must lie in a range where the polynomial approximation through the Taylor expansion is quantitatively valid: if the solution appears to be dubious, one must re-examine it by another polynomial approximation with a different center point of the Taylor expansion. One can also use polynomial of higher degrees for this purpose.
In summary, the present work shows that the concept of the "molecular orbital algebraic equation" by means of the "polynomial approximation to molecular integrals" is applicable to the realistic first principles electronic structure calculation, as well as its potentiality to several fundamental problems which are difficult to be handled by the conventional method. At the same time, we must recognize that the hardship to be surmounted is large. However, the improvement in the computer architecture is so rapid that we can expect the achievement of the sufficient accuracy by the present method and its application to complex and large material in future, spurred by the refinement of the symbolic computation theory.   Using these integrals, together with LCAO coefficient, we can construct the total energy functional for hydrogen molecule. By means of Taylor expansion at some inter-atomic distance R0, and setting z=1, we can make the polynomial approximation.
In the polynomial approximation, some of numerical coefficients take lengthy expression constructed of rational and irrational numbers. Thus we should re-express coefficients in decimal numbers, truncate them, and express them in rational numbers of moderate sizes, so that we can manipulate them in a modest cost in the symbolic computation. integer, expressed as powers of ten will suffice.) This integer-multiplied energy functional is hereafter processed by SINGULAR package, by which the polynomial equations are generated through symbolic differentiation of the energy functional and numerically solved. In the next section, an example of SINGULAR script is given. Overlapping: A | A =1. (E^(2*R*z)*(R^20*z^19)) + (121455523460201203410637500*