X-ray absorption Debye-Waller factors from ab initio molecular dynamics
An ab initio equation of motion method is introduced to calculate the temperature-dependent mean square vibrational amplitudes which appear in the Debye-Waller factors in x-ray absorption, x-ray scattering, and related spectra. The approach avoids explicit calculations of phonon-modes, and is based instead on calculations of the displacement-displacement time correlation function from ab initio density functional theory molecular dynamics simulations. The method also yields the vibrational density of states and thermal quantities such as the lattice free energy. Illustrations of the method are presented for a number of systems and compared with other methods and experiment.
Thermal vibrations give rise to exponentially damped Debye-Waller (DW) factors in x-ray absorption spectra (XAS), x-ray diffraction (XRD), and related spectra. For example, in x-ray absorption fine structure spectra (XAFS) where refers to the mean square relative displacement (MSRD) of a given bond , is the photo-electron wave number, and the absolute temperature. Crozier et al. (1988) In XRD, and similarly in neutron diffraction (ND) and the Mössbauer effect, is the mean-square displacement of an atom along the momentum transfer vector , being the instantaneous displacement vector. Due to their strong variation with temperature, energy, and the geometrical structure of a material, accurate DW factors are crucial to a quantitative analysis of XAS; conversely, the lack of precise Debye-Waller factors is one of the main limitations to accurate structure determinations from experiment, especially for coordination numbers.
Various methods have been developed for obtaining these DW factors. Phenomenological models, e.g., correlated Einstein and Debye models Sevillano et al. (1979); Van Hung and Rehr (1997) are widely used in fitting but are often only semi-quantitative. More generally, they can be calculated in terms of Debye-integrals over appropriate projected vibrational densities of states (VDOS). In small molecules, explicit sums over modes can be used to calculate the VDOS.Loeffen et al. (1996); Dimakis and Bunker (2006) Such sums can also be used for periodic solids, both for crystallographic Debye-Waller factors and other thermodynamic quantities.Baroni et al. (2001); Lee and Gonze (1995); Rignanese et al. (1996) For complex materials, however, calculating and summing over modes can be a computational bottleneck. As an alternative, a Lanczos algorithm can be used to evaluate the VDOS, starting from a dynamical matrix (or Hessian), that can be obtained either from force-field models,Poiarkova and Rehr (2001); Krappe and Rossner (2002) or first principles DFT calculations.Vila et al. (2007) At high temperatures, brute force classical MD (or DFT/MD) methods can also be used to obtain moments of vibrational distribution functions,Vila et al. (2008) but such methods can fail at low temperatures when quantum statistics dominate. First principles DFT methods can also be computationally demanding, especially in complex systems.
In an effort to speed up the calculations, we present here a first principles approach based on an ab initio equation of motion (AEM) approach using DFT molecular dynamics calculations of displacement-displacement time-correlation functions. The method is a generalization of the equation of motion method Rehr and Alben (1977) for calculations of the VDOS, which was adapted for calculations of Debye-Waller factors based on force-field models.Poiarkova and Rehr (1999) However since accurate force-field models are not generally available, especially for complex molecules and solids, DFT or other ab initio methods are needed. Because depends primarily on the vibrational structure in the local environment around a given bond , the calculations can be carried out using relatively small clusters of atoms, without the use of periodic boundary conditions or other symmetry considerations. Thus the approach is applicable to general aperiodic materials.
Ii Equation-of-motion method
The theory used in the present study is a first principles extension of the equation of motion approach Rehr and Alben (1977); Poiarkova and Rehr (1999) for calculations of the VDOS and thermodynamic quantities that can be expressed as Debye-integrals over the VDOS. Our ab initio equation of motion (AEM) extension builds in dynamical structure in terms of first principles DFT calculations for a general structure, but does not rely on explicit calculations of the dynamical matrix (DM). The technique builds in Bose-Einstein statistics, and allows one to calculate the DW factors and related thermal properties either in real-time or the frequency domain. The AEM method has a number of computational advantages. It can be efficient even for large systems, since the method is local and diagonalization of huge matrices is not necessary. Also the computational time scales linearly with the size of a cluster. Anharmonic effects such as lattice expansion can be added using a cumulant expansion.Vila et al. (2007)
Our AEM method is based on calculations of the displacement-displacement correlation function in real time, using solutions of the coupled Newton’s equations of motion with DFT/MD methods. Such correlation functions are Fourier transforms of projected vibrational densities of states (VDOS), which are defined uniquely by the initial conditions. Physically the VDOS can be interpreted as the “sound” of a lattice “plucked” along a given set of initial displacements. Here is the number of atoms in the system which is centered within the region of interest and typically a few near-neighbors in radius. Regarding the total lattice potential energy of the crystal lattice as a function of the local atomic displacements from their thermal-equilibrium positions , and making use of a quasi-harmonic approximation, the equations of motion can be written as
with given initial displacements , and zero initial velocities . Here denote reduced displacements at site where is the atomic mass, and is the dynamical matrix of order . The matrix consists of second derivatives of the potential energy with respect to the atomic displacements and , where , are atomic sites and . Formally, the reduced displacement vectors can be expanded in normal coordinates and eigenmodes as
Substituting this relation into Eq. (1), leads to a standard eigenvalue problem for the normal modes,
After evaluating the thermal averages using Bose-Einstein statistics, one obtains for the normal coordinates
In applying these results for calculations of interest here, it is convenient to define a normalized displacement state
For example, for the MSRD for a given near-neighbor bond ,Poiarkova and Rehr (1999) the initial displacement state has ; , and otherwise , where is the reduced mass. A frequency domain expression for the MSRD can then be obtained from Eqs. (2-4) and summing over all modes, i.e.,
is the projected VDOS contributing to relative vibrational motion along and . The maximum frequency in Eq. (7) can be estimated from the relation where is the coordination number and is the near-neighbor force constant.
In order to obtain an equivalent time-domain expression for the VDOS, we calculate the cosine-transform of the displacement-displacement time-correlation function with an ad hoc exponential damping factor that limits the maximum time of the integration and MD runs,
Thus as a consequence of the damping factor, the projected VDOS of Eq. (10) is broadened by narrow -like functions of width typically chosen to be about 5% of the bandwidth. This broadening also smooths the otherwise discrete spectrum of the finite system used, but has practically no effect on integrated quantities. The spectral width is determined by the cutoff parameters and . These cutoff parameters also focus on the local environment by cutting off long distance behavior. The time-correlation function in Eq. (9) is
where is the number of non-vanishing displacements in . Instead of using the 2nd order differential equations in Eq. (1), in our approach the displacement state vector is determined by integrating the equations of motion numerically using Velocity-VerletVerlet (1967a, b) molecular dynamics with initial conditions as in ,
where , and are, respectively, the instantaneous position, velocity and acceleration of atom . The acceleration and is the force on atom . The Hellmann-Feynman theorem ensures that the forces can be calculated as the expectation value of the analytical derivative of the Hamiltonian with respect to the nuclear positions. This algorithm is efficient since an explicit calculation of the dynamical matrix at each time-step is not necessary.
Eqs. (7), (9) and (14) are the key formulas used in our AEM calculations. Throughout this work results obtained with Eqs. (7) and (9) will be labeled AEM-FT, while those obtained with Eq. (14) will be labeled AEM-RT. The form of Eq. (14) shows that, it is not essential to determine the VDOS as an intermediate step, and hence that can be calculated directly from the corresponding displacement-displacement autocorrelation function. Note that in the time domain the analog of the Bose-Einstein weight factor is . At long times when the weight factor is negative and reduces to at high temperatures and to at low. Due to the exponential damping, the net time integration limit is usually several vibrational cycles and typically requires about 25–35 time-steps per cycle for accuracy to a few percent. In addition, the singular behavior of the integrands in Eq. (9) and (14) must be handled with care. This is especially important at low temperatures due to zero-point motion. Thus in the time-domain, we further stabilize the long time behavior by convolving the time-correlation function with the inverse Fourier transform of a smoothed, low-frequency cutoff function , where is an appropriate cutoff frequency. In the frequency domain in Eq. (7) we replace the very low-frequency region with a similar cutoff or a Debye-model chosen to fit the very low frequency behavior of . All the integrals in our implementation of the AEM method are evaluated using the trapezoidal rule, which is appropriate for highly oscillatory integrands.
ii.2 Maximum Entropy Method
Since the MSRDs are obtained from Debye-integrals over the VDOS, a precise determination of the spectra is not important, as long as the leading moments are accurate. Thus, as an alternative approach, the projected density of states can be obtained approximately using the Maximum Entropy Method (MEM).Press et al. (1992) In this approach the VDOS is approximated as
where is the sampling interval in the time domain and is the desired order of the approximation. The MEM approach is well suited to represent phonon densities with sharp resonances, due to the presence of poles in Eq. (15). The coefficients can be obtained by solving the system of linear equations
and is the number of MD evolution steps. Although the MEM method can be less efficient than the direct FT approach, we find that it can be more stable in the low frequency region, since it is less sensitive to non-periodic trends in the time evolution. The reduced efficiency arises from the high order of approximation () needed to achieve an accurate representation of at all frequencies. Throughout this work results obtained with Eqs. (7) and (15) will be labeled AEM-MEM.
ii.3 Multiple scattering
The above real-time AEM method can also be used to calculate the MSRD for a given XAFS multiple-scattering path with legs. This MSRD corresponds to the mean-square fluctuation in the effective MS path length Poiarkova and Rehr (1999)
Here , where represent the directional unit vectors between the site and the sites before and after, along the multiple-scattering path . In analogy with the single scattering results, we obtain expressions similar to Eq. (7) for and Eq. (8) for , but with the weights in mode given by
These weights can be interpreted as the normalized probability that an initial displacement state , corresponding to a multiple-scattering path stretch, is in vibrational mode . Thus the initial displacements in the state are , (). Here the inverse reduced mass is
which is defined so that and is normalized.
ii.4 Other Dynamical Properties
Other dynamical properties can be obtained similarly, by generalizing the seed-state appropriately.Vila et al. (2007) For example, when the seed state is defined as a single-atom displacement, the resulting correlation function yields the mean square atomic displacements in x-ray scattering DW factors. Also, when all symmetry unique Cartesian atomic displacements are added, one obtains the total VDOS per site . This permits calculations of thermodynamic functions such as the vibrational free energy per site,Vila et al. (2007)
where is the Boltzmann constant. Finally, if the seed state is initialized with atomic displacements perpendicular to instead of parallel to it, we can generate the mean-square transverse displacement , which provides a correction to the lattice expansion.Vila et al. (2007)
ii.5 Computational Details
The micro-canonical (i.e., NVE) ensemble MD simulations for the applications presented here were done using vaspKresse and Furthmüller (1996) for the crystalline systems and siesta E. Artacho, D. Sánchez-Portal, P. Ordejón, A. García, and J. M. Soler (1999); J. M Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón, and D. Sánchez-Portal (2002) for the Zn-imidazole complex. These codes were chosen on the basis of efficiency, although in principle, any program capable of NVE dynamics can be interfaced with the AEM codes used in this work. The vasp simulations used standard ultrasoft pseudopotentials, and were optimized for efficiency in MD runs. The Ge calculations used a 222 -point grid with a plane-wave cutoff of 105 eV, while for ZrWO the grid was 444 and the cutoff was 297 eV. The siesta calculations used Troullier-Martins norm-conserving pseudopotentialsTroullier and Martins (1991) and standard double- basis sets with a single polarization function (DZP). The confinement-energy shift defining the numerical atomic orbitals was 10 meV. Finally, the Hartree and exchange-correlation potentials were represented on a real-space grid with a plane-wave-equivalent cutoff of 120 Ry within a (18.4 Å) cell. Both crystalline and molecular simulations used the PBE functional.Perdew et al. (1996) We have previously shown that the choice of exchange-correlation functional plays an important role in obtaining accurate MSRDs for metallic systems.Vila et al. (2007) However, here we only focus on non-metallic and molecular systems, for which the PBE functional yields reasonable accuracy compared to experiment.Vila et al. (2007)
ii.6 Efficiency Considerations
The efficiency of the AEM method depends on three factors: 1) The number of individual MSRDs that need to be computed, 2) the minimum and maximum frequencies that contribute to the VDOS, and 3) the quality of the ab initio MD. First, if a large number of MSRDs is needed, the computation of the full DM may be preferable since it yields all necessary DW factors with minimal additional effort. However, in most XAFS analysis only a handful of local DW factors need to be known accurately while those for more distant shells can approximated roughly using correlated Debye or Einstein models. For example, in the case of the coordination shell around a metallic center in a complex biomolecule the AEM approach can provide an efficient alternative to the Lanczos DM approach. Second, if a given MSRD has similar contributions from low and high frequency modes, the MD must have a short enough time-step to accurately represent the high frequency (25–35 steps per cycle) and a total run time with sufficient cycles of the low frequency (4-8 cycles). Third, the AEM approach can take advantage of efficient implementations of DFT energies and forces such as those used here, without relying on analytic second derivatives needed in the Lanczos DM approach or the equations of motion in Eq. (1).
Of the applications presented here, results for Ge and Zn-tetraimidazole can be more efficiently treated using the Lanczos DM approach. In the case of Ge this is due to the simplicity of the unit cell. In the case of Zn-tetraimidazole, first there are a relatively small number of modes and second the modes cover a broad range of frequencies that would require small time-steps and a long total simulation time to represent accurately. On the other hand, the zirconium tungstate (ZrWO) system, illustrates the definite advantage of the AEM approach for complex systems, since only a handful of MSRDs are needed for XAFS, while the unit cell contains hundreds of atoms. Based on our experience with the DM Lanczos approach, we estimate that the AEM approach would be nearly two orders of magnitude faster than a dynamical matrix calculation.
As a relatively simple test case, the AEM was applied to a crystalline germanium system using an 64-atom supercell generated by repeating 222 times the diamond cubic cell, with the experimental lattice constant of 5.6575 Å. The MD simulations used a 2 fs timestep and a total simulation time of 4.5 ps. The initial structure was generated by introducing a 4.8% bond stretch to one of the nearest neighbor pairs in the cell.
The correlation function resulting from the Velocity-Verlet time evolution is shown in Fig. 1. As expected, the oscillations are dominated by a single mode with a period of about 117 fs, associated with the Ge-Ge optical mode stretch. This dominant behavior can also be observed in the VDOS shown in Fig. 2, where the optical modes are centered at about 8.5 THz. An integration time of about 2 ps is adequate to obtain phonon-spectra with a spectral broadening of about 5%. The centroid of the VDOS is located at about 8 THz, in good agreement with Einstein models for the nearest neighbor single-scattering path with an Einstein frequency of 7.55 THz.Dalba et al. (1999) It should be noted that although the integration time for optical mode is well above that needed for convergence, the net integration time for lower frequencies around 5 THz is just adequate. Due to the singular behavior in Eq. (6), an adequate time integration for the lower frequency components is essential, and is especially important at low temperatures for some of the systems discussed in the next section.
The MSRDs calculated for the nearest neighbor Ge-Ge bond are shown in Fig. 3. The agreement with experiment is quite good, with an average error of 4% for the AEM-FT approach and 2% for the AEM-MEM approach. For comparison, the DM-Lanczos approach has an average error of 2%. Fig. 3 also shows the results obtained with the real-time approach of Eq. (14) and a frequency cutoff of 1.7 THz as in the FT and MEM approaches. As expected, given the formal equivalence between Eq. (14) the FT approach with an intermediate calculation of the VDOS, the results are nearly identical.
To explore the accuracy and efficiency of the AEM-FT and AEM-MEM approaches, we have also integrated the correlation function both for shorter times and for larger time steps. For Ge we find that the total integration time can be reduced to about 1 ps without significant loss of accuracy. This corresponds to about 10 periods of the 8.5 THz dominant frequency. For integration times of about 500 fs the mean error for the MSRD increases to 8% for the FT approach and to 16% for MEM. From the point of view of the length of the time step, both the FT and MEM approaches are extremely resilient. In both cases the mean errors for the Ge MSRD remain constant with time steps up to 24 fs. This corresponds to approximately five samples per period of the 8.5 THz frequency. Such large time steps, however, might not be feasible within the MD simulation itself due to loss of energy conservation in the Verlet algorithm.
As an example of other dynamical quantities that can be obtained with the AEM approach, Fig. 4 shows the total phonon density of states for Ge calculated with the FT and MEM approaches. For comparison broadened dynamical matrix Lanczos and experimentalNelin and Nilsson (1972) results are also included. This VDOS was obtained by applying a single atomic displacement along the axis, as described in II.3, and by propagating as for for 4.5 ps. Overall, the centroid of the DOS is accurately reproduced by all methods: The centroid of the experimental DOS is located at 5.8 THz, while the FT and MEM approaches place it at 6.0 and 5.7 THz, respectively. The spread (i.e., 2nd moment) of the DOS is also well reproduced with the FT and MEM, giving 2.8 and 2.9 THz, respectively, versus 2.6 THz in the experiment. Finally, all methods reproduce the positions and weights of main features of the experimental VDOS quantitatively. On average the positions of the peaks deviate by at most 0.4 THz (i.e., about 4% of full bandwidth) and the relative weights are within 5% of those observed in experiment.
The accuracy of the total VDOS can also be gauged by comparing with the experimentally measured atomic MSD for Ge. Fig. 5 shows the MSD computed using the total VDOS shown in Fig. 4. The AEM results are in excellent agreement with those obtained with the full DM Lanczos approach and in good agreement with the available experimental results Peng et al. (1996) except at low temperatures.
As an example of a complex molecule, Zn-tetraimidazole was simulated using the full structure shown in Fig. 6.
This structure was optimized in siesta and one of the equivalent Zn-N bonds was distorted with a 3.4% bond stretch. The MD simulations used a 3 fs timestep and a total simulation time of 3.9 ps. Given its large number of degrees of freedom, the dynamics of Zn-tetraimidazole are significantly more complicated than those of Ge. This can be seen in the correlation function shown in Fig. 7, which exhibits a superposition of several modes.
The dominant contributions can be analyzed by examining the VDOS in Fig. 8.
The DM approach exhibits three dominant frequencies at 5, 13 and 25 THz, which contribute 32, 18 and 24%, respectively, of the MRSD value. It is interesting to note that the weight of the associated poles is 9, 13 and 31%, further highlighting the importance of the correct representation of the low frequency modes. In principle, the Zn-N path should be dominated by low frequency Zn-ligand tetrahedral modes. Loeffen et al.Loeffen et al. (1996) find that these modes appear at about 6.5 THz, in fair agreement with our principal contribution at 5 THz. Although the VDOS calculated with the AEM-FT and AEM-MEM approaches are in good agreement with each other, they have small differences with respect to the Lanczos DM VDOS. For instance, the mode at 25 THz is blueshifted about 2 THz in the real-time approaches. Fig. 9 shows that the agreement between the MSRDs calculated from the different VDOS is quite good. At 8% error, the theoretical results are less accurate than those obtained for Ge. They are, however, still within the error margins of the available experimental value at 20K. The larger error is likely due to the quality of the basis set used in the siesta calculations.
Our final example is zirconium tungstate (ZrWO), a ceramic that exhibits negative thermal expansion (NTE). This system is quite challenging, having a complex unit cell that puts the calculation of the DM for the Lanczos approach beyond the reach of our current implementation and computational capabilities. Here we have applied the AEM approach to a 352-atom supercell (Fig. 10) made of 222 repetitions of the unit cell. The simulations used the experimental unit cell lattice constant of 9.1546 Å and a timestep of 4 fs, for a total simulation time of 1.5 ps. ZrWO has several interactions of interest, including Zr-Zr, W-W, W-O and two inequivalent nearest-neighbor Zr-O bonds with distances 2.03 and 2.11 Å. In principle, any of these interactions can be studied using the AEM approach. As a proof of principle here we study the of the shortest of the Zr-O bonds by using an initial structure corresponding to a 3.8% bond stretch.
For a Zr-O distortion, the dynamics of ZrWO are not as complex as those observed for Zn-tetraimidazole. The correlation function (Fig. 11) is mostly dominated by a mode with a 40 fs period superposed on a mode with a period approximately three times longer. Visual inspection of the MD trajectory reveals that the 40 fs mode is associated principally with the longitudinal Zr-O stretch mode. These vibrational modes can be clearly seen at about 25 and 8 THz, respectively, in the VDOS shown in Fig. 12. As in the previous examples, the agreement between the AEM-MEM and AEM-FT approaches is very good. The agreement with the mode frequencies observed in the experimental Raman spectrum is also quite good. The FT and MEM VDOS show modes at approximately 7.7, 24.5, 27.7 and 31.7 THz, compared to the experimental peaks at 5.7-11.8, 23.8, 27.9 and 31.0 THz. The 5.7-11.8 THz peaks are associated mostly with modes located on the tungstate ion and with some low frequency WO modes.Ravindran et al. (2003) The 23.8, 27.9 and 31.0 THz peaks correspond exclusively to asymmetric WO modes. It is interesting to note that the dynamics of this system are quite complex. Since the WO units are very stiff, the ZrO units must rotate as the WO units translate. Cao et al. (2003) Thus, a simple distortion of the Zr-O bond is able to activate both the low and high frequency modes.
Fig. 13 shows the nearest-neighbor Zr-O MSRD as a function of temperature. As expected, given the similarity of their VDOS, the FT and MEM values are in very good agreement. The direct integration RT approach also agrees well with the FT approach, at least for higher temperatures. The agreement with experimentCao et al. (2003) is also quite good, with all theories falling within the experimental error bars for most of the temperature range. The largest disagreement occurs in the 80-140K region where other MSRDs (W-W and Zr-Zr) are known to have an anomaly that is likely related to the NTE.Cao et al. (2003)
We have introduced an ab initio equation of motion (AEM) method for calculations of the MSRDs , needed for Debye-Waller factors in x-ray absorption, x-ray scattering, and related spectra. The method is based on calculations of displacement-displacement time correlation functions from ab initio density functional theory molecular dynamics simulations, using the Velocity-Verlet time-evolution algorithm. Thus the approach avoids the need for explicit calculations of phonon-modes or the dynamical matrix. The AEM method builds in Bose-Einstein statistics and yields the vibrational density of states (VDOS) as either cosine Fourier Transforms of displacement-displacement correlation functions or through the Maximum Entropy Method. The MSRDs and other thermal quantities such as the lattice free energy, are obtained in terms of Debye-integrals over the VDOS. Alternatively, the MSRDs can be computed directly from the correlation functions by using the time-domain counterpart of the Bose-Einstein weight factor. Application of the method to a number of systems show that the approach is computationally advantageous for large, complex systems, and is in quantitative agreement with other methods and with experimental results.
Acknowledgements.We thank J. Kas, J. Vinson, S. Williams, and F. Bridges for comments and suggestions. This work is supported in part by NSF Grant PHY-0835543 (FDV and JJR). One of us (VEL) thanks the REU program at the University of Washington in summer 2010, which is supported by NSF REU Grant PHY-0754333, where part of this work was carried out.
- E. D. Crozier, J. J. Rehr, and R. Ingalls, in X-Ray Absorption: Principles, Applications, Techniques of EXAFS, SEXAFS, and XANES, edited by D. C. Koningsberger and R. Prins (Wiley, New York, 1988), p. 375.
- E. Sevillano, H. Meuth, and J. J. Rehr, Phys. Rev. B 20, 4908 (1979).
- N. Van Hung and J. J. Rehr, Phys. Rev. B 56, 43 (1997).
- P. W. Loeffen, R. F. Pettifer, and J. Tomkinson, Chemical Physics 208, 403 (1996).
- N. Dimakis and G. Bunker, Biophysical J. 91, L87 (2006).
- S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
- C. Lee and X. Gonze, Phys. Rev. B 51, 8610 (1995).
- G.-M. Rignanese, J.-P. Michenaud, and X. Gonze, Phys. Rev. B 53, 4488 (1996).
- A. Poiarkova and J. J. Rehr, J. Synchrotron Radiat. 8, 313 (2001).
- H. J. Krappe and H. H. Rossner, Phys. Rev. B 66, 184303 (2002).
- F. D. Vila, J. J. Rehr, H. H. Rossner, and H. J. Krappe, Phys. Rev. B 76, 014301 (2007).
- F. Vila, J. J. Rehr, J. Kas, R. G. Nuzzo, and A. I. Frenkel, Phys. Rev. B 78, 121404 (2008).
- J. J. Rehr and R. Alben, Phys. Rev. B 16, 2400 (1977).
- A. V. Poiarkova and J. J. Rehr, Phys. Rev. B 59, 948 (1999).
- L. Verlet, Phys. Rev. 159, 98 (1967a).
- L. Verlet, Phys. Rev. 165, 201 (1967b).
- W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in Fortran: The Art of Scientific Computing, Second Edition (Cambridge University Press, Cambridge, 1992), pp. 565–566.
- G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
- E. Artacho, D. Sánchez-Portal, P. Ordejón, A. García, and J. M. Soler, Phys. Status Solidi B 215, 809 (1999).
- J. M Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón, and D. Sánchez-Portal, J. Phys.: Condens. Matter 14, 2745 (2002).
- N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991).
- J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- G. Dalba, P. Fornasini, R. Grisenti, and J. Purans, Phys. Rev. Lett. 82, 4240 (1999).
- G. Nelin and G. Nilsson, Phys. Rev. B 5, 3151 (1972).
- L.-M. Peng, G. Ren, S. L. Dudarev, and M. J. Whelan, Acta Crystallographica Section A 52, 456 (1996).
- T. R. Ravindran, A. K. Arora, and T. A. Mary, Phys. Rev. B 67, 064301 (2003).
- D. Cao, F. Bridges, G. R. Kowach, and A. P. Ramirez, Phys. Rev. B 68, 014303 (2003).