Expeditious Stochastic Calculation of Random-Phase Approximation

Subscriber access provided by MCGILL UNIV

Letter

Expeditious Stochastic Calculation of Random-Phase Approximation Energies for Thousands of Electrons in 3 Dimensions Daniel Neuhauser, Eran Rabani, and Roi Baer J. Phys. Chem. Lett., Just Accepted Manuscript • DOI: 10.1021/jz3021606 • Publication Date (Web): 08 Mar 2013 Downloaded from http://pubs.acs.org on March 10, 2013

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry Letters is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 5

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

Expeditious Stochastic Calculation of Random-Phase Approximation Energies for Thousands of Electrons in 3 Dimensions Daniel Neuhauser†, Eran Rabani‡ and Roi Baer † Department

of Chemistry and Biochemistry, University of California, Los Angeles CA 90095, USA. of Chemistry, The Sackler Faculty of Exact Sciences, Tel Aviv University, Tel Aviv 69978, Israel. Fritz Haber Center for Molecular Dynamics, Institute of Chemistry, Hebrew University, Jerusalem 91904 Israel.

‡ School 

ABSTRACT: A fast method is developed for calculating the Random-Phase-Approximation (RPA) correlation energy for density functional theory. The correlation energy is given by a trace over a projected RPA response matrix and the trace is taken by a stochastic approach using random perturbation vectors. For a fixed statistical error in the total energy per electron the method scales at most, quadratically with the system size, however, in practice, due to self-averaging, requires less statistical sampling as the system grows and the performance is close to linear scaling. We demonstrate the method by calculating the RPA correlation energy for cadmium selenide and silicon nanocrystals with over 1500 electrons. We find that the RPA correlation energies per electron are largely independent of the nanocrystal size. In addition, we show that a correlated sampling technique enables calculation of the energy difference between two slightly distorted configurations with scaling and a statistical error similar to that of the total energy per electron.

Local and semi-local correlation functionals of Kohn-Sham (KS) density functional theory (DFT) fail to describe longrange van der Waals interactions and other types of dynamical screening effects. ADDIN EN.CITE One route for overcoming these deficiencies is RPA theory1, 3-6 based on the KS-DFT adiabatic connection formalism7-9 in combination with the fluctuation dissipation theorem.10 In recent years this approach, especially when combined with exact exchange, was used successfully for treating various ailments of KSDFT in molecular and condensed matter systems5, 6, 11-18. The greatest hurdle facing widespread use of RPA is its exceedingly high computational cost. Several approaches have been developed 5, 6, 13, 19, 20 for reducing the naïve ( ̃ ) RPA scaling to ( ̃ ), ( ̃ is a measure of system size); however, this is still expensive. The problem is aggravated when planewaves or real-space grids are used, suffering from the huge number of unoccupied states and the strong reliance of the RPA energy on the unoccupied energies.21-23 In the present letter, we develop a stochastic sampling method for estimating the RPA correlation energy. Related sampling techniques have been recently developed by us for estimating the rate of multiexciton generation in nanocrystals (NCs),24 for a linear scaling calculation of the exchange energy,25 and for overcoming the computational bottleneck in Møller-Plesset second order perturbation theory (MP2).26 RPA, applied on top of a grid or plane waves calculation, starts from the KS or generalized KS Hamiltonian ̂ which can be applied to any wavefunction in linear scaling.27 For a closed-shell system of electrons on a grid/basis of size , the lowest KS orbitals ( ) of ̂ are occupied and are unoccupied.28 The RPA correlation energy can be written ∑ ̃ as5 (̃ ), where ̃ are eigenvalues of “ ̃ ” defined by: ̃( ) and

(

)( )

̃

( )

(1)

(2)

̃ ( )

∫

where

( )

( )

( )

are the Cou-

| - |

lomb integrals, is a coupling strength parameter ( is full-strength Coulomb coupling and is the noninteracting limit), and - is the difference of eigenvalues of ̂ . Usually, all ̃ are real; however ̃ , although having pure imaginary eigenvalues, is a real matrix operating on real vectors ( ). Note that ̃ are also the eigenvalues of the matrix ( -

) appearing in standard RPA treatments.5

-

An alternative formulation starts from the expression:5 ̃

[̃( ) (̃( )

)

]

(3)

̃ ( ). However, the calculation of ∑ ̃ where ̃( ) ̃( ) is still prohibitive for large systems because of the high cost of diagonalization of the ( - ) ( - ) ̃ matrix (in grid representations and easily reach and respectively). Our formulation is based on linear-response time-dependent Hartree approach.29, 30 is still given by Eq. ‎(3) but ̃( ) is replaced by the following trace: 29, 30 ( )

[

( ̂ ( ))]

(4)

( ) ( ), where ( ) is the Heaviside step funcHere tion, which we approximate as ( ) (- ); ̂ is a linear operator, revealed when linearizing the time-dependent KS equations (see refs. 14 and 30 for details): ̂( ) ( ) and

(

[

(̂ ]

) (̂

)

)

(5)

( ) are functions, originally describing the time-

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

dependent Hartree response of the -th KS orbital ( ), but are used here as stochastic perturbations as detailed below. [

]( )

∫

( )

is the Hartree perturbation poten-

| - |

tial depending linearly on the ’s via ( )

∑

( ) ( )

(6)

One can expand ̂ in the basis of the KS orbitals and obtain its matrix having positive imaginary eigenvalues and an equivalent negative set. ( - ) eigenvalues can be classified as “occupied-unoccupied” transitions , and eigenvalues are occupied-occupied transitions . Obviously, the dimensions of and ̃ differ, as the latter describes only occupied-unoccupied transitions. Nonetheless, within the occupied-unoccupied space the matrices and eigenvalues are identical:31 ̃

̃

Page 2 of 5

( )

(

( )

( ))

( )

(

( ))

̂(

( ) ( ))

(13)

Note that ̂ is a real operator (Eq. ‎(5)) so all calculations are done on real functions. ( ) is half the eigenvalue range of ̂ . The are numerical coefficients obtained as follows: First, prepare a series of length , ( ( )), - and then set ̃ (

)

) where { ̃ } are the discrete Fou-

(for

rier transform of { }. The series length is chosen large enough so that the sum in Eq. ‎(10) converges, i.e. | | is smaller than a prescribed tolerance. We rely on correlated sampling to reduce the statistical error in computing . is computed for 3 values of and – (with ) using the same random number seeds, and then the RPA energy is estimated by:

(7) [ ( )

( ) in Eq. ‎(4) is computed in two principal steps:

∑ (

) ( )]

(14)

1) The trace is replaced by an average (denoted by curly brackets) over random perturbation vectors ( ): ( )

( ̂ ( ))| ( )⟩}

)|

{⟨(

(8)

( ) is shorthand notation for the entire set of ’s.

’s and

2) The selection of the random vector and is done by generating two vectors composed of random complex ( ) phases at each grid point : ( ) where is the grid spacing and ( ) is a random number between 0 and . Then one sets: ( )

(

( )⟨ |

)

⟩

( )⟨ |

⟩

(9) ̂ 3) The action of the operator ( ) on ( ) is performed using an iterative modified Chebyshev polynomial expansion approach, so that: ( )

∑

(10)

where )| (

{⟨(

( ) ( ) )⟩}

(11)

are the modified Chebyshev residues, and ( culated iteratively: (

(

)

(

))

̂(

( ) ( ))

(

(

(

)

( ))

( ) ( )

Figure 1: The error of the stochastic estimate for the RPA energy per electron with respect to the exact value (black curve) and the variance (red curve) as function of stochastic iterations for SiH4. Table 1: Parameters for the CdSe nanocrystals NCs. Shown are the number of Cd ( ) and Se ( ) atoms, electrons ( ), NC diameter ( ), the numerical effort involved in operating with in a perturbation vector , where is the number of gridpoints and the occupied-unoccupied energy gap .

) are cal-

D (nm)

(

)

19 )

(12)

with We now demonstrate the performance of the stochastic method by applying it to calculate the RPA correlation energies of 2 ACS Paragon Plus Environment

Page 3 of 5

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


spherical cadmium-selenide (CdSe) and hydrogen passivated silicon NCs, where the Hamiltonian ̂ is constructed from a semiempirical pseudopotential model.32, 33 The occupied states of the NCs were obtained using the filter diagonalization technique34 with the implementation described in Refs. 32, 35

-

. We used for approximating the step function ( ) and (see discussion of Figure 4), and , slightly larger than half the maximal eigenvalue range for both NCs. Various features of the NCs are summarized in Table 1 and 2. Table 2: Same as Table 1, but for hydrogen passivated silicon NCs.

D (nm)

(

)

As a test, we compared the stochastic estimate and a full summation calculation of the RPA energy for on a point grid (see Figure 1). The deviation of the stochastic approach from the exact value is within the statistical variance for nearly all stochastic iterations. We found that for stochastic iterations the stochastic estimate deviates by ~ meV from the full summation value. As explained below, the number of required stochastic iterations decreases considerably (for a fixed error per electron) as the size of the system increases.

Figure 2 shows the RPA correlation energies for CdSe and silicon NCs up to electrons along with a comparison to MP2 energies obtained using the Neuhauser-Rabani-Baer (NRB) method.26 The RPA correlation energy depends weakly on the NC size in contrast to that of MP2. This is because the NC gaps decrease with system size and MP2 energies are sensitive to small gaps (diverging for metals). The RPA energy of silicon is somewhat above that of CdSe, and is within the LDA bulk limit21 range of . The insets of Figure 2 show the corresponding statistical errors normalized to 1000 stochastic iterations. The errors decrease when the number of electrons in the system increases. This shows that the algorithm profits from statistical selfaveraging. The statistical error of the CdSe NCs is approximately twice smaller than that for silicon despite having similar gaps for the same NC size. This suggests that the statistical errors are not trivially correlated with the gap. Figure 3 shows the total CPU time for calculations that yield a statistical error of meV per electron. The method scales, at most, quadratically with system size but in practice, due to self-averaging, requires considerably less statistical sampling as the system grows and the resulting performance is close to linear scaling. Furthermore, in comparison, the RPA CPU time for the same statistical error is an order of magnitude smaller than the CPU time required for the MP2 calculations. Regarding memory requirements, for the RPA scales quadratically with system size (17GB for the largest silicon NC) and linearly for MP2.

Figure 3: The CPU times for achieving a statistical error of meV per electron for the RPA and MP2 calculations of CdSe NCs.

Figure 2: RPA (blue) and MP2 (red) correlation energies per electron vs. the number of electrons , for silicon (top) and CdSe (bottom) NCs. Insets show the statistical errors in the RPA energies, normalized to 1000 stochastic iterations.

In some cases, the Chebyshev interpolation suffers from instabilities. This is shown in the Figure 4 for the NC where we plot the RPA energy estimate as function of the length of the Chebyshev expansion. The results (dashed line) clearly diverge as the Chebyshev expansion length grows. To ( ) alleviate this problem we project out from and ( ) , after each operation of ̂ , the occupied orbitals , where is a small system-dependent integer. Figure 4 shows the RPA correlation energy as a function of



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the Chebyshev interpolation order ( ). For the RPA energy with ( ) and without projection are nearly identical (to within statistical error). For higher interpolation orders than shown the RPA correlation energy without projection diverges while the projected one remains stable. For the Chebyshev polynomial visibly diverges already when the order is larger than 300 but becomes stable when projection is used with . Here too, the difference between the projected and non-projected correlation energies (before divergence) is smaller than the statistical fluctuations, which is large in the diverged result due to the small number of stochastic iterations done.

approach developed here bloats, by orders of magnitude, the size of systems that can be treated using RPA theory. DN was supported by the DOE MEEM center, award DESC0001342. RB was supported by the US-Israel Binational Foundation (BSF). RB and ER gratefully thank the Israel Science Foundation, grant numbers 1020/10 and 611/11, respectively. 1

2 3 4 5 6 7 8 9 10 11 12

Figure 4: The RPA correlation energy as function of the length of the Chebyshev interpolation polynomial for silicon NCs. Dashed lines are the results without projection while the solid lines are for for and for .

The present method is not only useful to calculate the RPA energy per electron. In fact, it can also be used to obtain, for example, the RPA corrections to chemical reaction energy profiles or to adsorption energies. In this respect, we have developed an approach, based on a combination of adiabatic evolution and correlated sampling to calculate energy differences between two structural configurations.36 Preliminary results for (enabling comparison to the full summation result) and , moving one atom by 0.01 , indicate that the statistical error of the energy gradient (i.e., the force due the RPA energy) is of the comparable to the statistical error of the total RPA energy per electron (much smaller than the statistical error in the total energy). For and the statistical error per stochastic iteration37 in the RPA energy is and , respectively, compared to the statistical error in the RPA energy per electron and , respectively. Summarizing, we have presented a new stochastic approach for calculating RPA energies for large electronic systems of exceptional size. The method scales formally as ( ̃ ) in terms of memory and CPU time but due to self-averaging has a near-linear scaling CPU time performance. We calculated the RPA correlation energy for CdSe and silicon NCs up to diameters of with over electrons. The stochastic

Page 4 of 5

13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

J. P. Perdew and K. Schmidt, in Density functional theory and its application to materials : Antwerp, Belgium, 8-10 June 2000, edited by V. E. Van Doren, C. Van Alsenoy and P. Geerlings (American Institute of Physics, Melville, N.Y., 2001). N. Marom, A. Tkatchenko, M. Rossi, V. V. Gobre, O. Hod, M. Scheffler, and L. Kronik, J. Chem. Theor. Comp 7, 3944 (2011). B. G. Janesko, T. M. Henderson, and G. E. Scuseria, J. Chem. Phys. 131, 034110 (2009). M. Fuchs and X. Gonze, Phys. Rev. B 65, 235109 (2002). H. Eshuis, J. E. Bates, and F. Furche, Theor. Chem. Acc. 131 (2012). X. G. Ren, P. Rinke, C. Joas, and M. Scheffler, J. Mater. Sci. 47, 7447 (2012). D. C. Langreth and J. P. Perdew, Sol. Stat. Comm. 17, 1425 (1975). O. Gunnarsson and B. I. Lundqvist, Phys. Rev. B 13, 4274 (1976). D. C. Langreth and J. P. Perdew, Phys. Rev. B 15, 2884 (1977). H. B. Callen and T. A. Welton, Phys. Rev. 83, 34 (1951). Y. Andersson, D. C. Langreth, and B. I. Lundqvist, Phys. Rev. Lett. 76, 102 (1996). J. F. Dobson and J. Wang, Phys. Rev. Lett. 82, 2123 (1999). J. Paier, X. Ren, P. Rinke, G. E. Scuseria, A. Gruneis, G. Kresse, and M. Scheffler, New J. Phys. 14, 043002 (2012). R. Baer and D. Neuhauser, J. Chem. Phys. 121, 9803 (2004). G. E. Scuseria, T. M. Henderson, and D. C. Sorensen, J. Chem. Phys. 129, 231101 (2008). J. G. Angyan, R. F. Liu, J. Toulouse, and G. Jansen, J. Chem. Theor. Comp 7, 3116 (2011). J. Paier, B. G. Janesko, T. M. Henderson, G. E. Scuseria, A. Gruneis, and G. Kresse, J. Chem. Phys. 132, 094103 (2010). L. Schimka, J. Harl, A. Stroppa, A. Gruneis, M. Marsman, F. Mittendorfer, and G. Kresse, Nature Materials 9, 741 (2010). D. Rocca, Z. Bai, R.-C. Li, and G. Galli, J. Chem. Phys. 136, 034111 (2012). X. G. Ren, P. Rinke, and M. Scheffler, Phys. Rev. B 80, 045402 (2009). H. V. Nguyen and S. de Gironcoli, Phys. Rev. B 79, 205114 (2009). J. Harl and G. Kresse, Phys. Rev. B 77, 045136 (2008). D. Lu, H.-V. Nguyen, and G. Galli, J. Chem. Phys. 133, 154110 (2010). R. Baer and E. Rabani, Nano Lett. 12, 2123 (2012). R. Baer and D. Neuhauser, J. Chem. Phys. 137, 051103 (2012). D. Neuhauser, E. Rabani, and R. Baer, J. Chem. Theor. Comp 9, 24 (2012). The numerical effort of the stochastic RPA calculation scales quadratically with system size provided the numerical effort involved in applying the underlying Hamiltonian to a wave function scales linearly with system size. This happens when the underlying Hamiltonian is derived from KS DFT since the potential is local. Furthermore, it will also be true for a generalized KS DFT Hamiltonian having short range exchange. When long-range exchange is present in the Hamiltonian, it is considerably more difficult to apply the Hamiltonian to a wave


Page 5 of 5

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

28

29 30 31

32 33 34 35 36 37


function in a linear scaling way. We denote occupied states by indices i,j,k,l=1,..,N and unoccupied states by a,b,c,d=N+1,…,M; the indices p,q,s,t=1,…,M denote "unrestricted" states. D. Neuhauser and R. Baer, J. Chem. Phys. 123, 204105 (2005). S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001). The occupied-occupied eigenvalues, although affecting R will not affect the RPA correlation energy calculated by Eq. (3). A proof will be given in a future publication. E. Rabani, B. Hetenyi, B. J. Berne, and L. E. Brus, J. Chem. Phys. 110, 5355 (1999). L. W. Wang and A. Zunger, J. Phys. Chem. 98, 2158 (1994). M. R. Wall and D. Neuhauser, J. Chem. Phys. 102, 8011 (1995). S. Toledo and E. Rabani, J. Comp. Phys. 180, 256 (2002). Deatails of this approach will be described in a future publication. To estimate the statistical error for I stocahstic iterations divide by the square root of I.


Expeditious Stochastic Calculation of Random-Phase Approximation

Recommend Documents