Taisuke Ozaki (JAIST) and Hiori Kino (NIMS)
Atomic orbitals as a basis set have been used for a long time in the electronic structure calculations of molecules and bulks. Especially, in covalent molecular systems, oneparticle wave functions are well described by a linear combination of atomic orbitals (LCAO) because of the nature of localization in the electronic states, which is a reason why chemists prefer to use the atomic orbitals, e.g. Gaussian orbitals.[1,2,3,4] On the other hand, in the solid state physics, LCAO has been regarded as a somewhat empirical method such as a tool for an interpolation of electronic structure calculations with a high degree of accuracy.[5] However, during the last decade, LCAO has been attracting much interest from different points of views, since great efforts have been made not only for developing O() methods of the eigenvalue problem,[6,7,8,9,10,11] but also for making efficient and accurate localized orbitals[12,13,14,15,16] as a basis set being suitable for O() methods to extend the applicability of density functional theories (DFT) to realistic large systems. Most O() methods are formulated under an assumption that a basis set are localized in the real space.[11] Therefore, the locality of the atomic orbitals can be fully utilized in DFT calculations coupled with O() methods. In addition, even if a minimal basis set of atomic orbitals is employed for valence electrons, it has been reported that a considerable accuracy is achieved in many systems.[12,17,18,19,20,21] This fact suggests that the matrix size of the eigenvalue problem is notably reduced compared to other localized basis methods such as finite elements method.[22,23] These aspects of LCAO encourage us to employ the atomic orbitals in the large scale O() DFT calculations.
However, several important problems still remain in the applications of LCAO to DFT calculations in spite of these advantages. The most serious drawback in LCAO is the lack of a systematic improvement of atomic orbitals in terms of the computational accuracy and efficiency. One expects that a basis set such as double valence orbitals with polarization orbitals for valence electrons provides a way for balancing a relatively small computational effort and a considerable degree of accuracy. Nevertheless, when a higher degree of accuracy is required, we find the lack of a systematic and simple way for increasing the number of atomic orbitals and for improving the shape of atomic orbitals at a satisfactory level within our knowledge of the LCAO method. Therefore, to overcome the drawback and to benefit the advantages, desirable atomic orbitals as a basis set must be developed with the following features: (i) the computational accuracy and efficiency can be easily controlled by simple parameters as few as possible, (ii) once the number of basis orbitals are fixed, which means that an upper limit is imposed on the computational efforts, the accuracy can be maximized by optimizing the shape of the atomic orbitals. Along this line, recently, accurate basis sets have been constructed in several ways.[12,13,14,15] Kenny et al. constructed a basis set so that atomic orbitals can span the subspace defined by selected and occupied states of reference systems as much as possible.[12] Junquera et al. optimized the shape and cutoff radii of atomic orbitals for reference systems by using the downhill simplex method.[13] However, in these approaches, the serious problem in LCAO has not been solved at a satisfactory level, since the transferability of these optimized orbitals may be restricted to systems similar to the reference systems used for the optimization in terms of atomic environments and states such as the coordination number and the charge state. A more complete treatment for the optimization scheme is to variationally optimize atomic orbitals of each atom located on different environments in a given system.[15,16] Furthermore, it is difficult to control the computational accuracy by simple parameters since the procedure for generating more multiple orbitals than a minimal basis set is not unique in these approaches proposed previously.[12,13]
In this paper, we present the first systematic study of convergence properties for numerical atomic orbitals ranging from H to Kr, in which it is shown that the computational accuracy can be controlled by two simple parameters: a cutoff radius and the number of atomic orbitals. Our comprehensive study not only provides a solution for the above first criterion, but also reveals the limitation in the applicability of LCAO to metallic systems, especially, alkaline and alkaline earth metals. Moreover, starting from our primitive orbitals, a simple and practical method is presented for variationally optimizing numerical atomic orbitals of each atom in a given system based on the force theorem. The orbital optimization scheme enables us to maximize the computational accuracy within a given number of basis orbitals, which fulfills the above second criterion. This paper is organized as follows. In Secs. II and III, we present a method for generating numerical atomic orbitals, and show the convergence properties within DFT in dimers ranging from H to Kr and selected bulks with respect to a cutoff radius and the number of orbitals. In Sec. IV, a simple method is presented for variationally optimizing numerical atomic orbitals of each atom in a given system, and the convergence properties of the optimized orbitals is discussed. In Sec. V, we conclude together with discussing applicability and limitation of the LCAO method.
Let us expand a KohnSham (KS) orbital of a given system
using numerical atomic orbitals
in a form of LCAO:
(1) 
We generate the numerical primitive orbitals based on the following conditions: (i) the atomic orbitals must completely vanish within the computational precision beyond a cutoff radius, and must be continuous up to the third derivatives around the cutoff radius so that matrix elements for the kinetic operator are continuous up to the first derivatives. (ii) a set of atomic orbitals are generated by simple parameters as few as possible, which means that orbitals are systematically available as many as we want. The first condition is needed if the atomic orbitals are coupled with O() methods which suppose the locality of basis sets in the real space. Also the number of nonzero elements of Hamiltonian and overlap matrices is exactly proportional to the number of atomic orbitals if the first condition is satisfied. Once the geometrical structure is given, the structure of the sparce matrices can be predictable through a connectivity table which is prepared from the geometrical structure. Thus, both the computational efforts and the size of memory for evaluating and storing the matrix elements scale as O(). Another cutoff scheme which neglects small elements of Hamiltonian and overlap matrices should not be used because it violates the variational principle.[24] The continuity of atomic orbitals assumed in the first condition is necessary to realize a stable geometry optimization and molecular dynamics (MD) simulations. The second condition we assumed is indispensable in order to obtain a systematic convergence with respect to simple parameters as few as possible.
Our primitive orbitals as a basis set[25]
are orbitals of eigenstates, including excited states,
of an atomic KohnSham equation with
a confinement pseudo potential in a semilocal form for each
angular momentum quantum number .[17,13,16]
To vanish the radial wave function of the outside of the
confinement radius , we modify the atomic core potential
in the all electron calculation of an atom and
keep the modified core potential in the generation of pseudo
potential as follows:
(2) 

The continuity of the modified core potential up to the first derivatives provides that of the radial wave functions up to the third derivatives. Thus, the elements of Hamiltonian and overlap matrices are continuous up to the first derivatives, including that for the kinetic operator. The sharpness of rising edge can be easily controlled by tuning of the radius . If the height of the wall is large enough, the tail of radial function vanishes within the double precision beyond even for excited states unbound without the modified potential. In this study, we used 0.2 (a.u.) and 2.0 (Hartree) for and , respectively, for all elements we considered. It should be mentioned that different modified potentials have been also proposed to shorten the tail of radial wave functions.[17,12,13] However, to the best of our knowledge, there is no modified potentials which can vanish the tail of the excited states at a cutoff radius within the double precision, while keeping the continuity of radial functions. On the other hand, our modified potential with a large can generate a set of continuous numerical radial functions with the complete vanishing tail. In the atomic DFT calculations with the modified potential, there are technical details to generate excited states in a numerical stable way. When the radial differential equation is solved from a distance, the starting radius must be marginally larger than . We determined by a relation (a.u.) for all elements we considered. If is considerably larger than , then numerical instabilities appear due to the large . If the differential equation is solved from the origin, we follow the usual prescription in the atomic DFT calculations.[26] The choice of the matching point, at which two wave functions solved from the origin and a distance are merged, is also an important factor to obtain the excited states in a numerical stable way. In the all electron calculation, we adopt a slightly outside of the most outer peak as the matching point in the usual way.[26] On the other hand, we use a fixed matching point in the calculation of wave functions for the excited states of under pseudo potentials with the modified core potential. This is because the most outer peak often arises near the cutoff radius as the number of nodes increases, which causes numerical instabilities in the solution of the differential equation from the origin. In all elements, we used the fixed matching point at which the logarithmic radial mesh is divided in the ratio of 3 to 1, measured from the origin. Using the fixed matching point, a set of radial wave functions as primitive orbitals were generated without numerical instabilities.[25]


In the DFT calculations using the primitive orbitals, the computational accuracy and efficiency can be controlled by two simple parameters: the cutoff radius and the number of orbitals per atom. The systematic control by two parameters for the accuracy and efficiency is similar to that of spherical wave basis sets.[14] However, we guess that a relatively small number of orbitals may be needed to obtain the convergent result compared to the spherical wave basis sets, since the primitive orbitals are prepared for each element unlike the spherical wave basis sets.[14] This is one of reasons why we use the eigenstates of an atomic KohnSham equation with the confinement pseudo potentials as the primitive orbitals. In order to investigate the convergence properties with respect to the cutoff radius and the number of orbitals per atom, we first show the total energy and the equilibrium bond length of dimer molecules ranging from hydrogen to krypton atoms. The calculation of a dimer molecule could be a severe test for the convergence with respect to basis orbitals, since the neighbouring atom is only one for each atom. In this case a sufficient contribution from orbitals belonging to the other atoms is not anticipated to well express the peripheral region of each atom in KS orbitals, which means that the convergence rate of dimer molecules could be the slowest one. Therefore, the convergence properties of dimer molecule provide a practical guideline in an optimum choice of basis sets for each element.


In Figs. 27 the total energies and the equilibrium bond lengths for dimer molecules from H to Kr are shown as a function of the number of primitive orbitals for various cutoff radii , where the total energies were calculated at the experimental bond length and the equilibrium bond lengths were computed by a cubic spline interpolation for the energy curve as a function of bond length. In most cases, a homonuclear diatomic molecule for each element was investigated. However, if the homonuclear diatomic molecule is not well resolved experimentally, or is highly weak binding, the monoxide or the monohydride was calculated instead of the homonuclear diatomic molecule. In our all DFT calculations, factorized norm conserving pseudo potentials[27,28] were used with multiple projectors proposed by Blochl.[29] In addition to valence electrons, semicore electrons were also considered in making of pseudo potentials for several elements such as alkaline, alkaline earth, and transition elements. Moreover, the nonlinear partial core correction (NLPCC)[30] was considered in the evaluation of the exchangecorrelation terms except for a hydrogen and a helium atom. A relativistic correction was not included in the generation of pseudo potentials. The basis set superposition error (BSSE)[31] was not corrected, since the dissociation energy was beyond the scope of this paper. The electronic states and the cutoff radii used in the generation of these pseudo potentials[25] are shown in Table I, where the cutoff radii are given in parentheses. We limited the study within the local spin density approximation (LSDA)[32] to the exchangecorrelation interactions, since we would like to focus our attention on the convergence properties with respect to the basis orbitals. Also the real space grid techniques were used with the energy cutoff given in the caption of figures in numerical integration[13] and the solution of Poisson's equation using the fast Fourier transformation (FFT). All calculations in this study were performed using our DFT code, OpenMX[25], which is designed for the realization of large scale calculations. We discuss the convergence properties with respect to the basis orbitals in the following categorized elements in Figs. 27.


Representative elements
In representative elements such as H, B, C, N, O, and F, systematic convergence properties are observed as expected. As the cutoff radius and the number of primitive orbitals increase, the total energy converges systematically. Along with the energy convergence, the calculated equilibrium bond length converges at the experimental value within an error of a few percentages. When the energy convergence is carefully observed from the left to the right elements in the first row, we find a trend that primitive orbitals with a higher angular momentum are required in order to achieve an enough convergence. Note that the primitive orbitals with a higher angular momentum beyond valence orbitals allocated to valence electrons are referred to as polarization orbitals according to quantum chemistry in this paper. In H and C, the valence orbitals are almost enough to accomplish the energy convergence. The inclusion of polarization orbitals, and orbitals for hydrogen and carbon atoms, are not effective for the energy convergence. In N and F, the total energies are significantly reduced by the inclusion of orbitals. The overestimated equilibrium bond lengths in the calculations by the and valence orbitals shorten in accordance with the inclusion of orbitals. The same trend as that observed in the first row is found in the convergence properties of the second row elements with respect to angular momentum of primitive orbitals. In P, S, and Cl the polarization orbitals are relevant to obtain the convergent results, while the inclusion of the polarization orbitals insensibly decreases the total energy in Al and Si. In the third row, the polarization orbitals are required in all representative elements from Ga to Br for the energy and geometrical convergences. If the valence and orbitals are only considered, the equilibrium bond lengths tend to be overestimated. Based on the calculations, rough estimations of appropriate cutoff radii of basis orbitals might be given as 5.0, 6.5, and 6.5 a.u. for representative elements of the first, the second, and the third rows, respectively. Although these rough estimations provide a practical guideline in an optimum choice of cutoff radius for each element, it should be mentioned that there is often an exceptional case in which a larger cutoff radius is required for the accurate description. As such an exceptional case, we can point out that a relatively larger cutoff radius must be used when an atom is negatively charged up, since the electrons tend to be far from the atom due to the repulsive interaction between electrons. In contrast, a smaller cutoff radius could be enough to achieve a sufficient convergence when an atom is positively charged up, and when an atom has a high coordination number which means that the atom is surrounded by the other many atoms. We will again discuss the cutoff radius in the later discussion based on numerical results.

Moreover, the convergence with respect to the basis orbitals could be confirmed when the electronic state is carefully examined. In Table II, we find that the electronic configurations of the experimental ground state are well reproduced in the first row representative elements by only the inclusion of valence orbitals. However, the polarization orbitals are more important for the second and third row elements to obtain a convergent result in the electronic configuration. Actually, the use of only the valence orbitals fails to predict the ground state of Si. The basis orbitals, Si6.522, gives to be the ground state of Si, while the ground state is correctly predicted as by the inclusion of orbitals. Although the ground state of Al is determined as , even if the orbitals are included, experimental investigations reveal that the ground state is .[40] The discrepancy between our theoretical prediction and the experiments is attributed to LSDA to the exchangecorrelation interaction as reported by Martinez et al.[66] They also obtained as the ground state of Al within LSDA, while their GGA calculation correctly predicts that the ground state of Al is and the state is the lowest excited state.[66] For the third row representative elements, we did not observe the same kind of discrepancy as that of Si in prediction of the ground state. However, the inclusion of orbitals is recommended for the geometrical convergence as mentioned before. Table III enumerated the equilibrium bond lengths which are calculated using the orbitals with the largest cutoff radius and the greatest number of orbitals for each dimer in Figs. 27. The calculated bond lengths of dimer molecules by representative elements are consistent with both the experimental and the other theoretical values, which supports that our primitive basis orbitals provide the convergent results comparable to the other DFT calculations.
Transition elements
As well as the representative elements, we find systematic convergence properties with respect to the cutoff radius and the number of primitive orbitals for the first row transition elements. Both the total energy and equilibrium bond length converge with increasing of the cutoff radius of orbitals and the number of orbitals. Interestingly, the polarization orbitals do not play an important role in the energy convergence. Within the valence orbitals the energy convergence is almost achieved in all the transition elements that we considered. However, there would be a possibility that the polarization orbitals become more effective for transition elements in a higher row which lies downward in the periodic table we have not studied, like the representative elements. In Table II, we find that the predicted electronic configurations of the ground states are almost consistent with the experimental results within the use of only valence orbitals, except for V and Ni. In V the electronic configuration, , observed experimentally[44] is correctly reproduced by the inclusion of the polarization orbitals, while only the use of the valence orbitals predicts as the ground state of V. The nonmonotonic convergence of the equilibrium bond length of V found in Fig. 5 can be attributed to the transition of singlet to triplet in the calculated ground state. For Ni, recent experiments[48] by resonant twophoton ionization spectroscopy using argon carrier gas suggest that the ground state is or , which would be a mixture of and or and . In contrary, we obtained , which is a hole state, as the ground state within LSDA. Although the calculated ground state is not consistent with the experiments[48], but is equivalent to the other theoretical prediction[70] within generalized gradient approximations (GGA) to the exchangecorrelation interaction. Thus, the discrepancy between the theoretical prediction and the experiments must be attributed to the poor description to exchangecorrelation interaction or the limitation of single configuration method such as DFT rather than the quality of basis orbitals. Also, in Table III we see that the calculated bond lengths are comparable to both the experimental and the other theoretical values. For the first row transition elements, a rough estimation of appropriate cutoff radius of basis orbitals might be given as 7.0 a.u. based on the calculations of dimer molecules, which provides a tradeoff between the computational accuracy and efficiency. However, as discussed in the representative elements, it should be noted that the appropriate choice of cutoff radius depends on atomic environments and states such as the coordination number and the charge state.


Alkaline and alkaline earth elements
Intrinsically, it would be hard to describe KS orbitals of systems consisting of alkaline and alkaline earth elements by using LCAO with short range atomic orbitals. In Fig. 8 shows nodeless pseudo radial wave functions of state in K, Fe, and Br atoms, which are generated by the pseudo potential calculations. The comparison definitely reveals that the radial wave function of an alkaline element has a longer tail compared to those of a transition element and a representative element, which might make the use of basis orbitals with a short cutoff radius difficult. Actually, we find in Figs. 27 that a large cutoff radius is required for dimer molecules of alkaline elements to achieve a sufficient convergence in the total energy and the equilibrium bond length. At least 9.0, 9.0, and 10.0 a.u. of the cutoff radii are needed for Li, Na, and K atoms, respectively. It should be noted that the requirement of a larger cutoff radius for alkaline elements causes great computational demands in the evaluation of the Hamiltonian matrix elements, since the number of grids in the sphere defined with a cutoff radius scales as O(). In contrast, the energy dependency on the number of orbitals is not so large, which represents that the convergence is almost accomplished even if a small number of orbitals are employed. This energy independency with respect to number of orbitals helps us to reduce the computational costs in the application of LCAO to alkaline metal systems in spite of the requirement of a large cutoff radius. For alkaline earth elements, the monoxide molecules were investigated for ease of calculations, since accurate calculations require a considerably larger cutoff energy for the numerical integration and the solution of Poisson's equation due to highly weak binding of the homonuclear diatomic molecules of alkaline earth elements. From Figs. 27 we find that the total energies and the equilibrium bond length systematically converge as the cutoff radius and the number of primitive orbitals increase like representative elements. In these monoxide molecules, 7.0 a.u. of the cutoff radius gives a tradeoff between the computational accuracy and efficiency for all the alkaline earth elements. The inclusion of polarization oribtals is essential for the convergence in CaO, while the orbitals do not play an important role in the energy and geometrical convergences of BeO and MgO. This is because the state of a calcium atom exists near the valence state. Therefore, the orbitals of CaO considerably contribute in and KS orbitals. In fact, the LCAO coefficient of the orbital with a nodeless radial function of the calcium atom is about 0.3 in the occupied orbital of the monoxide molecule along axis, in which Ca7.0222 and O5.022 were used as the basis orbitals. In Table II gives the electronic configurations for the ground state of the alkaline element dimers and the monoxide molecules of alkaline earth elements. The calculated electronic configurations for the ground state are wholly consistent with those determined experimentally. Also, we find in Table III that the calculated equilibrium bond lengths are comparable to both the experimental and the other theoretical values for the alkaline element dimers and the monoxide molecules of alkaline earth elements. These results support that our primitive basis orbitals give a complete basis set even for the alkaline and alkaline earth elements, while the cutoff radius required for the convergence is larger than those of representative and transition elements.
Rare gas elements
It has been reported that LDA and GGA fail to predict the equilibrium bond lengths and the dissociation energies of dimer molecules consisting of rare gas elements weakly binding by Van der Waals interactions.[62] However, our attention in this study is to know the convergence properties with respect to basis orbitals. Therefore, we investigated homonuclear diatomic molecules of rare gas elements, He, Ne, Ar, and Kr within LDA, and compared the calculated equilibrium bond lengths with the other theoretical values calculated by LDA. The dissociation energies of the rare gas dimers are significantly small unlike dimers of representative and transition elements. Therefore, we had to use higher cutoff energies, which are 262, 343, 290, and 290 (Ryd) for He, Ne, Ar, and Kr, for the numerical integration and the Poisson's equation so that the computations are not buried in numerical errors. It seems that the convergence in the total energy of rare gas dimers depends on only the cutoff radius of primitive orbitals. For all the rare gas dimers 7.0 a.u. of the cutoff radius are enough to accomplish a sufficient convergence. On the other hand, the calculated equilibrium bond lengths have a dependency on the number of orbitals, especially in Ar and Kr, however, the double valence with single polarization orbitals provide almost convergence results even for Ar and Kr. In Table III we see that the calculated equilibrium bond lengths are comparable to the other theoretical values calculated by all electron calculations within LDA. In addition, the calculated electronic configurations for the ground states are consistent with those reported experimentally. Thus, our primitive orbitals could be a systematic basis set for rare gas dimers as well as representative, transition, and alkaline and alkaline earth elements, while the Van der Waals interactions in the rare gas dimers is not correctly taken into account in LDA calculations.
Bulks
For selected bulk systems, we investigated the convergence properties of total energy, lattice constant, bulk modulus, and magnetic moment with respect to basis orbitals as shown in Fig. 9. The bulk modulus was calculated by a least square fitting of the total energy curve to Murnaghan's equation of state.[75] The magnetic moment in the bcc iron are evaluated from the difference in Mulliken charges of up and down spins. In the graphite carbon, the lattice constant of only the plane was varied with a lattice constant of the axis fixed at the experimental value. For all bulk systems the same systematic convergence as that for dimer molecules is achieved as the cutoff radius and the number of primitive orbitals increase. When the shorter cutoff radius of primitive orbitals was used, the calculated lattice constant (bulk modulus) tend to be shorter (greater) than the experimental values. As the cutoff radius increases, these calculated values converge at experimental values within an error of a few percentages. For carbon in the diamond and graphite, 4.0 and 4.5 a.u. of the cutoff radius might be regarded as a tradeoff between the computational accuracy and efficiency, respectively. The comparison in an appropriate cutoff radius for various systems suggests that the coordination number of each atom is an important factor which determines an adequate values of the cutoff radius. Our rough estimations of an appropriate cutoff radius are 5.0, 4.5, and 4.0 a.u. for the carbon dimer C, the diamond carbon, and the graphite carbon of which the coordination numbers are one, three, and four. Thus, an appropriate cutoff radius could be in inverse proportion to the coordination numbers. This observation is also confirmed in a comparison of the convergence with respect to a cutoff radius for iron, an iron dimer Fe and the bcc iron. In the bcc iron, we can obtain almost convergent results using 4.5 a.u. of the cutoff radius, while 7.0 a.u. are needed to obtain the convergence in Fe.
Here, it must be stressed that the calculation of the bcc iron illustrates a limitation in applicability of LCAO method. We found that it was difficult to perform reliable calculations of the bcc iron using Fe4.5333 or more. In particular, it was unable to perform meaningful calculations due to numerical errors for the bcc iron with a shorter lattice parameter. The problem comes from the overcompleteness of basis orbitals. In the bcc iron with a shorter lattice parameter and the use of a longer cutoff radius of basis orbitals, eigenvalues of the overlap matrix can be negative, which means that the basis orbitals is not linearly independent. A remedy in this case is to reduce matrices by removing eigenvectors corresponding to the negative eigenvalues. However, the eigenvalues of the overlap matrix in bulk systems distribute continuously as a function of energy, thus, many illconditioned eigenvalues, which are positive but almost zero, appear in the eigenvalue spectrum. The division by the illconditioned eigenvalues brings about the numerical instability we met. If basis orbitals with a larger cutoff radius is used for a dense system with atoms of large coordination numbers, a small number of optimized orbitals should be used to avoid the numerical instabilities.
In this Section, we present a simple and practical method for
variationally optimizing numerical basis orbitals of each atom located
on different environments in a given system.[16]
Starting from the primitive basis orbitals discussed in the Sec. II,
the shape of radial wave functions of each atom are variationally
optimized within a given cutoff radius so that the total energy
is minimized based on the force theorem. The orbital optimization scheme
promises us to reduce the computational cost, while keeping
a high degree of accuracy.
Although the primitive radial wave function ,
the eigenstate of atomic KS equation with a confinement potential,
was used as radial basis orbital of the KS orbital in
a form of LCAO in Sec. II, here, we reconsider a different expression
for and thus
.
To give a variational degree of freedom of
,
we furthermore expand
using primitive orbitals
as follows:
(3) 

The contraction coefficients can be easily optimited by the two step optimization scheme. The details of the two step optimization scheme has already been described in elsewhere [16]. In this two step optimization scheme, the atomic orbitals are optimized variationally in the same two step procedure as that of the geometry optimization in terms of instead of atomic positions. The radial parts of basis orbitals in each atom located on different environment are automatically varied so that the total energy is minimized, which is a quite important advantage of our scheme compared to the other optimization method [12,13]. In the later part of this section, we demonstrate capability of our method based on numerical results. Figure 10 shows the convergence properties of total energies for a carbon dimer C, a methane molecule CH, carbon and silicon in the diamond structure, an ethane molecule, and a hexafluoro ethane molecule as a function of the number of primitive and optimized orbitals. The orbital optimization was done by ten iterative steps according to Eq. (7), in which each step includes ten SCF loops on the average. We see that the unoptimized orbitals provide systematic convergent results for not only molecules, but also bulk systems as the number of orbitals increase as discussed in Sec. II. Moreover, remarkable convergent results are obtained using the optimized orbitals for all systems. The small set of optimized orbitals rapidly converge to the total energies calculated by a larger number of the primitive orbitals, which implies that the computational effort can be reduced significantly with a high degree of accuracy. Only the restricted contraction was investigated in this study, since we found in the previous study that the unrestricted optimization is not effective to minimize the total energy.[16] Also, the restricted optimization guarantees the rotational invariance of the total energy. Therefore, our study was limited within the restricted optimization.
In Fig. 4 the radial parts of the minimal orbitals obtained by the restricted optimization for CH and CF are shown with those of the lowest primitive orbitals of an carbon atom for comparison. It is observed that the tails of both the optimized and orbitals shrink compared to the primitive orbitals in CH and CF. In addition, we find that the orbitals of the carbon atom in CF more shorten than that of CH. The considerable shortening tail of the orbital is related to change in the effective charge of carbon atom. Decomposed Mulliken populations of the carbon atom are 1.05 and 2.67 in CH, and 0.86 and 2.00 in CF for and orbitals, respectively. So, we see that the deviation for the orbital is larger than that for the orbital in a comparison of the decomposed Mulliken population. Therefore, the large shortening tail of orbital in CF is due to increase of effective core potential for electrons. The comparison between CH and CF clearly reveals that the shape of the basis orbital can automatically vary within the cutoff radius in order to respond to the change of the environment of each atom, while minimizing the total energy.
First, we performed the geometry optimization with the orbital
optimization as a preconditioning for a C molecule.
Before doing the geometry optimization, the orbital optimization was
performed by ten iterative steps, which includes ten SCF loops per
step on the average.

Second, we show that it is significantly effective for the realization of a high degree of accuracy and efficiency to construct a database of optimized orbitals for a specific group such as proteins. Proteins are constructed from twenty kinds of amino acid residues. So, we categorized atoms in the residues as eleven, eighteen, four, three, and two kinds of hydrogen, carbon, nitrogen, oxygen, and sulfur atoms from a chemical point of view in consideration of chemical environment and connectivity. To construct a database of optimized orbitals for the categorized atoms, the structures of tripeptides, GlyXGly, are considered for the orbital optimization, where X could be one of twenty kinds of amino acid residues. The structure of a GlyXGly was optimized by a molecular mechanics (MM) using a software TINKER[85] with a AMBER98 force field[86] before the orbital optimization. Then, for the optimized GlyXGly the restricted orbital optimization was performed by ten iterative steps with nine SCF loops per step on the average, in which LDA was employed to exchangecorrelation interaction and the cutoff energy of 130 Ryd was used for numerical integration and the solution of Poisson's equation. In a series of optimizations, the basis specifications were given as H4.55251, C5.0525251, N4.5525251, O4.5525251, and S6.0525251. Because of the basis specifications, the orbitals stored in the database are well optimized for the use of double valence plus polarization orbitals. However, the basis sets maybe provide a better performance even for the other specifications of basis sets compared to the original primitive basis sets, when the basis sets are used for calculations of proteins. Following the construction of the database, we performed SCF calculations of small polypeptides to investigate the transferability of the optimized orbitals for proteins. In Table V shows total energies of small polypeptides, Metenkephalin[82], Valorphin[83], and dynorphin A[84], calculated using unoptimized primitive orbitals and the optimized orbitals stored in the database, where the geometrical structures of the small polypeptides are generated by a simulated annealing method using a software TINKER[85] with a AMBER98 force field[86] to apply the optimized orbitals to an arbitrary structural conformation. We see that a set of the basis orbitals optimized for proteins give a lower energy than the primitive orbitals in all the polypeptides, which shows that the optimized orbitals well span the occupied states of proteins beyond the primitive orbitals. This illustrates that the database construction for a specific system promises us a substantial improvement in the basis convergence, while keeping the same computational efforts as that of the primitive orbitals. The details of the database construction for proteins will be presented elsewhere.

To conclude, we have presented the first systematic study for numerical atomic basis orbitals ranging from H to Kr. Our primitive orbitals as a basis set are eigenstates, including excited states, of an atomic KohnSham equation with a confinement pseudo potential in a semilocal form for each angular momentum quantum number . The scheme has been discussed for generating the systematic basis orbitals in a numerical stable way. The comprehensive investigation of convergence properties shows that our primitive basis orbitals could be one of practical solutions as a systematic basis set in large scale DFT calculations for a wide variety of systems. Using the primitive orbitals, the computational accuracy and efficiency are systematically controlled by two simple parameters: the cutoff radius and the number of orbitals per atom. As the cutoff radius and the number of primitive orbitals increase, the total energy and the physical quantities converge systematically. The convergence properties of total energy and equilibrium bond length for dimer molecules with respect to basis orbitals provide a practical guideline in an optimum choice of a cutoff radius, the number and the maximum angular momentum of basis orbitals for each element. In addition, our widespread study shows limitations of the LCAO method to metallic systems and dense structures with a large coordination number. In alkaline and alkaline earth elements, valence orbitals tend to have a longer tail, which makes applications of the LCAO method to the systems difficult due to increase of computational costs. In dense structures such as bcc, fcc, and hcp, the primitive basis orbitals often become overcomplete. Owing to the overcompleteness, we have difficulty in the systematic improvement of basis orbitals for systems with a dense structure. Therefore, careful treatments are required in the applications of LCAO method to a such kind of systems. In spite of the difficulty, we believe that the primitive orbitals can be a systematic basis set in a wide variety of materials, especially for highly covalent systems such as organic molecules.
Furthermore, we have developed a simple and practical method, based on the force theorem, for variationally optimizing the radial shape of numerical atomic orbitals. The optimization algorithm similar to the geometry optimization allows us to fully optimize atomic orbitals within a cutoff radius for each atom in a given system. The optimized orbitals well reproduce convergent results calculated by a larger number of primitive orbitals. The comparison between CH and CF demonstrates that the basis orbital can automatically vary within the cutoff radius in order to respond to the change of the environment of each atom, while minimizing the total energy. As practical applications of the orbital optimization, we have demonstrated two examples: the geometry optimization coupled with the orbital optimization of a C molecule and the preorbital optimization for a set of amino acid residues. The former shows that the small set of optimized orbitals promises to greatly reduce the computational effort with a high degree of accuracy. The latter demonstrates that it is significantly effective for the realization of a high degree of accuracy and efficiency to construct a database of optimized orbitals for a specific group such as proteins. The scheme also could be a remedy for the problem of the overcompleteness. Thus, we conclude that the optimized orbitals are well suited for large scale DFT calculations.
We would like to thank K. Kitaura and T. Miyazaki for helpful suggestions on the orbital optimization. This work was partly supported by NEDO under the Nanotechnology Materials Program, and Research and Development for Applying Advanced Computational Science and Technology of Japan Science and Technology Corporation (ACTJST).