- Email: [email protected]
Towards a field-free quantum Monte Carlo approach to polarizabilities of excited states: Application to the n = 2 hydrogen atom
Chemical Physics Letters 445 (2007) 345–349 www.elsevier.com/locate/cplett
Towards a field-free quantum Monte Carlo approach to polarizabilities of excited states: Application to the n = 2 hydrogen atom Yu Li b
a,1
, Jan Vrbik b, Stuart M. Rothstein
a,c,*
a Department of Physics, Brock University, St. Catharines, Ontario, Canada L2S 3A1 Department of Mathematics, Brock University, St. Catharines, Ontario, Canada L2S 3A1 c Department of Chemistry, Brock University, St. Catharines, Ontario, Canada L2S 3A1
Received 3 February 2007; in final form 31 July 2007 Available online 7 August 2007
Abstract Accurate calculation of excited state polarizabilities by diffusion Monte Carlo remains elusive. Here, the straight-forward approach based on the finite-field approximation is problematic due to the shifting of the physical nodes upon imposition of an external field. In this Letter we take the first steps towards developing a field-free approach to polarizabilities of excited states. We explain why the physical nodes introduce surprising sampling difficulties for some of the estimators and how to overcome these for the simplest case, namely the n = 2 hydrogen atom. We point out confusion in the literature in the definition and calculation of polarizabilities for this system. 2007 Elsevier B.V. All rights reserved.
1. Introduction Decades of research devoted to developing diffusion Monte Carlo (DMC) have firmly established the advantages of this approach to electronic structure calculations: rapid convergence with basis-set size, favorable scaling with number of electrons, no calculation and storage of large numbers of integrals, and codes that are naturally suited for parallel computation. Monographs [1,2] and recent reviews, e.g. [3], have been devoted to these issues. Accurate calculation of polarizabilities, defined as the energy change upon imposition of an external field (or field gradient) in the zero field limit, remains elusive for excited states. Here the commonly used finite field approximation, e.g. [4,5], is problematic in DMC due to the shifting of the physical nodes on application of the field.
* Corresponding author. Address: Department of Chemistry, Brock University, 500 Glenridge Ave, St. Catharines, Ontario, Canada L2S 3A1. Fax: +1 905 682 9020. E-mail address: [email protected] (S.M. Rothstein). 1 Present address: Department of Electrical and Computer Engineering, McMaster University, Hamilton, Ontario, Canada L8S 4L8.
0009-2614/$ - see front matter 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2007.07.100
We are developing a field-free approach to polarizabilities of excited states; here two or more physical nodes are present. Estimators for field derivatives of the local energy within each nodal region are derived and sampled in a straight-forward manner. Also required are field derivatives of the nodal populations of walkers, estimated by the correlation of the dipole moment with nodal-region membership for walkers occupying the individual nodal regions. Sampling these quantities is non-trivial: a perfect DMC algorithm does not allow node crossing, and the correlation is unphysically large, in principle, infinite. It is the imperfections in the algorithm that allow node crossing, and this breaks the correlation, but not necessarily sufficiently fast to avoid a false covariance at small values of the time step. We show how to overcome this surprisingly difficult sampling problem for the simplest case, the n = 2 hydrogen atom, whose hybrid orbitals partition into two nodal regions. Our computed polarizabilities are in agreement with some, but not all of the values reported in the literature. We attribute their disagreements to confusion in the definition and calculation of polarizabilities for this system.
346
Y. Li et al. / Chemical Physics Letters 445 (2007) 345–349
2. Two nodal regions in an applied field Following Buckingham [6], when an atom or molecule is perturbed by a non-uniform, static electric field, F, the Hamiltonian is given by X 1X^ bF ¼ H b0 ^a Fa H ð1Þ l Hab F ab ; 3 ab a b 0 is the zero-field, unperturbed Hamiltonian, l ^ is where H ^ the quadrapole moment of the the dipole moment, and H atom or molecule. The Greek subscripts denote Cartesian variables, and Fa and Fab represent the electric field and field gradient, respectively. Our approach is to explicitly differentiate the estimator of the expected value of the perturbed local energy in the mixed distribution: R R /F b P ðpÞ b H Fw 2 w dR /F HwF w w dR w ij wij;F E L;F ðRij Þ w Ee;F ¼ R ; ¼ R/ 2 P ðpÞ F /F w dR w dR w ij wij;F ð2Þ
Restricted sampling is performed as follows: after a large number of standard drift-diffusion moves, ensuring the correct proportion of walkers are distributed among the nodal regions, we reject any drift-diffusion move that attempts to move a walker across a nodal boundary. The usual (but this time unrestricted) variational simulation with accumulated, weights provides a good estimate of WA via: R /F 2 P ðpÞ R I A w w dR I A w/F dR ij wij;F I A ðRij Þ W A;F ¼ R ; ð8Þ ¼ R/ 2 P ðpÞ F w/F dR w dR w ij wij;F where IA is the indicator function of being in Region A: it equals to 1 when we are inside A, becomes zero otherwise. This mode of sampling allows occasional, and effectively accidental, node crossing. Consider the special case of an electric field applied along the z-axis, Fz. The weight assigned to the ith walker on the jth iteration is given as follows: j Y
ðpÞ
wij;F z ¼
^zi;k F z EL;F z Þ; exp½sðEL;0 ðRi;k Þ l
ð9Þ
k¼jL
and then take its zero field limit. Here and below, w is an inputted, field-independent trial function, /F is the perturbed exact wave function, and separates a random variable on its left from its Monte Carlo estimator on its right. Rij is the jth generation of the ith walker, and the quantities ðpÞ wF are weighting factors, expfs½EL;F ðRÞ EL;F g, at time step s, accumulated over a specific number of past iterations, used to weigh the perturbed local energy, EL;F b H Fw . It is averaging with these weights that converts the simw ulated drift-diffusion (variational) distribution to the mixed one. For the special case of two nodal regions (say, A and B) Ee;F ¼ W A;F EAL;F þ W B;F EBL;F ; where, formally, R w/F dR A R W A;F ¼ R ; w/ dR þ B w/F dR F A R ^ F w=wÞw/F dR ðH A ; EL;F ¼ A R w/F dR A
ð3Þ
ð4Þ
ð10Þ
where f^ lz gij
j X
ð11Þ
^zi;k : l
k¼jL
P
P ðpÞ P ðpÞ ðpÞ ij ðowij;F z =oF z ÞI A ij wij;F z I A ij owij;F z =oF z P ðpÞ P ðpÞ P ðpÞ ij wij;F z ij wij;F z ij wij;F z
oW A =oF z jF z ¼0
P
ð5Þ
s
P ðpÞ P ðpÞ ðpÞ lz gij wij;0 I A lz gij wij;0 ij f^ ij wij;0 I A ij f^ P ðpÞ P ðpÞ P ðpÞ ij wij;0 ij wij;0 ij wij;0
! jF z ¼0
!
:
ð12Þ
ð6Þ
Note that in the zero-field limit the weights w(p) (Eq. (9)) involve only the unperturbed local energies. Thus, the above estimator and those quoted below are evaluated by way of a field-free simulation, without recourse to the so-called ‘finite-field approximation’. Rewriting Eq. (12) in a compact notation, where Æ æ denotes the expected value yields:
ð7Þ
oW A =oF z jF z ¼0 ¼ sðhI A f^ lz gi hI A i hf^ lz giÞ
and E00e;F ¼ W 00A;F ðEAL;F EBL;F Þ þ 2W 0A;F ðEAL;F 0 EBL;F 0 Þ þ W A;F EAL;F 00 þ W B;F EBL;F 00 ;
ðpÞ
owij;F z ðpÞ ¼ sf^ lz gij wij;F z ; oF z
Thus, from Eq. (8)
and WB,F = 1 WA,F. The first and second derivatives of the exact energy, Ee,F, with respect to field strength are given, respectively, as follows: E0e;F ¼ W 0A;F ðEAL;F EBL;F Þ þ W A;F EAL;F 0 þ W B;F EBL;F 0 ;
^ is the dipole moment, k the iteration label, L a where l chosen number of past iterations. and j the current iteration. It is straight-forward to derive an estimator for W0A , the derivative of WA with respect to Fz, and then take its zerofield limit.
where EAL;F 0 , and EAL;F 00 are estimates of derivatives of the energy, with sampling restricted to Region A or B. In the zero-field limit, these quantities are the dipole moment and a polarizability, respectively, for Region A.
¼ sCovðI A ; f^ lz gÞ;
ð13Þ
and oW A =oF z jF z ¼0 ¼ oW B =oF z jF z ¼0 :
ð14Þ
Y. Li et al. / Chemical Physics Letters 445 (2007) 345–349
By virtue of Eq. (11), the above covariance is the sum of covariances of the indicator function and dipole moment for the present and past iterations. For the second derivative of WA, we have W 00A jF z ¼0 ¼ s2 ðhI A f^ lz g2 i 2hI A f^ lz gi hf^ lz gi þ 2hI A i 2
l2z giÞ hf^ lz gi hI A i hf^ lz g; f^ lz gÞ; ¼ s2 CumðI A ; f^
ð15Þ
where the cumulant is defined as follows: CumðX ; Y ; ZÞ hðX hX iÞ ðY hY iÞ ðZ hZiÞi: ð16Þ Returning now to the energy derivatives, estimators are found in a similar manner, but done by restricted sampling: in Region A, implied by subscript A: 0
EAL jF z ¼0 ¼ h^ liA sCovðEL ; f^ lgÞA ;
ð17Þ
and 00
l; f^ lgÞA þ s2 CumðEL ; f^ lg; f^ lgÞA : EAL jF z ¼0 ¼ 2sCovð^
347
value: a[(1 log10(s)/b)/s], where a and b are determined from analysis of correlograms between iterations of the electron dipole moment. For 2px and 2py, the azz polarizabilities equal to 154.9(2.8) a.u. and 155.5(3.0) a.u., respectively (Fig. 1). Here, there are no Monte Carlo sampling difficulties. For these states, we can use unrestricted versions of Eqs. (17) and (18). Assigning Region ‘A’ where sp+ is positive, and ‘B’ where 0 it is negative, the zero time step limits for EAL (Eq. (17)) and B0 EL are equal to 1.580 6(38) a.u. and 3.397 1(87) a.u., respectively (Fig. 2). Also in that limit, WA (Eq. (8)) and WB are equal to 0.083 032(22) and 0.916 967(96), respectively (0.0803 and 0.9197, analytically) (Fig. 3). Thus, the dipole moment (Eq. (6)) for sp+ is equal to 2.98(1) a.u., in agreement with its analytical value of 3 a.u. Estimation of second derivatives is also straight-forward. 00 00 The zero time step limits for EAL (Eq. (18)) and EBL are equal to 24.70(64) and 137.5(2.8) a.u., respectively (Fig. 4).
ð18Þ 3. a polarizability for the n = 2 hydrogen atom It is well known [7] that the n = 2 hydrogen atom splits into four states in the presence of a static electric field along the z-axis: sp+ ” 2s + 2pz, degenerate 2px and 2py, and sp– ” 2s 2pz. The linear combinations of 2s and 2pz orbitals, the so-called ‘sp’ hybrid orbitals, have the non-trivial two nodal regions that are the focus of this Letter. Several simplifications arise by virtue of w being an exact solution to the unperturbed Schroedinger equation: all weights w(p) are unity, and the first terms in Eqs. (6) and (7) as well as the second terms in Eqs. (17) and (18) vanish. Algorithmic parameters were set as follows: The selected time steps are s = 0.05, . . . , 0.20 [0.05] a.u. Recognizing difficulties of Monte Carlo sampling of the nodal regions, we performed 20 000 iterations, a relatively large number, to equilibrate our ensembles to generate the correct distribution of configurations, 1000 at each time step. A configuration can not only accidentally cross a node, but can also land so close to it that the drift becomes extremely large, pushing an electron too far from a region of reasonable probability. To alleviate this, we truncated the velocity components such that the truncation vanishes in the s ! 0 limit [8]: Fi if jF i j 6 1=s; Fi ¼ ð19Þ 1=ssignðF i Þ otherwise: This standard notation for a component of the drift velocity should not be confused with the same symbol used above to define the electric field direction! To capture the serial correlation of the dipole moment sequence, we set the L parameter to have the following
Fig. 1. azz polarizability for 2px(h) and 2py (O). The data was fitted by linear regression.
Fig. 2. E0L for positive (bottom) and negative (top) nodal regions for sp+. The data was fitted by linear regression.
348
Y. Li et al. / Chemical Physics Letters 445 (2007) 345–349
To remedy this difficulty, we devised the following ‘double-extrapolation’ scheme: First, we truncate the drift velocity more severely than in Eq. (19) by F if jFj 6 F lim ; ð20Þ F¼ ðF lim =jFjÞF otherwise;
Fig. 3. Number of walkers in positive (bottom) and negative (top) nodal regions for sp+. There are 1000 walkers in total. The data was fitted by linear regression.
changing Flim so as to truncate 0.5, . . . , 2.5 [0.5] % of the velocities at s = 0.025 a.u. Next, we repeat this by doubling s and adjusting Flim to truncate twice as many velocities (1.0, . . . ,5.0 [1.0]%). We now have five pairs of estimates of W 0A , each with a time step bias. To remove this bias, we extrapolated each pair to zero time step (Table 1). We expect these zero-time step estimates of W 0A to be independent of truncation scheme; clearly they are not (Fig. 5). This is due to insufficient speed of nodal leakage, more pronounced for less severe truncation of F. The correct value is thus obtained in the limit of unrestricted migration between the two regions. To find this limit, we fit the data to an exponential model: W 0A ¼ A þ B expðx=CÞ;
ð21Þ 0 A.
where constant A(=3.98(57)) is our estimate of W Substituting this and other results quoted above into Eq. (7), the sp polarizability is equal to 168(3) a.u. It is surprising that ‘leakage’ of walkers across the nodes, due to imperfections in the DMC algorithm, is
Table 1 Estimation of W 0A as a function of drift truncation scheme, where Flim (Eq. (20)) is chosen to truncate a given % of drift velocities at each time step
Fig. 4. E00L for positive (bottom) and negative (top) nodal regions for sp+. The data was fitted by linear regression.
The main difficulty lies with estimating W0A (Eq. (13)). Covariance between any two quantities ‘normally’ tends to zero as their iteration lag increases to infinity. However, ^Þ ¼ given a perfect algorithm (no node crossing), Cov ðI A ; l W A ðh^ liA h^ liAþB Þ is non-zero, even after an infinite numP ^Þ is ber of iterations. Thus CovðI A ; f^ lgÞ ¼ CovðI A ; l infinite at all time steps. This ‘false covariance’ is broken by ‘nodal leakage’: accidental node crossing inherent in the standard, imperfect DMC algorithm. We discovered empirically that increasing the truncation of the drift velocity increases the probability of walkers migrating between the two regions, exploited below. For large time steps, leakage is so common that this false covariance does not interfere with the estimation of CovðI A ; f^ lgÞ, albeit at the cost of time step bias. The reverse is true for small time steps, when leakage becomes so slow that the false covariance starts to dominate.
Truncation scheme
s = 0.050 % s = 0.025 % s=0 Truncation Truncation
1 2 3 4 5
5.27(26) 3.91(13) 3.47(12) 3.32(11) 3.40(9)
1.0 2.0 3.0 4.0 5.0
7.59(29) 5.34(23) 4.60(17) 4.14(16) 3.86(16)
0.5 1.0 1.5 2.0 2.5
Zero time step estimate obtained by linear regression.
Fig. 5. Fit to exponential model for W 0A for sp+.
9.91(53) 6.78(41) 5.73(34) 4.97(31) 4.31(29)
Y. Li et al. / Chemical Physics Letters 445 (2007) 345–349
349
Table 2 Static a-polarizability for the second excited state of hydrogen atom
Acknowledgements
Method
This work was supported, in part, by grants from the Natural Sciences and Engineering Research Council of Canada (NSERC). We acknowledge and thank Harris J. Silverstone for his guidance in interpreting the conflicting literature values that appear in Table 2. We also thank the referees for their helpful comments.
a
Parabolic coordinate separation Green’s functionb Perturbation theoryc Dalgarno and Lewis method for perturbation theoryd B-splinese Coulomb Green’s functionf Partial sum approximationsg Mapped-Fourier grid methodh This work
State cited
a (a.u)
sp 2s 2p sp
168.0 120 176 167.5
2p 2s 2pjmj=1 2s/2pjmj=1/2p0 sp/2pjmj=1
216 120 156 120/156/216 168(3)/155(3)
˚ 3. a:u: ¼ e2 a20 ¼ 1.648778 · 1041 C2 m2 J1 = 0.148185 A a Ref. [9]. b Ref. [10]. c Ref. [11]. d Ref. [12]. e Ref. [13]. f Ref. [14]. g Ref. [15]. h Ref. [16].
sufficient to accurately estimate the values of nodal populations WA and WB. This is only with the s1 drift truncation, (Eq. (19)), but not with the % truncation scheme used here to estimate the population field derivates. Our calculations are in good agreement with previous ones for the sp hybrid states, and for the 2p states orthogonal to them (Table 2). For papers addressing polarizabilities of the ‘2s’ and ‘2p’ or ‘2p0’ states, we emphasize there are no such states under an external electric field applied along the z-axis, and that these two states should be combined to form the correct physical states. We believe in most cases discrepancies have arisen from the authors concentrating on the verification of their methodology rather than on the physical quantity itself.
References [1] B.L. Hammond, W.A. Lester Jr., P.J. Reynolds, Monte Carlo Methods in Ab Initio Quantum Chemistry, World Scientific, Singapore, 1994. [2] J.B. Anderson, Quantum Monte Carlo: Origins, Development, Applications Oxford, 2007. [3] A. Aspuru-Guzik, W.A. Lester Jr., in: P.G. Ciarlet, C. Le Bris (Eds.), Handbook of Numerical Analysis: Computational Chemistry, vol. 10, North Holland, Amsterdam, 2003, p. 485. [4] J. Kobus, D. Moncrieff, S. Wilson, J. Phys. B 40 (2007) 877. [5] S.A. Alexander, R.L. Coldwell, G. Aissing, A.J. Thakkar, Int. J. Quant. Chem. S26 (1992) 213. [6] A.D. Buckingham, Adv. Chem. Phys. 12 (1967) 107. [7] J.P. Lowe, Quantum Chemistry, Academic Press, 1993, p. 411. [8] M.F. DePasquale, S.M. Rothstein, J. Vrbik, J. Chem. Phys. 89 (1988) 3629. [9] P.S. Epstein, Phys. Rev. 28 (1926) 695. [10] B.J. Laurenzi, D.G. Williams, G.S. Bhatia, J. Chem. Phys. 61 (1974) 2077; K. McDowell, J. Chem. Phys. 65 (1976) 2518. [11] A.A. Radzig, B.M. Smirnov, Reference Data on Atoms, Molecules and Ions, Springer, Berlin, 1985, p. 119. [12] G.K. Sapra, V.S. Bhasin, L.S. Kothari, Int. J. Mod. Phys. 8 (1994) 727. [13] M.I. Bhatti, K.D. Coleman, W.F. Perger, Phys. Rev. A 68 (2003) 044503. [14] A.A. Krylovetsky, N.L. Manakov, S.I. Marmo, J. Phys. B 38 (2005) 311. [15] V. Koch, D. Andrae, Theor. Chem. Acc. 114 (2005) 380. [16] S. Cohen, S.I. Themelis, J. Chem. Phys. 124 (2006) 134106.
Towards a field-free quantum Monte Carlo approach to polarizabilities of excited states: Application to the n = 2 hydrogen atom Yu Li b
a,1
, Jan Vrbik b, Stuart M. Rothstein
a,c,*
a Department of Physics, Brock University, St. Catharines, Ontario, Canada L2S 3A1 Department of Mathematics, Brock University, St. Catharines, Ontario, Canada L2S 3A1 c Department of Chemistry, Brock University, St. Catharines, Ontario, Canada L2S 3A1
Received 3 February 2007; in final form 31 July 2007 Available online 7 August 2007
Abstract Accurate calculation of excited state polarizabilities by diffusion Monte Carlo remains elusive. Here, the straight-forward approach based on the finite-field approximation is problematic due to the shifting of the physical nodes upon imposition of an external field. In this Letter we take the first steps towards developing a field-free approach to polarizabilities of excited states. We explain why the physical nodes introduce surprising sampling difficulties for some of the estimators and how to overcome these for the simplest case, namely the n = 2 hydrogen atom. We point out confusion in the literature in the definition and calculation of polarizabilities for this system. 2007 Elsevier B.V. All rights reserved.
1. Introduction Decades of research devoted to developing diffusion Monte Carlo (DMC) have firmly established the advantages of this approach to electronic structure calculations: rapid convergence with basis-set size, favorable scaling with number of electrons, no calculation and storage of large numbers of integrals, and codes that are naturally suited for parallel computation. Monographs [1,2] and recent reviews, e.g. [3], have been devoted to these issues. Accurate calculation of polarizabilities, defined as the energy change upon imposition of an external field (or field gradient) in the zero field limit, remains elusive for excited states. Here the commonly used finite field approximation, e.g. [4,5], is problematic in DMC due to the shifting of the physical nodes on application of the field.
* Corresponding author. Address: Department of Chemistry, Brock University, 500 Glenridge Ave, St. Catharines, Ontario, Canada L2S 3A1. Fax: +1 905 682 9020. E-mail address: [email protected] (S.M. Rothstein). 1 Present address: Department of Electrical and Computer Engineering, McMaster University, Hamilton, Ontario, Canada L8S 4L8.
0009-2614/$ - see front matter 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2007.07.100
We are developing a field-free approach to polarizabilities of excited states; here two or more physical nodes are present. Estimators for field derivatives of the local energy within each nodal region are derived and sampled in a straight-forward manner. Also required are field derivatives of the nodal populations of walkers, estimated by the correlation of the dipole moment with nodal-region membership for walkers occupying the individual nodal regions. Sampling these quantities is non-trivial: a perfect DMC algorithm does not allow node crossing, and the correlation is unphysically large, in principle, infinite. It is the imperfections in the algorithm that allow node crossing, and this breaks the correlation, but not necessarily sufficiently fast to avoid a false covariance at small values of the time step. We show how to overcome this surprisingly difficult sampling problem for the simplest case, the n = 2 hydrogen atom, whose hybrid orbitals partition into two nodal regions. Our computed polarizabilities are in agreement with some, but not all of the values reported in the literature. We attribute their disagreements to confusion in the definition and calculation of polarizabilities for this system.
346
Y. Li et al. / Chemical Physics Letters 445 (2007) 345–349
2. Two nodal regions in an applied field Following Buckingham [6], when an atom or molecule is perturbed by a non-uniform, static electric field, F, the Hamiltonian is given by X 1X^ bF ¼ H b0 ^a Fa H ð1Þ l Hab F ab ; 3 ab a b 0 is the zero-field, unperturbed Hamiltonian, l ^ is where H ^ the quadrapole moment of the the dipole moment, and H atom or molecule. The Greek subscripts denote Cartesian variables, and Fa and Fab represent the electric field and field gradient, respectively. Our approach is to explicitly differentiate the estimator of the expected value of the perturbed local energy in the mixed distribution: R R /F b P ðpÞ b H Fw 2 w dR /F HwF w w dR w ij wij;F E L;F ðRij Þ w Ee;F ¼ R ; ¼ R/ 2 P ðpÞ F /F w dR w dR w ij wij;F ð2Þ
Restricted sampling is performed as follows: after a large number of standard drift-diffusion moves, ensuring the correct proportion of walkers are distributed among the nodal regions, we reject any drift-diffusion move that attempts to move a walker across a nodal boundary. The usual (but this time unrestricted) variational simulation with accumulated, weights provides a good estimate of WA via: R /F 2 P ðpÞ R I A w w dR I A w/F dR ij wij;F I A ðRij Þ W A;F ¼ R ; ð8Þ ¼ R/ 2 P ðpÞ F w/F dR w dR w ij wij;F where IA is the indicator function of being in Region A: it equals to 1 when we are inside A, becomes zero otherwise. This mode of sampling allows occasional, and effectively accidental, node crossing. Consider the special case of an electric field applied along the z-axis, Fz. The weight assigned to the ith walker on the jth iteration is given as follows: j Y
ðpÞ
wij;F z ¼
^zi;k F z EL;F z Þ; exp½sðEL;0 ðRi;k Þ l
ð9Þ
k¼jL
and then take its zero field limit. Here and below, w is an inputted, field-independent trial function, /F is the perturbed exact wave function, and separates a random variable on its left from its Monte Carlo estimator on its right. Rij is the jth generation of the ith walker, and the quantities ðpÞ wF are weighting factors, expfs½EL;F ðRÞ EL;F g, at time step s, accumulated over a specific number of past iterations, used to weigh the perturbed local energy, EL;F b H Fw . It is averaging with these weights that converts the simw ulated drift-diffusion (variational) distribution to the mixed one. For the special case of two nodal regions (say, A and B) Ee;F ¼ W A;F EAL;F þ W B;F EBL;F ; where, formally, R w/F dR A R W A;F ¼ R ; w/ dR þ B w/F dR F A R ^ F w=wÞw/F dR ðH A ; EL;F ¼ A R w/F dR A
ð3Þ
ð4Þ
ð10Þ
where f^ lz gij
j X
ð11Þ
^zi;k : l
k¼jL
P
P ðpÞ P ðpÞ ðpÞ ij ðowij;F z =oF z ÞI A ij wij;F z I A ij owij;F z =oF z P ðpÞ P ðpÞ P ðpÞ ij wij;F z ij wij;F z ij wij;F z
oW A =oF z jF z ¼0
P
ð5Þ
s
P ðpÞ P ðpÞ ðpÞ lz gij wij;0 I A lz gij wij;0 ij f^ ij wij;0 I A ij f^ P ðpÞ P ðpÞ P ðpÞ ij wij;0 ij wij;0 ij wij;0
! jF z ¼0
!
:
ð12Þ
ð6Þ
Note that in the zero-field limit the weights w(p) (Eq. (9)) involve only the unperturbed local energies. Thus, the above estimator and those quoted below are evaluated by way of a field-free simulation, without recourse to the so-called ‘finite-field approximation’. Rewriting Eq. (12) in a compact notation, where Æ æ denotes the expected value yields:
ð7Þ
oW A =oF z jF z ¼0 ¼ sðhI A f^ lz gi hI A i hf^ lz giÞ
and E00e;F ¼ W 00A;F ðEAL;F EBL;F Þ þ 2W 0A;F ðEAL;F 0 EBL;F 0 Þ þ W A;F EAL;F 00 þ W B;F EBL;F 00 ;
ðpÞ
owij;F z ðpÞ ¼ sf^ lz gij wij;F z ; oF z
Thus, from Eq. (8)
and WB,F = 1 WA,F. The first and second derivatives of the exact energy, Ee,F, with respect to field strength are given, respectively, as follows: E0e;F ¼ W 0A;F ðEAL;F EBL;F Þ þ W A;F EAL;F 0 þ W B;F EBL;F 0 ;
^ is the dipole moment, k the iteration label, L a where l chosen number of past iterations. and j the current iteration. It is straight-forward to derive an estimator for W0A , the derivative of WA with respect to Fz, and then take its zerofield limit.
where EAL;F 0 , and EAL;F 00 are estimates of derivatives of the energy, with sampling restricted to Region A or B. In the zero-field limit, these quantities are the dipole moment and a polarizability, respectively, for Region A.
¼ sCovðI A ; f^ lz gÞ;
ð13Þ
and oW A =oF z jF z ¼0 ¼ oW B =oF z jF z ¼0 :
ð14Þ
Y. Li et al. / Chemical Physics Letters 445 (2007) 345–349
By virtue of Eq. (11), the above covariance is the sum of covariances of the indicator function and dipole moment for the present and past iterations. For the second derivative of WA, we have W 00A jF z ¼0 ¼ s2 ðhI A f^ lz g2 i 2hI A f^ lz gi hf^ lz gi þ 2hI A i 2
l2z giÞ hf^ lz gi hI A i hf^ lz g; f^ lz gÞ; ¼ s2 CumðI A ; f^
ð15Þ
where the cumulant is defined as follows: CumðX ; Y ; ZÞ hðX hX iÞ ðY hY iÞ ðZ hZiÞi: ð16Þ Returning now to the energy derivatives, estimators are found in a similar manner, but done by restricted sampling: in Region A, implied by subscript A: 0
EAL jF z ¼0 ¼ h^ liA sCovðEL ; f^ lgÞA ;
ð17Þ
and 00
l; f^ lgÞA þ s2 CumðEL ; f^ lg; f^ lgÞA : EAL jF z ¼0 ¼ 2sCovð^
347
value: a[(1 log10(s)/b)/s], where a and b are determined from analysis of correlograms between iterations of the electron dipole moment. For 2px and 2py, the azz polarizabilities equal to 154.9(2.8) a.u. and 155.5(3.0) a.u., respectively (Fig. 1). Here, there are no Monte Carlo sampling difficulties. For these states, we can use unrestricted versions of Eqs. (17) and (18). Assigning Region ‘A’ where sp+ is positive, and ‘B’ where 0 it is negative, the zero time step limits for EAL (Eq. (17)) and B0 EL are equal to 1.580 6(38) a.u. and 3.397 1(87) a.u., respectively (Fig. 2). Also in that limit, WA (Eq. (8)) and WB are equal to 0.083 032(22) and 0.916 967(96), respectively (0.0803 and 0.9197, analytically) (Fig. 3). Thus, the dipole moment (Eq. (6)) for sp+ is equal to 2.98(1) a.u., in agreement with its analytical value of 3 a.u. Estimation of second derivatives is also straight-forward. 00 00 The zero time step limits for EAL (Eq. (18)) and EBL are equal to 24.70(64) and 137.5(2.8) a.u., respectively (Fig. 4).
ð18Þ 3. a polarizability for the n = 2 hydrogen atom It is well known [7] that the n = 2 hydrogen atom splits into four states in the presence of a static electric field along the z-axis: sp+ ” 2s + 2pz, degenerate 2px and 2py, and sp– ” 2s 2pz. The linear combinations of 2s and 2pz orbitals, the so-called ‘sp’ hybrid orbitals, have the non-trivial two nodal regions that are the focus of this Letter. Several simplifications arise by virtue of w being an exact solution to the unperturbed Schroedinger equation: all weights w(p) are unity, and the first terms in Eqs. (6) and (7) as well as the second terms in Eqs. (17) and (18) vanish. Algorithmic parameters were set as follows: The selected time steps are s = 0.05, . . . , 0.20 [0.05] a.u. Recognizing difficulties of Monte Carlo sampling of the nodal regions, we performed 20 000 iterations, a relatively large number, to equilibrate our ensembles to generate the correct distribution of configurations, 1000 at each time step. A configuration can not only accidentally cross a node, but can also land so close to it that the drift becomes extremely large, pushing an electron too far from a region of reasonable probability. To alleviate this, we truncated the velocity components such that the truncation vanishes in the s ! 0 limit [8]: Fi if jF i j 6 1=s; Fi ¼ ð19Þ 1=ssignðF i Þ otherwise: This standard notation for a component of the drift velocity should not be confused with the same symbol used above to define the electric field direction! To capture the serial correlation of the dipole moment sequence, we set the L parameter to have the following
Fig. 1. azz polarizability for 2px(h) and 2py (O). The data was fitted by linear regression.
Fig. 2. E0L for positive (bottom) and negative (top) nodal regions for sp+. The data was fitted by linear regression.
348
Y. Li et al. / Chemical Physics Letters 445 (2007) 345–349
To remedy this difficulty, we devised the following ‘double-extrapolation’ scheme: First, we truncate the drift velocity more severely than in Eq. (19) by F if jFj 6 F lim ; ð20Þ F¼ ðF lim =jFjÞF otherwise;
Fig. 3. Number of walkers in positive (bottom) and negative (top) nodal regions for sp+. There are 1000 walkers in total. The data was fitted by linear regression.
changing Flim so as to truncate 0.5, . . . , 2.5 [0.5] % of the velocities at s = 0.025 a.u. Next, we repeat this by doubling s and adjusting Flim to truncate twice as many velocities (1.0, . . . ,5.0 [1.0]%). We now have five pairs of estimates of W 0A , each with a time step bias. To remove this bias, we extrapolated each pair to zero time step (Table 1). We expect these zero-time step estimates of W 0A to be independent of truncation scheme; clearly they are not (Fig. 5). This is due to insufficient speed of nodal leakage, more pronounced for less severe truncation of F. The correct value is thus obtained in the limit of unrestricted migration between the two regions. To find this limit, we fit the data to an exponential model: W 0A ¼ A þ B expðx=CÞ;
ð21Þ 0 A.
where constant A(=3.98(57)) is our estimate of W Substituting this and other results quoted above into Eq. (7), the sp polarizability is equal to 168(3) a.u. It is surprising that ‘leakage’ of walkers across the nodes, due to imperfections in the DMC algorithm, is
Table 1 Estimation of W 0A as a function of drift truncation scheme, where Flim (Eq. (20)) is chosen to truncate a given % of drift velocities at each time step
Fig. 4. E00L for positive (bottom) and negative (top) nodal regions for sp+. The data was fitted by linear regression.
The main difficulty lies with estimating W0A (Eq. (13)). Covariance between any two quantities ‘normally’ tends to zero as their iteration lag increases to infinity. However, ^Þ ¼ given a perfect algorithm (no node crossing), Cov ðI A ; l W A ðh^ liA h^ liAþB Þ is non-zero, even after an infinite numP ^Þ is ber of iterations. Thus CovðI A ; f^ lgÞ ¼ CovðI A ; l infinite at all time steps. This ‘false covariance’ is broken by ‘nodal leakage’: accidental node crossing inherent in the standard, imperfect DMC algorithm. We discovered empirically that increasing the truncation of the drift velocity increases the probability of walkers migrating between the two regions, exploited below. For large time steps, leakage is so common that this false covariance does not interfere with the estimation of CovðI A ; f^ lgÞ, albeit at the cost of time step bias. The reverse is true for small time steps, when leakage becomes so slow that the false covariance starts to dominate.
Truncation scheme
s = 0.050 % s = 0.025 % s=0 Truncation Truncation
1 2 3 4 5
5.27(26) 3.91(13) 3.47(12) 3.32(11) 3.40(9)
1.0 2.0 3.0 4.0 5.0
7.59(29) 5.34(23) 4.60(17) 4.14(16) 3.86(16)
0.5 1.0 1.5 2.0 2.5
Zero time step estimate obtained by linear regression.
Fig. 5. Fit to exponential model for W 0A for sp+.
9.91(53) 6.78(41) 5.73(34) 4.97(31) 4.31(29)
Y. Li et al. / Chemical Physics Letters 445 (2007) 345–349
349
Table 2 Static a-polarizability for the second excited state of hydrogen atom
Acknowledgements
Method
This work was supported, in part, by grants from the Natural Sciences and Engineering Research Council of Canada (NSERC). We acknowledge and thank Harris J. Silverstone for his guidance in interpreting the conflicting literature values that appear in Table 2. We also thank the referees for their helpful comments.
a
Parabolic coordinate separation Green’s functionb Perturbation theoryc Dalgarno and Lewis method for perturbation theoryd B-splinese Coulomb Green’s functionf Partial sum approximationsg Mapped-Fourier grid methodh This work
State cited
a (a.u)
sp 2s 2p sp
168.0 120 176 167.5
2p 2s 2pjmj=1 2s/2pjmj=1/2p0 sp/2pjmj=1
216 120 156 120/156/216 168(3)/155(3)
˚ 3. a:u: ¼ e2 a20 ¼ 1.648778 · 1041 C2 m2 J1 = 0.148185 A a Ref. [9]. b Ref. [10]. c Ref. [11]. d Ref. [12]. e Ref. [13]. f Ref. [14]. g Ref. [15]. h Ref. [16].
sufficient to accurately estimate the values of nodal populations WA and WB. This is only with the s1 drift truncation, (Eq. (19)), but not with the % truncation scheme used here to estimate the population field derivates. Our calculations are in good agreement with previous ones for the sp hybrid states, and for the 2p states orthogonal to them (Table 2). For papers addressing polarizabilities of the ‘2s’ and ‘2p’ or ‘2p0’ states, we emphasize there are no such states under an external electric field applied along the z-axis, and that these two states should be combined to form the correct physical states. We believe in most cases discrepancies have arisen from the authors concentrating on the verification of their methodology rather than on the physical quantity itself.
References [1] B.L. Hammond, W.A. Lester Jr., P.J. Reynolds, Monte Carlo Methods in Ab Initio Quantum Chemistry, World Scientific, Singapore, 1994. [2] J.B. Anderson, Quantum Monte Carlo: Origins, Development, Applications Oxford, 2007. [3] A. Aspuru-Guzik, W.A. Lester Jr., in: P.G. Ciarlet, C. Le Bris (Eds.), Handbook of Numerical Analysis: Computational Chemistry, vol. 10, North Holland, Amsterdam, 2003, p. 485. [4] J. Kobus, D. Moncrieff, S. Wilson, J. Phys. B 40 (2007) 877. [5] S.A. Alexander, R.L. Coldwell, G. Aissing, A.J. Thakkar, Int. J. Quant. Chem. S26 (1992) 213. [6] A.D. Buckingham, Adv. Chem. Phys. 12 (1967) 107. [7] J.P. Lowe, Quantum Chemistry, Academic Press, 1993, p. 411. [8] M.F. DePasquale, S.M. Rothstein, J. Vrbik, J. Chem. Phys. 89 (1988) 3629. [9] P.S. Epstein, Phys. Rev. 28 (1926) 695. [10] B.J. Laurenzi, D.G. Williams, G.S. Bhatia, J. Chem. Phys. 61 (1974) 2077; K. McDowell, J. Chem. Phys. 65 (1976) 2518. [11] A.A. Radzig, B.M. Smirnov, Reference Data on Atoms, Molecules and Ions, Springer, Berlin, 1985, p. 119. [12] G.K. Sapra, V.S. Bhasin, L.S. Kothari, Int. J. Mod. Phys. 8 (1994) 727. [13] M.I. Bhatti, K.D. Coleman, W.F. Perger, Phys. Rev. A 68 (2003) 044503. [14] A.A. Krylovetsky, N.L. Manakov, S.I. Marmo, J. Phys. B 38 (2005) 311. [15] V. Koch, D. Andrae, Theor. Chem. Acc. 114 (2005) 380. [16] S. Cohen, S.I. Themelis, J. Chem. Phys. 124 (2006) 134106.