CC BY
15
Development of the orbital-free approach for hetero-atomic systems
V. G. Zavodinsky1, O.A. Gorkusha2 institute for Material Science, Khabarovsk, 680042, 153 Tikhookeanskaya str., Russia 2Institute of Applied Mathematics, Khabarovsk, 680000, 54 Dzerzhinskogo str., Russia [email protected], [email protected]
PACS 03.65.-w DOI 10.17586/2220-8054-2016-7-6-1010-1016
The key problem of the orbital-free approach is calculation of kinetic energy, especially for hetero-atomic systems. In this work, we used the mono-atomic functionals of kinetic energy to construct the kinetic functionals of complicated systems. We constructed some atomic weights associated with densities of single atoms and then calculated kinetic functions for some atomic complexes. For the examples of SiC, SiAl, AlC, SiO and CO dimers we have demonstrated possibility of our approach to find equilibrium interatomic distances and dissociation energies for hetero-atomic systems.
Keywords: orbital-free, density functional, hetero-atomic systems, interatomic distances, dissociation energies. Received: 4 April 2016 Revised: 10 May 2016
1. Introduction
Nanotechnology requires simulation methods, which could operate with huge numbers of atoms - up to millions. The most effective quantum methods (for example the Kohn-Sham (KS) method [1]) can work with only hundreds. Therefore, researchers are obliged to use for large nanosystems some less accurate methods with empiric potentials (for example [2, 3]).
The Kohn-Sham method is based on density functional theory (DFT) [4]. The orbital-free (OF) approach also follows this theory, however, it operates with the electron density only (without wave functions or orbitals) and if properly developed, can be applied for the simulation of very large systems: up to millions atoms [5]. Several groups [5-13] are working in this area with different success, and the calculation of the kinetic energy is noted as a main problem. In our previous papers [14-16], we suggested that there is no universal way to describe the kinetic energy of different atoms and compounds. We proposed some simple formulas for systems containing atoms of identical types and simulated the dimers and trimers with metallic and covalent bonds. For examples of Al, Si, and C, we obtained equilibrium interatomic distances, binding energies and interbonding angles in good accordance with published data. Now we try to describe how it is possible to extend our approach to systems with different types of atoms.
2. A general description of the OF approach
As it is known, DFT claims that the energy E of the ground state of any quantum system can be found by minimization of the some functional depending only on the electronic density of this system p(r):
E[p] = J e(p)dr = J V(r)p(r)dr + 1 J ^(r)p(r)dr + J eex-c(p)dr + J efci„(p)dr, (1)
/// )
1--dr/ is the electrostatic electron potential Hartree, eex-c and
£kin are exchange-correlation and kinetic energies (per electron). Minimization of (1) means solution the following equation:
F [p] = ^ = V (r) + ^(r) + Mex-c(p) + Mfcin(p) = 0, (2)
Sp
where p have to satisfy the condition f p(r)dr = N, N is the number of electrons in the system, pex-c(p) = SSex-c (p) , \ S£kin(p) Sp , №n(p) = "ST".
There are some realistic approximations for exchange-correlation potential pex-c(p) there; the potential Hartree y(r) may be calculated using Fourier transformations or Poisson equations; the external potential V(r) usually consists of atomic potentials or pseudopotentials. The only real problem is the kinetic potential ^kin.
3. Pseudopotential approach
In practice, the DFT calculations are simpler if one uses pseudopotentials instead of full electron potentials. Therefore, let us rewrite the above equations in the pseudopotential approach, and, for simplicity, let us limit ourselves by s- and p-components of pseudopotentials and a diatomic system. Their distribution on more complicated cases is possible without any trouble. Thus, we will present the total density p12 as a sum of partial densities:
P12 = Pl2-s + P12-p + ...
The electron energy of this system E12 = /e12(p12-s, p12-p)dr must be minimal with the condition I (p12-s + P12-P)dr = N12, where e12 is the electron energy per electron for the two-atomic system with the total number of electrons N12. In the other words, we have to find the density p12 that satisfies the system of two equations:
F12-S = 0, F12-P = 0.
Here
F12-s = = V1-s(r) + V2-s(r) + ^12(r) + M?X-c(P12) + Mk2"-s(P12-s), (3a)
¿P12-S
F12-P = = V1-p(r) + V2-p(r) + ^12(r) + M?2-C(P12) + Ml2-p(P12-p), (3b)
OP12-p
where V1-s(r), V2-s(r), V1-p(r) and V2-p(r) are s and p components of pseudopotentials of the first and second atoms, y12(r) and M?2-c(p12) are the electrostatic and exchange-correlation potentials calculated for the total electron density p12 of a dimer, Mi2"-s(p12-s) and Mi2"-P(p12-P) are partial kinetic potentials depending on corresponding partial densities p12-s and p12-p.
Thus we can write equations for finding p12-s and p12-p:
V1-s(r) + V2-s(r) + ^12 (r) + M?rc(P12) + Mk2"-s(Ps) = 0, (4a)
V1-p(r) + V2-p(r) + ^12(r) + M?2-c(P12) + m22-p(Pp) = 0, (4b)
Obviously, for two isolated atoms we can write equations similar to (4a) and (4b):
V1-s(r) + ¿?(r) + M?x-c(p?) + m2-" (p?-s) = 0, V1-p(r) + ^?(r) + M?x-c(p?) + Mi-np(Pi-p) = 0 (5a)
V2-S(r) + ^l(r) + m2x-c(p2) + M2kins(p2-s) = 0, V2-p(r) + ^l(r) + m2x-c(p2) + M2fc->0-p) = 0 (5b) As p?-s, p?-P, p2-s, and are equilibrium atomic densities taken from DFT calculations, we can write for V1-S(r), V1-p(r), V--s(r) and Vs-p(r):
V1-S(r) = -^°(r) - M?x-c(p?) - M2-ns(P?-s), V1-p(r) = -^°(r) - M?x-c(p?) - M2-np(P?-p), (6a)
V2-s(r) = -^2(r) -M2x-c(p0) -m2*-s(p2-s), V2-p(r) = -^°(r) - M2x-c(p2) -m2->2-p). (6b) Putting (6a) and (6b) in (4) we obtain:
^12(r)-^a(r)-^a(r)+Mf2-c(P12)-Mfx-c(p?)-M2x-c(p2)+M22n-s(P12-s)-M2-ns(P?-s)-M22ins(p2-s)=0, (7a)
^12(r)-^a(r)-^a(r)+M?2-c(P12)-M?x-c(p?)-Mex-c(p0)+M22n-p(p12-p)-M2il>tp)-M2->2-p) = 0. (7b) The kinetic dimer functionals m22"-s(p12-s) and m22"-p(p12-p) may be presented as follows:
m22-s(P12-s) = M2-" (P?-s) + s(p2-s) + Am22-s(P12-s ), (8a)
Mk2-p(p12-p) = Mk-"p(p?-p) + Mk-np(p2-p) + am22-p(p12-p), (8b)
where Am22-s(p12-s) and A^22-p(p12-p) are unknown functions of partial densities of the two-atomic system. These functions must approach zero if the interatomic distance approaches to infinity. Thus we can take the following simple approximation for them:
Am22-S(P12-S « vA<ns(p)12-s - v£-s(p?_s) - vA-s(p°-s), (9a)
AMk2n-p(p12-p « vA-p(p12-p - vA-p(p?-p) - vAypl-p), (9b)
where vA-^Ps) and vA®np(pp) are some functions having the same kind for single atoms and for dimers formed from atoms of the same type A.
We can solve these equations with some fitting functions vA- and vA- and then calculate the total energy. We find the test functions from the simple request: they must lead to the equilibrium interatomic distances and binding energy for dimers. We hope that these functions will be suitable for more complicated systems in future.
The electron energy of a dimer El'im contains the electrostatic energy:
ElVtat = / [Vl-s(r) + V2-s (r)] Pl2-sdr + J [Vl-P(r) + V2-p(r)] Pi2-pdr ^12(r)p12dr, (10)
the exchange-correlation energy El2-c = / e!2-c(p12)dr and the kinetic energy:
Ei2n = E?in + E2kin + JJ [v1fc2-s(pi2-s)dpi2-s + V1fc2-s(pi2-p)dpi2-s]dr. (11)
To find the equilibrium distance and the total energy Etd°tm, we need to add the repelling energy E[Hp = Z Z
where Z1 and Z2 are positive charges of atomic ions with coordinates R1 and R2. Thus E^°m =
R- R2I'
El-stat+ El2x-c + E1™ + E[lp. The binding energy for a dimer (per one atom) would be calculated as follows: Eb = 2 (Edm - 2Ea), where Ea is the atomic energy.
4. Dimers with identical atoms
We took Al, Si, and C as test elements. We used the FHI98pp [17] package as a generator of pseudo-potentials and equilibrium partial electron densities. We calculated exchange and correlation potentials in the local density approach [18,19]. Studied atoms were located in a cubic cell of the L size (L=30 a.u.; 1 a.u. = 0.529 A). The cell was divided on 150 x 150 x 150 elementary sub-cells for the integration with the step AL of 0.2 a.u. The results of these calculations were compared with published data.
We used the same types of kinetic functions ^?in and for isolated atoms and dimers and trimers, however they were found different for different types of atoms. Namely, we used for Al: v?in = 1.0ps1/4'5; vpin = 22.0ppA'5; for Si: v?in = 8.0ps1/1'5; v1in = 1.6pp/3;
for C: v? in = 1.75ps1/3; v1in = 1.8pp/3.
Calculated values of interatomic distances and binding energies for the Al2, Si2, and C2 dimers are collected in Table 1 in comparison with other data. Agreement is rather satisfactory, when one considers that other calculated data are often differing from experimental results and each other.
Table 1. Equilibrium distances d and binding energies Eb (absolute values, per atom) for Si2, Al2 and C2 in comparison with known data
Dimer Source d, A Eb ,eV
Si2 Our OF method 2.2 1.8
Other calculations 2.21a 2.23b 1.599a 1.97b
Experiment 2.24c 3.0c
Al2 Our OF method 2.8 1.4
Other calculations 2.95d 2.511 1.23d 1.551
Experiment 2.56f 1.56f
C2 Our OF method 1.4 3.0
Other calculations 1.24g 1.36h 2.6g 2.7h
Experiment 1.24® 3.1®
Notations: a[20], b[21], c[22], d[23],e[24], f [25], g[26], h[27], ®[28].
5. Dimers with different atoms
Let us rewrite equations (8a) and (8b) for a dimer contained atoms of types A and B:
(12a) (12b)
The functions AMAiB-s(pAB-s) and A^Ai_n-p(PAB-P) have to be approximately equal to atomic functions
MaB-s(Pab-s) = MAms(pA-s) + MBms(pB-s) + AMAi^-s(PAB-s), MAB-p^AB-p) = MA-p^A-pp) + MB-p^B-p) + AMab-p(PAB-p).
kin
AMa-s(p) or AMBi"s(p) near atoms A or B, but they have to be mixtures of the atomic functions in the whole
and Mab-s(p) is to summarize
,kin
space. It seems to us that the simplest way to construct the functions Amab-s(p) the atomic functions with some weights:
Amab-s(pab-s) = wa-sAma-s(pab-s) + wb-sA^-sGpab-s);
AMaB-p(PAB-p) = WA-p AMa-"p(PAB-p) + WB -pAMBi-p(PAB-p).
For AmA®!!s(p) and AmB®-s(p) we propose the following approximations:
AMA-s(PAB-s) « vA-s(PAB-s) - V?-s (pA-s) - vA-s(pB-s), AmB®-s(PAB-s ) « vB-s(PAB-s) - vB-s(pA-s) - vB-s(pB-s),
ApAmp(pAB-p) ~ vAi"p(pAB-p) - vAi"p(pA-p) - ^A-ypB-p^
AmB®-p(pab-p) « Vki-p(pAB-p) - v
p
kin B-p
p/ "A-p(pB-p)
(pA-p) - vB%(pB-p
(13a) (13b)
(14a) (14b) (14c) (14d)
where
,kin 'A-s
yA-p> vB-s
vB- and vB- are functions related to atoms A and B. For Al, Si and C they are found in
the previous Section.
The weights WA-s, WA-p, WB-s and WB-P may be determined through Gauss functions fitted to atomic densities:
Wa-s
Wb_s
W^p
W^p
aA-s exp
(r-Ra)2
ftA —s
aA-s exp
(r-Ra)2
ftA — s
+ aB-s exp -
(r-RB )2
ftB — s
aB-s exp
(r-Rb )2
ftB — s
(r-Ra)2
ftA — s
aA-s exp "A/ ) + aB-s
aA-p exp
exp
(r-Rb )2
ft B — s
(r-Ra)2
, (r-Ra)2 \ i / (r-Rb )2
aA-p exp i-a^J + aB-p exp (- pBB—pr
aB-pexp
(r-Rb )2
ftB — p
(r-RA)^ + aR ex^_ (r-RB )2
ftA—p i + aB-p exp I ftB —p
aA-p exp
(15a)
(15b)
(15c)
(15d)
An example of fitting of the densities and weights is demonstrated in Figure 1. Values for a and P for Si, Al, C and O are presented in Table 2.
Table 2. Parameters of weight functions (a and P) for Si, Al, C and O atoms
k
Type of atom as ap Ps Pp
Si 0.065 0.040 3.5 4.5
Al 0.065 0.005 3.5 4.5
C 0.200 0.160 1.5 1.8
O 0.300 0.450 1.0 1.5
Fig. 1. The s-densities and the weights Ws for oxygen (dashed) and silicon atoms (solid)
We fulfilled calculations for the SiC, SiAl, AlC, SiO, CO, and AlO dimers with parameters shown in Table 2. The kinetic functions for oxygen = Up^3 — 1.0ps and v^®" = 1.5pp/3'5 — 1.0pp have been found through simulation of the SiO dimer and then they were used for other oxygen contained dimers. Results of calculations are presented in Table 3. Unfortunately, we did not find published data for the all studied heteroatomic dimers. Therefore, we compared our OF results also with results calculated by us in the framework of the KS DFT approach using the well-known package FHI96md [17].
Table 3. Equilibrium distances d and energies of dissociation (absolute values) for SiC, SiAl, AlC, SiO, CO, and AlO
Dimer Source d, A Ed,eV
SiC Our OF calculations 1.9 6.9
Published calculations 1.8° 7.7°
Our KS FHI96md calculations 1.69 6.66
SiAl Our OF calculations 2.5 3.8
Our KS FHI96md calculations 2.30 3.10
AlC Our OF calculations 2.0 6.1
Our KS FHI96md calculations 1.83 4.32
SiO Our OF calculations 1.6 7.0
Our KS FHI96md calculations 1.51 12.06
Experiment 1.52b 7.2b
CO Our OF calculations 1.0 9.5
Our KS FHI96md calculations 1.11 15.96
Experiment 1.13c 9.6c
AlO Our OF calculations 1.8 3.0
Our KS FHI96md calculations 1.55 9.0
Experiment 1.62c 5.27c
Notations: °[30], b[31], c[32].
One can see from Table 3 that our OF equilibrium distances slightly exceed the experimental ones as well as the KS calculation results, except the CO dimer. As for the dissociation energy, the OF results for the SiO and CO dimers are closer to experimental data than the KS ones. The OF result for SiC correlates with our KS and known calculated values (there is no experimental data). OF calculated energies for the Al contained dimers are rather far from experimental values and from results of KS calculations. The reason for this discrepancy requires future investigation; however, it is remarkable that the KS results are also far from experimental energies in many cases. As a whole, we can conclude that the OF method is able to give us a rather satisfactory information on interatomic distances and energies of systems containing different atoms.
6. Conclusion
We showed the possibility for simulating the interactions of atoms of non-identical types in the framework of the orbital-free version of the density functional theory. For this purpose, we used a rather simple technique, namely: first, the atomic kinetic functions were found for homo-atomic dimers Si2, Al2, C2 and for the SiO dimer; second, some atomic weights were proposed using Gaussians associated with atomic densities; third, kinetic functions for hetero-atomic dimers were constructed. Equilibrium interatomic distances and dissociation energies for the SiC, SiAl, AlC, SiO and CO dimers were found to be in satisfactory agreement with the Kohn-Shem calculations and experimental data.
As the calculation of the kinetic energy is a key point in the modeling of polyatomic systems in the orbital-free approach, it is possible to consider that our work opens a direct way to design an effective modeling method for complicated nanosystems and macromolecules with a large number of atoms.
References
[1] Kohn W., Sham J.L. Self-Consistent Equations including Exchange and Correlation Effects. Phys. Rev., 1965, 140, P. A1133-A1138.
[2] Konovalenko I.S., Zolnikov K.P., Psakhie S.G. Molecular dynamics investigation of deformation response of thin-film metallic nanostruc-tures under heating. Nanosystems: Physics, Chemistry, Mathematics, 2011, 2(2), P. 76-83.
[3] Smolin A.Yu., Roman N.V., Zolnikov K.P., Psakhie S.G., Kedrinskii V.K. Simulation of structural transformations in copper nanoparticles under collision. Nanosystems: Physics, Chemistry, Mathematics, 2011, 2(2), P. 98-101.
[4] Hohenberg H., Kohn W. Inhomogeneous Electron Gas. Physical Review, 1964, 136, P. B864-B871.
[5] Hung L., Carter E.A. Accurate Simulations of Metals at the Mesoscale: Explicit Treatment of 1 Million Atoms with Quantum Mechanics. Chemical Physics Letters, 2009, 475, P. 163-170.
[6] Wang Y.A., Carter E.A. Orbital-free kinetic-energy density functional theory. Progress in Theoretical Chemistry and Physics, Kluwer, Dordrecht, 2000, 117 p.
[7] Huajie Chen, Aihui Zhou. Orbital-Free Density Functional Theory for Molecular Structure Calculations. Numerical Mathematics: Theory, Methods and Applications, 2008, 1, P. 1-28.
[8] Baojing Zhou, Ligneres V.L., Carter E.A. Improving the orbital-free density functional theory description of covalent materials. Journal Chemical Physics, 2005, 122, P. 044103-044113.
[9] Karasiev V.V., Trickey S.B. Issues and challenges in orbital-free density functional calculations. Computational Physics Communications, 2012, 183, P. 2519-2527.
[10] Karasiev V.V., Chakraborty D., Shukruto O.A., Trickey S.B. Nonempirical generalized gradient approximation free-energy functional for orbital-free simulations. Physical Review B, 88, P. 161108-161113(R).
[11] Wesolowski T.A. Approximating the kinetic energy functional Ts[p]: lessons from four-electron systems. Molecular Physics, 2005, 103, P. 1165-1167.
[12] Junchao Xia, Chen Huang, Ilgyou Shin, Carter E.A. Can orbital-free density functional theory simulate molecules? The Journal of Chemical Physics, 2012, 136, P. 084102(13).
[13] Lehtomaki J., Makkonen I., Caro M.A., Haiju A. and Lopez-Acevedo O. Orbital-free density functional theory implementation with the projector augmented wave method. Journal Chemical Physics, 2014, 141, P. 234102(7).
[14] Zavodinsky V.G., Gorkusha O.A. A practical way to develop the orbital-free density functional calculations. Physical Science International Journal, 2014, 4(6), P. 880-891.
[15] Zavodinsky V.G., Gorkusha O.A. Quantum-Mechanical Modeling without Wave Functions. Physics of the Solid States, 2014, 56(11), P. 2329-2335.
[16] Zavodinsky V.G., Gorkusha O.A. New Orbital-Free Approach for Density Functional Modeling of Large Molecules and Nanoparticles. Modeling and Numerical Simulation of Material Science, 2015, 5, P. 39-46.
[17] Fuchs M., Scheffler M. Ab initio pseudopotentials for electronic structure calculations of poly-atomic systems using density-functional theory. Computational Physics Communications, 1999,119, P. 67-98.
[18] Perdew J.P., Zunger A. Self-interaction correction to density functional approximation for many-electron systems. Physical Review B, 1981, 23, P. 5048-5079.
[19] Ceperley D.M., Alder B.J. Ground state of the electron gas by a stochastic method. Physical Review Letters, 1980, 45, P. 566-569.
[20] Mukhtarov A.P., Normurodov A.B., Sulaymonov N.T. Umarova F.T. Charge States of Bare Silicon Clusters up to Si8 by Non-Conventional Tight-Binding Method. Journal of nano- and electronic physics, 2015, 7, P. 01012(7).
[21] Raghavachari K., Logovinsky V. Structure and bonding in small silicon clusters. Phys. Rev. Lett., 1985, 55, P. 2853-2856.
[22] Tomanek D., Schluter M.A. Structure and bonding of small semiconductor clusters. Phys Rev B, 1987, 36, P. 1208-1217.
[23] Martinez A., Vela A. Stability of charged aluminum clusters. Physical Review B, 1994, 49, P. 17464(4).
[24] Kumar V., Sundararajan V. Ab initio molecular-dynamics studies of doped magic clusters and their interaction with atoms, Physical Review b, 1998, 57, P. 4939-4942.
[25] Herzberg G. Spectra of Diatomic Molecules. Van Nostrand, New York, 1950.
[26] Chelikowsky J.R., Chou M.Y. Ab initio pseudopotential-local density description of the structural properties of small carbon clusters. Physical Review B, 37, P. 6504-6507.
[27] Fougere P.F., Nesbet R.K. Electronic Structure of C2. J. Chem. Phys, 1966, 44, P. 285-297.
[28] Huber K.P., Herzberg G. Molecular Spectra and Molecular Structure. IV. Constants of Diatomic Molecules. Reinhold, New York, 1979.
[29] Beckstedte M., Kley A., Neugebauer J., Scheffler M. Density functional theory calculations for poly-atomic systems: electronic structure, static and elastic properties and ab initio molecular dynamics. Computational Physics Communications, 1997, 107, P. 187-205.
[30] Patrick A.D., Xiao Dong, Allison T.C., Blaisten-Barojas E. Silicon carbide nanostructures: A tight binding approach. Journal of Chemical Physics, 2009, 130, P. 244704. http://dx.doi.org/10.106371.3157282
[31] Deskbook of the Chemist. Moscow, Chemistry, 1982, P. 336-341. http://wwwchemway.ru/bd_chem/tbl_mol/w_tbl_r_m_08.php.
[32] Andersen T., Haugen H.K., Hotop H. Binding energies in atomic negative ions: III. J. Phys. Chem. Ref., 1999, 28, P. 1511-1534.
