Orbital Optimization Method: Ver. 1.0
The original has been published in Phys. Rev. B 67, 155108 (2003).

Abstract

A simple and practical method for variationally optimizing numerical atomic orbitals used in density functional calculations is presented based on the force theorem. The derived equation provides the same procedure for the optimization of atomic orbitals as that for the geometry optimization. The optimized orbitals well reproduce convergent results calculated by a larger number of unoptimized orbitals. In addition, we demonstrate that the optimized orbitals significantly reduce the computational effort in the geometry optimization, while keeping a high degree of accuracy.

1 INTRODUCTION

During the last decade, to extend the applicability of density functional theories (DFT) to realistic large systems, great efforts have been done for developing O($N$) methods of the eigenvalue problem [1, 2, 3, 4, 5, 6] and for making efficient and accurate localized orbitals [7, 8, 9, 10] as a basis set being suitable for O($N$) methods. Among these studies, one of important and unresolved problems is how atomic orbitals as a basis set are constructed to maximize both the computational efficiency and accuracy. One expects that a basis set such as double valence orbitals with polarized orbitals for valence electrons provides a way for balancing a relatively small computational effort and a considerable degree of accuracy. Along this line, accurate basis sets were constructed in several ways [7, 8, 9, 10]. Kenny et al. constructed a basis set, so that atomic orbitals span the subspace defined by selected and occupied states of reference systems as much as possible [7]. Junquera et al. optimized the shape and cutoff radii of atomic orbitals for reference systems by using the downhill simplex method [8]. However, the transferability of these optimized orbitals might 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 is to optimize atomic orbitals of each atom located on different environments in a given system [10]. In addition, the complete optimization procedure should be simple and efficient practically. In this paper, to overcome the difficulty, we present a simple and practical method, based on the force theorem, for variationally optimizing numerical atomic orbitals of each atom in a given system.

2 VARIATIONAL OPTIMIZATION

Let us expand a Kohn-Sham (KS) orbital ${\psi }_{\mu }$ of a given system using numerical atomic orbitals ${\varphi }_{i\alpha }$ in a form of linear combination of atomic orbitals (LCAO):

${\psi }_{\mu }(𝐫)$ $=$ ${\displaystyle \sum _{i\alpha }}{c}_{\mu ,i\alpha }{\varphi }_{i\alpha }(𝐫-{𝐫}_i),$ (1)

where $i$ is a site index, $\alpha \equiv (plm)$ an organized orbital index, and ${\varphi }_{i\alpha }\equiv Y_{lm}{R}_{iplm}$. For simplicity we consider only nonspin-polarized systems and an non-Bloch expression of the one-particle wave functions, but the extentions of the below description to spin-polarized systems and Bloch wave functions are straightforward. Note that a radial wave function $R_{iplm}$ depends on not only an angular momentum quantum number $l$, but also a site index $i$, a multiplicity index $p$, and a magnetic quantum number $m$. To give a variational degree of freedom of ${\varphi }_{i\alpha }$, we furthermore expand ${\varphi }_{i\alpha }$ using primitive orbitals ${\chi }_{i\eta }$ as follows:

${\varphi }_{i\alpha }(𝐫)$ $=$ ${\displaystyle \sum _q}{a}_{i\alpha q}{\chi }_{i\eta }(𝐫),$ (2)

where $\eta \equiv (qlm)$, in which the indices $l$ and $m$ denote the same as those of the index $\alpha$, and ${\chi }_{i\eta }\equiv Y_{lm}{R}_{iql}^{\prime }$. Note that a primitive radial wave function $R_{iql}^{\prime }$, which is discussed later on, is independent on $m$, and that the coefficients $a_{i\alpha q}$ are independent variables on the eigenstate $\mu$. Substituting Eq. (2) into Eq. (1), we have

${\psi }_{\mu }(𝐫)$ $=$ ${\displaystyle \sum _{i\alpha }}{\displaystyle \sum _q}{c}_{\mu ,i\alpha }{a}_{i\alpha q}{\chi }_{i\eta }(𝐫-{𝐫}_i).$ (3)

Although the expansion of a KS orbital by Eq. (3) is linearized for each variable $c_{\mu ,i\alpha }$ or $a_{i\alpha q}$, however, the primitive orbitals ${\chi }_{i\eta }$ are expanded by the product of two variables $c_{\mu ,i\alpha }$ and $a_{i\alpha q}$ in a nonlinearized form. Therefore, it is difficult to directly find the minimum of the KS total energy $E_{\mathrm{tot}}$ for the ground state with respect to $c_{\mu ,i\alpha }$ and $a_{i\alpha q}$ as a generalized eigenvalue problem which can be derived in the usual LCAO. To avoid the difficulty, here, we propose a two step optimization scheme, in which the coefficients $a_{i\alpha q}$ are optimized after $c_{\mu ,i\alpha }$ are determined with a set of fixed $a_{i\alpha q}$. Considering $\partial E_{\mathrm{tot}}/\partial c_{\mu ,i\alpha }=0$ for the KS total energy $E_{\mathrm{tot}}$ with the orthonormalization relation $⟨{\psi }_{\mu }|{\psi }_{\nu }⟩={\delta }_{\mu \nu }$ among one-particle wave functions ${\psi }_{\mu }$ and fixed contraction coefficients $a_{i\alpha q}$, we have a well-known KS matrix equation with respect to the coefficient $c_{\mu ,i\alpha }$ as follows:

${\displaystyle \sum _{j\beta }}⟨{\varphi }_{i\alpha }|\widehat{H}|{\varphi }_{j\beta }⟩c_{\mu ,j\beta }$ $=$ ${\epsilon }_{\mu }{\displaystyle \sum _{j\beta }}⟨{\varphi }_{i\alpha }|{\varphi }_{j\beta }⟩c_{\mu ,j\beta },$ (4)

where $\widehat{H}$ is a KS Hamiltonian, and ${\epsilon }_{\mu }$ is a KS eigenvalue of the system. There is no restriction to solve Eq. (4) or to find $\partial E_{\mathrm{tot}}/\partial c_{\mu ,i\alpha }=0$. So, we can use any solution method which could be an exact diagonalization method, iterative methods such as Car-Parinello (CP) method [11] and conjugate gradient (CG) method [12], and O($N$) methods. On the other hand, regarding $c_{\mu ,i\alpha }$ as dependent variables on $a_{i\alpha q}$ and assuming that the Kohn-Sham equation is solved self-consistently with respect to $c_{\mu ,i\alpha }$, we can derive the following equation based on the force derivation of non-orthogonal orbitals as follows:

${\displaystyle \frac{\partial E_{\mathrm{tot}}}{\partial a_{i\alpha q}}}$ $=$ ${\displaystyle \frac{\delta E_{\mathrm{tot}}}{\delta \rho (𝐫)}}{\displaystyle \frac{\delta \rho (𝐫)}{\delta a_{i\alpha q}}}$ (5)

$=$ $4{\displaystyle \sum _{\mu }}{n}_{\mu }{\displaystyle \sum _{i\alpha ,j\beta }}{c}_{\mu ,i\alpha }{c}_{\mu ,j\beta }⟨{\displaystyle \frac{\partial {\varphi }_{i\alpha }}{\partial a_{i\alpha q}}}|\widehat{H}|{\varphi }_{j\beta }⟩+4{\displaystyle \sum _{\mu }}{n}_{\mu }{\displaystyle \sum _{i\alpha ,j\beta }}{\displaystyle \frac{\partial c_{\mu ,i\alpha }}{\partial a_{i\alpha q}}}{c}_{\mu ,j\beta }⟨{\varphi }_{i\alpha }|\widehat{H}|{\varphi }_{j\beta }⟩$

$=$ $2{\displaystyle \sum _{j\beta }}\left({\mathrm{\Theta }}_{i\alpha ,j\beta }⟨{\chi }_{i\eta }|\widehat{H}|{\varphi }_{j\beta }⟩-E_{i\alpha ,j\beta }⟨{\chi }_{i\eta }|{\varphi }_{j\beta }⟩\right),$

where $n_{\mu }$ is an occupancy number for the eigenstate $\mu$, ${\mathrm{\Theta }}_{i\alpha ,j\beta }$ a bond-order, and $E_{i\alpha ,j\beta }$ an energy bond-order. The final equation in Eq. (5) is derived by taking into account Eq. (4) and the orthonormalization relation $⟨{\psi }_{\mu }|{\psi }_{\nu }⟩={\delta }_{\mu \nu }$. It should be noted that Eq. (5) excludes any derivative, and that ${\mathrm{\Theta }}_{i\alpha ,j\beta }$ and $E_{i\alpha ,j\beta }$ only have to be evaluated for the contracted atomic orbital ${\varphi }_{i\alpha }$ in Eq. (5), which implies that additional computational costs are not required to evaluate Eq. (5). Once we obtain a self-consistent solution of the Kohn-Sham equation, Eq. (4), with a set of given coefficients $a_{i\alpha q}$, then Eq. (5) gives the gradient of $E_{\mathrm{tot}}$ with respect to $a_{i\alpha q}$ within small computational costs. This fact shows apparently that the atomic orbitals can be optimized variationally in the same two step procedure as that of the geometry optimization in terms of $a_{i\alpha q}$ instead of atomic positions. Therefore, the contraction coefficients $a_{i\alpha q}$ are optimized iteratively, coupled with the self-consistent solution of Eq. (4), as follows:

Step 1 $\text{Self-consistently solving of Eq. (}\text{<a href="\#S2.E4" title="(4) ‣ 2 VARIATIONAL OPTIMIZATION ‣ Orbital Optimization Method: Ver. 1.0 The original has been published in Phys. Rev. B 67, 155108 (2003)." class="ltx\_ref"><span class="ltx\_text ltx\_ref\_tag">4</span></a>}\text{)},$

Step 2 $a_{i\alpha q}^{(n+1)}=a_{i\alpha q}^{(n)}-\lambda {\left({\displaystyle \frac{\partial E_{\mathrm{tot}}}{\partial a_{i\alpha q}}}\right)}_{a^{(n)}},$

$n:=n+1,$

where an optimum $\lambda$ is determined, so that the norm of the gradients at $a=a^{(n+1)}$,

$N_{\mathrm{G}}^{(n+1)}$ $=$ ${\displaystyle \sum _{i\alpha ,q}}{\left({\displaystyle \frac{\partial E_{\mathrm{tot}}}{\partial a_{i\alpha q}}}\right)}_{a^{(n+1)}}^2,$ (7)

becomes a minimum with respect to $\lambda$ under the fixed ${\mathrm{\Theta }}_{i\alpha ,j\beta }$ and $E_{i\alpha ,j\beta }$. Substituting Eqs. (5) and (2) into Eq. (7), and considering $\partial N_{\mathrm{G}}^{(n+1)}/\partial \lambda =0$ with the fixed ${\mathrm{\Theta }}_{i\alpha ,j\beta }$ and $E_{i\alpha ,j\beta }$, we have $\lambda =B/A$ for the minimum of $N_{\mathrm{G}}^{(n+1)}$ with

$A$ $=$ ${\displaystyle \sum _{i\alpha q}}{\left({\displaystyle \sum _{j\beta q^{\prime }}}{D}_{i\alpha q,j\beta q^{\prime }}{\left({\displaystyle \frac{\partial E_{\mathrm{tot}}}{\partial a_{j\beta q^{\prime }}}}\right)}_{a^{(n)}}\right)}^2,$ (8)

$B$ $=$ ${\displaystyle \sum _{i\alpha q,j\beta q^{\prime }}}{D}_{i\alpha q,j\beta q^{\prime }}{\left({\displaystyle \frac{\partial E_{\mathrm{tot}}}{\partial a_{i\alpha q}}}\right)}_{a^{(n)}}{\left({\displaystyle \frac{\partial E_{\mathrm{tot}}}{\partial a_{j\beta q^{\prime }}}}\right)}_{a^{(n)}},$ (9)

where

$D_{i\alpha q,j\beta q^{\prime }}$ $=$ ${\mathrm{\Theta }}_{i\alpha ,j\beta }{H}_{i\eta ,j{\eta }^{\prime }}-E_{i\alpha ,j\beta }{S}_{i\eta ,j{\eta }^{\prime }}$ (10)

with $H_{i\eta ,j{\eta }^{\prime }}\equiv ⟨{\chi }_{i\eta }|\widehat{H}|{\chi }_{j{\eta }^{\prime }}⟩$ and $S_{i\eta ,j{\eta }^{\prime }}\equiv ⟨{\chi }_{i\eta }|{\chi }_{j{\eta }^{\prime }}⟩$. Taking into account the sparseness of both the Hamitonian and overlap matrices in the real space, we find that the computational efforts to evaluate Eqs. (8) and (9) scale linearly. Therefore, we can easily evaluate $\lambda$, leading no intensive computational demands. After achieving the self-consistent field (SCF) for Eq. (4) with respect to the coefficients $c_{\mu ,i\alpha }$ by an usual SCF procedure, the contraction coefficients $a_{i\alpha q}$ are updated by Eq. (2) and renormalized so that ${\varphi }_{i\alpha }$ is normalized. Thus, the two step optimization scheme enables us to optimize the contraction coefficients $a_{i\alpha q}$ along the stationary minimum line of the KS total energy functional with respect to $c_{\mu ,i\alpha }$, while keeping $⟨{\psi }_{\mu }|{\psi }_{\nu }⟩={\delta }_{\mu \nu }$. It is found that about five iterative procedures of the two step optimization, which includes the solution of Eq. (4) and the optimization of $a_{i\alpha q}$ by Eq. (2), are enough to accomplish a sufficient convergence of $a_{i\alpha q}$ for our test systems. Again it should be mentioned that the bond-order and the energy bond-order are required for only the contracted atomic orbital ${\varphi }_{i\alpha }$ in Eq. (5). This is a crucial point to make our optimization procedure efficient, since the Kohn-Sham equation based on Eq. (4) can be solved for not large ${\chi }_{i\eta }$, but small ${\varphi }_{i\alpha }$. Once the contraction coefficients $a_{i\alpha q}$ are fixed after the orbital optimization, the Hamiltonian and overlap matrices are directly constructed for the small ${\varphi }_{i\alpha }$ without constructing the elements for the larger primitive orbitals, since we can directly utilize the contracted orbital ${\varphi }_{i\alpha }$ as a numerical table because of the use of numerical orbitals, which is also a reason why the optimization scheme could be totally efficient.

A rather technical but important problem still remains in the application of the two step optimization method. To avoid the orbital optimization to a local minimum, we have developed the following careful procedure to provide a good initial guess for a set of coefficients $a_{i\alpha q}$:

(I) The partition of the system. A cluster is constructed including the nearest neighboring atoms $j$ for each atom $i$.

(II) Solving of the Kohn-Sham equation of each cluster. By non self-consistently solving the Kohn-Sham equation of each cluster which is centered on an atom $i$, we obtain the coefficients $d_{\nu ,iqlm}$ for the atom $i$ of an one-particle wave function

${\phi }_{\nu }^{(i)}$ $=$ ${\displaystyle \sum _{jqlm}}{d}_{\nu ,jqlm}{\chi }_{jqlm}$ (11)

with an eigenvalue ${ϵ}_{\nu }^{(i)}$.

(III) Construction of $a$ from $d$. Then, $a_{i0lmq}$ is given as

$a_{i0lmq}$ $=$ $\mathrm{N}{\displaystyle \sum _{\nu }}\mathrm{sgn}(d_{\nu ,i0lm})d_{\nu ,iqlm}f(({ϵ}_{\nu }^{(i)}-{\mu }_i)/k_{\mathrm{B}}T),$ (12)

where $f$ is the Fermi function, ${\mu }_i$ a local chemical potential for the cluster $i$, and N a normalization factor. For , the coefficients $a_{iplmq}$ are generated by the Gram-Schmidt orthonormalization from $a_{i0lmq}$ and $d_{\nu ,iqlm}$ in order of the magnitude of $f(({ϵ}_{\nu }^{(i)}-{\mu }_i)/k_{\mathrm{B}}T){\sum }_q|d_{\nu ,iqlm}|$.

To find a good initial guess for the coefficients $a_{i\alpha q}$, we tried to estimate the ratio of coefficients $a_{i\alpha q}$ of the eigenstates of the whole system expanded by ${\chi }_{i\eta }$ from the local cluster for each atom $i$ by the above treatment. Then, the eigenstates of the cluster are weighted by the Fermi function, so that the contribution of the lower states is taken into account as much as possible. Also, the contraction coefficients $a_{iplmq}$ for are generated using the Gram-Schmidt method to avoid the overcompleteness of contracted basis orbitals. We found that the procedure provides good initial coefficients $a_{i\alpha q}$ in all of our test systems, and did not observe that the orbital optimization is trapped to any serious local minimum, while the other trial was trapped to a local minimum often. The additional cost for the above procedure (I)-(III) is almost negligible, when the orbital optimization method is applied as a preconditioning of the geometry optimization as discussed later on.

3 PRIMITIVE ORBITALS

The primitive orbitals ${\chi }_{i\eta }$ we used are eigenstates of an atomic Kohn-Sham equation with confinement pseudo potentials [8, 13]. To vanish the radial wave function $R_{iql}^{\prime }$ of the outside of the confinement radius $r_{\mathrm{c}}$, we modify the atomic core potential $V_{\mathrm{core}}(r)$ in the all electron calculation of an atom and the generation of pseudo potential as follows:

$V_{\mathrm{core}}(r)$ $=$ (13)

where $b_0$, $b_1$, $b_2$, and $b_3$ are determined, so that the value and the first derivative are continuous at both $r_1$ and $r_{\mathrm{c}}$. Figure 1 shows radial wave functions for $l=0$ of a carbon atom under the confinement pseudo potential. The eigenstates construct an orthonormal basis set at the same atomic position and vanish beyond $r_{\mathrm{c}}$ within the double precision. Because of the complete vanishing tail of numerical orbitals, we find that non-zero elements of Hamiltonian and overlap matrices can be exactly proportional to the number of atoms.

Figure 1: The radial wave function for $l=0$ of a carbon atom under the confinement pseudo potential obtained from Eq. (13), where 4.5 (a.u.), 4.3 (a.u.), and 3.0$\times$ ${10}^4$ (Hartree) are used for $r_{\mathrm{c}}$, $r_1$ and $h$, respectively.

In Fig. 2 the total energy for a carbon dimer calculated using the eigenstates as a basis set is shown as a function of the number of orbitals for various cutoff radii $r_{\mathrm{c}}$. Factorized norm conserving pseudo potentials [13] and the local density approximation (LDA) [14] to the exchange-correlation interactions were used in our all DFT calculations. Also the real space grid techniques were used with the energy cutoff of 113 (Ryd) for numerical integrations [8]. As the cutoff radius and the number of orbitals increase, the total energy converges systematically. Thus, we see that the primitive orbitals ${\chi }_{i\eta }$ itself are systematic basis sets controlled by two simple parameters, the cutoff radius and the number of orbitals, in the same manner as spherical wave basis sets [9]. In addition, a relatively small number of orbitals may be needed to obtain the convergent result compared to the spherical wave basis sets, since the primitive basis set is prepared for each element, unlike the spherical wave basis sets [9]. These are reasons why we use the eigenstates of an atomic Kohn-Sham equation with the confinement pseudo potentials as the primitive orbitals. A systematic study for convergence properties as a function of the cutoff radius and the number of orbitals will be presented for several elements including first row elements, alkaline metals and transition metals elsewhere.

Figure 2: The total energy for a carbon dimer calculated using the eigenstates as a basis set as a function of the number of basis orbitals per atom for the cutoff radius $r_{\mathrm{c}}$ of 3.5, 4.0, 4.5, 5.0, 5.5 and 6.0 (a.u.).

For the later discussion, here, we introduce an abbreviation of the basis orbital as C4.5-s62${}^{*}$p62, where C indicates the atomic symbol, 4.5 is the cutoff radius $r_{\mathrm{c}}$ (a.u.) used in the generation, s62 means that two optimized orbitals are constructed from six primitive orbitals for the s-orbital, and the symbol ${}^{*}$ signifies the restricted optimization that the radial wave function $R$ is independent on the index $m$. In case of snn such as s66, corresponding to no optimization, snn can be simplified as sn.

4 NUMERICAL RESULTS

Figure 3 shows the convergence properties of total energies for a carbon dimer C${}_2$, a methane molecule CH${}_4$, and the diamond as a function of the number of unoptimized and optimized orbitals. The orbital optimization was done by five iterative steps according to Eq. (2), in which each step includes ten SCF loops. We see that the unoptimized orbitals provide systematic and rapid convergent results for not only molecules C${}_2$ and CH${}_4$, but also a bulk system diamond, as the number of orbitals increase. Moreover, remarkable convergent results are obtained using the optimized orbitals for all systems. The small optimized orbitals rapidly converge to the total energies calculated by a larger number of unoptimized orbitals, which implies that the computational effort can be reduced significantly with a high degree of accuracy. For three systems the effect of the restriction for the orbital optimization is almost negligible, which encourages us to use the restriction, since the restricted optimization guarantees the rotational invariance of the total energy. In Fig. 4 the radial parts of the minimal orbitals obtained by the restricted optimization for the diamond are shown with those of the lowest primitive orbitals of a carbon atom for comparison. It is observed that the tails of both the optimized s- and p-orbitals shrink compared to the primitive orbitals, which clearly reveals that the basis orbital can automatically vary within the cutoff radius to minimize the total energy.

Figure 3: The total energy for a carbon dimer C${}_2$, a methane CH${}_4$, and the diamond as a function of the number of unoptimized (unopt) orbitals and optimized orbitals with (rest) and without (unrest) the restriction. The total energy and the number of orbitals are defined as those per atom for C${}_2$ and the diamond, and as those per molecule for CH${}_4$. The energy cutoff of 113, 113, and 222 (Ryd) were used for the numerical integrations in C${}_2$, CH${}_4$, and the diamond, respectively. The two step convergence of C${}_2$ is due to the inclusion of d-orbitals.

Figure 4: The radial wave function of the minimal orbitals obtained by the restricted optimization for the diamond and the lowest primitive orbitals of a carbon atom. The optimization was done in the same conditions as those in Fig. 3.

Finally, as an illustration of the orbital optimization, we performed the geometry optimization with the orbital optimization as a preconditioning for the most stable conformer of a neutral glycine molecule [15, 16] which is the smallest amino acid. Before doing the geometry optimization, the orbital optimization was performed by five iterative steps, which includes ten SCF loops per step, for an initial structure optimized by a molecular mechanics (MM). Then, the geometry optimization was done using the optimized orbitals by fifty steepest decent (SD) steps with a variable prefactor for accelerating the convergence, which includes twenty SCF loops per step. The optimized geometrical parameters are given in Table 1 together with the total energy and the computational time per MD step. In case of the unoptimized orbitals SN, TN, and TNDP, as the number of orbitals increase, we find the decrease of the total energy and the convergent geometrical parameters comparable to the experimental [16] and the other theoretical values [15]. Although there are some deviations in the optimized parameters calculated using TNDP from the other theoretical values [15], the deviations may be attributed to the pseudo potentials rather than the basis orbitals, since we verified that the optimized parameters of the glycine depend on the cutoff radii in the pseudo potential generation. Comparing to the unoptimized and optimized minimal orbitals SN and SN${}^{\prime }$, it is found that the geometrical parameters are significantly improved without increasing the computational time. In case of the optimized orbitals SNP${}^{\prime }$ remarkable improvements are obtained in both the geometrical parameters and the computational time. The optimized orbitals SNP${}^{\prime }$ provide a convergent result comparable to TNDP with a great reduction of the computational time. The computational time required for the orbital optimization of SN${}^{\prime }$ occupies only 3 % of that of the whole calculation. So the orbital optimization can be regarded as a preconditioning before doing the geometry optimization or the molecular dynamics. Of course, it is possible to perform the orbital optimization during the geometry optimization. It is worth mentioning that the orbital optimization can be combined with an O($N$) method [1, 2, 3, 4, 5, 6], since only ${\mathrm{\Theta }}_{i\alpha ,j\beta }$ and $E_{i\alpha ,j\beta }$, which are calculated by the O($N$) method, are required in Eq. (5). Therefore, the orbital optimization can be applied to large-scale systems in O($N$) operations.

Table 1: Optimized geometrical parameters (Å and degrees) of the most stable conformer of a neutral glycine molecule. A denotes the atomic symbol C, N, or O. The computational time per MD step was measured using one CPU on a Sharp Mebius PC-GP1-C7U. The energy cutoff of 113 (Ryd) was used for the numerical integrations in all calculations. The results by the other theory were taken from Ref. [15], and the experimental values from Ref. [16].

	SN	TN	TNDP	SN${}^{\prime }$	SNP${}^{\prime }$	Other theory	Expt.
	H4.0-s1	H4.0-s3	H4.0-s3p2	H4.0-s31${}^{*}$	H4.0-s31${}^{*}$p21${}^{*}$	LDA/DZP
	A4.5-s1p1	A4.5-s3p3	A4.5-s3p3d2	A4.5-s31${}^{*}$p31${}^{*}$	A4.5-s32${}^{*}$p31${}^{*}$d21${}^{*}$	Full potential
$r$(C$-$C)	1.555	1.535	1.528	1.515	1.528	1.510	1.532
$r$(N$-$C)	1.530	1.480	1.490	1.502	1.444	1.439	1.469
$r$(C$=$O)	1.353	1.231	1.235	1.281	1.238	1.218	1.207
$r$(C$-$O)	1.498	1.365	1.349	1.416	1.350	1.348	1.357
$r$(O$-$H)	1.144	0.998	0.987	1.010	0.995	0.988	-
$\mathrm{\angle }$(NCC)	104.8	106.5	108.2	108.4	108.8	114.8	112.1
$\mathrm{\angle }$(CC$=$O)	136.8	127.9	128.3	128.9	125.9	124.9	125.1
$\mathrm{\angle }$(COH)	96.8	107.9	105.7	105.8	106.8	105.6	-
C=O$\mathrm{\cdots }$N	3.132	2.998	2.905	2.998	2.882	2.827	-
Energy (Hartree)	-55.662	-55.981	-56.106	-55.818	-56.036	-	-
Time(s)/MD step	32	86	217	34	84	-	-

5 CONCLUSIONS

To conclude, we have developed a simple and practical method based on the force theorem for variationally optimizing numerical atomic orbitals used in density functional calculations. 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 illustration of geometry optimization with the orbital optimization for a small molecule clearly shows that the small optimized orbitals promise to greatly reduce the computational effort with a high degree of accuracy.

ACKNOWLEDGMENTS

We would like to thank H. Kino for helpful suggestions and discussions. This work was partly supported by NEDO under the Nanotechnology Materials Program.