A polynomial Ansatz for Normconserving Pseudopotentials
Abstract
We show that efficient normconserving pseudopotentials for electronic structure calculations can be obtained from a polynomial Ansatz for the potential. Our pseudopotential is a polynomial of degree ten in the radial variable and fulfils the same smoothness conditions imposed by the TroullierMartins method [Phys. Rev. B 43, 1993 (1991)] where pseudopotentials are represented by a polynomial of degree twentytwo. We compare our method to the TroullierMartins approach in electronic structure calculations for diamond and iron in the bcc structure and find that the two methods perform equally well in calculations of the total energy. However, first and second derivatives of the total energy with respect to atomic coordinates converge significantly faster with the plane wave cutoff if the standard TroullierMartins potentials are replaced by the pseudopotentials introduced here.
I Introduction
Density Functional Theory (DFT) Hohenberg and Kohn (1964); Kohn and Sham (1965) is nowadays the most common computational method for calculating electronic structure properties of molecules and solids from first principles. The standard formulation of DFT is based on local, semilocal and hybrid exchangecorrelation functionals Martin (2004), and has proven to be very successful in describing the electronic structure of weaklycorrelated systems Kohn (1999); Pople (1999); Bowler (2016); Lejaeghere et al. (2016).
A second pillar of DFT calculations is the concept of pseudopotentials Martin (2004). They allow one to reduce the numerical cost of electronic structure calculations significantly by eliminating atomic core states that do not participate in chemical bonding. Modern DFT calculations are based on normconserving pseudopotentials (NCPPs) Hamann et al. (1979); Bachelet et al. (1982); Vanderbilt (1985); Kerker (1980); Troullier and Martins (1991); Rappe et al. (1990); Reis et al. (2003), ultrasoft pseudopotentials (USPPs) Vanderbilt (1990); Laasonen et al. (1993); Kresse and Hafner (1994); Moroni et al. (1997); Garrity et al. (2014) or projector augmented wave (PAW) methods Blöchl (1994); Holzwarth et al. (1998); Kresse and Joubert (1999).
Calculations based on NCPPs are typically carried out in the KleinmanBylander approach Kleinman and Bylander (1982) where a set of pseudopotentials for different angular momenta are represented by one local potential and a set of separable, nonlocal projectors. This method enables computationally efficient electronic structure calculations based on simple planewave representations of the Hamiltonian P. Giannozzi et al. (2009). The USPP and PAW methods have been developed to further improve (i) the transferability of pseudopotentials among different chemical environments and oxidation states and (ii) the numerical efficiency of ab initio calculations. However, NCPPs are significantly easier to implement than USPPs and PAWs due to the reconstruction/augmentation term in these methods Miwa (2011). Hence NCPPs remain the method of choice in more advanced calculations like, e.g., densityfunctional perturbation theory Baroni et al. (2001) or manybody perturbation theory Hedin (1965); Hybertsen and Louie (1985), and further improving their efficiency is a highly desirable goal Hamann (2013).
One of the most successful methods for generating NCPPs is the TroullierMartins (TM) method Troullier and Martins (1991). At the heart of this approach is an Ansatz for the pseudowavefunction (PWF), and its parameters are chosen such that normconservation and a set of smoothness conditions are met Troullier and Martins (1991). Like in other NCPP methods Hamann et al. (1979); Bachelet et al. (1982); Vanderbilt (1985); Kerker (1980); Rappe et al. (1990), the pseudopotential itself is obtained from the PWF by an inversion of the radial Schrödinger equation. In the case of the TM method, the pseudopotential turns out to be a polynomial of degree twentytwo in the radial variable due to the specific Ansatz for the PWF.
Here we pursue a different approach and show that efficient NCPPs can be obtained from a polynomial Ansatz for the pseudopotential. In this way, the final step of inverting the radial Schrödinger equation is not required. Our pseudopotentials are represented by a polynomial of degree ten in the radial variable and we impose the same smoothness conditions as in the TM method.
We carry out electronic structure calculations for diamond and iron with the Quantum Espresso code P. Giannozzi et al. (2009) and analyse the performance of our pseudopotentials compared to the TM potentials. While both pseudopotential methods perform equally well in calculations of the total energy and pressure, we find that our pseudopotentials speed up the convergence of calculations for atomic forces and phonon frequencies. Our pseudopotential scheme thus promises big computational savings when calculating structural relaxations and phonon frequencies in large systems.
Note that our pseudopotential method is related to unpublished work by von Barth and Car who suggested the generation of NCPPs by varying the parameters in an Ansatz for the pseudopotential. Their proposed parametrisation of the pseudopotential is different from ours and has been used, e.g., in Dal Corso et al. (1993).
This paper is organised as follows. In the Methods section II we briefly review the theory of NCPPs (Sec. II.1) and the TM method (Sec. II.2). These two sections set the stage for our novel approach for constructing NCPPs which is described in Sec. II.3. We then compare the performance of our NCPPs with the standard TM approach in electronic structure calculations for diamond and iron. These results are presented in Sec. III, and a brief summary is given in Sec. IV.
Ii Methods
In Sec. II.1 we briefly review the general theory of NCPPs. The TM method Troullier and Martins (1991) and our novel approach for generating NCPPs are described in Secs. II.2 and II.3, respectively.
ii.1 Normconserving pseudopotentials
The generation of pseudopotentials starts with an atomic DFT calculation which results in the selfconsistent allelectron (AE) potential
(1) 
where is the the Coulomb potential of the ion core, is the Hartree energy of all electrons, is the exchangecorrelation potential and is the electron density Martin (2004). For nonrelativistic calculations the AE radial wavefunctions with principal quantum number , orbital angular momentum and energy are determined by the radial Schrödinger equation,
(2) 
where is the radial coordinate of the electron. We introduce a shorthand notation for states corresponding to valence states and denote their energies and radial wavefunctions by and , respectively. For each of these valence states one introduces an dependent pseudopotential (we denote all pseudopotential quantities by caligraphic letters), and the nodeless pseudowavefunction (PWF) is the solution to the nonrelativistic radial Schrödinger equation with potential and energy ,
(3) 
In order to replace the AE potential by the set of dependent pseudopotentials in electronic structure calculations, the following conditions have to be met Hamann et al. (1979); Martin (2004):

The ground state energy of the pseudopotential has to coincide with the AE valence energy ,
(4) 
The PWF must coincide with the AE wavefunction outside a cutoff region ,
(5) and we have suppressed the dependence of for simplicity.

The norm of the PWF has to equal the norm of the AE wavefunction inside the core region,
(6) This normconservation condition ensures that the first energy derivative of the logarithmic derivatives of the AE and pseudo wavefunctions agrees at Hamann et al. (1979).
Since is nodeless by construction, the pseudopotential can be obtained from the PWF by inverting the radial Schrödinger equation (3),
(7) 
Conditions (NC1) and (NC2) thus enforce that coincides with outside the cutoff region . Conversely, for implies together with conditions (NC1), (NC3) and the normalisation of the radial wavefunction that condition (NC2) is met. An equivalent set of conditions is thus given by (NC1), (NC3) and (NC2’), where (NC2’)is defined as follows:

must coincide with outside the cutoff region,
(8) where is the potential inside the core.
In order to use the pseudopotential in electronic structure calculations one has to generate ionic pseudopotentials,
(9) 
where and are the Hartree and exchangecorrelation potentials calculated from the PWFs, respectively. The semilocal pseudopotential operator is then given by
(10) 
where projects the wavefunction onto the th angular momentum component.
ii.2 Review of the TroullierMartins method
The TM method Troullier and Martins (1991) for generating NCPPs makes the following Ansatz for the PWF,
(11) 
Outside the core region coincides with the AE wavefunction , and hence (NC2) is automatically fulfilled. Inside the core region , where is a polynomial of degree twelve which contains only even powers of ,
(12) 
The seven coefficients in Eq. (12) are chosen such that the PWF and the potential exhibit the following smoothness conditions at and :

should have zero curvature at the origin, i.e., , where the double prime denotes the second derivative with respect to . This condition can be cast into a nonlinear equation for the coefficients and and reads Troullier and Martins (1991),
(13) 
and its first four derivatives should be continuous at . These five conditions result in a system of linear equations for five of the coefficients and can be solved with standard methods.
(S1) and (S2) constrain six out of the seven coefficients in Eq. (12). The last free parameter is determined by condition (NC3) in Sec. II.1,
(14) 
This is a highly nonlinear equation, and a solution is not guaranteed for all cutoff parameters .
The above procedure fully determines the PWF. Finally, is obtained from the PWF via Eq. (7) by setting such that condition (NC1) is fulfilled. With the help of Eq. (11), the pseudopotential in the core region can be written as
(15) 
The definition of in Eq. (12) allows us to write as a polynomial in ,
(16) 
We thus find that the pseudopotential in the core region is a polynomial of degree twentytwo, and contains only even powers of . The explicit expressions for the coefficients in terms of are given in Appendix A.
ii.3 Variation of the potential
Here we introduce a novel approach for generating NCPPs via the direct variation of a polynomial Ansatz (PA) for the potential. We label this potential by and associated quantities with the superscript PA. The pseudopotential outside the core region must coincide with the AE potential such that condition (NC2’) is satisfied. Inside the core region we make the following Ansatz,
(17) 
in Eq. (17) contains only even powers of like of the TM method in Eq. (16). However, our Ansatz for is not equivalent to the TM method since the degree of the polynomial is only ten and thus much lower than the degree of .
Next we show that our Ansatz for is sufficient for creating a NCPP with the same smoothness conditions (S1) and (S2) of the TM method. First, condition (S1) can be satisfied by setting in Eq. (17). The Ansatz in Eq. (17) then contains only a constant term and polynomials that are at least of order 4, and hence the second derivative of vanishes at .
Second, the five conditions in (S2) concerning the continuity of the PWF and its derivatives at can be cast into conditions on the potential via Eq. (7). We find that and its first two derivatives must be continuous at , and these three conditions result in a linear system of equations for three of the coefficients . Here we choose to express , and in terms of and , and the explicit expressions for these coefficients are given in Appendix A.
By imposing the smoothness conditions (S1) and (S2) as described above, only depends on the two coefficients and . These parameters must be chosen such that conditions (NC1) and (NC3) are satisfied. To this end we find the ground state energy and wavefunction of by solving the eigenvalue problem in Eq. (3) on a discrete grid. The optimal parameters and for satisfying conditions (NC1) and (NC3) are found by standard Newton methods, and we find fast convergence for suitable starting values. In particular, we find that the eigenvalues of our PA pseudopotentials match those of the AE calculation very accurately. For all systems considered in Sec. III, the maximal difference between PA and AE eigenvalues is 0.5 meV for s orbitals, and for p and d orbitals the deviation is even smaller by at least one order of magnitude.
Iii Results
In order to test the method introduced in Sec. II.3, we compare the performance of the PA and TM pseudopotentials in electronic structure calculations for Carbon and Iron. All DFT calculations based on pseudopotentials are performed with the Quantum Espresso code P. Giannozzi et al. (2009) using the PZ functional Perdew and Zunger (1981), and AE calculations are carried out with the same exchangecorrelation functional and the ELK code elk ().
In a first step, we compare the screened pseudopotentials for the TM and PA methods. To this end we perform nonrelativistic (semirelativistic) DFT calculations for a single C atom (Fe ion). We then generate the TM and PA pseudopotentials for the valence states as described in Sec. II, and the results are shown in Fig. 1. Note that we choose the same cutoff parameters for the TM and PA potentials for a fair comparison of the two methods. The 2S and 2P pseudopotentials for C are shown in Fig. 1(a), and we find that the 2S potentials of the two methods differ only very little. The differences are more pronounced for the 2P pseudopotentials, in particular near where and differ by about . Next we discuss the pseudopotentials of the TM and PA methods for shown in Fig. 1(b). We find that the 3S and 3P potentials are very similar for both methods. On the contrary, the 3D potential of the PA and TM methods differ by near . Furthermore, the TM potential for the 3D state “oscillates” around the corresponding PA potential, i.e., the potentials cross twice near and . This behaviour is consistent with the fact that the degree of the polynomial describing in the core region is more than twice as large as in the case of .
In the case of iron we performed tests of pseudopotential transferability for different oxidation states and different occupations of the 3d orbitals. In Tab. 1 we report the energy difference between the AE and pseudopotential eigenvalues for a set of spinpolarised configurations. We find that the energy differences for the TM and PA pseudopotentials are of the same order of magnitude, and no clear trend distinguishing the two methods can be established. It follows that the TM and PA pseudopotentials are both transferable to a similar extent.
Fe TM  Fe PA  
Configuration  
Fe: 3d 4s  0.0166  0.0072  0.0223  0.0016 
Fe: 3d 4s  0.0132  0.0012  0.0187  0.0067 
Fe: 3d 4s  0.0083  0.0166  0.0086  0.0161 
Fe: 3d 4s  0.0081  0.0132  0.0117  0.0087 
In order to investigate the transferability of the pseudopotentials generated via the TM and PA methods further, we introduce the dimensionless logarithmic derivative
(18) 
where is the Bohr radius and is the radial wavefunction corresponding to energy and orbital angular momentum . can be either an AE wavefunction or a PWF, and for an ideal pseudopotential the logarithmic derivatives of the PWFs and the AE wavefunctions match over a wide range of energies. A comparison of for AE, TM and PA wavefunctions in the case of Fe is shown in Fig. 2 for four different angular momentum channels. We find that the logarithmic derivatives of the TM and PA wavefunctions are very similar, and both of them follow the results of the AE calculation closely. This result is consistent with our findings in Tab. 1 and shows that the PA method results in pseudopotentials that are as transferable as their TM counterparts.
Next we consider bulk electronic structure calculations for carbon and iron. In all calculations we use the KleinmanBylander representation Kleinman and Bylander (1982) of the pseudopotentials and choose the pseudopotential with the largest angular momentum as the local potential. First, we calculate the equilibrium lattice constants of diamond and bcc iron. To this end we evaluate the energy vs. volume curve and fit it to the Vinet equation of state Vinet et al. (1987). The results are reported in Tab. 2 and show that the two PP methods result in practically equivalent results for carbon. In the case of iron, differences between the two pseudopotential methods are also small.
Atom  (Å)  (GPa)  
C TM  3.5402  458.1  3.558 
C PA  3.5434  456.3  3.609 
Fe TM  2.7956  189.8  6.769 
Fe PA  2.8104  184.7  8.458 
Fe FPLAPW  2.7479  193.1  4.170 
Results from an AE calculation for iron are shown in the last row of Tab. 2 and demonstrate that both pseudopotential methods overestimate the equilibrium lattice spacing by about 2%. The lattice spacing obtained from the TM method is 0.5% smaller than the one obtained with the PA method and thus slightly closer to the AE result. Note that the pressure derivative of the bulk modulus in Tab. 2 takes on large values for all three methods in the case of iron. This problem has been reported in Zhang et al. (2010) and is due to an incipient magnetic transition at volumes larger than the equilibrium, and is an artefact of the fitting.
In order to better distinguish between the two pseudopotential methods, we systematically investigate the convergence of the total energy, the pressure and higher derivatives of the total energy with the plane wave cutoff energy . For carbon we use the structure of cubic diamond at the experimental lattice spacing of 3.567 Å, while for iron we use the ferromagnetic bcc structure with a lattice spacing of 2.870 Å. In order to rule out any possible bias from the experimental lattice spacing, we also perform calculations at other values and obtain qualitatively identical results. The convergence of the total energy and pressure with E for carbon is shown in Figs. 3(a) and (b), respectively.
We find that the rate of convergence for these two quantities is similar for both methods. The corresponding results for iron are shown in Figs. 4(a) and (b), and here the PA pseudopotential gives rise to a well converged pressure at smaller values of compared to the TM method.
The convergence of the first and second derivatives of the total energy with respect to atomic coordinates is investigated as follows. For analysing the first energy derivative, we displace one atom by of the lattice constant in the direction and extract the value of the restoring force as a function of . Information about the convergence of the second energy derivatives is obtained by calculating the frequency of the zonecenter optical phonon mode in the undistorted lattice. The results for carbon are shown in Figs. 3(c) and (d), and illustrate that the PA potentials give better results for forces and phonon frequencies at small cutoff energies compared to the TM method. This trend is even more pronounced for iron as can be seen in Figs. 4(c) and (d). In particular, the phonon frequency calculation converges much faster for the PA pseudopotentials compared to the TM method.
Finally, for completeness we provide the converged values of the total energy, pressure, force and phonon frequencies for C and Fe in Tab. 3. We find that both pseudopotential methods result in very similar values for most quantities. The largest deviation occurs for pressure in the case of Fe, where the TM and PA values differ by . This discrepancy is related to the fact that the TM and PA pseudopotentials give rise to slightly different equilibrium lattice constants as shown in Tab. 2. On the other hand, the differences between TM and PA values for total energy, force and phonon frequency are less than for both C and Fe.
Atom  (Ry)  (kPa)  (Ry/au)  () 

C TM  22.870  101.5  0.2356  1286.47 
C PA  22.858  89.2  0.2362  1288.14 
Fe TM  245.207  182.4  0.0365  221.04 
Fe PA  244.487  145.0  0.0379  222.80 
Iv Summary and discussion
In this paper we have introduced a novel method for generating NCPPs. We represent the pseudopotential by a polynomial of degree ten in the radial variable and impose the same set of smoothness conditions as the TM method Troullier and Martins (1991). The latter approach makes an Ansatz for the PWF and obtains the pseudopotential by an inversion of the radial Schrödinger equation. The resulting TM potential is also a polynomial in the radial variable, but its degree of twentytwo is significantly higher than the polynomial representing our pseudopotentials. A direct Ansatz for the pseudopotential instead of an Ansatz for the wavefunction thus allows us to impose the same set of smoothness conditions as the original TM method, while at the same time the degree of the polynomial representing the pseudopotential can be significantly reduced.
This reduced polynomial degree comes at an increased numerical cost for the generation of our pseudopotentials compared to the TM method. The numerically most expensive step for generating a TM potential is in finding a solution to the norm conservation condition, see Eq. (14). Each integral requires operations, where is the number of points in the discrete grid. On the other hand, our method requires to solve an eigenvalue problem on the same grid and thus needs operations. However, we find that the search for the optimal coefficients resulting in a normconserving potential with the correct ground state energy converges quickly with Newton methods.
In Sec. III we compared our PA pseudopotentials with the TM method in electronic structure calculations for carbon and iron in the bcc structure. A comparison of the screened PA and TM pseudopotentials for C and reveals that the two methods result in similar potentials, and the largest differences occur for the states with highest angular momentum. Moreover, we find that the TM potentials of higher angular momentum states cross the PA potentials several times as a function of the radial variable. These ”oscillations“ are a signature of the higherorder polynomials in the TM potentials. The main result of this work is that our PA potentials speed up the convergence of the first and second derivatives of the total energy with respect to atomic coordinates. More specifically, we find our new scheme produces accurate forces and Hessians at lower cutoff energies than the TM method. This could constitute a big computational saving when calculating structural relaxations and phonon frequencies in large systems. It follows that the increased cost in generating our PA pseudopotentials is by far outweighed by the promised savings in largescale electronic structure calculations.
Acknowledgements.
The authors acknowledge financial support from the National Research Foundation and the Ministry of Education, Singapore. DJ acknowledges funding from the European Research Council under the European UnionÊ¼s Seventh Framework Programme (FP7/20072013)/ERC Grant Agreement no. 319286, QMAC.Appendix A Coefficients
References
 Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
 Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
 Martin (2004) R. M. Martin, Electronic Structure (Cambridge University Press, Cambridge, 2004).
 Kohn (1999) W. Kohn, Rev. Mod. Phys. 71, 1253 (1999).
 Pople (1999) J. A. Pople, Rev. Mod. Phys. 71, 1267 (1999).
 Bowler (2016) D. R. Bowler, Journal of Physics: Condensed Matter 28, 421001 (2016).
 Lejaeghere et al. (2016) K. Lejaeghere, G. Bihlmayer, T. Björkman, P. Blaha, S. Blügel, V. Blum, D. Caliste, I. E. Castelli, S. J. Clark, A. Dal Corso, S. de Gironcoli, T. Deutsch, J. K. Dewhurst, I. Di Marco, C. Draxl, M. Dułak, O. Eriksson, J. A. FloresLivas, K. F. Garrity, L. Genovese, P. Giannozzi, M. Giantomassi, S. Goedecker, X. Gonze, O. Grånäs, E. K. U. Gross, A. Gulans, F. Gygi, D. R. Hamann, P. J. Hasnip, N. A. W. Holzwarth, D. Iuşan, D. B. Jochym, F. Jollet, D. Jones, G. Kresse, K. Koepernik, E. Küçükbenli, Y. O. Kvashnin, I. L. M. Locht, S. Lubeck, M. Marsman, N. Marzari, U. Nitzsche, L. Nordström, T. Ozaki, L. Paulatto, C. J. Pickard, W. Poelmans, M. I. J. Probert, K. Refson, M. Richter, G.M. Rignanese, S. Saha, M. Scheffler, M. Schlipf, K. Schwarz, S. Sharma, F. Tavazza, P. Thunström, A. Tkatchenko, M. Torrent, D. Vanderbilt, M. J. van Setten, V. Van Speybroeck, J. M. Wills, J. R. Yates, G.X. Zhang, and S. Cottenier, Science 351 (2016).
 Hamann et al. (1979) D. R. Hamann, M. Schlüter, and C. Chiang, Phys. Rev. Lett. 43, 1494 (1979).
 Bachelet et al. (1982) G. B. Bachelet, D. R. Hamann, and M. Schlüter, Phys. Rev. B 26, 4199 (1982).
 Vanderbilt (1985) D. Vanderbilt, Phys. Rev. B 32, 8412 (1985).
 Kerker (1980) G. P. Kerker, J. Phys. C: Solid State Phys. 13, L189 (1980).
 Troullier and Martins (1991) N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991).
 Rappe et al. (1990) A. M. Rappe, K. M. Rabe, E. Kaxiras, and J. D. Joannopoulos, Phys. Rev. B 41, 1227 (1990).
 Reis et al. (2003) C. L. Reis, J. M. Pacheco, and J. L. Martins, Phys. Rev. B 68, 155111 (2003).
 Vanderbilt (1990) D. Vanderbilt, Phys. Rev. B 41, 7892 (1990).
 Laasonen et al. (1993) K. Laasonen, A. Pasquarello, R. Car, C. Lee, and D. Vanderbilt, Phys. Rev. B 47, 10142 (1993).
 Kresse and Hafner (1994) G. Kresse and J. Hafner, Journal of Physics: Condensed Matter 6, 8245 (1994).
 Moroni et al. (1997) E. G. Moroni, G. Kresse, J. Hafner, and J. Furthmüller, Phys. Rev. B 56, 15629 (1997).
 Garrity et al. (2014) K. F. Garrity, J. W. Bennett, K. M. Rabe, and D. Vanderbilt, Computational Materials Science 81, 446 (2014).
 Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
 Holzwarth et al. (1998) N. A. W. Holzwarth, G. E. Matthews, A. R. Tackett, and R. B. Dunning, Phys. Rev. B 57, 11827 (1998).
 Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
 Kleinman and Bylander (1982) L. Kleinman and D. M. Bylander, Phys. Rev. Lett. 48, 1425 (1982).
 P. Giannozzi et al. (2009) P. Giannozzi et al., J. Phys.: Condens. Matter 21, 395502 (2009).
 Miwa (2011) K. Miwa, Phys. Rev. B 84, 094304 (2011).
 Baroni et al. (2001) S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
 Hedin (1965) L. Hedin, Phys. Rev. 139, A796 (1965).
 Hybertsen and Louie (1985) M. S. Hybertsen and S. G. Louie, Phys. Rev. Lett. 55, 1418 (1985).
 Hamann (2013) D. R. Hamann, Phys. Rev. B 88, 085117 (2013).
 Dal Corso et al. (1993) A. Dal Corso, S. Baroni, and R. Resta, Phys. Rev. B 47, 3588 (1993).
 Perdew and Zunger (1981) J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
 (32) Http://elk.sourceforge.net/.
 Vinet et al. (1987) P. Vinet, J. R. Smith, J. Ferrante, and J. H. Rose, Phys. Rev. B 35, 1945 (1987).
 Zhang et al. (2010) H. L. Zhang, S. Lu, M. P. J. Punkkinen, Q.M. Hu, B. Johansson, and L. Vitos, Phys. Rev. B 82, 132409 (2010).