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, one-particle 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 Kohn-Sham (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 non-zero 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 Kohn-Sham 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 Kohn-Sham 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. 2-7 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
exchange-correlation 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 exchange-correlation 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. 2-7.
![]() |
![]() |
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.5-
2
2, 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 exchange-correlation
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. 2-7.
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 non-monotonic 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 two-photon
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 exchange-correlation interaction. Thus, the discrepancy
between the theoretical prediction and the experiments must be
attributed to the poor description to exchange-correlation 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 trade-off 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. 2-7 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. 2-7 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 trade-off 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.0-
2
2
2 and O5.0-
2
2 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 trade-off 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.5-3
3
3 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
ill-conditioned eigenvalues, which are positive but almost zero,
appear in the eigenvalue spectrum. The division by the ill-conditioned
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 C
F
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 C
H
and C
F
.
In addition, we find that the
orbitals of the carbon atom
in C
F
more shorten than that of C
H
.
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 C
H
, and 0.86 and 2.00 in C
F
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 C
F
is due to increase of effective core potential for
electrons. The comparison between C
H
and C
F
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, Gly-X-Gly, are considered for the orbital
optimization, where X could be one of twenty
kinds of amino acid residues. The structure of a Gly-X-Gly 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 Gly-X-Gly 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 exchange-correlation 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.5-52
51
,
C5.0-
52
52
51
,
N4.5-
52
52
51
,
O4.5-
52
52
51
,
and
S6.0-
52
52
51
.
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,
Met-enkephalin[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 Kohn-Sham 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 C
F
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 pre-orbital
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 (ACT-JST).
. Paz,
D. S
nchez-Portal, and
E. Artacho, Phys. Rev. B 64, 235111 (2001);
J. M. Soler et al., J. Phys.:Condens. Matter 14, 2745 (2002)
and references therein.