On a possibility to develop a full-potential orbital-free modeling approach Текст научной статьи по специальности «Физика»

On a possibility to develop a full-potential orbital-free modeling approach

V. G. Zavodinsky1, O. A. Gorkusha2 institute for Material Science, Tikhookeanskaya str., 153, Khabarovsk, 680042, Russia 2Institute of Applied Mathematics, Khabarovsk Division, Dzerzhinsky, 54, Khabarovsk, 680000, Russia

vzavod@mail.ru

PACS 03.65.-w DOI 10.17586/2220-8054-2019-10-4-402-409

We studied a principal opportunity to develop a full-potential orbital-free method for modeling of multi-atomic systems using results of Kohn-Sham calculations for single atoms. We have obtained equilibrium bond lengths and binding energies for homoatomic dimers Li2, Be2, B2, C2, N2, O2, F2, Na2, Mg2, Al2, Si2, P2, S2, and Cl2 as well as for heteroatomic dimers CSi, CB, CN, CO, SiO, NO, AlO, AlC, and NaCl. We analyzed our results and concluded that they are coordinated with experimental data not worse, than the results received by means of full-electrons calculations by the Kohn-Sham method.

Keywords: orbital-free, density functional, full-potential.

Received: 25 May 2019

Revised: 17 July 2019, 17 August 2019

1. Introduction

Modern technologies - nanotechnologies, bio-engineering, drug production and so on - need powerful tools to predict properties of systems containing hundreds of thousands and millions of atoms. Traditional quantum-mechanical approaches, such as Hartree-Fock theory, Quantum Chemistry methods and the Kohn-Sham version of Density Functional Theory (DFT), do not provide an opportunity to operate with large amounts of atoms; their limits do not exceed thousand atoms, even using pseudopotentials. Methods of empirical potentials and Monte-Carlo allow operation with large systems but they do not provide reliability of results. On the other hand, opportunities for increa-seing computer speed are near their physical limit, thus it is useless to hope for solution of the problem by this way. An imperative need is a new modeling method, which would combine quantum-mechanical accuracy with a possibility of operating with enormous number of atoms. The idea of such method has arisen in 1964 when Hohenberg and Kohn formulated the theorem [1] that the ground state energy of any quantum system is completely defined by its electronic density. In the same work, they have declared that there is a certain universal functional E[p], which minimization has to lead to the equilibrium the electronic density p and total electron energy E. The functional E[p] has been written in the following form:

E[p] = J e(p)dr = J V(r)p(r)dr + 1 J p(r)p(r)dr + J £ex_c(p)dr + J ekin{p)d,r, (1)

///)

--- dr/ is the electro-

|r - r/|

static electron Hartree potential, eex-c(p) and ekin(p) are exchange-correlation and kinetic energies (per electron). Minimization of (1) means solution has the following equation:

5E[p] = J ^Sp(r) = 0, (2)

from where, according to the general variation theory, it follows:

de(p)

gp = V (r)+ <^(r) + pex-c(p)+ Pkin(p) = 0 (3)

with the condition J p(r)dr = N, where N is the number of electrons in the system,

3ekin(p) 8eex-c(p)

Pkin(P) = —dp—, Pex-c(P) =-dp--(4)

are so called kinetic and exchange-correlation potentials.

The Hartree potential y(r) may be calculated using Fourier transformations or Poisson equations, the external potential V (r) usually consists of atomic potentials or pseudopotentials. There are some realistic approximations

for the exchange-correlation potential pex-dp) (for example [2-4]). The only real problem is the kinetic potential [¿kin (p), or the kinetic energy £kin(p).

There were attempts to use the Thomas and Fermi (TF) approximation [5,6] based on the free electron approach:

€Fn(.p) = M2/3 p5/3, plFn(p) = 2 (3n2)2/3 p1/3. (5)

This functional was found to be absolutely inadequate (all molecules were unstable); the von Weizsacker (vW) correction AeYin(p) [7] added to the TF:

AsWWn (P) = 8 / ^(f dr (6)

also did not solve the problem (binding energies were incorrect).

There were no other serious innovations in this area therefore Kohn and Sham have offered a compromise approach [8]. They have proposed to find the kinetic energy Ekin by solving some one-electron equation, the Hamilton-ian of which depended only on the electronic density:

- ^A^i(r) + Veff (r)Vi(r) = £iV>i(r), (7)

Veff (r) = V (r) + ^(r) + pex-c(p), (8)

2 ^

Ekin = -1 J ^ ^i(r)A^i(r)dr = ^ £i -J ^ ^i(r)Veff (r)^(r)dr, (9)

where p(r) = ^ |^i(r|2, ^i(r) are wave functions, or orbitals, ei - energy of ¿-state. Equation (7) is called the

i

Kohn-Sham (KS) equation.

The KS equation has become widely known, on its basis, many effective computing programs were created, many tasks of modeling of polyatomic systems were solved; however, as has been mentioned above, its opportunities are currently exhausted. The orbital-free (OF) approach as a version of the density functional theory could present an alternative to the KS method. It is a consecutive development of the idea of Hohenberg-Kohn [1] that the basic state of a quantum system can be completely described by means of electronic density. The advantage of this approach is obvious: operating only with an electronic density instead of numerous wave functions, it allows to increase sharply the speed of calculations and to include in consideration huge number of atoms.

The first attempts to develop an OF method of modeling began about 20 years ago. Those were simulations of liquid metals in the jellium approximation [9,10]. Then works of other researchers (see for example reviews and original articles [11-14]) were applied to some simple molecules and solids. All of these were based on the use of special pseudopotentials, with most of them trying to use the TF and vW approximations for the kinetic energy in different combinations. However, these attempts did not have a large degree of success and were not widely adopted. It seems to us that the main reason for their inadequate efficiency is that they try to use some universal functional of kinetic energy for all systems. However, it was recently shown [15,16], the Hohenberg-Kohn idea about existence of the universal density functional leading to the energy minimum was not strictly proved. Therefore, the search of specific types of kinetic potentials is necessary and justified.

In our recent works [17-21], we described an orbital-free pseudopotential approach for modeling of nanosystems containing atoms with s,p and d electrons. The key point of the approach was finding of the kinetic energy using some functions special for each type of atoms. This approach was tested on clusters contained C, Al, Si, O, Ti, and Cu atoms and demonstrated good agreement with Kohn-Sham method and experimental data. However, construction of pseudopotentials is a rather ambiguous operation and brings this method closer to the category of semi-empirical ones. In our present work, we describe an attempt of the orbital-free full-potential approach working without pseudopotentials.

2. General points of the approach

First, let us consider a single atom of any type A, which equilibrium total electron density p(r) may be calculated easily by the KS method. According to (3), we can write an equation for finding the single-atom kinetic potential

Mkin :

^¿(r; A) = -^ARf - ^(r) - [ex-c(r), (10)

where Z (A) is the total nuclear charge, R is an atomic position. If we know pk^ as a function of r and p as a function of r we can find r(p), substitute it into m^M and obtain [^(p). Using p(r) obtained from total potential KS calculations (FHI98pp [22]), we constructed kinetic potentials for B, C, N, O. They are plotted in Fig. 1.

c

zs

_Q

s_

03

M

c CD

1

o

CD c

2

10000 8000 6000 4000 2000 0

14 12 10 8 6 4 2 0

Li; O I N 1

■ 1 \ ■ 1 ■ 1 0 \ 50 100 1 1 1 1 1 1 1 1 < 150 200 250 300

. 1 . 1 . 1 . 1 . 1 . 1 . 1 \ \ r O N --

0,00 0,25 0,50 0,75 1,00 1,25

Density, arb. units

1,50 1,75

2,00

Fig. 1. Kinetic potentials calculated for equilibrium single-atomic full-electrons densities. The top panel demonstrates the total view of (p) for different atoms; the down panel shows ^¿i (p) for small densities

One can see that plotted curves differ from each other's very much. In accordance with works [15,16], there is no universal kinetic potential for different quantum systems. However, a question arises: If we have the single-atomic kinetic potentials for an each kind of atoms why we cannot use them for polyatomic systems? What is the difference between a single atom and for example a dimer?

Let us consider a dimer consisted of two boron atoms using the full-potential KS code Elk [23]. The equilibrium electron density of this dimer is plotted in Fig. 2.

Distance, at. units

Fig. 2. The equilibrium electron density of the boron dimer. The interatomic distance is 3.0 A

The kinetic potential isn't calculated in the KS approach. But we can determine it. Because we have the equilibrium KS density, we can find the two-atomic kinetic potential (r) for the boron dimer, according to (3):

№kin

(r) = -

Z (boron) Z (boron)

- ^(r) - ^ex-c(r)

(11)

,r - Ri| |r - R2I

where Z(boron) = 5; Ri, R2 are coordinates of 1 and 2 atoms of a dimer; ^(r), fiex-c(r) are calculated using the equilibrium KS density.

It is clear from Fig. 2 that in the case of the dimer the two different types of regions exist. In the first one (I) the density grows from zero at the large distance from the atom to the maximum value at the atomic nuclear. The second region type (II) means the space between the atoms. Our calculations showed (Fig. 3) that in the I region the two-atomic kinetic potential behaves identically the one-atomic potential. It is naturally because an electron in this

region is far from the atom 2 and interacts with the atom 1 only. In the II region an electron interacts with the both atoms; the kinetic potential decreases, and the level of the decreasing depends on the interatomic distance: less the distance less the kinetic potential.

Density, arb. units

Fig. 3. Behavior of the kinetic potential in different regions of the boron dimer

Obviously, the direct way to provide correct results for orbital-free calculations is to fulfill KS calculations for different atomic positions, to find the kinetic potential in the each space points, and to use it for OF calculations. However, such way is very expensive and has no sense. It seems to us, the optimal method is to find some regularities from the KS calculations on the dimer level and extend them to polyatomic systems in the OF approach. The simplest way is to introduce some median kinetic potential acting in the multi-atomic system, leading to the correct binding energy and atomic configuration. We constructed the following fitting expression for the kinetic potentials for any homoatomic dimers :

phi

homo dimer (

r; A) = Phi (r; A)

1.0 - a exp(-ß • Zlal(A))

(12)

where A is the type atoms, d is the dimer length, a and p are fitting constants, Zvai (A) is the number of valence electrons in the dimer atom. Parameters a and p control values of the binding energy and the bond length. If they are fitted for one system (for B2 for example) and are satisfied then for other dimers, this stage of our approach will be successful and we can pass to modeling more complicated systems.

As different atoms have different functions for kinetic energy, some procedure to calculate the total kinetic functions in the space of the atomic system has to be developed. Near each atom it has to be approximately equal to its specific function, but it has to be equal to mixture of the specific atomic functions between atoms. It seems that, in the case of a heteroatomic dimer, the simple way to construct the total function plth^0 dimer (r) is to summarize the specific atomic kinetic functions with some weights:

p hetero dimer phin

(r) = Wi(r) • Ph

homo dimer

(r; Ai) + W2(r) • ph

homo dimer

(r; A2),

Wi(r)

exp ( - ^)

exp ( - ^) + exp ( - ^)'

W2(r) =

exp ( - ^)

exp ( - ^) + exp ( - ^)'

(13)

(14)

(15)

where Ai ,A2 are the types of atoms in the dimer, Ri, R2 are coordinates of 1 and 2 atoms, Bi, B2 are parameters, with which the Gauss functions limit areas of valence electrons, and the equilibrium values of d are used for

Ph

homo dimer homo dimer

(r, Ai) and pkin (r, A2) here.

3. Details of calculations

The full-electrons consideration of atoms meets some peculiarities, which makes it rather difficult. One of them is a sharp intensive peak of the electron density centered on the atomic nucleus. These peaks correspond to localized core states, which do not participate in interatomic interactions and usually are considered as "frozen" (for example in the package Elk [23]. Such dividing of electron density helps to avoid awful operations with intensive sharp peaks and to construct realistic computer codes. Here, we also follow that technique and divide the atomic density to core

and valent components. In particularly, we consider that the B atom has two core electrons and three valence ones and their densities distributions are shown in Fig. 4.

U,U I | -.-1-.-1-.-1-.-1-.-1-.-1-.-1-

0,0 0,2 0,4 0,6 0,8 1,0 1,2 1,4 1,6

Distance, at. units

Fig. 4. Electron density of a single boron atom. The solid line (Total) demonstrates the total atomic density, the dashed line presents the density of the 2s2 2p1 electrons, the Valent-OF line shows the valent density used in our OF approach

Let us introduce a function F12 (r) for a homoatomic dimer:

Fi2 (r) = iT-A^ + jZ-lbl + ^12 (r) + ^(P12) + Pex-c(Pi2), (16)

where Z(A) is an atomic nuclear charge, pi2(r) = piore(r) + p2ore(r) + pV!l(r),

/// )

r12—t, pex-c(pi2) is calculated using some well-known approaches (LDA in our case).

|r - r/1

Because we started with equilibrium states for single atoms 1 and 2 we can, according (10), write

_HAL -Ta(r) - Mni(pa) - Mex-c(pa), (17)

|r - Ri|

-^M - - Pex-M, (18)

r R2

where pa and pa are equilibrium electron densities of the 1 and 2 atoms, t1 and t1 are electrostatic potentials formed by these densities.

Substituting (17) and (18) into (16) we obtain

Fi2(r) = T12 (r) - t1 (r) - t1 (r) + Mk2ii(pi2) - (p1) - rfl(pa) + Pex-c(pi2) - Hex-M - Hex-M- (19)

Our purpose is to find such density p12, which would turn (19) into zero.

Taking into account that the core densities do not change due to the interatomic interactions, we can write the iteration equation for valent density pH1 in any point r.

pial(r; i) = pial(r; i - 1) + Kiter ■ Fi2(r; i - 1)pial(r; i - 1), (20)

where Kiter is an iteration parameter controlling the procedure convergence. The start step (i = 0) means that the density pVI1 (r; 0) is a sum of equilibrium atomic valence densities: p\ ^ (r; 0) = (pVaI(|r-R 1 |))a + (pi al (|r-R2|))a.

If we have found the equilibrium electron density we can calculate all components of the total energy. Let us remind that the total energy of the dimer (Et0t)12 is a sum of the nuclear-nuclear repulsive energy (Erep)12, the

Coulomb energy (EC)12, the Hartree energy (EH)12, the exchange-correlation energy (Eex-c) 12 and the kinetic energy (Ekin) 12 :

Z(A) ■ Z(A)

(Erep) 12

|Ri - R2I '

(Ec) 12 — — J ^

Z (A) + Z (A)

|r — Ri| |r — R2I

• P12 (r)dr,

(EH) 12 = 1 J ^(P12) ■ Pl2(r)dr, (Eex-c)12 = / = / pex-c(p12)dp1

(Ekin) 12 = J e^P^M^ ek2n(P12) = J Pkin (P12 )dP12 .

(21) (22)

(23)

(24)

(25)

4. Calculations and discussions 4.1. Homoatomic dimers

Let us not foorget that we chose the B2 dimer as a test object. We have found for it the binding energy Eb = 1.8 eV and the equilibrium length d = 1.59 A taking parameters a = 1.08, and ft = 0.13. Then we used these values for all dimers of the Li-F and the Na-Cl rows. Results for binding energies are shown in Fig. 5; the dimer lengths are collected in Table 1. Experimental data are taken from the book [24]. The cited book contains experimental data taken from different sources; they often differ from each other very significantly. In these cases we averaged them and gave average values with two numbers after a decimal point. Our results are yielded with one number after a decimal point, because now we cannot guarantee higher precision for technical reasons. First of all, it is connected with the fact that the equilibrium states of the modeled systems are found by "manual" change of distance between atoms. In this work, we changed interatomic distances by the increments of 0.01 A and we observed that we cannot guarantee the energy accuracy more than 0.1 eV. We supposed that such accuracy is sufficient for comparing our results with experimental data and the KS results.

We calculated mean absolute energy and distance deviations of OF and KS approaches with respect to reference experimental data. They are 0.6 eV, 1.2 eV and 0.09 A, 0.03 A for OF and KS, respectively. Thus we see that our approach provides better energy accuracy than the KS approach, however it is less accurate when finding interatomic distances.

We did not attract the published data to comparison with our results. The matter is that these data, as a rule, significantly differ from each other because they were obtained by different methods within different approaches and approximations. Instead of this, we took a modern code (Elk [23]) based on the full-electrons KS-DFT approach and fulfilled calculations for systems of interest to us using the same approach for the exchange-correlation interaction that we used in our OF calculations.

Table 1. Equilibrium bond lengths d(A) for studied dimers

Method Li2 Be2 B2 C2 N2 O2 F2 Na2 Mg2 Al2 Si2 P2 S2 Cl2

OF 3.02 2.43 1.59 1.16 1.00 1.16 1.59 3.02 4.02 2.33 2.16 1.91 1.91 1.85

KS 3.07 2.46 1.63 1.23 1.11 1.18 1.37 3.07 3.89 2.54 2.22 1.83 1.90 1.96

Exp. [24] 3.09 2.47 1.59 1.24 1.10 1.15 1.417 3.07 3.891 2.56 2.32 1.90 1.88 1.99

Notations : "OF" are our OF calculations,

"KS" means full-potential calculations using the Elk [23]

It is clear from Fig. 5 that our OF energy results are in good agreement with the experimental data in many cases, they demonstrate even better agreement than KS calculations. The OF dimer lengths (see Table 1) also agree well with experimental ones, although this agreement is sometimes worse than in the case of the KS calculations. In general, we can say that our OF approach demonstrates a rather good ability to describe the interaction of atoms, especially taking into consideration the fact that parameters a and ft were fitted only for one type of dimers (B2).

Fig. 5. Binding energies for studied dimers. "Experiment" [24], "KS-full-potential" means calculations using the Elk code [23], "OF" are our OF calculations

4.2. Heteroatomic dimers

The dimers CSi, CB, CN, CO, SiO, NO, AlO, AlC, NaCl were taken as heteroatomic pairs interesting for technical and chemical applications. The Gaussian weight parameters B for corresponding atoms are shown in Table 2. Calculated equilibrium bond lengths and dissociation energies are collected in Table 3. As in the case of homoatomic dimers, we compare our results with experimental data and results of KS full-electrons calculations (Elk code [23]).

Table 2. The values of the parameters B for weight functions

C Si B N O Al Na Cl

B (a.u.)2 1.5 3.0 1.6 1.3 1.3 4.0 7.0 2.6

Table 3. Dissociation energies Ed and bond lenghts d for heteroatomic dimers

Source C-Si C-B C-N C-O Si-O N-O Al-O Al-C Na-Cl

Ed, eV OF 5.7 6.4 12.9 12.7 6.7 7.9 4.9 3.9 4.2

KS 5.9 6.32 13.35 14.7 15.55 10.95 2.3 1.69 5.07

Exp. [24] 4.6 4.6 7.7 11.09 8.26 6.5 5.2 2.8 4.23

d, A OF 1.8 1.4 1.2 1.3 1.7 1.1 1.7 1.9 2.4

KS 1.74 1.43 1.11 1.15 1.53 1.08 1.64 1.90 2.41

Exp. [24] 1.7 1.49 1.15 1.11 1.60 1.15 1.62 1.96 2.36

Analysis of Table 3 shows that our OF approach describes interactions of studied atoms not worse than KS Elk-code.

5. Conclusions

In this work, we demonstrated that it is possible to design a full-potential orbital-free approach for modeling of atomic systems using one-atomic kinetic potentials obtained from Kohn-Sham calculations. We proposed a practical way to construct the two-atomic kinetic potential for the B2 dimer and used successfully this way for dimers from

Li2 up to Cl2. We have generalized our OF method on heteroatomic systems. For this purpose, we used some weight functions and calculated equilibrium energies and bond lengths for dimers CSi, CB, CN, CO, SiO, NO, AlO, AlC, NaCl. The analysis of our results shows that they will be coordinated with experiments not worse than results of full-electrons KS calculations. Certainly, for development of this method in full, it is necessary to overcome a set of difficulties: namely, to construct kinetic potentials for multi-atomic systems, to develop a relaxation procedure, to develop a spin dependent version of the orbital free approach and the OF approach atoms with d-electrons. However, these problems seems to us surmountable, and an opportunity to model very big nanosystems - of hundreds thousands atoms and more - will be as a result had as an award.

References

[1] Hohenberg H., Kohn W. Inhomogeneous Electron Gas. Physical Review, 1964,136, P. B864-B871.

[2] Perdew J. P., Zunger A. S. Self-interaction correction to density functional approximation for many-electron systems. Physical Review, 1981, 23, P. 5048-5079.

[3] Ceperley D. M., Alder B. J. Ground state of the electron gas by a stochastic method. Physical Review, 1980, 45, P. 566-569.

[4] [4] Perdew J. P., Wang Y. Accurate snd simple density functional for the electronic exchange energy. Physical Review, 1986, 33, P. 88008802.

[5] Thmas L. H. The calculation of atomic field. Proc. Cambr. Phil. Soc., 1927, 23, P. 542-548.

[6] Fermi E. Un metodo statistic per la determinazione di alcune priorieta dell'atomo. Rend. Accad. Lincei., 1927, 6, P. 602607.

[7] v. Weizsacker C. F. Theorie de Kernmassen. Z. Physik, 1935, 96, P. 431458.

[8] Kohn W., Sham J. L. Self-consistent equations including exchange and correlation effects. Phys. Rev., 1965,140, P. A1133-A1138.

[9] Garcia,- Gonzalez P., Alvarellos J. E., Chacon E. Nonlocal symmetrized kinetic- energy density functional: Application to simple surfaces. Phys. Rev., 1998, 57, P. 48574862.

[10] Gomez S., Gonzalez L.E., Gonzalez D. J., Stott M. J., Dalgic S., Silbert M.J. Orbital free ab initio molecular dynamic study of expanded liquid Cs. Non-Cryst. Solids, 1999, 250-252, P. 163-167.

[11] Wang Y. A, Carter E. A. Orbital- free kinetic- energy density functional theory. In: Theoretical Methods in Condensed Phase Chemistry. Schwartz, S.D., Ed. Springer, Dordrecht.: 2002, P. 117-184.

[12] Huajie Chen, Aihui Zhou. Orbital- free density functional theory for molecular structure calculations. Numerical Mathematics: Theory, Methods and Applications, 2008, 1, P. 1-28.

[13] 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.

[14] Karasiev V. V., Chakraborty D., Trickey S. B. Progress on New Approaches to Old Ideas:Orbital-Free Density Functionals. In: Many-Electron Approaches in Physics, Chemistry and Mathematics. Mathematical Physics Studies. Eds: Bach V, Delle S. L. Ed. Springer, Dordrecht.: 2014, P. 113-135.

[15] Sarry A. M., Sarry M. F. To the density functional theory. Physics of Solid State, 2012, 54(6), P. 1315-1322.

[16] Bobrov V. B., Trigger S. A. The problem of the universal density functional and the density matrix functional theory. Journal of Experimental and Theoretical Physics, 2013, 116(4), P. 635-640.

[17] Zavodinsky V. G., Gorkusha O. A. 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-47.

[18] Zavodinsky V. G., Gorkusha O. A. Development of an orbital free approach for simulation of multiatomic nanosystems with covalent bonds. NANOSYSTEMS: PHYSICS, CHEMISTRY, MATHEMATICS, 2016, 7(3), P. 427-432.

[19] Zavodinsky V. G., Gorkusha O. A. Development of the orbital free approach for heteroatomic systems.NANOSYSTEMS: PHYSICS, CHEMISTRY, MATHEMATICS, 2016, 7(6), P. 1010-1016.

[20] Zavodinsky V. G., Gorkusha O. A. New Orbital Free Simulation Method Based on the Density Functional Theory. Applied and Computational Mathematics, 2017, 6(4), P. 189-195.

[21] Zavodinsky V. G., Gorkusha O. A. Orbital- free modelling method for materials contained atoms with d- electrons. International Journal of Scientific Research in Computer Science, Engineering and Information Technology, 2018, 3(7), P. 57-62.

[22] 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.

[23] URL: http://elk.sourceforge.net.

[24] Huber K. R., Herzberg G. Molecular Spectra and Molecular Structure. IV. Constants of Diatomic Molecules. Litton Educational Publishing, N.Y.: 1979. 732 p.