Chemical Physics Letters 635 (2015) 201–204
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Improving DIIS convergence for metallic systems using Gaussian basis set David Sulzer a,b , Satoru Iuchi b , Koji Yasuda b,c,∗ a b c
Institute for Molecular Science, 38 Nishigo-Naka, Myodaiji, Okazaki 444-8585, Japan Graduate School of Information Science, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan EcoTopia Science Institute, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan
a r t i c l e
i n f o
Article history: Received 16 April 2015 In final form 25 June 2015 Available online 6 July 2015
a b s t r a c t Self-consistent calculations on metallic systems in Gaussian basis sets are likely to encounter convergence problems because of the long-wavelength charge sloshing. Using a simple formula for the charge response of the Fock matrix, a correction to DIIS have been derived to suppress sloshing effects. The proposed method is similar in spirit to the Kerker preconditioner in band-structure calculations. The new method is applied to systems like Ru4 (CO), Pt13 , Pt55 and (TiO2 )24 and the convergence is compared to previously available DIIS methods. For systems with small HOMO–LUMO gap, this new method manages to converge while previous DIIS methods do not. © 2015 Elsevier B.V. All rights reserved.
1. Introduction Solving the self-consistent field (SCF) problem in the Hartree–Fock or in the density functional theory is a necessary step for quantum chemistry calculations. Among various methods, the combination of the energy DIIS (EDIIS) [1] and the commutator DIIS (CDIIS) [2,3] is regarded as the best choice for Gaussian basis sets [4], because it works satisfactory for small molecules and insulators. However, for large metal clusters with very small HOMO–LUMO gap or metallic systems, it shows the slow convergence and sometimes fails, because the charge oscillation cannot be sufficiently damped. In this case, a damping technique with a small mixing parameter to smooth the oscillation is often used, but this considerably reduces the convergence speed. A quadratic convergent SCF (QCSCF) method [5] is robust and considered as more reliable for these systems, but this method is more expensive in computations. To overcome these problems, we propose a simple yet efficient method to deal with the convergence difficulty in metallic systems for Gaussian basis set. This new method is inspired by approaches often used in periodic systems with plane wave basis set; the combination of DIIS that minimizes the density residual, the special preconditioner
∗ Corresponding author at: Graduate School of Information Science, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan. E-mail addresses:
[email protected] (D. Sulzer),
[email protected] (S. Iuchi),
[email protected] (K. Yasuda). http://dx.doi.org/10.1016/j.cplett.2015.06.063 0009-2614/© 2015 Elsevier B.V. All rights reserved.
to damp the long wavelength charge oscillations [6] and the smearing of the sharp Fermi level. These techniques are known to greatly enhance the SCF convergence in metallic systems [7]. However, this preconditioner is closely tied to the plane wave basis set and cannot be used for Gaussian basis set. Thus, we propose an adaptation of these damping techniques to Gaussian basis set. We found that this new method was not superior to EDIIS + CDIIS for small molecules as expected, but greatly improved convergence for metal clusters, in which EDIIS + CDIIS fails to converge. This new method should allow the ab initio study of these systems with a reasonable computational cost.
2. Method In the usual Gaussian basis calculation [5] the error of the self-consistency R i at the ith iteration can be measured by the commutator of the density matrix P i and the Fock matrix F i . Assuming the linearity of R with respect to P, the DIIS method guesses the solution as a linear combination of previous Fock matri-
i
˛ F . Then, the matrix F is diagonalized, and the ces as F = j=1 j j new density matrix is generated by selecting the occupied orbitals: P i+1 = ( F), where is the distribution function. We propose to calculate the correction ı P, to the DIIS extrapolated density matrix: P = holds at the first-order as P + ıP = (F + ıF).
i
˛P, j=1 j j
so that the self-consistency
(1)
202
D. Sulzer et al. / Chemical Physics Letters 635 (2015) 201–204
The correction ı P induces the Fock matrix response ı F = F( P + ı P) − F( P). This ı F can be approximated using the orthonormal basis as ıF = W1 ıP + W2 1,
(2)
where Wi are constants determined below, and 1 is the unit matrix. We found that the correction based on this drastic approximation greatly improves the convergence for metallic systems. The off-diagonal elements of ı P are determined with the commutator equation, [ F + ı F, ( F) + P − ( F) + ı P] = 0. Since F was already diagonalized, it was solved in molecular orbital (MO) basis. Then, keeping only the nonzero lowest-order terms, we obtain [(i − j ) − W1 (ni − nj )]ıPij + (i − j )[Pij − (F)ij ] = 0,
(3)
where i and ni are the eigenvalues of F and ( F), respectively. The off-diagonal elements of ı P are then given as ıPij = Wij [(F)ij − Pij ], with
Wij = 1 − W1
ni − nj
(4)
.
(5)
To reduce singularity, the smooth Fermi–Dirac distribution function was used for at a temperature T as (F) = [1 + eˇ(F−1) ]
−1
,
(6)
T = ˇ−1
with kB being the Boltzmann constant, and the where kB chemical potential is determined to satisfy the electron number conservation: Tr( F) = N. Hereafter, kB is dropped in equations for simplicity. When |i − j | in Eq. (5) is small (less than 10−5 a.u.), the weight factor Wij is replaced with its limit [1 + ˇni (1 − nj )]−1 . The diagonal elements of ı P are determined from the diagonal part of Eq. (1) as ıPkk = Wkk [nk − Pkk − ˇW2 nk (1 − nk )],
(7)
with Wkk = [1 + ˇW1 nk (1 − nk )]
−1
.
(8)
The constant W2 , which is the correction to the chemical potential, is determined for ı P in order to satisfy the normalization: k ıPkk = 0. When there is no fractional occupation, ıPkk becomes nk − Pkk . The correction Eq. (4) is ( F) − P multiplied by the orbitaldependent damping factor Wij and satisfies 0 < Wij ≤ 1 when W1 ≥ 0. The damping factor Wij is considered as a small quantity for partially occupied orbitals, while other elements in the occupiedoccupied or the empty-empty block are close to 1. When nk = 1/2 the diagonal element Wkk reaches the minimum [1 + ˇW1 /4]−1 . If terms containing W1 is zero, then Wij = 1 and ı P = ( F) − P. The parameter W1 expresses the Fock response caused by the change of the density matrix. It can be also understood as the averaged second derivative of the energy with respect to the density matrix and is used to calculate the change of the density matrix in the current step. In the actual implementation, we applied a simple formula for W1 as W1 = 4(1/W3 − 1)T,
(9)
where W3 is an adjustable parameter. This formulation ensures that the weight of Wkk is at least W3 . W3 = 0.1 was chosen because it works well for examples presented below. For the temperature, we used T = max({1/M
ij
2.1. Comparison to Kerker preconditioner The Kerker preconditioner was developed to enhance the convergence of metallic systems [6,7] in the field of solid state plane wave calculations. Usually, a SCF cycle consists of three steps. First, the input density i determines the Fock matrix Fi . Then, by diagonalizing it and distributing electrons in the orbitals, the new density new is calculated. Self-consistency requires that the change of the density R() = new − i is zero. Otherwise, the DIIS technique extrapolates the optimum density from densities of previous itera-
i
−1
i − j
where M is the number of atomic orbitals (AOs). The T is an intensive quantity and decreases as the SCF error becomes small until reaching the predefined value Tgoal . Finally, the density matrix for the next iteration is given as P + Aı P, where A = 0.8 is the second adjustable parameter. If T is close to zero compared to the HOMO–LUMO gap, the next density matrix becomes A( F) + (1 − A) P. Thus, the present method reduces to the standard CDIIS for A = 1.
[F, P]2ij }
1/2
, Tgoal ),
(10)
˛ by minimization of the 2-norm of R(opt ). tions as opt = j=1 j j Because of the linearity assumption, R(opt ) approximates the selfconsistency error when opt is the input density. The best correction to the input density (ı) which leads to the self-consistency is searched as follows. The density correction ı induces the Fock matrix change, ıF, in which the major contribution comes from the electrostatic potential. It is given by solving the Poisson equation for a Coulomb potential in momentum space, ıF(k) =
4 ı(k), k2
(11)
where k is the wave vector. This ıF is considered as a perturbation that changes the electron density slightly. By using the Thomas–Fermi approximation, slowly varying potential ıF induces the density response as, ı = −
k02 4
ıF,
(12)
where k0 is the local Fermi vector. The self-consistency requires opt + ı = new + ı ,
(13)
where new is the density after diagonalization when we use opt as the input density. Under linearity assumption we can estimate it as new − opt ≈ R(opt ).
(14)
Hence the self-consistency can be achieved for ı =
k2 k2
+ k02
R(opt ).
(15)
It is interesting to notice that a momentum-dependent factor k2 /(k2 + k02 ) damps the long-wavelength charge oscillation. This factor is called the Kerker preconditioner. The final formula is deduced: i+1 = opt + A
k2 k2
+ k02
R(opt ),
(16)
in which A is an adjustable parameter. In the actual implementation, the local Fermi vector k0 is also regarded as an adjustable parameter. The proposed new method is comparable with the Kerker preconditioner, because both of them correct the DIIS extrapolated density or the density matrix. This correction induces ıF, which then changes the density or the density matrix. In plane wave basis, the induced ıF is calculated accurately by solving the Poisson
D. Sulzer et al. / Chemical Physics Letters 635 (2015) 201–204
203
(a) Ru4(CO)
(b) (TiO2)24
(c) Pt13
(d) Pt55
Figure 1. Comparison of SCF energy convergence between EDIIS + CDIIS and the present method for various molecules. The initial guess and the final energy are also shown.
equation, while the density response ı should be approximated. The self-consistency requirement leads to a damping factor of the density in momentum space. On the other hand, in Gaussian basis, the density matrix response to ıF can be solved exactly in MO basis while ıF is approximated with Eq. (2), which results in a damping factor of the density matrix in the MO basis. The validity of the drastic approximation used will be examined in the next section.
Refs. [17,18]. These are distorted Ih symmetry structures corresponding to local minima, and coordinates of Pt13 and Pt55 are given in supplementary materials. The density fitting technique [19] with the SVP density basis [20,21] was used to accelerate the Coulomb matrix computation of Pt55 and (TiO2 )24 .
4. Results and discussion 3. Computational details The performance of the new method has been examined by comparing it with that of the standard EDIIS + CDIIS method. Three types of molecules and/or clusters were investigated; i.e., small molecules, a semi-conductor and metal clusters. All the calculations were performed using the modified version of Gaussian 09 (G09) program [8], and the default options of G09 were used unless otherwise indicated. All the calculations retained 20 Fock and density matrices as the DIIS subspace. Point group symmetry were turned off in all the calculations. A Tgoal = 298 K was chosen in this study unless otherwise stated. Considering small molecules, Cd2+ imidazole complex ([Cd(Im)]2+ ), tetrahedral UF4 , one-bond stretched SiH4 , and Ru4 (CO) were calculated at the B3LYP/3-21G, UB3LYP/LANL2DZ, LDA/6-31G*, and UB3LYP/LANL2DZ, respectively. These choices were employed previously in Refs. [4,9,10] to discuss the performance of several DIIS methods including EDIIS + CDIIS and thus suitable as a test case for the newly proposed method in this study. For the semi-conductor, a (TiO2 )24 cluster defining one layer of the anatase crystals [11] cut along the 1 0 1 direction was used. The geometry was generated from ideal crystallographic structure and the PBE functional [12] was used with a 6-31G basis set [13,14]. For metal clusters, icosahedral (Ih) Pt13 and Pt55 molecules were calculated using the PBE functional [12] with the double-zeta Stuttgart/Dresden ECP [15,16]. Geometries were extracted from
First, in order to check whether the convergence problem can be solved by various choices of methods and relevant parameters, various combinations of DIIS and smearing techniques as implemented in the G09 program [22] were examined. The DIIS subspace was increased up to 50, and the damping factor was reduced to 0.1. These choices showed almost the same convergence properties, and the default EDIIS + CDIIS seems to be the most relevant. Hence, a comparison between EDIIS + CDIIS with standard parameter and the present methods is reported in terms of the SCF convergence. Figure 1 displays the most interesting cases: Ru4 (CO), (TiO2 )24 , Pt13 , and Pt55 . Results for other molecules are found in Figure S1 in the supplementary materials. For a small molecule with large HOMO–LUMO gap, little improvement is expected over CDIIS, and this was confirmed by the calculations for [Cd(Im)]2+ , UF4 , and SiH4 . The EDIIS + CDIIS method was better suited than the present method for these systems, because EDIIS + CDIIS shows the similar or the faster convergence (Figure S1 in the supplementary materials). Interestingly for UF4 the present method resulted in slightly lower energy (6 × 10−4 a.u.) at the expense of longer iterations. Because HOMO–LUMO gap is much larger than the room temperature, no visible fractional occupation was observed in these molecules. Also, if the temperature is small enough compared to the HOMO–LUMO gap, the present method reduces to the usual CDIIS with damping. For Ru4 (CO) standard methods failed to converge, whereas the new proposed method achieved convergence (Figure 1a). Both calculations started
204
D. Sulzer et al. / Chemical Physics Letters 635 (2015) 201–204
from the initial guess generated with G09 options, guess = (harris, mixed). However, it was reported in Ref. [4] that EDIIS + CDIIS with 20 DIIS matrices rapidly led to the convergence from a brokensymmetry triplet solution as the initial guess. Our solution seems to be a local minimum, because the energy was higher than the value reported [4]. Thus, both methods should not be applied as a black-box for this system. By contrast, the present method successfully overcame convergence problems of EDIIS + CDIIS for metallic and semi-conductor clusters. In the case of the semi-conductor cluster, (TiO2 )24 , the convergence could not be achieved even after 300 iterations by EDIIS + CDIIS, whereas the new method converges after 102 iterations with a standard temperature of 298 K (Figure 1b). The orbital energy is characterized by a small HOMO–LUMO gap (0.01 eV) hindering the convergence for the EDIIS + CDIIS method. Both methods show a similar behavior for the first 50 iterations. Then, the proposed method allows to obtain the converged density and the energy after about 100 iterations, whereas the energy of EDIIS + CDIIS oscillates and is not improved. In the metallic test case, the convergence stagnated with the usual EDIIS + CDIIS method as seen for the Pt13 cluster (Figure 1c) and Pt55 cluster (Figure 1d), and the calculation was terminated after 300 iterations without observing convergence. Because of the narrow HOMO–LUMO gap, the level-shift method is difficult to apply without precise knowledge of the electronic structure. Yet, the new method quickly leads to the convergence for the Pt13 cluster as seen in Figure 1c. In the case of the Pt55 cluster, Tgoal = 580 K (0.05 eV) was used to improve the convergence. The same smearing parameter is often used in the plane-wave calculations. The density matrix slowly approached the solution (Figure 1d) and the convergence is observed after 266 iterations for the new method. We found that the convergence behavior was influenced by the selection of the DIIS vectors. We kept 10 latest vectors and replaced the least stable one in the rest with a new one. To see the temperature dependence on the SCF convergence, the calculations with Tgoal = 298 K and 1500 K were also performed. It is commonly accepted that the smearing with higher temperature leads to faster convergence for difficult cases in plane-wave calculations. We also observed the same trend: after 100 iterations an energy difference of 10−8 a.u. was observed for 1500 K while 10−6 a.u. for 580 and 298 K. Finally, the standard CDIIS with finite temperature smearing have been tried for Pt13 . This method led the convergence after 44 iterations if we keep the temperature at T = 298 K. Note that we had to modify the program to do so, since G09 always turns off the smearing near the convergence. Nevertheless, this simple method did not achieve convergence for Pt55 and (TiO2 )24 at T = 298 nor at T = 580 K. Thus, the method proposed in this letter is a good alternative for studying the properties of large cluster with a small HOMO–LUMO gap. In these cases, the present method was able to converge to the solution, while EDIIS + CDIIS or similar methods with a similar computational cost encounter convergence difficulties.
Fermi–Dirac distribution of electrons and a gradual cooling to the target temperature was used to smooth the convergence. The proposed method was applied to test molecules including Ru4 (CO), Pt13 , Pt55 , and (TiO2 )24 . The convergence properties were compared to those of DIIS (EDIIS + CDIIS) or similar methods and we found that the present method was not better than EDIIS + CDIIS for small molecules. Nevertheless, for metallic and semi-conductor clusters, the present method was superior to EDIIS + CDIIS, as it managed to converge to the solution while EDIIS + CDIIS or similar methods did not. Thus, the proposed method provides another route to overcome convergence difficulties in metallic systems. With a more elaborated formula of the Fock matrix response the method is expected to attain even better convergence property. Acknowledgements This work was partially supported by the Computational Materials Science Initiative (CMSI), Japan. Some computation were performed using Research Center for Computational Science, Okazaki, Japan. Appendix A. Supplementary Data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.cplett.2015.06. 063 References [1] [2] [3] [4] [5] [6] [7] [8]
[9] [10] [11] [12] [13] [14] [15] [16]
5. Conclusion The origin of the slow SCF convergence in a metallic system is the huge charge response and the resultant long wavelength charge sloshing. Inspired by the Kerker preconditioner used in the bandstructure calculations, a simple model for the charge response of the Fock matrix was proposed, and the correction to the commutator DIIS (CDIIS) of Pulay was derived. This new correction can be understood as an orbital dependent damping for CDIIS. Then, a
[17] [18] [19] [20] [21] [22]
K.N. Kudin, G.E. Scuseria, E. Cancès, J. Chem. Phys. 116 (19) (2002) 8255. P. Pulay, Chem. Phys. Lett. 73 (2) (1980) 393. P. Pulay, J. Comput. Chem. 3 (4) (1982) 556. A.J. Garza, G.E. Scuseria, J. Chem. Phys. 137 (5) (2012) 054110. T. Helgaker, P. Jorgensen, J. Olsen, Molecular Electronic-Structure Theory, John Wiley & Sons, 2009. G.P. Kerker, Phys. Rev. B 23 (6) (1981) 3082. G. Kresse, J. Furthmüller, Comput. Mater. Sci. 6 (1) (1996) 15. M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G.A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H.P. Hratchian, A.F. Izmaylov, J. Bloino, G. Zheng, J.L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J.A. Montgomery Jr., J.E. Peralta, F. Ogliaro, M. Bearpark, J.J. Heyd, E. Brothers, K.N. Kudin, V.N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J.C. Burant, S.S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J.M. Millam, M. Klene, J.E. Knox, J.B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R.E. Stratmann, O. Yazyev, A.J. Austin, R. Cammi, C. Pomelli, J.W. Ochterski, R.L. Martin, K. Morokuma, V.G. Zakrzewski, G.A. Voth, P. Salvador, J.J. Dannenberg, S. Dapprich, A.D. Daniels, Ö. Farkas, J.B. Foresman, J.V. Ortiz, J. Cioslowski, D.J. Fox, Gaussian 09, Revision C.01, Gaussian Inc., Wallingford, CT, 2009. X. Hu, W. Yang, J. Chem. Phys. 132 (5) (2010) 054109. Y.A. Wang, C.Y. Yam, Y.K. Chen, G. Chen, J. Chem. Phys. 134 (24) (2011) 241103. R. Wyckoff, Miscellaneous Inorganic Compounds, Silicates, and Basic Structural Information, vol. 4, Wiley Interscience, New York, 1968. J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. W.J. Hehre, R. Ditchfield, J.A. Pople, J. Chem. Phys. 56 (5) (1972) 2257. V.A. Rassolov, J.A. Pople, M.A. Ratner, T.L. Windus, J. Chem. Phys. 109 (4) (1998) 1223. T.H. Dunning Jr., P.J. Hay, in: H.F. Schaefer III (Ed.), Modern Theoretical Chemistry, vol. 3, Plenum Press, New York, 1976. D. Andrae, U. Häußermann, M. Dolg, H. Stoll, H. Preuß, Theor. Chim. Acta 77 (2) (1990) 123. T. Futschek, J. Hafner, M. Marsman, J. Phys. Condens. Matter 18 (42) (2006) 9703. J.L.F. Da Silva, H.G. Kim, M.J. Piotrowski, M.J. Prieto, G. Tremiliosi-Filho, Phys. Rev. B 82 (20) (2010) 205424. B.I. Dunlap, N. Rösch, S. Trickey, Mol. Phys. 108 (21-23) (2010) 3167. K. Eichkorn, O. Treutler, H. Öhm, M. Häser, R. Ahlrichs, Chem. Phys. Lett. 240 (4) (1995) 283. K. Eichkorn, F. Weigend, O. Treutler, R. Ahlrichs, Theor. Chem. Acc. 97 (1–4) (1997) 119. A.D. Rabuck, G.E. Scuseria, J. Chem. Phys. 110 (2) (1999) 695.