XANES calculation
Contents |
Abstract
We have implemented an original procedure for calculating XANES spectra within the BigDFT code. Our approach consists firstly in projecting the photoelectron wave-function onto the pseudopotential functions basis with the help of a reversed PAW projector^{[1]} and, then, propagating this initial state by iterative applications of the Hamiltonian. The Hamiltonian is the one provided by the BigDFT code. The obtained spectra therefore correspond exactly to the electronic structure of the self consistent ground state.The spectra are built by obtaining at each new application of the Hamiltonian a new coefficient of the spectral decomposition in Chebychev polynomials. The method therefore systematically improves the calculated spectra resolution at each step and can be stopped once the required resolution is reached. Moreover, the Chebychev polynomials have a higher nodes’ density at the extrema of the spectral range than in the middle. For practical applications this makes convergence even faster because the broadening due to lifetime is narrow near the edge and large elsewhere. The execution times represent a break-through for this kind of calculations and reflect the exceptional features of the BigDFT code. For a cluster contained in a 25Å side cube, the spectra are obtained within half an hour, running on a single processor.
Introduction
The interaction between matter and a photon described by a wave vector Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): { k}
and polarisation Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): { \epsilon}
, is written, discarding elastic Thomson scattering (the Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): A^2(r)
term in the Schr{ö}dinger equation which becomes dominant off-resonance), discarding also the spin-magnetic field interaction^{[2]}, and using CGS units, as : Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): H_{int} = \left( \displaystyle \frac{2 \pi \hbar c^2 }{ \omega V}\right)^{1/2} ({ a^\dagger }_{ k,\epsilon} exp(-i {\textbf{ k} \cdot \textbf{ r}} ) + a_{k ,\epsilon} exp(i { \textbf{ k} \cdot \textbf{ r}} ))\frac{e}{ m c} { \textbf{ p} \cdot \textbf{ \epsilon} }
where Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): { r}
is the electron coordinate, Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): { p} the kinetic moment, Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): { a^\dagger }_{ k,\epsilon}
, Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): a_{ k ,\epsilon}
are creation annihilation photon operators, Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): \omega= kc
, and Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): V
is the space volume.
The photon absorption is a first order perturbative process whose probability per unit of time Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): w(\omega) , for photon with polarization Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): \epsilon
state, is obtained from Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): H_{int} matrix elements by the Fermi golden rule. For one photon in the Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): V volume: Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): w(\hbar \omega) = \frac{(2 \pi)^2 e^2 \hbar }{\hbar \omega V m^2 } \sum_n \left|< n | exp(i { \textbf{ k} \cdot \textbf{ r}} ) { \textbf{ p} \cdot \textbf{ \epsilon} } | 0 >\right|^2 \delta( \hbar \omega-E_n)
where Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): <n|
is a complete set of eigenstates of matter, Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): <0| being the initial state, Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): E_n is the energy difference between Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): <n| and Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): <0| and we have retained only the resonating denominator. The absorption cross section is the quantity directly observable in experiment. It is obtained from the above equation dividing by the photon flux Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): c/V
- Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): \sigma(\hbar \omega) = \frac{(2 \pi)^2 e^2 \hbar }{\hbar \omega m^2 c } \sum_n \left|< n | exp(i { \textbf{ k} \cdot \textbf{ r}} ) { \textbf{ p} \cdot \textbf{ \epsilon} } | 0 >\right|^2 \delta( \hbar \omega-E_n)
The matrix element can be simplified in terms of Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): \textbf{ r}
by using the following identity for Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): \textbf{ p}
- Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): exp(i { k \cdot r} ) { p \cdot \epsilon } \simeq ( 1 + i k \cdot r + ... ) \frac{i m}{\hbar} [ H, { r \cdot \epsilon }] = \displaystyle\frac{i m}{\hbar} ([ H, { r \cdot \epsilon }] + i [ H, { k \cdot r} { r \cdot \epsilon }] /2 + ...)
Such substitution leads to:
Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): \sigma(\hbar \omega) = (2 \pi)^2 \alpha_0 \hbar \omega f(\hbar \omega)
with Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): f(\hbar \omega) = \sum_n \left| < n | i > \right|^2 \delta( \hbar \omega-E_n)
where
Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): |i> = { \textbf{ r} \cdot \textbf{ \epsilon} } + i\frac{( \textbf{ k \cdot r})( \textbf{ r \cdot \epsilon })}{2} +... | 0 >
is the state resulting from the application of the interaction operator on the ground state.
The Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): f(\hbar \omega)
is formally determined by its distribution moments which can be computed by iterative applications of the Hamiltonian : Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): \int f(E) E^n d E = \left| < i | H^n| i > \right|
To avoid numerical ill-conditioning coming from the numerical linear dependence of high order powers we proceed as explained in {weiss}. The Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): H
Hamiltonian is rescaled and shifted to a new Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): H' whose spectra is comprised within Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): ]-1,1[
- Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): < n | H' | n > \in \,\, ]-1,1[
The overlap integrals between the spectra and Chebychev polynomials Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): T_n(x) is calculated by recurrence. The Chebychev polynomials Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): T_n(x) are: Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): T_n(x) = cos( n * arccos(x)) and satisfy the recurrence relation:
Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): T_{m+1}(x) = 2 x T_m(x) - T_{m-1} (x)
The overlap integral are: Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): \mu_n = \int f'(x) T_n(x) = < i | T_n(H') | i > |
and can be obtained by applying the recurrence relation to the left side of Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): |i>
state:
Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): T_{m+1}(H')| i > = 2 H' T_m(H')| i > - T_{m-1} (H')| i >
and then obtaining the scalar product. The spectra is finally recovered as:
Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): f( b+a x) = \frac{ (\mu_0+ 2 \sum_{n=1}^{\inf} \mu_n T_n(x)) }{a \pi \sqrt{1 - x^2} }
where Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): b
and Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): a are the shift and the scaling factor respectively, involved in the Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): H' \rightarrow H transformation. The spectra can be obtained up to an arbitrary degree of resolution by applying iteratively the recurrence relation for a sufficient number of times. To avoid oscillation due to the truncation of the summation, a convolution is applied multiplying the Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): \mu_n components by the Jackson kernel^{[3]}.
Use of the code
The binary executable for absorption spectra calculation is named abscalc. Besides the usual files which are needed for a basic BigDFT calculation, abscalc needs to find, in the current work directory, at least the file ’{input.abscalc}’. This file is specific to absorption calculations and its structure will be explained here. The usual files needed for a basic BigDFT calculation and by abscalc are, we recall it here: ’{input.dft}’, ’{posinp.xyz}’, and the pseudopotential files ’{psppar.XX}’ for all the elements entering in the structure.
Others files may be needed when one selects the option of importing the self consistent potential from a previous SCF calculation, as will be explained later.
Concerning the ’{posinp.xyz}’ file we suggest, for a good result, to set a large box diameter of about Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): 25 Å or more, and periodic boundary conditions.
The structure of the file ’{input.abscalc}’ is the following:
The input parameter {in%iabscalc_type} which can be equal to Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): 1 or Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): 2 . In the first case the spectra is calculated as a series of Chebychev polynomials, while in the second case the Lanczos tridiagonalisation is applied. We suggest the first option because the Lanczos method requires, in order to keep orthogonality, that all the vectors of the tridiagonal base be stored in memory and this limits the number of steps that one can perform.
The parameter {in%iat_absorber}, an integer number, is the position (starting from one) in the {posinp.xyz} file of the absorber atom.
{in%L_absorber} : this parameter is the angular moment of the multipole component of the interaction operator that one wants to consider. It is Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): 1 for dipole interaction, Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): 2 for quadrupolar ...
{in%Gabs_coeffs(i), i=1,2*in%L_absorber +1}.
This are the components of the interaction operator, in Cartesian coordinates. These correspond to the BigDFT internal representation of tensors, which is based on spherical tensor.
For dipolar interaction these will be therefore the three Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): x,y,z components. In this present version the input are real numbers. Entry for the imaginary parts will have to be added in the future versions for magnetic dichroism.
{in%potshortcut } : This parameter determines how the local potential is obtained. If it is set to Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): 1 the potential is obtained from the superposition of atomic charges. If it is set to Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): 2 the potential is read from a previous SCF calculation. In this case the program searches, in the current work directory, the file {b2B_xanes.cube}. Such file contains the atomic positions and the local potential of a previous SCF run and can be created setting the parameter {in%output_grid} to Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): -2 in the file {input.dft}. The SCF potential is interpolated, using interpolant wavelets, to get the potential on the target grid. The SCF grid of the previous run must have less points than the target grid. This is easily the case, because in XANES calculations one usually sets a large diameter for the box. The fact of choosing a large box limits the spurious oscillations due to bands effects and to interference with the replica of the absorber atom in the other Bravais cells.
{in%nsteps } : this is the total number of Hamiltonian applications. The larger this number, the better the resolution you get in the calculated spectra.
Note about the absorber pseudopotential
In the actual version of the code (Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): 1.4 ) for the absorber pseudopotential one must use a Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): Z+1
approximation. For future versions of the code the development of on-the fly pseudopotentials, corresponding to the excited electronic structure, is foreseen.
Description of the output
In case of {in%iabscalc_type} equal to Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): 1 , the spectra is written, in the current work directory, with the name : {cheb_spectra_NUM } where {NUM} is the number of Chebychev components. The number of components increases at each Hamiltonian application and ,during the calculation, several partial results, with increasing number of components, are written.
The reversed PAW method
The BigDFT code, in the actual version as well as in the previous ones, uses a norm-conserving pseudo-potential. From this pseudo-potential we extract the corresponding PAW projector. This is the reverse of the PAW method. In the PAW method the projector determines the pseudopotential. In our case instead, the PAW projector is a function of the pseudo-potential. The implementation of such function is detailed here. We need to construct the projector for the absorber atom only because the core orbitals don’t overlap other atoms. For the absorber atom we proceed by calculating the SCF solution in the isolated atom case, once using pseudo-potentials, and another using the all-electrons potential. Then, for these potentials, we calculate an almost-complete set of eigen-functions which will be the basis for the projector. We need to do this only for selected angular moments. For example for the Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): K
edge and dipolar interaction only Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): l=1 eigen-functions are calculated.
To calculated the almost-complete set of eigen-functions we solve the radial Schrödinger equation in the interval Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): [0,R_c]
for the lowest Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): N_c eigenvalues. The radius Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): R_c must be large enough to contain the core orbitals. We have hard-coded it to the value of Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): 5au
. We have set Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): N_c=200 . With this choice of parameters the maximum eigenvalue of the almost-complete basis is of the order of Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): XXX Hartree . We can compare this value with the maximum eigenvalue representable with wavelets on a grid with spacing Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): d=0.3au , for a typical BigDFT calculation, which is Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): Emax = XXXX . Our almost-complete basis is indeed complete, after this comparison, because it covers completely the BigDFT spectra. Naming Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): \tilde \psi_{nlm}
and Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): \psi_{nlm} the solutions for the pseudo-potential and the all-electrons case respectively, the projector is thus written : Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): P= \sum_{nlm} |\tilde \psi_{nlm}><\psi_{n+\Delta_l,lm}| where Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): \Delta_l is the number of core orbitals contained in the pseudopotential. We show in fig xxx the radial parts Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): \tilde \psi_{nl}(r) and Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): \psi_{n+\Delta_l,l}(r) for Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): l=1 in the case of Iron, with a HGH pseudopotential
with Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): Zion=16 , for Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): n=1,2,xx,xxx . In figure we show the corresponding pseudo-eigenvalues Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): \tilde E_{n,l=1}
and the all-electrons ones Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): E_{n+\Delta_l,l=1}
. We can see that the agreement is within XXX for the first XXx eigenvalues and then increase still remaining better than XXX per cent, while the agreement between the psuedo-wavefunctions and the all-electrons ones remains excellent.
Test against an exactly solvable case.
In progress
- ↑ This reversed PAW procedure, that we describe and document here, looks similar to the one which seems to be used in Taillefumier et al., Phys. Rev. B 66, 195107 (2002)
- ↑ M. Blume, in Resonant Anomalous X-Ray Scattering, edited by G. Materlik, J. Sparks and K. Fisher (Elsevier, Amsterdam, 1994), p. 495.
- ↑ A. Weiss, G. Wellein, A. Alvermann and H. Fehske Rev. Mod. Phys. 78, 275 (2006)