- Email: [email protected]
Effective numerical method for theoretical studies of small atomic clusters
Computational Materials Science 35 (2006) 359–365 www.elsevier.com/locate/commatsci
Effective numerical method for theoretical studies of small atomic clusters B. Gervais a, E. Giglio
a,*
, A. Ipatov
a,b,1
, J. Douady
a
a
b
CIRIL, unite´ mixte CEA-CNRS-ENSICAEN, BP 5133, F-14070 Caen Cedex 05, France St. Petersburg Polytechnical University, Polytechnicheskaja 29, 195251 St. Petersburg, Russian Federation Received 24 April 2004; accepted 27 July 2004
Abstract We present an effective numerical method to study the electronic and ionic dynamics of small atomic systems, like alkali-metal clusters. Our approach is based on the density functional theory (DFT) combined with molecular dynamics (MD). The time-dependent Kohn–Sham (KS) equations describing the electronic subsystem are solved self-consistently in the coordinate space. Their numerical solution is performed using the time-splitting method on a three-dimensional coordinate grid by means of the generalized pseudospectral method (GPS). The KS orbitals are expanded over a spherical harmonics basis, and their radial part are represented on a non-uniform radial grid. We employ an optimal radial mapping of the collocation points consistent with the GPS method, which allows achieving a higher density of radial points inside the cluster where electronic density is larger and the potential varies rapidly. This method allows evaluating both the Coulomb potential and the kinetic energy operator with high accuracy. Unlike FFT based methods, the GPS method with non-uniform radial mapping, allows to use large grids at lower cost and thus to study efficiently negatively charged cluster. As an example, we present the number of emitted electrons of a negatively charged sodium cluster Na 7 , irradiated by a femtosecond laser pulse, as a function of the grid size. We show that the value of the electron emission converges only for large grids. 2005 Elsevier B.V. All rights reserved. PACS: 31.15.Ev; 31.15.p; 31.15.Qg; 36.40.c Keywords: Local density approximation; Time splitting; Harmonic expansion; Grid mapping; Spectral methods; 3D TDDFT; Molecular dynamics
1. Introduction The electronic and structural properties of small clusters, in particular alkali clusters, have been the subject of numerous recent investigations, from both experimental and theoretical point of view (for overview see [1–6]). The theoretical methods employing the real-time approach [7–9] are known as a powerful tool to study the optical properties of finite many-electron systems,
*
1
Corresponding author. Fax: +33 2 3145 4714. E-mail address: [email protected] (E. Giglio). Acknowledgments to CEA and CNRS for financial support.
0927-0256/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2004.07.015
their excitation processes, photoabsorption and electron emission. For instance, the surface dipole plasmon resonance is one of the most important collective excitations observed in metallic clusters [3,4,10]. The recent development of laser technologies and the use of femtosecond lasers has led to experimental studies [11–13] of clusters in strong laser field and stimulated the development of theoretical methods that allow for an adequate description of optical response beyond the linear regime [7, 14–19]. The real-time methods describe the dynamical properties in highly excited clusters in a natural way and can be effectively employed in studies of clusters ionization, thermal relaxation and fission for free, deposited and embedded systems.
360
B. Gervais et al. / Computational Materials Science 35 (2006) 359–365
The main goal of this work is to present a numerical method based on the real-time technique that possesses the simplicity of the direct discretization methods with high accuracy and computational effectiveness. Our method is based on the density functional theory (DFT) combined with molecular dynamics (MD). Assuming that a proper theoretical description of the system requires to account for the coupled dynamics of both electronic and ionic degrees of freedom, we treat the real-time motion of the electronic system using the time-dependent local density approximation (TDLDA) [7,9] method with average-density self-interaction correction ADSIC [20] in full coordinate-space grid, while the ionic motion of the cluster core is described using the standard methods of MD. We evaluate the cluster electronic system in real time, that allows to go beyond the adiabatic description of the system dynamics [4,5] and to consider, for instance, relatively long time scales, i.e., typically from a few hundreds fs to a few ps, accounting for ionic motion. We generalize the GPS Method introduced by Chu et al. mainly for atomic systems [21–24]. The latter approach was developed for atomic systems with cylindrical symmetry that naturally reduces to 2D problems. We employ for the first time this method for a full 3D computation and for atomic clusters. It allows to investigate the dynamics of such systems in providing self-consistent solution of the time-dependent Kohn–Sham equations. The proposed approach allows to improve a lot the efficiency of the numerical algorithms in comparison with alternative methods based on the FFT approach [25,8]. Especially for small and middle range clusters, for which larger grids can be used and longer time scales computed. Atomic system of units (h = e = me = 1) is used throughout the paper.
2. Theoretical approach The valence electrons of the Na cluster are described using the density functional theory methods [28]. The single-electron wave functions wa(r, t) follow the timedependent Kohn–Sham [30] equations d D i wa ðr; tÞ ¼ þ V KS ½r; q" ðr; tÞ; q# ðr; tÞ wa ðr; tÞ; dt 2 ð1Þ where " and # stands for the spin of electrons, with the corresponding electronic density defined as qr ðr; tÞ ¼ PN r 2 jw ðr; tÞj where r = ", # and N" and N# are the i i number of spin up and spin down electrons, respectively. The total density is q(r, t) = q"(r, t) + q#(r, t), and the number of electrons Ne = N" + N#. For such a system, the Kohn–Sham potential for spin up VKS" reads: VKS"(r, t) = Vext(r, t) + VCoul(r, t) +
Vxc"(r, t), where the external potential Vext takes into account the interaction with external field and ionic cores. The Coulomb potential VCoul is determined by solving the Poisson equation without retardation. Within ADSIC approximation [20], the exchange-correlation potential Vxc"(r, t) reads: Z 1 qðr0 Þ dr0 þ V LSD V xc" ¼ xc" ½q" ðr; tÞ; q# ðr; tÞ N jr r0 j q" ðr; tÞ LSD ;0 ; ð2Þ V xc" N" where V LSD xc" is a local density exchange and correlation potential for spin up. Similar equations are derived for spin down. The interaction among the ions is modeled as a simple Coulomb interaction of point charges. The molecular dynamics equations of ionic motion reduces to the following: " Z d Pi ¼ rRi V ext ðRi ; tÞ qðr; tÞV ps ðjRi rjÞ dr dt # X 1 þ ; ð3Þ jRi Rj j j6¼i d Pi Ri ¼ ; dt M ion
ð4Þ
where Mion is the ionic mass, Ri and Pi the position and momentum of the ion respectively.
3. Computational approach 3.1. Time propagation of the KS orbitals The computation of the Kohn–Sham orbitals is performed on a 3D spherical grid. We define va(r, t) = rwa(r, t), so that the KS equation (1) reads: " # 2 b L o 1 d2 þ þ Vb r ðtÞ va ðr; tÞ; ð5Þ i va ðr; tÞ ¼ ot 2 dr2 2r2 where Vb r ðtÞ is the time-dependent Kohn–Sham potential for spin r. We use a second order time-splitting algorithm [7] for the time propagation of the KS orbitals: bU b 2 ðtÞfa ðr; tÞ þ Oðdt3 Þ; fa ðr; t þ dtÞ ¼ Q r
ð6Þ
to the KS orbitals by the relation where fa is connected b y ðr; tÞva ðr; tÞ, with fr ðr; tÞ ¼ U r dt b b U r ðr; tÞ ¼ exp ið V r ðr; tÞ V 0 ðrÞÞ ; ð7Þ 2 b ¼ exp i H b 0 dt ; Q ð8Þ
B. Gervais et al. / Computational Materials Science 35 (2006) 359–365 2
b 0 ¼ 1 d22 þ bL 2 þ V 0 ðrÞ and V0(r) is a spheriwhere H 2 dr 2r cally symmetric potential. It is usually chosen as the spherical average of the total potential at the initial time V(r, t = 0). However, it can be updated to account for possible significant modifications of the potential in the course of the system dynamics. This partitioning of the potential part of the propagator between the timedependent and time-independent operators permits to reduce the error resulting from the time-splitting b does not depend on time t so that scheme. Note that Q it can be computed once for each given V0(r). The b numerical evaluation of the action of the operator Q on va(r, t) depends on the radial grid choice and is detailed in the next sections. For local potentials, the numerical evaluation of b r ðtÞva ðr; tÞ is straightforward. Indeed, for any given U set of collocation points {rg} it requires only the evalua b r ðtÞva ðr; tÞ tion the potential V(rg, t) to obtain U ¼ r¼rg b r ðrg ; tÞva ðrg ; tÞ. For non-local potentials of integral U form, the method is applicable as well. Indeed, the action of the operator Vb r ðtÞ on va(r, t) can be determined accurately by Gauss quadrature: Z Vb r ðtÞva ðtÞ r¼rg ¼ dr0 V r ðrg ; r0 ; tÞva ðr0 ; tÞ X wg0 V r ðrg ; rg0 ; tÞva ðrg0 ; tÞ. ð9Þ ¼ g0
The algebraic equation (9) can be easily put into a symmetric form and then diagonalized to evaluate the corresponding exponential matrix U r ðrg ; rg0 ; tÞ. The method is efficient provided the non-local part extent over a limited set of grid points so that: (i) the matrix V r ðrg ; rg0 Þ is sparse and (ii) it is block diagonal. This is usually the case for standard pseudopotentials like for example those proposed by Goedecker et al. [26] and for interatomic distances encountered in normal situation where the ions move in a limited range around their equilibrium positions, so that the bond length does not become excessively short. 3.2. Spherical harmonic expansion The grid points we use are distributed on concentric spheres of radius ri at solid angles Xq = (ha, /b). The N/ azimuthal angle coordinate points /b are uniformly spaced in the interval [0, 2p] so that /b ¼ N2p/ . The Nh polar angles ha are distributed non-uniformly with ya = cos(ha) taken as the roots of the Legendre polynomials P N h ðyÞ. Since we solve the 3D problem on a spherical grid, it is natural to split the operator into the radial and the angular parts, using the expansion of the full 3D functions va(r, X) into spherical harmonics Yp(h, /), with p l, m. Moreover, the use of spherical harmonics allows eliminating the so-called pole problem [27]. We define the radial functions Rap as the projection over the spherical harmonics Yp:
va ðr; X; tÞ ¼ Rap ðr; tÞ
¼
Z
X
361
Rap ðr; tÞY p ðh; /Þ;
ð10Þ
p
ðva ðr; XÞÞY p ðh; /Þ dX.
ð11Þ
Within this approach, the operator (8) now reads: X Y p Q b p Y p ; b ¼ Q ð12Þ p
b p i ¼ expði H b p ¼ hY p j QjY b 0;l dtÞ, is a sum of rawhere Q dial operators acting in 1-dimensional-subspace, for which 2 b 0;l ¼ 1 d þ lðl þ 1Þ þ V 0 ðrÞ. H 2 dr2 2r2
ð13Þ
The projections (11) are performed efficiently at each time step using a method based on the FFT algorithm for integration over the azimuthal angle /, and Gauss–Legendre quadrature for integration over the azimuthal angle h: Z nm ðri ; hj Þ ¼ vðri ; hj ; /Þ expðim/Þ d/ N
/ 2p X vðri ; hj ; /k Þ expðim/k Þ; N / k¼1 Z jmj Rlm ðri Þ ¼ C lm nm ðri ; hÞP l ðcos hÞ d cos h
¼
¼ C lm
Nh X
ð14Þ
jmj
wj nm ðri ; hj ÞP l ðcos hj Þ;
j¼1
where Clm is the normalization constant of the spherical harmonics, P ml are the associated Legendre polynomials and {wk} are the weights for Gauss–Legendre quadrature. We truncate the spherical harmonics expansion (10) at l = Lmax. The number of azimuthal points N/ is conventionally chosen to N/ = 2(Lmax + 1). For the polar coordinate we use Nh = Lmax + 1. This particular choice for the grid allows achieving a good numerical accuracy with a moderate number of points. We use typically Nr = 80 radial points, Lmax = 19, so that the total number of the grid points is N = 2Nr(Lmax + 1)2 = 64,000, which is rather modest in comparison with (2Nr)3 necessary for an equivalent cubic grid with the same extension [7]. 3.3. Radial components: radial grid mapping Using the above expansion, the time propagation of each Rp is obtained independently by application of b p . We need therefore to construct the corresponding Q these evolution operators through the discretization b 0;l , which is done by means and diagonalization of H of the pseudospectral method as proposed by Chu et al. [23,29]. However, we introduce a new radial mapping, which allows for more flexibility than the mapping
362
B. Gervais et al. / Computational Materials Science 35 (2006) 359–365
proposed for Coulomb problems [21,23]. Our mapping is better suited for clusters and molecules. We shall map the radial domain [0, rmax], where rmax is the largest radius we use, onto the new domain [1, 1], so that the Legendre pseudospectral method can be applied [31]. Since we need a higher density of points inside the cluster where the potential is rapidly varying and where the electronic density is relatively large, we use the following mapping r(x), with x 2 [1, 1]: rðxÞ ¼ rmax by þ ð1 bÞy 2 ; ð15Þ 1 arcsinðcxÞ þ1 . y yðxÞ ¼ 2 arcsin c Various radial coordinates distributions obtained using the mapping function (15) with different set of parameters are presented in Fig. 1. Following [21] we introduce the new function S p ðxÞ ¼ Rp ðrðxÞÞ 1 . Using Sp, u and the variable x, , with uðxÞ ¼ pffiffiffiffiffiffi uðxÞ 0 r ðxÞ
the transformed problem reads: b 0;p S p b 0;p Rp ¼ 1 K H u 1 1 4 d2 u3 u00 lðl þ 1Þ u þ ¼ þ þ V 0 ðrÞ u2 S p ; u 2 dx2 2r2 2 u00 ¼
d2 u dx2
u2 S p ðxÞ
Nr X
gi ðxÞu2i sp ðxi Þ;
ð18Þ
i¼1
where gi(x) is the ith cardinal function for Gauss–Lobatto method [23], possessing the property gi(xj) = dij, xi is the ith root of the first derivative of the Legendre polynomial P 0N r ðxÞ, and Nr + 1 is the total number of points including the two extremities. The eigenproblem b 0;p S p ¼ S p is formulated on our grid as, using the K following notation ui = u(xi), Spi = Sp(xi), 3 00 Nr X 1 4 00 ui ui lðl þ 1Þ 2 þ ui gj ðxi Þuj spj þ þ V 0 ðri Þ u2i spi 2 2 2r 2 i j¼1 ¼ u2i spi ;
ð19Þ
which can be put into a symmetric form in a straightforward way. In Eq. (19), g00j ðxÞ is the second derivative of the cardinal function. For the sake of completeness and since the formulae given in [29] contains errors, we give here the explicit expression: P N r ðxi Þ ; i 6¼ j; P ðxi xj Þ N r ðxj Þ N r ðN r þ 1Þ ; i ¼ j. g00j xi ¼ ð1 x2i Þ g00j xi ¼
ð16Þ
so that the time propagation reads: b p Rp ¼ exp i H b 0;p dt Rp Q 1 b 0;p dt u2 S p ; ¼ exp i K ð17Þ u b has been which can be performed after the operator K discretized on our grid and diagonalized. According to the pseudospectral method, and for functions Sp that fulfill the boundary conditions Sp(x = 1) = Sp(x = 1) = 0, we expand Sp as
2
2
ð20Þ
With the smooth pseudopotentials used in this work, we can safely use a time step dt = 0.1 au. The method is fast enough to achieve time propagation up to a few ps, i.e., typically 106 iterations. Our experience shows that, using the method described here, we propagate the wave function with an error less 109 for the chosen dt = 0.1 and the grid parameters reported in Section 3.2. 3.4. Poisson equation In the course of the time propagation, we need to evaluate the Coulomb potential. This is done also with
60 1.4
grid mesh R(R) [a0]
=1 =1 =1 =0.5 =0.9 =0.7
R(xi) [a0]
40
20
=1 =1 =1 =0.5 =0.9 =0.7
1.2 1 0.8 0.6 0.4 0.2
0
–1
–0.5
0 xi
0.5
1
0
0
20
40
60
R [a0]
Fig. 1. Left: Points on the radial grid using the mapping function (15) as a function of the {xi} (i = 1, . . . , 79) roots of the Legendre polynomial of order 80. Right: Grid mesh DR at point Ri on the radial grid using the mapping function (15). The dashed line shows a uniform distribution, whereas the dotted line gives small meshes close to the center of the grid and larger at the border.
B. Gervais et al. / Computational Materials Science 35 (2006) 359–365
363
the help of the GPS method and the spherical harmonics expansion applied to the electronic density q(r, t) and to the Coulomb potential VCoul(r, t). This approach allows for a precise and fast computation of the Coulomb potential as well as for the kinetic energy operator evaluation. At a given time t we project the electronic density q(r, t) onto the spherical harmonics basis, obtain qlm(r, t) similarly to (11) and solve the corresponding Poisson equation for each (lm): 2 o lðl þ 1Þ þ ð21Þ rV 0lm ðr; tÞ ¼ 4prqlm ðr; tÞ r2 or2 with the zero boundary conditions rV 0lm ðr ¼ 0; tÞ ¼ rmax V 0lm ðr ¼ rmax ; tÞ ¼ 0. The complete Coulomb potential is thus obtained as X V lm ðr; tÞY lm ðXÞ; ð22Þ V ðr; tÞ ¼ lm
V lm ðr; tÞ ¼ V 0lm ðr; tÞ þ
rl ð2l þ 1Þr2lþ1 max
Z r0
dr0 qlm ðr0 ; tÞr0lþ2 ; max
ð23Þ where the integration is performed in a straightforward way by Gauss–Legendre integration associated with our grid. The relative accuracy of the radial second derivative is of the order of 109 like for the evaluation of the kinetic energy in the previous section, which ensures an accurate evaluation of the Coulomb potential. 4. Numerical example 4.1. Ground state relaxation As an example of application of our numerical algorithm, we present a study of a negatively charged sodium cluster Na 7 . Due to the large radial expansion of the first unoccupied KS wave functions, it is necessary, for a proper numerical representation, to consider a large radial grid. We take here a radial grid of rmax = 60a0 with the following grid parameters c = 0.9, b = 0.7. The ionic configuration of the ground state isomer, represented in the inset of Fig. 2, consists of a flat pentagonal ring with two extra ions symmetrically placed on the C5 axis. The equilibrium ionic geometry of the cluster was obtained from damped ab initio molecular dynamics by solving equations of motion (3) and (4) with the Beeman algorithm [32]. The ionic motion was damped during relaxation by means of a small extra friction force experienced by each atom. The electronic ground state configuration is found by standard imaginary time propagation associated with damping of the electronic kinetic energy [33], self-consistently with relaxation of the ionic subsystem. We stop the minimum energy search when relative variation of the electronic density and energy is less than 106 and all the forces
Fig. 2. Radial electronic density of the highest occupied and 3 lowest unoccupied KS orbitals of Na 7 . The bending of the third unoccupied orbital close to the border of the grid is due to the size effect, because the KS orbital is constraint to be 0 at the border. The inset represents ground state ionic configuration of Na 7.
in the system become less than 107, what in our case corresponds to a final ionic temperature lower than 1 K. With respect to other methods, this relaxation scheme has an advantage to be fully consistent, and does not require any further assumption. 4.2. Physical observables Real-time propagation of the system is an elegant way to obtain the dipole response of the system S(x) [7]. For a given ground state obtained using the method described above, we perturb the electronic system by applying a small boost q, so that w0a ðrÞ ¼ expðiqrÞwa ðrÞ. Applying the same boost to all single-electron wave functions causes a dipole moment d(t) to develop in _ time. We record the time derivative of the dipole dðtÞ _ and simply take its Fourier transform dðxÞ, which is the response we seek. The absorption energies of the system correspond to the peaks in SðxÞ ¼ 23 R hP i d_ m ðxÞ=pqm . Since we use a finite time step, the m¼x;y;z
peaks have a small finite width, and the dipole oscillator strength of a given transition n 0 is deduced from S(x) by the following relation: Z fn0 ¼ dx SðxÞ; ð24Þ where the integral is taken over a small interval around the energy of the transition. To prove the effectiveness of our numerical algorithm for a large time scale calculations, the monitoring of the dipole moment d(t) was done by real-time evolution of the electron wave functions during the interval of 2.5 ps with time step dt = 0.1. The calculations was performed with ionic motion taken into account and using absorbing boundary conditions at r = rmax [7,8].
364
B. Gervais et al. / Computational Materials Science 35 (2006) 359–365 1 -
Na7 Na2
0.9 0.8
Nesc
0.7 0.6 0.5 0.4 0.3 0.2
0
50
100
150
200
250
300
350
Rmax [a0]
Fig. 3. The dipole spectrum of the Na 7 cluster excited by the small boost q, with qx = 0, qy = qz = 0.001. The full and dashed lines differ in the grid size used for the calculation.
In Fig. 3 we present the dipole spectra obtained using an initial small boost excitation qx = 0, qy = qz q = 0.001, what corresponds to an energy supplied to the system of about 3.4 leV, i.e., to a purely linear regime of excitation. The discrete lines below the ionization potential (marked on the plot as IP) are slightly broadened due to the cluster ionic motion. We note that the spectral structure above the ionization potential has no direct physical meaning and results from the finite size of the radial cell we use. To illustrate the effect of the cell radius Rmax on the result of the calculations we present the dipole spectra obtained by our method for two different values of Rmax. One can see that both calculations give practically the same results for the two main discrete peaks that lie below the IP, while the structure of the pseudo-continual part of the spectrum is sensitive to the cell size. It is important to note that the dipole spectra S(x) obtained reproduce the sum rule S(x) dx = Ne [10] with good accuracy, what is the criterion of sufficient completeness of our basis we keep throughout all the time domain, despite a small electronic charge loss (about 107) due to the absorption of the electron density at the border. In Fig. 4 we plot the total emitted electron charge, defined as the difference between the initial charge and the final charge (here after 500 fs), as a function of the grid size Rmax. Test cases are Na 7 and Na2 excited by a femtosecond laser pulse, tuned to an eigenfrequency of the system, which does not depend on the grid size. The parameters of the laser pulse are defined in the caption of Fig. 4. We obtain that the total emitted charge decreases strongly with increasing grid size. For Na2 we note that the emission reaches its asymptotic value for Rmax = 300a0 whereas for Na 7 the asymptotic value is still not attained. This shows that very large grids are
Fig. 4. Total emitted electronic charge of Na 7 and Na2 after a short laser pulse as a function of the grid size Rmax. We used the following laser parameters (FWHM [fs], intensity [W/cm2] and frequency [eV]) = (56, 5.67 · 1011, 0.727) and (28, 7.9 · 1010, 2.06) for Na 7 and Na2, respectively.
necessary to handle properly the electronic emission, especially for negatively charged clusters.
5. Conclusions A numerical method to study small atomic systems in full 3D coordinate space and real time was presented. Our approach is based on the TDLDA and MD methods. The non-adiabatic propagation of three-dimensional wave functions is performed using the GPS method combined with spherical harmonic expansion to solve the time-dependent Kohn–Sham equations. The radial mapping is used in the course of the time propagation to evaluate both the Coulomb potential and the kinetic energy operator with high accuracy. The dipole spectra of negatively charged sodium clusters Na 7 are calculated, both with frozen ionic core and with full real-time ionic dynamics. We checked the effect of the grid size on physical observables like electronic emission and dipole spectra, and showed that very large grids are necessary to study excited negatively charged clusters.
References [1] W.A. de Heer, Rev. Mod. Phys. 65 (1993) 611. [2] M. Brack, Rev. Mod. Phys. 65 (1993) 677. [3] H. Haberland (Ed.), Clusters of Atoms and Molecules, Springer Series in Chemical Physics, vol. 52, Springer, Berlin, Heidelberg, New York, 1994. [4] U. Kreibig, M. Vollmer, Optical Properties of Metal Clusters, Springer-Verlag, Berlin, Heidelberg, 1995. [5] W. Ekardt (Ed.), Metal Clusters, Wiley, New York, 1999. [6] V.K. Ivanov, A.N. Ipatov, in: J.P. Connerade (Ed.), Correlations in Clusters and Related Systems, World Scientific, Singapore, 1996.
B. Gervais et al. / Computational Materials Science 35 (2006) 359–365 [7] F. Calvayrac, P.-G. Reinhard, E. Suraud, C.A. Ullrich, Phys. Rep. 337 (2000) 493. [8] F. Calvayrac, P.-G. Reinhard, E. Suraud, Ann. Phys. (NY) 255 (1997) 125. [9] K. Yabana, G.F. Bertsch, Phys. Rev. B 54 (1996) 4484. [10] G.F. Bertsch, R.A. Broglia, Oscillations in Finite Quantum Systems, Cambridge University Press, Cambridge, 1994. [11] R. Schlipper, R. Kusche, B.V. Issendorff, H. Haberland, Phys. Rev. Lett. 80 (1998) 1194. [12] R. Schlipper, R. Kusche, B.V. Issendorff, H. Haberland, Appl. Phys. A 72 (2001) 255. [13] T. Ditmire et al., Phys. Rev. Lett. 78 (1997) 2732. [14] F. Catara, Ph. Chomaz, N. Van Giai, Phys. Rev. B 48 (1993) 18207. [15] K. Hagino, Phys. Rev. B 60 (1999) R2197. [16] G.B. Bertsch, N. Van Giai, N. Vinh Mau, Phys. Rev. A 61 (2000) 33202. [17] L. Gerchikov, C. Guet, A. Ipatov, Phys. Rev. A 66 (2003) 53202. [18] L. Gerchikov, A. Ipatov, J. Phys. B 36 (2003) 1193. [19] C. Kohl, E. Suraud, P.-G. Reinhard, Eur. Phys. J. D 11 (2000) 115. [20] C. Legrand, E. Suraud, P.-G. Reinhard, J. Phys. B. 35 (2002) 1115.
[21] [22] [23] [24] [25] [26] [27]
[28]
[29] [30] [31]
[32] [33]
365
G. Yao, S.I. Chu, Chem. Phys. Lett. 197 (1992) 413. J. Wang, S.I. Chu, Chem. Phys. Lett. 227 (1994) 663. X.M. Tong, S.-I. Chu, Phys. Rev. A 57 (1998) 452. X. Chu, S.-I. Chu, Phys. Rev. A 63 (2001) 23411. G. Lauritsch, P.-G. Reinhard, Int. J. Mod. Phys. C 5 (1994) 65. S. Goedecker, M. Teter, J. Hutter, Phys. Rev. B 54 (1998) 1703. D. Gottlieb, S.A. Orszag, Numerical Analysis of Spectral Methods: Theory and Applications, SIAM Publications, Philadelphia, 1977. R.M. Dreizler, E.U.K. Gross, Density Functional Theory: An Approach to the Quantum Many-body Problem, Springer, Berlin, 1990. J. Wang, S.-I. Chu, C. Laughlin, Phys. Rev. A 50 (1994) 3208. W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) 1133. D. Gottlieb, M. Yousuff, S.A. Orszag, in: R. Voigt, D. Gottlieb, M.Y. Hussaini (Eds.), Spectral Methods for Spatial Differential Equations, SIAM, Philadelphia, 1984. D. Frenkel, B. Smit, Understanding Molecular Simulation, Academic Press, London, 1996. V. Blum, G. Lauritsch, J.A. Maruhn, P.-G. Reinhard, J. Comp. Phys. 100 (2001) 364.
Effective numerical method for theoretical studies of small atomic clusters B. Gervais a, E. Giglio
a,*
, A. Ipatov
a,b,1
, J. Douady
a
a
b
CIRIL, unite´ mixte CEA-CNRS-ENSICAEN, BP 5133, F-14070 Caen Cedex 05, France St. Petersburg Polytechnical University, Polytechnicheskaja 29, 195251 St. Petersburg, Russian Federation Received 24 April 2004; accepted 27 July 2004
Abstract We present an effective numerical method to study the electronic and ionic dynamics of small atomic systems, like alkali-metal clusters. Our approach is based on the density functional theory (DFT) combined with molecular dynamics (MD). The time-dependent Kohn–Sham (KS) equations describing the electronic subsystem are solved self-consistently in the coordinate space. Their numerical solution is performed using the time-splitting method on a three-dimensional coordinate grid by means of the generalized pseudospectral method (GPS). The KS orbitals are expanded over a spherical harmonics basis, and their radial part are represented on a non-uniform radial grid. We employ an optimal radial mapping of the collocation points consistent with the GPS method, which allows achieving a higher density of radial points inside the cluster where electronic density is larger and the potential varies rapidly. This method allows evaluating both the Coulomb potential and the kinetic energy operator with high accuracy. Unlike FFT based methods, the GPS method with non-uniform radial mapping, allows to use large grids at lower cost and thus to study efficiently negatively charged cluster. As an example, we present the number of emitted electrons of a negatively charged sodium cluster Na 7 , irradiated by a femtosecond laser pulse, as a function of the grid size. We show that the value of the electron emission converges only for large grids. 2005 Elsevier B.V. All rights reserved. PACS: 31.15.Ev; 31.15.p; 31.15.Qg; 36.40.c Keywords: Local density approximation; Time splitting; Harmonic expansion; Grid mapping; Spectral methods; 3D TDDFT; Molecular dynamics
1. Introduction The electronic and structural properties of small clusters, in particular alkali clusters, have been the subject of numerous recent investigations, from both experimental and theoretical point of view (for overview see [1–6]). The theoretical methods employing the real-time approach [7–9] are known as a powerful tool to study the optical properties of finite many-electron systems,
*
1
Corresponding author. Fax: +33 2 3145 4714. E-mail address: [email protected] (E. Giglio). Acknowledgments to CEA and CNRS for financial support.
0927-0256/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2004.07.015
their excitation processes, photoabsorption and electron emission. For instance, the surface dipole plasmon resonance is one of the most important collective excitations observed in metallic clusters [3,4,10]. The recent development of laser technologies and the use of femtosecond lasers has led to experimental studies [11–13] of clusters in strong laser field and stimulated the development of theoretical methods that allow for an adequate description of optical response beyond the linear regime [7, 14–19]. The real-time methods describe the dynamical properties in highly excited clusters in a natural way and can be effectively employed in studies of clusters ionization, thermal relaxation and fission for free, deposited and embedded systems.
360
B. Gervais et al. / Computational Materials Science 35 (2006) 359–365
The main goal of this work is to present a numerical method based on the real-time technique that possesses the simplicity of the direct discretization methods with high accuracy and computational effectiveness. Our method is based on the density functional theory (DFT) combined with molecular dynamics (MD). Assuming that a proper theoretical description of the system requires to account for the coupled dynamics of both electronic and ionic degrees of freedom, we treat the real-time motion of the electronic system using the time-dependent local density approximation (TDLDA) [7,9] method with average-density self-interaction correction ADSIC [20] in full coordinate-space grid, while the ionic motion of the cluster core is described using the standard methods of MD. We evaluate the cluster electronic system in real time, that allows to go beyond the adiabatic description of the system dynamics [4,5] and to consider, for instance, relatively long time scales, i.e., typically from a few hundreds fs to a few ps, accounting for ionic motion. We generalize the GPS Method introduced by Chu et al. mainly for atomic systems [21–24]. The latter approach was developed for atomic systems with cylindrical symmetry that naturally reduces to 2D problems. We employ for the first time this method for a full 3D computation and for atomic clusters. It allows to investigate the dynamics of such systems in providing self-consistent solution of the time-dependent Kohn–Sham equations. The proposed approach allows to improve a lot the efficiency of the numerical algorithms in comparison with alternative methods based on the FFT approach [25,8]. Especially for small and middle range clusters, for which larger grids can be used and longer time scales computed. Atomic system of units (h = e = me = 1) is used throughout the paper.
2. Theoretical approach The valence electrons of the Na cluster are described using the density functional theory methods [28]. The single-electron wave functions wa(r, t) follow the timedependent Kohn–Sham [30] equations d D i wa ðr; tÞ ¼ þ V KS ½r; q" ðr; tÞ; q# ðr; tÞ wa ðr; tÞ; dt 2 ð1Þ where " and # stands for the spin of electrons, with the corresponding electronic density defined as qr ðr; tÞ ¼ PN r 2 jw ðr; tÞj where r = ", # and N" and N# are the i i number of spin up and spin down electrons, respectively. The total density is q(r, t) = q"(r, t) + q#(r, t), and the number of electrons Ne = N" + N#. For such a system, the Kohn–Sham potential for spin up VKS" reads: VKS"(r, t) = Vext(r, t) + VCoul(r, t) +
Vxc"(r, t), where the external potential Vext takes into account the interaction with external field and ionic cores. The Coulomb potential VCoul is determined by solving the Poisson equation without retardation. Within ADSIC approximation [20], the exchange-correlation potential Vxc"(r, t) reads: Z 1 qðr0 Þ dr0 þ V LSD V xc" ¼ xc" ½q" ðr; tÞ; q# ðr; tÞ N jr r0 j q" ðr; tÞ LSD ;0 ; ð2Þ V xc" N" where V LSD xc" is a local density exchange and correlation potential for spin up. Similar equations are derived for spin down. The interaction among the ions is modeled as a simple Coulomb interaction of point charges. The molecular dynamics equations of ionic motion reduces to the following: " Z d Pi ¼ rRi V ext ðRi ; tÞ qðr; tÞV ps ðjRi rjÞ dr dt # X 1 þ ; ð3Þ jRi Rj j j6¼i d Pi Ri ¼ ; dt M ion
ð4Þ
where Mion is the ionic mass, Ri and Pi the position and momentum of the ion respectively.
3. Computational approach 3.1. Time propagation of the KS orbitals The computation of the Kohn–Sham orbitals is performed on a 3D spherical grid. We define va(r, t) = rwa(r, t), so that the KS equation (1) reads: " # 2 b L o 1 d2 þ þ Vb r ðtÞ va ðr; tÞ; ð5Þ i va ðr; tÞ ¼ ot 2 dr2 2r2 where Vb r ðtÞ is the time-dependent Kohn–Sham potential for spin r. We use a second order time-splitting algorithm [7] for the time propagation of the KS orbitals: bU b 2 ðtÞfa ðr; tÞ þ Oðdt3 Þ; fa ðr; t þ dtÞ ¼ Q r
ð6Þ
to the KS orbitals by the relation where fa is connected b y ðr; tÞva ðr; tÞ, with fr ðr; tÞ ¼ U r dt b b U r ðr; tÞ ¼ exp ið V r ðr; tÞ V 0 ðrÞÞ ; ð7Þ 2 b ¼ exp i H b 0 dt ; Q ð8Þ
B. Gervais et al. / Computational Materials Science 35 (2006) 359–365 2
b 0 ¼ 1 d22 þ bL 2 þ V 0 ðrÞ and V0(r) is a spheriwhere H 2 dr 2r cally symmetric potential. It is usually chosen as the spherical average of the total potential at the initial time V(r, t = 0). However, it can be updated to account for possible significant modifications of the potential in the course of the system dynamics. This partitioning of the potential part of the propagator between the timedependent and time-independent operators permits to reduce the error resulting from the time-splitting b does not depend on time t so that scheme. Note that Q it can be computed once for each given V0(r). The b numerical evaluation of the action of the operator Q on va(r, t) depends on the radial grid choice and is detailed in the next sections. For local potentials, the numerical evaluation of b r ðtÞva ðr; tÞ is straightforward. Indeed, for any given U set of collocation points {rg} it requires only the evalua b r ðtÞva ðr; tÞ tion the potential V(rg, t) to obtain U ¼ r¼rg b r ðrg ; tÞva ðrg ; tÞ. For non-local potentials of integral U form, the method is applicable as well. Indeed, the action of the operator Vb r ðtÞ on va(r, t) can be determined accurately by Gauss quadrature: Z Vb r ðtÞva ðtÞ r¼rg ¼ dr0 V r ðrg ; r0 ; tÞva ðr0 ; tÞ X wg0 V r ðrg ; rg0 ; tÞva ðrg0 ; tÞ. ð9Þ ¼ g0
The algebraic equation (9) can be easily put into a symmetric form and then diagonalized to evaluate the corresponding exponential matrix U r ðrg ; rg0 ; tÞ. The method is efficient provided the non-local part extent over a limited set of grid points so that: (i) the matrix V r ðrg ; rg0 Þ is sparse and (ii) it is block diagonal. This is usually the case for standard pseudopotentials like for example those proposed by Goedecker et al. [26] and for interatomic distances encountered in normal situation where the ions move in a limited range around their equilibrium positions, so that the bond length does not become excessively short. 3.2. Spherical harmonic expansion The grid points we use are distributed on concentric spheres of radius ri at solid angles Xq = (ha, /b). The N/ azimuthal angle coordinate points /b are uniformly spaced in the interval [0, 2p] so that /b ¼ N2p/ . The Nh polar angles ha are distributed non-uniformly with ya = cos(ha) taken as the roots of the Legendre polynomials P N h ðyÞ. Since we solve the 3D problem on a spherical grid, it is natural to split the operator into the radial and the angular parts, using the expansion of the full 3D functions va(r, X) into spherical harmonics Yp(h, /), with p l, m. Moreover, the use of spherical harmonics allows eliminating the so-called pole problem [27]. We define the radial functions Rap as the projection over the spherical harmonics Yp:
va ðr; X; tÞ ¼ Rap ðr; tÞ
¼
Z
X
361
Rap ðr; tÞY p ðh; /Þ;
ð10Þ
p
ðva ðr; XÞÞY p ðh; /Þ dX.
ð11Þ
Within this approach, the operator (8) now reads: X Y p Q b p Y p ; b ¼ Q ð12Þ p
b p i ¼ expði H b p ¼ hY p j QjY b 0;l dtÞ, is a sum of rawhere Q dial operators acting in 1-dimensional-subspace, for which 2 b 0;l ¼ 1 d þ lðl þ 1Þ þ V 0 ðrÞ. H 2 dr2 2r2
ð13Þ
The projections (11) are performed efficiently at each time step using a method based on the FFT algorithm for integration over the azimuthal angle /, and Gauss–Legendre quadrature for integration over the azimuthal angle h: Z nm ðri ; hj Þ ¼ vðri ; hj ; /Þ expðim/Þ d/ N
/ 2p X vðri ; hj ; /k Þ expðim/k Þ; N / k¼1 Z jmj Rlm ðri Þ ¼ C lm nm ðri ; hÞP l ðcos hÞ d cos h
¼
¼ C lm
Nh X
ð14Þ
jmj
wj nm ðri ; hj ÞP l ðcos hj Þ;
j¼1
where Clm is the normalization constant of the spherical harmonics, P ml are the associated Legendre polynomials and {wk} are the weights for Gauss–Legendre quadrature. We truncate the spherical harmonics expansion (10) at l = Lmax. The number of azimuthal points N/ is conventionally chosen to N/ = 2(Lmax + 1). For the polar coordinate we use Nh = Lmax + 1. This particular choice for the grid allows achieving a good numerical accuracy with a moderate number of points. We use typically Nr = 80 radial points, Lmax = 19, so that the total number of the grid points is N = 2Nr(Lmax + 1)2 = 64,000, which is rather modest in comparison with (2Nr)3 necessary for an equivalent cubic grid with the same extension [7]. 3.3. Radial components: radial grid mapping Using the above expansion, the time propagation of each Rp is obtained independently by application of b p . We need therefore to construct the corresponding Q these evolution operators through the discretization b 0;l , which is done by means and diagonalization of H of the pseudospectral method as proposed by Chu et al. [23,29]. However, we introduce a new radial mapping, which allows for more flexibility than the mapping
362
B. Gervais et al. / Computational Materials Science 35 (2006) 359–365
proposed for Coulomb problems [21,23]. Our mapping is better suited for clusters and molecules. We shall map the radial domain [0, rmax], where rmax is the largest radius we use, onto the new domain [1, 1], so that the Legendre pseudospectral method can be applied [31]. Since we need a higher density of points inside the cluster where the potential is rapidly varying and where the electronic density is relatively large, we use the following mapping r(x), with x 2 [1, 1]: rðxÞ ¼ rmax by þ ð1 bÞy 2 ; ð15Þ 1 arcsinðcxÞ þ1 . y yðxÞ ¼ 2 arcsin c Various radial coordinates distributions obtained using the mapping function (15) with different set of parameters are presented in Fig. 1. Following [21] we introduce the new function S p ðxÞ ¼ Rp ðrðxÞÞ 1 . Using Sp, u and the variable x, , with uðxÞ ¼ pffiffiffiffiffiffi uðxÞ 0 r ðxÞ
the transformed problem reads: b 0;p S p b 0;p Rp ¼ 1 K H u 1 1 4 d2 u3 u00 lðl þ 1Þ u þ ¼ þ þ V 0 ðrÞ u2 S p ; u 2 dx2 2r2 2 u00 ¼
d2 u dx2
u2 S p ðxÞ
Nr X
gi ðxÞu2i sp ðxi Þ;
ð18Þ
i¼1
where gi(x) is the ith cardinal function for Gauss–Lobatto method [23], possessing the property gi(xj) = dij, xi is the ith root of the first derivative of the Legendre polynomial P 0N r ðxÞ, and Nr + 1 is the total number of points including the two extremities. The eigenproblem b 0;p S p ¼ S p is formulated on our grid as, using the K following notation ui = u(xi), Spi = Sp(xi), 3 00 Nr X 1 4 00 ui ui lðl þ 1Þ 2 þ ui gj ðxi Þuj spj þ þ V 0 ðri Þ u2i spi 2 2 2r 2 i j¼1 ¼ u2i spi ;
ð19Þ
which can be put into a symmetric form in a straightforward way. In Eq. (19), g00j ðxÞ is the second derivative of the cardinal function. For the sake of completeness and since the formulae given in [29] contains errors, we give here the explicit expression: P N r ðxi Þ ; i 6¼ j; P ðxi xj Þ N r ðxj Þ N r ðN r þ 1Þ ; i ¼ j. g00j xi ¼ ð1 x2i Þ g00j xi ¼
ð16Þ
so that the time propagation reads: b p Rp ¼ exp i H b 0;p dt Rp Q 1 b 0;p dt u2 S p ; ¼ exp i K ð17Þ u b has been which can be performed after the operator K discretized on our grid and diagonalized. According to the pseudospectral method, and for functions Sp that fulfill the boundary conditions Sp(x = 1) = Sp(x = 1) = 0, we expand Sp as
2
2
ð20Þ
With the smooth pseudopotentials used in this work, we can safely use a time step dt = 0.1 au. The method is fast enough to achieve time propagation up to a few ps, i.e., typically 106 iterations. Our experience shows that, using the method described here, we propagate the wave function with an error less 109 for the chosen dt = 0.1 and the grid parameters reported in Section 3.2. 3.4. Poisson equation In the course of the time propagation, we need to evaluate the Coulomb potential. This is done also with
60 1.4
grid mesh R(R) [a0]
=1 =1 =1 =0.5 =0.9 =0.7
R(xi) [a0]
40
20
=1 =1 =1 =0.5 =0.9 =0.7
1.2 1 0.8 0.6 0.4 0.2
0
–1
–0.5
0 xi
0.5
1
0
0
20
40
60
R [a0]
Fig. 1. Left: Points on the radial grid using the mapping function (15) as a function of the {xi} (i = 1, . . . , 79) roots of the Legendre polynomial of order 80. Right: Grid mesh DR at point Ri on the radial grid using the mapping function (15). The dashed line shows a uniform distribution, whereas the dotted line gives small meshes close to the center of the grid and larger at the border.
B. Gervais et al. / Computational Materials Science 35 (2006) 359–365
363
the help of the GPS method and the spherical harmonics expansion applied to the electronic density q(r, t) and to the Coulomb potential VCoul(r, t). This approach allows for a precise and fast computation of the Coulomb potential as well as for the kinetic energy operator evaluation. At a given time t we project the electronic density q(r, t) onto the spherical harmonics basis, obtain qlm(r, t) similarly to (11) and solve the corresponding Poisson equation for each (lm): 2 o lðl þ 1Þ þ ð21Þ rV 0lm ðr; tÞ ¼ 4prqlm ðr; tÞ r2 or2 with the zero boundary conditions rV 0lm ðr ¼ 0; tÞ ¼ rmax V 0lm ðr ¼ rmax ; tÞ ¼ 0. The complete Coulomb potential is thus obtained as X V lm ðr; tÞY lm ðXÞ; ð22Þ V ðr; tÞ ¼ lm
V lm ðr; tÞ ¼ V 0lm ðr; tÞ þ
rl ð2l þ 1Þr2lþ1 max
Z r0
dr0 qlm ðr0 ; tÞr0lþ2 ; max
ð23Þ where the integration is performed in a straightforward way by Gauss–Legendre integration associated with our grid. The relative accuracy of the radial second derivative is of the order of 109 like for the evaluation of the kinetic energy in the previous section, which ensures an accurate evaluation of the Coulomb potential. 4. Numerical example 4.1. Ground state relaxation As an example of application of our numerical algorithm, we present a study of a negatively charged sodium cluster Na 7 . Due to the large radial expansion of the first unoccupied KS wave functions, it is necessary, for a proper numerical representation, to consider a large radial grid. We take here a radial grid of rmax = 60a0 with the following grid parameters c = 0.9, b = 0.7. The ionic configuration of the ground state isomer, represented in the inset of Fig. 2, consists of a flat pentagonal ring with two extra ions symmetrically placed on the C5 axis. The equilibrium ionic geometry of the cluster was obtained from damped ab initio molecular dynamics by solving equations of motion (3) and (4) with the Beeman algorithm [32]. The ionic motion was damped during relaxation by means of a small extra friction force experienced by each atom. The electronic ground state configuration is found by standard imaginary time propagation associated with damping of the electronic kinetic energy [33], self-consistently with relaxation of the ionic subsystem. We stop the minimum energy search when relative variation of the electronic density and energy is less than 106 and all the forces
Fig. 2. Radial electronic density of the highest occupied and 3 lowest unoccupied KS orbitals of Na 7 . The bending of the third unoccupied orbital close to the border of the grid is due to the size effect, because the KS orbital is constraint to be 0 at the border. The inset represents ground state ionic configuration of Na 7.
in the system become less than 107, what in our case corresponds to a final ionic temperature lower than 1 K. With respect to other methods, this relaxation scheme has an advantage to be fully consistent, and does not require any further assumption. 4.2. Physical observables Real-time propagation of the system is an elegant way to obtain the dipole response of the system S(x) [7]. For a given ground state obtained using the method described above, we perturb the electronic system by applying a small boost q, so that w0a ðrÞ ¼ expðiqrÞwa ðrÞ. Applying the same boost to all single-electron wave functions causes a dipole moment d(t) to develop in _ time. We record the time derivative of the dipole dðtÞ _ and simply take its Fourier transform dðxÞ, which is the response we seek. The absorption energies of the system correspond to the peaks in SðxÞ ¼ 23 R hP i d_ m ðxÞ=pqm . Since we use a finite time step, the m¼x;y;z
peaks have a small finite width, and the dipole oscillator strength of a given transition n 0 is deduced from S(x) by the following relation: Z fn0 ¼ dx SðxÞ; ð24Þ where the integral is taken over a small interval around the energy of the transition. To prove the effectiveness of our numerical algorithm for a large time scale calculations, the monitoring of the dipole moment d(t) was done by real-time evolution of the electron wave functions during the interval of 2.5 ps with time step dt = 0.1. The calculations was performed with ionic motion taken into account and using absorbing boundary conditions at r = rmax [7,8].
364
B. Gervais et al. / Computational Materials Science 35 (2006) 359–365 1 -
Na7 Na2
0.9 0.8
Nesc
0.7 0.6 0.5 0.4 0.3 0.2
0
50
100
150
200
250
300
350
Rmax [a0]
Fig. 3. The dipole spectrum of the Na 7 cluster excited by the small boost q, with qx = 0, qy = qz = 0.001. The full and dashed lines differ in the grid size used for the calculation.
In Fig. 3 we present the dipole spectra obtained using an initial small boost excitation qx = 0, qy = qz q = 0.001, what corresponds to an energy supplied to the system of about 3.4 leV, i.e., to a purely linear regime of excitation. The discrete lines below the ionization potential (marked on the plot as IP) are slightly broadened due to the cluster ionic motion. We note that the spectral structure above the ionization potential has no direct physical meaning and results from the finite size of the radial cell we use. To illustrate the effect of the cell radius Rmax on the result of the calculations we present the dipole spectra obtained by our method for two different values of Rmax. One can see that both calculations give practically the same results for the two main discrete peaks that lie below the IP, while the structure of the pseudo-continual part of the spectrum is sensitive to the cell size. It is important to note that the dipole spectra S(x) obtained reproduce the sum rule S(x) dx = Ne [10] with good accuracy, what is the criterion of sufficient completeness of our basis we keep throughout all the time domain, despite a small electronic charge loss (about 107) due to the absorption of the electron density at the border. In Fig. 4 we plot the total emitted electron charge, defined as the difference between the initial charge and the final charge (here after 500 fs), as a function of the grid size Rmax. Test cases are Na 7 and Na2 excited by a femtosecond laser pulse, tuned to an eigenfrequency of the system, which does not depend on the grid size. The parameters of the laser pulse are defined in the caption of Fig. 4. We obtain that the total emitted charge decreases strongly with increasing grid size. For Na2 we note that the emission reaches its asymptotic value for Rmax = 300a0 whereas for Na 7 the asymptotic value is still not attained. This shows that very large grids are
Fig. 4. Total emitted electronic charge of Na 7 and Na2 after a short laser pulse as a function of the grid size Rmax. We used the following laser parameters (FWHM [fs], intensity [W/cm2] and frequency [eV]) = (56, 5.67 · 1011, 0.727) and (28, 7.9 · 1010, 2.06) for Na 7 and Na2, respectively.
necessary to handle properly the electronic emission, especially for negatively charged clusters.
5. Conclusions A numerical method to study small atomic systems in full 3D coordinate space and real time was presented. Our approach is based on the TDLDA and MD methods. The non-adiabatic propagation of three-dimensional wave functions is performed using the GPS method combined with spherical harmonic expansion to solve the time-dependent Kohn–Sham equations. The radial mapping is used in the course of the time propagation to evaluate both the Coulomb potential and the kinetic energy operator with high accuracy. The dipole spectra of negatively charged sodium clusters Na 7 are calculated, both with frozen ionic core and with full real-time ionic dynamics. We checked the effect of the grid size on physical observables like electronic emission and dipole spectra, and showed that very large grids are necessary to study excited negatively charged clusters.
References [1] W.A. de Heer, Rev. Mod. Phys. 65 (1993) 611. [2] M. Brack, Rev. Mod. Phys. 65 (1993) 677. [3] H. Haberland (Ed.), Clusters of Atoms and Molecules, Springer Series in Chemical Physics, vol. 52, Springer, Berlin, Heidelberg, New York, 1994. [4] U. Kreibig, M. Vollmer, Optical Properties of Metal Clusters, Springer-Verlag, Berlin, Heidelberg, 1995. [5] W. Ekardt (Ed.), Metal Clusters, Wiley, New York, 1999. [6] V.K. Ivanov, A.N. Ipatov, in: J.P. Connerade (Ed.), Correlations in Clusters and Related Systems, World Scientific, Singapore, 1996.
B. Gervais et al. / Computational Materials Science 35 (2006) 359–365 [7] F. Calvayrac, P.-G. Reinhard, E. Suraud, C.A. Ullrich, Phys. Rep. 337 (2000) 493. [8] F. Calvayrac, P.-G. Reinhard, E. Suraud, Ann. Phys. (NY) 255 (1997) 125. [9] K. Yabana, G.F. Bertsch, Phys. Rev. B 54 (1996) 4484. [10] G.F. Bertsch, R.A. Broglia, Oscillations in Finite Quantum Systems, Cambridge University Press, Cambridge, 1994. [11] R. Schlipper, R. Kusche, B.V. Issendorff, H. Haberland, Phys. Rev. Lett. 80 (1998) 1194. [12] R. Schlipper, R. Kusche, B.V. Issendorff, H. Haberland, Appl. Phys. A 72 (2001) 255. [13] T. Ditmire et al., Phys. Rev. Lett. 78 (1997) 2732. [14] F. Catara, Ph. Chomaz, N. Van Giai, Phys. Rev. B 48 (1993) 18207. [15] K. Hagino, Phys. Rev. B 60 (1999) R2197. [16] G.B. Bertsch, N. Van Giai, N. Vinh Mau, Phys. Rev. A 61 (2000) 33202. [17] L. Gerchikov, C. Guet, A. Ipatov, Phys. Rev. A 66 (2003) 53202. [18] L. Gerchikov, A. Ipatov, J. Phys. B 36 (2003) 1193. [19] C. Kohl, E. Suraud, P.-G. Reinhard, Eur. Phys. J. D 11 (2000) 115. [20] C. Legrand, E. Suraud, P.-G. Reinhard, J. Phys. B. 35 (2002) 1115.
[21] [22] [23] [24] [25] [26] [27]
[28]
[29] [30] [31]
[32] [33]
365
G. Yao, S.I. Chu, Chem. Phys. Lett. 197 (1992) 413. J. Wang, S.I. Chu, Chem. Phys. Lett. 227 (1994) 663. X.M. Tong, S.-I. Chu, Phys. Rev. A 57 (1998) 452. X. Chu, S.-I. Chu, Phys. Rev. A 63 (2001) 23411. G. Lauritsch, P.-G. Reinhard, Int. J. Mod. Phys. C 5 (1994) 65. S. Goedecker, M. Teter, J. Hutter, Phys. Rev. B 54 (1998) 1703. D. Gottlieb, S.A. Orszag, Numerical Analysis of Spectral Methods: Theory and Applications, SIAM Publications, Philadelphia, 1977. R.M. Dreizler, E.U.K. Gross, Density Functional Theory: An Approach to the Quantum Many-body Problem, Springer, Berlin, 1990. J. Wang, S.-I. Chu, C. Laughlin, Phys. Rev. A 50 (1994) 3208. W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) 1133. D. Gottlieb, M. Yousuff, S.A. Orszag, in: R. Voigt, D. Gottlieb, M.Y. Hussaini (Eds.), Spectral Methods for Spatial Differential Equations, SIAM, Philadelphia, 1984. D. Frenkel, B. Smit, Understanding Molecular Simulation, Academic Press, London, 1996. V. Blum, G. Lauritsch, J.A. Maruhn, P.-G. Reinhard, J. Comp. Phys. 100 (2001) 364.