CC BY
19
Development of an orbital-free approach for simulation of multi-atomic nanosystems with covalent bonds
V. G. Zavodinsky1, O.A. Gorkusha2
institute for Material Science, 153 Tikhookeanskaya str., 680042 Khabarovsk, Russia 2Institute of Applied Mathematics, 54 Dzerzhinskogo str., 680000 Khabarovsk, Russia
vzavod@mail.ru, o_garok@rambler.ru PACS 03.65.-w DOI 10.17586/2220-8054-2016-7-3-427-432
On the example of the three-atomic clusters Al3, Si3, and C3, it is shown that an orbital-free version of the density functional theory may be used for finding equilibrium configurations of multi-atomic systems with both metallic and covalent bonding. The equilibrium interatomic distances, interbonding angles and binding energies are found to be in good agreement with known data.
Keywords: Orbital-free, density functional, covalent bonding, angle dependence.
Received: 9 November 2015 Revised: 17 December 2015
1. Introduction
A number of works (for instance [1-8]) have been devoted to the development of the orbital-free (OF) version of the density functional theory (DFT) [9-13] over the last few years. This approach, in contrast to the Kohn-Sham (KS) method, does not use the wave functions (orbitals) and operates with the only electron density. The OF approach is the result of further development of the ideas of Hohenberg-Kohn [10]; i.e., that the basic state of a quantum system can be completely described by means of electronic density. There has been significant progress in the description of diatomic systems [8,14,15] and simple crystals [16]. However, an essential hindrance to further development of the OF method is the fact that the electronic density of a single (isolated) atom is spherical, i.e. the "orbital-free" atom has a shape of a ball, but balls form the close packed structures. For example, three identical atoms are obliged to form an equilateral triangle with corners of 60 degrees. At the same time, it is known that three atoms of silicon form an isosceles triangle with the main corner of about 80 degrees [17], atoms of carbon build a linear chain [18], and atoms of aluminum really behave like balls - they form a correct equilateral triangle [19]. The present work is an attempt to develop a technique which would allow us to adequately describe the geometry of interatomic bonds in polyatomic systems within the OF approach.
2. A general description of the OF approach
As is well-known the DFT claims that the energy Eel of the ground state of any quantum system can be found by minimization of some function which depends only on the electronic density of this system p:
Eei [p] = Ekin[p] + Eex[p] + Ec[p] + Eh [p] -J Vext(r)p(r)d3r, (1)
where Vext is an external potential, Ekin is kinetic energy, Eex is exchange energy, Ec is correlation energy, and EH is Hartree energy:
1 i p(r)p(r') jS.JS. |r — r'|
The total energy Etot is given by integral:
F r 1 1 f p(r)p(r>) ,,S jS ,
EH [p]=2j |r — r/| d rd r •
Etot = J Eei[p(r)]dSr.
(2)
Minimization of (1) with the condition J p(r)dSr = N means solving the following equation:
dEel [p] n (3)
7--M = 0, (3)
dp
where p is the Lagrange parameter gives one a sense of the electron chemical potential.
Introducing F [p] = SEe [p] - m, we obtain the equation:
Sp
F [p] = -Vext(r) + ^(r) + Mfcin(p) + Mex(p) + Mc(p) - M = 0, (4)
U , N i p(r') ,3 / , N SEkin[p] , s SEex[p] . , SEc[p]
where ^(r) = ^ TT| d r ' Mkin(p) = Sp , Mex(p) = ~¡j""' Mc(p) = .
There are some realistic approximations for exchange Mex(p) and correlation Mc(p) potentials; the potential of the electron-electron repulsion y>(r) may be calculated using Fourier transformations or Poisson equations; the external potential Vext(r) usually consists of atomic potentials or of pseudopotentials. The key problem is to find the potential of kinetic energy - Mfcm(p). In the Kohn-Sham approach, this problem is absent because the kinetic energy is calculated using electron orbitals (wave functions).
Quantum mechanical pseudopotentials are usually constructed for different angular states. Thus, we have to present the total density as a sum of partial densities:
p = ps + pp + pd + .... (5)
For the s-p case, we may write the equations:
Fs[ps,p] = -Vs(r) + ^(r)+ Mfcin(ps)+ Mex(p)+ Mc(p) - Ms =0,
Fp[pp, p] = -Vp(r) + ^(r) + Mkin(pp) + Mex(p) + Mc(p) - Mp = 0
where Vs(r) and Vp(r) are the s, p components of atomic pseudopotential. The electrostatic potential y(r), exchange and correlation potentials Mex(p) and Mc(p) are calculated through the total density p while partial kinetic potentials Mkin(ps) and Mpjn(pP) depend on corresponding partial densities ps and pp.
As designing of pseudopotentials is followed by the finding of equilibrium pseudo-wave functions, it is always possible to calculate partial densities of the isolated atom ps(r) and pp(r), as well as its full equilibrium electronic density p(r) = ps(r) + pp(r). If we know the type of these functions we could calculate energy Ekin and find the total energy of a single atom using Eq. (2). However, there is no basis to believe that there are some functions for kinetic functional and energy, which could be used for atoms of any type and any quantity. Moreover, as it was recently shown [20], the Hohenberg-Kohn idea about the existence of a universal density functional leading to the energy minimum was not strictly proved, and it is possible to say only about approximate solutions for the problem. However, in our case, the problem consists not in finding of the energy of a single atom (this can be done in more traditional ways), but in calculation of the energy of the interatomic interaction, and in finding the geometry of the polyatomic system corresponding to this interaction.
3. From dimers to trimers: a role of quantum rules and restrictions 3.1. Dimers
For simplicity, we will consider dimers consisting of atoms of one type. The elementary approach for the electronic density of such dimer pdim is the sum of densities of the atoms:
pdim(r) = poi(r - Ra)+ pai(r - RB ), (7)
where RA and RB are coordinates of points in which the A and B atoms with densities pat are situated. Then, the binding energy (per one atom) would be calculated as follows:
Eb =1 (Edim - 2Eat) , (8)
where Eat = E^f = EeVt(r)]d3r,
Edim = J Ee1[pdim(r)]d3r +
ZaZb (9)
|RA - RB 1 '
ZA and ZB are the positive charges of atomic cores equal to absolute values of charges of valence electrons.
We took Al, Si, and C as test elements (for the reasons stated above: trimers of these elements have essentially different geometric configurations). For all three elements we accepted the following function Mfcin(p) as the universal function for s-and p-states and:
Mkin(p) = 0.9p1/3 - 15p2. (10)
We used the FHI98pp [21] package as a generator of pseudo-potentials. We calculated exchange and correlation potentials in the local density approach [19,20]. The studied atoms were located in a cubic cell of the L size
Table 1. Equilibrium distances d and binding energies Eb (absolute values, per atom) for Si2, Al2 and C2 in comparison to known calculated data
Dimer Source of data d,A Eb, eV
Our method 2.1 2.0
Si2 Other calculations 2.23° 1.97°
2.21b 1.599b
Our method 2.4 1.2
Al2 Other calculations 2.46c 1.0d
2.95e 1.23e
Our method 1.3 4.8
C2 Other calculations 1.247 - 1.367f 4.7f
1.316s 3.5s
Notations: ° [17], b [25], c [23], d [19], e [27], f [28], s [29]
(L = 52 a.u.; 1 a.u. = 0.529 A). The cell was divided into 128x18x128 elementary sub-cells for the integration with the integration step AL of 0.406 a.u. Results of calculations were compared to the published data.
Calculated values of interatomic distances and binding energies for the Al2, Si2, and C2 dimers are collected in Table 1 in comparison with known published data. It is clear that agreement is rather good.
3.2. Trimers
To describe the angle dependence of interatomic bonding, we must analyze the reasons for this dependence in the standard quantum-mechanical approach, which uses wave functions and electronic states. For example, it is specified in the work [24] that the angle peculiarities of a cluster Si3 are defined by the Yang-Teller effect, which is caused by the existence of an energy gap between the occupied and empty states. In other words, the differences in semiconductor and metal small cluster structures are connected with the difference of their bond wave functions: namely, covalent atoms have localized functions oriented between nearest atoms, while metallic atoms have dispersed functions without spatial orientation.
In our case, both wave functions and electronic states are absent, and therefore, we cannot speak about any energy gap. In the OF approach, we deal only with the electronic density which defines all energies and structures of the polyatomic system. However, the main quantum-mechanical rules still remain fair in this case. Besides the Schrodinger's (or Kohn-Sham) equations, out of which wave functions and electron states arise, we must not forget Pauli's principle which specifies that in one quantum state there can be only two electrons (without taking into account spin). In our case, this principle may be paraphrased in the following way: a covalent bond is formed by two electrons, the common wave function of which is localized in the space between two nearest atoms. It is obvious that the quantity of the electrons which are responsible for this bond doesn't change as the distance between atoms changes (if, of course, the bond isn't broken at all and the electronic structure isn't completely reconstructed). In the case of metals, the conduction states are close each other and electrons can easily "flow" from one state to another during the changing of atomic geometry.
The above-mentioned concepts may be reformulated in the language of the electronic density: the density integral (nint) between atoms with covalent bonding, has to remain its value with a change of distance between atoms; in the case of metal bonding, the integral nint can have any possible value.
Certainly, there is a question: on what space do we have to provide integration, and what do we have to do with intermediate cases, with atoms of different types? We will leave these questions for the future, and now we will try to explain the differences in the structures of covalent and metal systems using the following homogeneous clusters Al3, Si3 and C3 as examples.
It is obvious that the space of integration has to be rather local, and at the same time, it has to give us the information on quantity of the electrons included in a covalent bond. In the present work we used the space having a shape of the slab situated between two nearest atoms (Fig. 1) and oriented perpendicularly to the plane in which the trimer triangle is placed. The thickness of a slab was taken as 2AL. As the number of the integration points can be changed with a change in the trimer configuration, the value of nint was normalized to one point.
Fig. 1. The scheme for the arrangement of atoms in a trimer. Dashed lines show the space on which the electronic density is integrated for definition the number of the electrons involved in the covalent bond; a is a corner between bonds with identical lengths of d
The electronic density of the trimer pirim may be found as follows:
Pirim(r) = Pai(r - Ra) + Pai(r - RB ) + pai(r - RC ), (11)
where RA, RB, and RC are coordinates of points in which the A, B, and C atoms with densities pai are situated. The binding energy is:
Eb = 1 (Eirim — 3Eai), (12)
where
n /neh / M ,3 , ZAZB ZAZC ZB ZC
Eirim = Eel [ptrim(r)]d3r + —-— + —-— + —-—, (13)
J |RA - RB | |RA - RC | |RB - RC |
ZA and ZB are the positive charges of atomic cores equal to absolute values to charges of valence electrons.
"trimer
Let us take the interatomic bond in a dimer as a standard. Let us call the value P = — (the relation
^dimer v "ini
of nini of trimer to nini of dimer) as "the bonding strength" and accept that for covalent bonds, the value of P shouldn't exceed 1.0 with a change of distance between atoms. For metallic bonds, P can have any possible value.
Calculated values P for clusters Al3, Si3, and C3 are presented in Fig. 2 as functions of the angle between interatomic bonds in the case when restrictions on these values are absent. For each angle, we found the values of interatomic distances, which corresponded to the minimum of the total energy of the cluster. One can see that P is approximately 1.0 at a = 180 ° and increases when a approaches 60 The maximum value (P = 1.40) is observed for carbon, which has the smallest interatomic distances. Interatomic distances in aluminum and silicon are approximately the same, therefore it should come as no surprise that the "bonding strengths" for the Al3 and Si3 trimers are approximately equal.
We repeated the calculations using thicker slabs (4AL and 8AL) for integration of nini(dimer) and nini(trimer) values between atoms in dimers and trimers and we found that results changed by no more than by 2 per cent.
In Figure 3 (curve A), we present the results for calculations of binding energies for the Al3, Si3 and C3 trimers without restriction of the "bonding strengths". From these curves it is clear that in all three cases, the maxima of the binding energy (on the absolute value) correspond to triangular clusters, "bonding strengths" in which significantly exceed the corresponding values, characteristic for linear chains. This result looks natural for aluminum as its states have the metallic, not localized character; but for the clusters of Si3 and C3 having covalent bonds it is necessary to introduce restrictions on values of P stipulated above. We have taken into account this condition (P = 1) and found dependences of the binding energy on the angle between bonds in the Si3 and C3 clusters (Fig. 3, curve B). One can see that atoms of carbon seek to form linear chains, while for silicon, neither a linear chain, nor an equilateral triangle is energetically favorable; atoms of silicon prefer to form an isosceles triangle with the angle a of 80 degrees.
Equilibrium values for interatomic distances d, angles a, and binding energies Eb for the trimers Al3, Si3, and C3 are collected in Table 2 (calculated for the condition P = 1) in comparison with known data. One can see that comparison is good. Thus, we showed that an orbital-free approach is capable of correctly describing the orientations of interatomic bonds in atomic clusters, as well as values of interatomic distances and binding energies.
Certainly, it is interesting to compare our results for interatomic densities with results of standart DFT-KS calculations. For this purpose, we calculated "bonding strentghs" P for Si3 and Al3 using the popular DFT-KS code FHI96md [31] for the same triangles as were studied above. We have found that P were equal 1 (±0.02) for all cases for Si3, but it increased up to 1.3 for the Al3 equilateral triangle. These results are in excellent agreement with ours.
Fig. 2. "The bonding strengths" in the Al3, Si3 and C3 trimers as functions of the corner between interatomic bonds in the case when there are no restrictions on the interatomic electron density
Fig. 3. Calculated dependence of binding energy (on atom) on the angle between interatomic bonds for the Al3, Si3 and C3 trimers. A) Values are obtained without restrictions on electronic density in interatomic bonds; B) values calculated with the condition P = 1.0
Table 2. Equilibrium distances d, angles a and binding energies Eb (absolute values, per atom) for Si3, Al3 and C3 in comparison with known calculated data
Trimer Source of data a, deg d,A Eb, eV
Our method «80 2.1 3.7
Si3 Other calculations 77.8a 2.26b 2.51b
78.10c 2.177c 2.93d
79.6e
Our method 60 2.3 2.4
Al3 Other calculations 60f 2.50f 1.74g
60h 2.55h 1.96f
Our method 180 1.2 7.0
C3 Other calculations 180* 1.29 j 6.8*
180k 1.3k 5.01
1.3161
Notations: ° [17], b [25], c [30], d [31], e [32], f [27], g [19], h [33], 4 [18], j [30], k [28], 1 [29] 4. Conclusion
We showed that the use of the restriction principle for the interatomic density (following from Pauli's principle) allows us to correctly describe the angular dependences of the interatomic bonding in polyatomic clusters within the orbital-free version of the density functional theory. In particular, it is possible to show that for the Al3 cluster, the equilateral triangle is favorable; the Si3 trimer is characterized by an isosceles triangle with angles of 80 and 50 degrees, and the three atoms of carbon are present as a linear chain. Calculated equilibrium interatomic distances and binding energy values are in fair agreement with known data. As the problem of correctly describing the angles between interatomic bonds is a key point in the modeling of polyatomic systems, it is possible to consider that our work opens a direct way to design an effective method for modeling of big nanosystems and supermolecules.
We have to note that the consideration which is carried out above is directly applicable only to the systems consisting of identical atoms and requires a special development for application to more complicated systems.
References
[1] Wang Y.A., Carter E.A. Orbital-free kinetic-energy density functional theory. In: Progress in Theoretical Chemistry and Physics, Kluwer, Dordrecht, 2000, 117 p.
[2] Huajie Chen, Aihui Zhou. Orbital-Free Density Functional Theory for Molecular Structure Calculations. Numerical Mathematics: Theory, Methods and Applications, 2008, 1, P. 1-28.
[3] 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-13.
[4] Karasiev V.V., Trickey S.B. Issues and challenges in orbital-free density functional calculations. Computational Physics Communications, 2012, 183, P. 2519-2527.
[5] 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, 2013, 88, P. 161108-13(R).
[6] Wesolowski T.A. Approximating the kinetic energy functional Ts[p]: lessons from four-electron systems. Molecular Physics, 2005, 103, P. 1165-67.
[7] 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.
[8] Zavodinsky V.G., Gorkusha O.A. Quantum-Mechanical Modeling without Wave Functions. Physics of the Solid States, 2014, 56 (11), P. 2329-35.
[9] Kohn W., Sham J.L. Self-Consistent Equations including Exchange and Correlation Effects. Phys. Rev. A, 1965, 140, P. 1133-38.
[10] Hohenbeg H., Kohn W. Inhomogeneous Electron Gas. Physical Review B, 1964, 136, P. 864-871.
[11] Sivasathya S., Thiruvadigal D.J. The effects of defects on electron transport in metallic single wall carbon nanotubes. Nanosystems: physics, chemistry, mathematics, 2013, 4 (3), P. 405-408.
[12] Hariharan R.M., Thiruvadigal D.J. Effect of anchoring atoms on transport properties of a carbon-dimer based molecular junctions: a first principles study. Nanosystems: physics, chemistry, mathematics, 2013, 4 (2), P. 294-298.
[13] Polyakov E.A., Vorontsov-Velyaminov P.N. Exact classical stochastic representations of the many-body quantum dynamics. Nanosystems: physics, chemistry, mathematics, 2015, 6 (4), P. 501-512.
[14] Junchao Xia, Chen Huang, Ilgyou Shin, Carter E.A. Can orbital-free density functional theory simulate molecules? The Journal of Chemical Physics, 2012, 136, 084102(13).
[15] 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.
[16] Carling K.M., Carter E.A. Orbital-free density functional theory calculations of the properties of Al, Mg and Al-Mg crystalline phases. Modelling and simulation in materials science and engineering, 2003, 11, P. 339-348.
[17] Raghavachari K., Logovinsky V. Structure and bonding in small silicon clusters. Phys. Rev. Lett., 1985, 55, P. 2853-2856.
[18] Van Orden A., Saykally R.J. Small carbon clusters: spectroscopy, structure, and energetics. Chemical Review, 1998, 98, P. 2313-57.
[19] Feng-Chuan Chuang, Wang C.Z., Ho K.H. Structure of neutral aluminum clusters Aln(2 < n < 23): Genetic algorithm tight-binding calculations. Phys. Rev. B, 2006, 73, 125431(7).
[20] Sarry A.M., Sarry M.F. To the density functional theory. Physics of Solid State, 2012, 54 (6), P. 1315-22.
[21] 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.
[22] Perdew J.P., Zunger A. Self-interaction correction to density functional approximation for many-electron systems. Physical Review B, 1981, 23, P. 5048-79.
[23] Ceperley D.M., Alder B.J. Ground state of the electron gas by a stochastic method. Physical Review Letters, 1980, 45, P. 566-569.
[24] Tomanek D., Schluter M.A. Structure and bonding of small semiconductor clusters. Phys. Rev. B, 1987, 36, P. 1208-17.
[25] 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, 01012(7).
[26] Nayak S.K., Khanna S.N., Jena P.J. Evolution of bonding in AlnN clusters: A transition from nonmetallic to metallic character. Physical Review B, 1998, 57, P. 3787-90.
[27] Matrnez A., Vela A. Stability of charged aluminum clusters. Physical Review B, 1994, 49, 17464(4).
[28] Karton A., Tarnopolsky A., Martin J.M.L. Atomization energies of the carbon clusters Cn (n = 2 — 10) revisited by means of W4 theory as well as density functional, Gn, and CBS methods. International Journal of Interface between Chemistry and Physics, 2009, 107, P. 977-1003.
[29] Afshar M., Babaei M., Kordbacheh A.H. First principles study on structural and magnetic properties of small and pure carbon clusters (Cn, n = 2 — 12). Journal of Theoretical and Applied Physics, 2014, 8, P. 103-108.
[30] McCarthy M.C., Thaddeus P. Rotational spectrum and structure of Si3. Physical Review Letters, 2003, 90, 213003(4).
[31] Liu B., Lu Z.Y., et al. Ionization of medium-sized silicon clusters and the geometries of the cations. Journal of Chemical Physics, 1998, 109, P. 9401-09.
[32] Raghavachari K., Rohlfing C.M. Bonding and stabilities of small silicon clusters: A theoretical study of Si7-Siio. Journal of Chemical Physics, 1988, 89, P. 2219-34.
[33] Tse J.S. Electronic structure of the dimer and trimer of aluminium. Theoretical Chemistry (Journal of Molecular Structures), 1988, 165, P. 21-24.
[34] Beckstedte M., Kley A., Neugebauer J., Scheffler M. Density functional theory calculation for poly-atomic systems: electronic structure, static and elastic properties and ab initio molecular dynamics. Computational Physics Communications, 1997, 107, P. 187-205.
