A Fast Density-Functional Method for Chemistry - American

Chapter 26

A Fast Density-Functional Method for Chemistry Xiao-Ping Li , Jan Andzelm , John Harris , and Anne M . Chaka Downloaded by STANFORD UNIV GREEN LIBR on October 8, 2012 | http://pubs.acs.org Publication Date: May 5, 1996 | doi: 10.1021/bk-1996-0629.ch026

1

1

1

2

Biosym/Molecular Simulation, 9685 Scranton Road, San Diego, CA 92121-3752 Lubrizol Corporation, 29400 Lakeland Boulevard, Wickliffe, OH 44092-2298

1

2

Standard methods for performing Kohn-Sham calculations involve self-consistency cycling a n d require either the solution of Poisson's equation for a density given as a quadratic sum over the orbital basis, or the use of a separate density basis and a fitting procedure. A n alternate approach using a functional different from but closely related to that of K o h n a n d Sham avoids both self-consistency cycling and density fitting and allows accurate results w i t h comparatively less effort. The main elements of this approach a n d some applications that illustrate its potential i n chemistry are discussed.

Over the past few years, the usefulness of the Kohn-Sham density f u n c t i o n a l / l o c a l density approximation ( D F T - L D A ) approach [1] i n chemistry has been w i d e l y recognized. Recently, the method was v a l i d a t e d i n systematic calculations for a set of transition-metal compounds i n v o l v i n g first- a n d second-row transition metals [2,3]. Determining structural information for transition-metal compounds has presented a significant challenge for theoretical methods. Hartree-Fock theory [4], so successful when applied to organic systems, fails to predict accurately the structures of organometallics and more accurate methods that i n c l u d e electron c o r r e l a t i o n e x p l i c i t l y are c o m p u t a t i o n a l l y demanding. Recently, a successful parametrization of a semi-empirical method has been accomplished for several transition metals [5]. The semi– empirical approach is very fast, but it is limited to systems for w h i c h parameters are available. There are many possible hybridization schemes and oxidation states of transition metals bonded to various ligands and this makes determination of semi-empirical or force-field parameters difficult and not easily transferable. It was found [2,3] that the local density approximation provides accurate structures i n v o l v i n g covalent metal– ligand bonds. In the case of dative bonds such as occur i n metal carbonyls, 0097-6156/96/0629-0388$15.00/0 © 1996 American Chemical Society In Chemical Applications of Density-Functional Theory; Laird, B., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 1996.

Downloaded by STANFORD UNIV GREEN LIBR on October 8, 2012 | http://pubs.acs.org Publication Date: May 5, 1996 | doi: 10.1021/bk-1996-0629.ch026

26. LI ET AL.

A Fast Density-Functional Method for Chemistry

389

the L D A gives somewhat shortened metal-ligand bond lengths. This is due to the tendency of the L D A to overestimate attractive nonbonded interactions and is corrected when gradient corrections are included. The availability of a first principles method that predicts accurate structures at a lesser computational cost than correlated ab initio methods is therefore rather important. Nevertheless, the standard implementation of the D F T method, via self-consistent solution of the K o h n Sham equations, remains a major computational problem a n d alternative methods of implementation that can reduce this cost are highly desirable. The main purpose of this paper is to outline the elements of such a method and to validate it via explicit calculation, particularly for organometallic molecules. The crucial advantage of the K o h n - S h a m D F T approach over conventional quantum chemical methods is the replacement of the manyelectron wavefunction by the density as the quantity that is varied, and the resulting reduction to a "one-particle" problem. In carrying out K o h n Sham calculations some computational difficulties arise that are due ultimately to the definition of the Kohn-Sham density functional, Eks/ in terms of V-representable densities. Each trial density must be generated by a Schrôdinger equation and cannot be chosen freely. In practice, this means that E k s cannot be minimized directly by optimizing within a restricted density basis. Furthermore, if atomic orbitals are used to solve the one-particle problem, the evaluation of Eks requires, in principle, fourcenter integrals. T o circumvent this, density fit (or density-projection) procedures are commonly used i n conjunction with self-consistency cycling (see, eg., Ref. [6]). A n alternate procedure is to exploit the properties of a different density functional, E , that has come to be known as the Harris Functional [7-9]. This possesses the same stationary points as the K o h n - S h a m functional but is defined on function space. The crucial advantage of Ε as distinct from Eks is that it can be evaluated using a density basis chosen for calculational convenience, for example, a sum over spherically symmetric site densities. The disadvantage is that, unlike E k s / does not obey a m i n i m u m principle. Finnis first pointed out that Ε displays a maximum for "reasonable" density variations [10]. Robinson and Farid showed that, within the local density approximation ( L D A ) , Ε displays a saddle point, but with positive curvature occurring only for density fluctuations composed solely of spatially rapidly varying components [11]. Since it is quite straightforward to eliminate this possibility, a maximum property of Ε can be guaranteed. Maximizing Ε is then the equivalent in the new approach to m i n i m i z i n g Eks (commonly achieved via "self-consistency cycling"). Other authors have exploited the favourable properties of E , notably Sankey and Niklewski [12], who restricted the trial density to a sum of overlapped neutral atom densities and made a series of approximations that resulted in an extremely fast tight-binding-like scheme. A series of applications, notably to large carbon systems, has shown that this approach E

In Chemical Applications of Density-Functional Theory; Laird, B., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 1996.


390

CHEMICAL APPLICATIONS OF DENSITY-FUNCTIONAL THEORY

can work rather well [13]. The scheme we describe here is more conservative and focuses in the first instance on universality (i.e. full coverage of the periodic table) rather than on efficiency for a limited class of elements. The scheme is a further development of a method suggested previously [14] and described in detail by Lin and Harris [15]. The most important innovation is provision for charge transfer which is crucial in general applications. Recently, the Sankey-Niklewski scheme has also been modified to allow for charge transfer [16]. The organization of this paper is as follows: In Section 2, we give a brief description of the scheme and detail those calculational features that differ from the proposal of Lin and Harris [15]. These refer in the main to the orbital and density bases used, the manner of performing integrals and the provision for charge transfer. In Section 3, we then quote results on the structure of a series of molecules, which illustrate that the present scheme gives an excellent description of structural properties over a wide class of bonding environments. Theoretical Discussion The theoretical elements underlying the basic scheme are described in detail in Ref. [15] and we restrict ourselves here to a brief summary. The energy functional, for Ν electrons in external (usually nuclear) field Vext(x) is:

Ε[η]=Σα ε -ί^){^ ( )-εχο( )+μ ^)} η

φ

η

χ

χ

χ€

Eq.l

where ε are the eigenvalues of the one-particle Schrodinger equation η

Eq.2 with potential construction

ν ί£(^)=ν χ ( ) ^) μ ^) 6

6

ί

5ί +φ

+

χο

Eq.3

Here, Φ is the Coulomb potential corresponding to n(x) and Eq.4 We assume the local density approximation (LDA) for the exchangecorrelation energy density


26. LI ET AL.

Λ Fast Density-Functional Method for Chemistry

ε ο(*)=ε«Η*))

391

Eq-s

Χ

with £x (n) the energy of a homogeneous electron gas having density n. A notable feature of the energy functional, Eq. 1, is that its evaluation does not involve explicitly the "out" density from the Schrôdinger equation, C

«(*) = Σ„Λ„|ψ„(ί)|

Eq.6


2

where the a are appropriate occupation numbers. The "in-density" appearing in E[n] in Eq. 1 is independent of the orbitals that are generated via the solution of Eq. 2 and it is this freedom that gives Eq. 1 its flexibility with regard to approximation. n

The evaluation of Eq. 1 and its derivatives is achieved by first choosing a convenient density basis. This defines the potential appearing in Eq. 2, which is then solved using an orbital basis. As an appropriate trial density we use the form

η(χ) = ΣΣυΖ ν(|χ-χ|) ι

Eq.7

where ι^^(|χ —X||) is an atom centered, spherically-symmetric site density and 2X is a parameter that may be fixed or can vary. In Eq. 7, the sum over " i " refers to a sum over sites, and the sum over "nu" allows for two or more density basis functions per site. Since only spherically symmetric site densities are used, the solution of Poisson's equation for any density of r

the form in Eq. 7 is trivial. If the fixed, the energy and forces are determined by a single evaluation of Eq. 1 and its derivatives with respect to the nuclear positions (formulae for which are given in Ref. [15]). The time limiting step in the calculation is the solution of the one-particle Schrôdinger equation as in Eq. 2, but only one solution is needed for a given set of nuclear locations so the speed up as compared with a traditional DFT code (other things being equal) is roughly the number of iterations the latter requires to become "self-consistent." This can be as small as 3 and as large as 50 depending on the application. The individual density basis functions in Eq. 7 are numerical and are generated using an atom program. For light atoms, two functions per site, one representing the core and one the valence, are often sufficient. This is because s- and p-orbitals in the first row have similar extents. Supplementary functions calculated using ionic configurations to improve accuracy can be added. For heavier atoms, in particular for transition elements it is important to allow for weighting of individual a r e


392

CHEMICAL APPLICATIONS OF DENSITY-FUNCTIONAL THEORY

shells. This is because transfer between shells that accompanies chemical bonding can cause quite drastic changes i n the density. Thus, i n the third row separate functions are used to represent the s- and d - shells. In general, it w i l l be necessary to allow for some variation i n the


multiplicative parameters, Z?, that influence the charge distribution. The variation may need to describe a charge transfer between sites or an effective change of intra-site electron configuration (eg., s => d transfer). This can be done conveniently using the maximum principle obeyed by Ε w i t h respect to density fluctuations. The "forces" acting on the charge variables 2^ are given by E q . 9 of Ref. 15 and all that is necessary is to maximize Ε subject to the constraint on the overall norm

ΣΣΖ^Ν i

Eq.8

υ

The m a x i m u m principle can be guaranteed since a l l density fluctuations corresponding to a given have the same spatial form. It is therefore easy to ensure that the electrostatic contribution to the curvature outweighs the exchange-correlation component, w h i c h is the condition for a m a x i m u m . (Another w a y of saying this is that a l l density fluctuations corresponding to a given site density have the same decomposition i n Fourier space. It is necessary only to ensure that each function of the density basis has sufficient weight i n low-k components.) In practice, this w i l l be true for any atom-derived density of the type considered. A potential problem w o u l d arise here only if fine tuning of the density were to be attempted by, eg., including density basis functions w i t h strong angular variations. Since the Harris functional is quite flat near to its saddle point, it is likely that acceptable accuracy can invariably be achieved without such fine timing. The one-particle Schrôdinger e q u a t i o n , E q . 2, is s o l v e d approximately by expanding the one-electron wavefunctons, Ψ„, i n a basis of atomic orbitals (AO's), χ., i n the usual way,

i

and solving the resulting secular problem

[H,-£A]c; = 0

Eq.10

Here Hij and Sij are Hamiltonian and overlap matrices, respectively. In a departure from the method i n Ref. [15], the orbital basis is numerical rather than analytic, comprising A O s generated using an atom program with finite range boundary conditions, so that the orbitals vanish


26.

L I ET AL.

A Fast Density-Functional

Method for Chemistry

393


continuously at a cutoff radius, r?. These radii can be chosen to combine accuracy w i t h efficiency and numerical stability. The r? determine the range beyond which overlap of AO's centered on different sites vanishes and, therefore, control the number of integrals that contribute to the H a m i l t o n i a n and overlap matrices. As Sankey and N i k l e w s k i , i n particular, have noted [12], finite range orbitals have a critical advantage over conventional AOs when calculations are carried out for extended systems. A l l AOs are generated using an atom program w i t h neutral atomic, and ionic configurations. The latter include "polarization functions" or "atom-unoccupied orbitals" that are important whenever the corresponding levels are low-lying. (Occasionally, truncated Slater orbitals are preferred for polarization functions.) The overall time-limiting feature of the scheme is the evaluation of the three-center integrals that contribute to the matrix elements in Eq. 10. Each of these is performed using the efficient weight-function method of Delley [6]. Energy derivatives are calculated as detailed in Ref, [15] (cf. Eq. 10. This is equivalent to the corresponding expression, also Eq. 10, in Ref. [17]), w i t h additional contributions that arise because of the dependence of the numerical grids on the nuclear positions. The energy gradients and energy are then consistent to high accuracy even when the numerical meshes are sparse. I n contrast to Sankey and Niklewski, we have not invoked approximating procedures for the three-center integrals. This is because, initially, our goal is to retain universality at the same level as conventional D F T - L D A schemes (some of the S a n k e y / N i k l e w s k i approximations breakdown for transition elements, for example). We believe, nevertheless, that there are many possibilities for approximation that can reduce the CPU requirement of the code very significantly. Geometry optimization can be carried out either directly, using a BFGS scheme [18], w i t h or without optimization of the charge parameters, or by dynamical methods. In the former case, it should be possible in principle to optimize the nuclear positions and the charge density variables simultaneously. However, we have not yet found a procedure for this that was sufficiently stable. Currently, the charge variables are optimized separately for each set of nuclear positions. This is a stable, but not efficient procedure. Dynamical methods (simulated annealing) employ M D algorithms and can, i n principle, propagate the charge variables along w i t h the nuclear positions. As detailed elsewhere [10], the density parameters can be assigned a (negative) mass such that they are continuously and adiabatically driven to the maximum of E[n]. This method, which is a real-space version of the propagation method of Car and Parrinello for plane-wave coefficients [19], achieves continuous updating of the trial density so this is o p t i m i z e d at all times at very little additional computational cost. Tliis is certainly the method of choice for simulated annealing and molecular dynamics if it is relatively straightforward to


394

C H E M I C A L APPLICATIONS O F DENSITY-FUNCTIONAL T H E O R Y

find appropriate masses for the density variables. This is a matter of building up experience.


Results A n important aspect of the calculational scheme described above is the limited nature of the density basis. Since the actual density of a system displaying covalent bonding cannot be represented i n the form given i n Eq. 7, it might seem that the use of this ansatz may involve unacceptable error whatever basis functions are used. In this section, we w i l l quote results for a series of molecules w h i c h demonstrate, so far as structural information is concerned, that the errors made are actually quite small. Table 1 shows structural parameters for a number of dimers and small organic molecules. These all-electron calculations were performed using the local density approximation w i t h the V W N parameterization, and the "standard basis" o p t i o n for density a n d orbitals of the F a s t S t r u c t / S i m A n n module (see below). The final column i n Table 1 shows the percentage error as compared w i t h experiment. For the most part, the errors are of the order of typical L D A errors. That is, the use of the functional i n E q . 1 i n conjunction with a trial density constructed from spherically symmetric site densities does not i n v o l v e m u c h loss of significance i n the results. N e x t we consider some transition-metal compounds i n v o l v i n g first- a n d second-row transition metals. It has p r e v i o u s l y been demonstrated [2,3] that D F T - L D A calculations (using D M o l ) give rather accurate structures for these compounds. In Table 2, we compared these earlier data w i t h corresponding results obtained u s i n g the current approach (as implemented i n the FastStruct/SimAnn module) and w i t h experimental data. The comparison does not yield a definitive measure of the difference between the simplified calculation and the K o h n - S h a m limit because slightly different exchange correlation functionals were used (von Barth-Hedin for D M o l and V W N for F S / S A ) , and the orbital basis functions were also somewhat different (the "standard" option of the respective modules). Nevertheless, the close agreement between the two sets of calculations allows the conclusion that the relative simplicity of the trial density i n E q . 7 does not obviate rather accurate calculation of structural parameters. F i n a l l y , we consider applications to larger organometallic molecules. The structure of bis(N-methyl-5-nitrosalicylideneaminato) nickel(II) was determined recently by Kamenar et al [20]. This is a square planar N i compound with 41 atoms, as displayed i n Figure 1. The F S / S A optimized structure of this molecule reproduces all the qualitative features of the experimental structure. To give a more quantitative comparison, selected bond distances and bond angles are presented i n Table 3. A s previously, the F S / S A results are close to the D M o l results, with a m a x i m u m difference of about 2% (in the C - N bond of the nitro


26.

L I ET AL.

A Fast Density-Functional Method for Chemistry

Table 1. Structures of dimers and small molecules System

Percent

Exp. [4]

r(LiLi) r(BeBe) r(BB) r(CC) r(NN) r(00) r(FF) r(NaNa) r(AlAl) r(SiSi) r(PP) r(SS) r(ClCl)

2.638 2.535 1.630 1.232 1.117 1.249 1.409 2.957 2.478 2.301 1.910 1.982 2.089

2.673 2.49 1.60 1.24 1.098 1.208 1.412 3.078 2.466 2.245 1.893 1.889 1.988

1.8% 1.9% -0.7% 1.8% 3.4% -0.2% -4.1% 0.5% 2.5% 0.9% 4.9% 5.1%

Cu LiF NaCl CuF CuCl TiCl

r(CuCu) r(LiF) r(NaCl) r(CuF) r(CuCl) r(TiCl)

2.178 1.594 2.346 1.760 2.059 2.148

2.220 1.564 2.361 1.745 2.051 2.170

-1.9% 1.9% -0.6% 0.8% 0.4% -1.0%

CH CH

2

CH

3

r(CH) 1.156 r(CH) 1.059

A Fast Density-Functional Method for Chemistry - American

Recommend Documents