Generation of pseudopotentials and optimization of basis sets for use with SeqQuest code

 

Final report

 

 

 

Introduction

 

Ab initio computations of periodical (1D, 2D, and 3D) systems can be made much more efficient with application of pseudopotentials. The pseudopotentials replace nuclei and inner electrons. This reduces the number of electrons, which have to be described in simulations, and reduces size of basis sets needed to expand one-electron wave functions leading to significant savings of computational time.

One of codes used in Materials and Process Simulation Center (MSC) is SeqQuest [1]. Pseudopotentials cannot be inputted to SeqQuest in usually applied parameterized form. SeqQuest code requires input of pseudopotentials in a specific format as numerical arrays. The last requirement makes more difficult application of pseudopotentials, described in literature. However, it is still possible to generate input data for pseudopotentials from published pseudopotentials, or from newly created pseudopotentials. The second approach has an advantage that one can generate pseudopotentials with the same computational approach, which then will be used for actual simulations. The most of published pseudopotentials were generated within Hartree-Fock approach. However, most of ab initio calculations are made with different flavors of Density Functional Theory. In particular, SeqQuest allows to produce computations in Local Density Approximation (LDA) [2] and with Perdew-Berke-Erzenkorf (PBE) exchange-correlation functional [3]. Pseudopotentials based on these functionals can be generated with the code written in Fritz-Haber Institute [4].

SeqQuest code employs local basis sets from Gaussian-type atomic orbitals. One-electron wave functions are built as linear combinations of the atomic orbitals (LCAO). Before ab initio calculations with local basis sets can be performed, one has to optimize the basis sets for all atoms in a system under study.

Before this project has started Hamman type [5] pseudopotentials were generated and basis sets were optimized for majority of s-, p- and d- elements using LDA functional. For some elements the same kind of pseudopotentials and basis sets were prepared with PBE functional as well. The goal of this project was to generate pseudopotentials and to optimize basis sets for the rest of s-, p- and d- atoms. Simultaneously, some of obtained pseudopotentials and basis sets were tested by calculations of Equations of State for different phases of several metals ( Sr, Zr, Pb, Y, Ru, Co ), metal oxides ( SrO, BaO, PbO ) and lead titanate perovskite.

 

Applied Methodology

 

Pseudopotential generation

 

To generate Hamman type [5] norm-conserving pseudopotentials I applied the package of computer programs written in Fritz-Haber Institute [4]. I generated pseudopotentials for all s-, p-, and d- elements (excluding noble gases), which have not been generated before this job started. Details, how to use this package at MSC can be found at [6]. Parameters used to generate the pseudopotentials are collected in Table 1 and also can be found at [7].

As a reference state I used neutral atoms in the ground state for halogen atoms and positively charged ions for all other cases. Charges of the ions sometimes were fractional. Such a choice was usually made for alkali and alkali earth elements. The choice to use positive ions as reference states for metals was made because in the most of chemical compounds atoms of these elements have positive charge. It is obvious, that the made choice is closer to situations, in which we need a good description of metal atoms and will not harm computations of metallic solids, where atoms of these elements are neutral.

Leaving outer core electron shells outside pseudopotentials allows significantly better description of bond lengths and binding energies. One-electron wavefunctions of these electrons are still sensitive to chemical environment of atoms because valence d-orbitals have large magnitude in the same space, where wavefunction of the outer core electrons are. Valence electrons of alkali and alkali earth atoms usually are nearly completely ionized in chemical compounds. In such situations p-orbitals of the outer core electron shells are in direct contact with neighboring atoms. Therefore I tried, where it was possible, to leave p-orbitals from outer core electronic shells of atoms unfrozen. This approach was applied for transition elements and alkali and alkali earth elements, except of the first row atoms, where there are not such orbitals. It is also impossible to apply this approach for p-elements, because the employed pseudopotential generation code does not allow for more than one free orbital in each orbital moment component of the pseudopotential. For example, generated pseudopotentials for atoms from Sr till Cd have free 4s, 3p, and 3d orbitals, and pseudopotential for atoms from Tl till Bi have free 5s, 4d, and 5p orbitals.

Matching radii determine distances from nuclei where pseudo-orbitals and their logarifmic derivatives have to coincide with atomic orbitals. I applied usual choice for the matching radii. Their values have to be between the last nodes of the valence orbitals and its last maximums. Also, the radii cannot be chosen too close to the last node of valence one-electron orbitals to avoid irregular behavior of pseudopotentials. Radii, used for generation of pseudopotentials are collected in Table 1 and at web-site [7].

Exchange-correlation interactions are not linear in respect to an electron density. This non-linearity is more pronounced in open-shell systems and for alkali and alkali earth metal atoms [8]. The regular semi-local pseudopotentials have only linear interaction with pseudo-wave-functions. The non-linear part of interactions between core and valence electrons is not accounted for. To correct the situation, a model core electron density ρm.core can be introduced [8] with the condition, that it coincides with core electron density at distances from nuclei beyond of cut-off radius rnlc. Then density ρ used for calculation of exchange-correlation energy is

 

ρ = ρm.core + ρpseudo ,

 

where ρpseudo is electron density calculated with pseudo-wavefunctions. Applied cut-off radii are shown in the Table 1 and at [7].

 

Basis set optimization

 

Application of new pseudopotentials in ab initio codes employing local Gaussian-type basis sets requires a suitable basis sets. Such basis sets can be obtained by fitting atomic pseudoorbitals, built during generation of pseudopotentials. Another approach is directly to optimize basis sets in a representative system. The technique utilized in this project combined both these approaches. At the first, atomic pseudoorbitals were fitted with several Gaussian atomic orbitals. Then the obtained expansions into Gaussian orbitals were put to input for SeqQuest as contracted basis functions. One or two most diffuse Gaussian exponents and their expansion coefficients in atomic orbitals were optimized in elementary solids. At the final step the most diffuse Gaussian was set as a separate basis function as well. In some cases, especially when the most diffuse Gaussian function has very small exponent parameter, two most diffuse Gaussians were splitted as a separate basis function and exponent and expansion coefficients were optimized.

Basis sets were optimized in elementary solids. In the majority of cases I used the most stable phase of the elementary for basis set optimization. However, there were several exceptions, where more symmetric phases with smaller number of atoms in unit cell were used. Basis sets for magnetic metals ( Co .. Fe) were optimized in ferroelectric states. Also, basis sets for halide atoms were optimized in potassium halides KX ( X= Cl, Br, I ). Lattice constants in all employed crystals had experimental values taken from [9].

Basis set optimization was performed with help a separate optimization code, written previously for different purposes. I modified interface of the code to make the code to work with SeqQuest. This code performs Conjugated Gradients optimizations. It basically contains two drivers. The one drives computations of gradients with numerical differentiation over 3 or 5 points for each variable. The second driver controls Conjugated Gradient optimization. The code creates necessary inputs for another code ( for example, for SeqQuest ), using prepared templates. Then it runs a single point calculation of another code and extracts final energy from output of this code. The extracted energy is used for calculation of gradients and in energy optimization procedure.

 

Optimization of occupation numbers in reference state

 

SeqQuest uses difference between actual electron density and electron density calculated as a sum of reference atomic densities to reach significant acceleration of computations. The reference atomic densities are given in SeqQuest input through occupation numbers of atomic orbitals. The criterion for reference state quality is the non-spherical part of electrostatic energy of interaction of electron density ( two-electron interaction). When this contribution is small, approximations made in energy computation in SeqQuest work better. Therefore, to improve performance of computations with SeqQuest it was necessary to minimize the non-spherical electrostatic energy with respect to the occupation numbers. This optimization was done with help of the same optimization code, which was used for optimization of basis sets. A few modifications needed for this job were made to the code. The first, modification was introduced into interface of the code to extract non-spherical electrostatic energy from the code output instead of the total energy. The second modification was necessary to maintain sum of the occupation numbers constant. This was needed to keep atom charge constant during the optimization.

 

Results of pseudopotential generation

 

Results of generation of pseudopotentials can be found at web-site [7] and in the directory /source/seqquest/atoms_pbe in MSCs computer system . These web-sites contain links to .atm files with pseudopotentials, optimized basis sets, and occupation numbers. These files can be also found in *_pbe_h subdirectories of /source/seqquest/atoms_pbe directory. The web-site [7] contain also detail information, how pseudopotentials were generated.

 

Testing of pseudopotentials and basis sets

 

Method of testing

 

Several pseudopotentials ( Sr, Zr, Pb, Y, Ru, Co ) were tested in more details by computations of Equations of States ( EoS) for Sr, Zr, Pb, Y, Ru, and Co metals, SrO and ZrO2 metal oxides, and PbTiO3 perovskite. EoS were calculated for bcc, fcc, simple cubic, diamond, a15, and hcp phases of Sr, Zr, Pb, Y, and Ru metals, for bcc, fcc, simple cubic, diamond, and hcp phases of ferromagnetic ( fully spin-polarized ) Co, for rock-salt phase of SrO, for tetragonal phase of PbO, for fluorite cubic, tetragonal, monoclinic and two orthorhombic phases of ZrO2, cubic, tetragonal and rhombohedral perovskite phases of PbTiO3. EoS computations were performed with SeqQuest code with newly generated pseudopotentials for Sr, Zr, Pb, Y, Ru, Co and previously generated pseudopotentials for Ti [10] and O [11] . Integration grids and Monkhorst- net applied for in the EoS calculations for each material and phase is given at the web site [7]. Sizes of integration grids were made to reach good convergence of energy, stresses, and forces, and keep volume per a node of the grid approximately the same in the different phases of the same materials. To calculate EoS I first calculated equilibrium lattice constants and atomic positions in unit cell at zero pressure. Then lattice constants and atomic positions were optimized at several different volumes of crystal unit cell. At the end EoS were written in form of energy as function of Wigner-Seitz radius, where Wigner-Seitz radius is a radius of a sphere with volume equal to volume per formula unit in the respective crystals ( volume per atom in metals, volume per SrO, ZrO2 in the oxides, volume per PbTiO3 in the studied perovskite).

On the next stage obtained EoS were fitted with analytical expression in form of Roses EoS [12] :

 

(1)

 

where E0 is the atomization energy, aws is Wigner-Seitz radius, aws(0) is its equilibrium value at zero pressure, s is scale factor, and k is a parameter. All parameters of the Roses EoS were optimized with the Least Square Deviation technique to reproduce crystal energies at sets of volumes where crystal energies were calculated with SeqQuest code. Representation of EoS in the analytical form allows for further analysis of material properties. Bulk modulus B and its derivative B over pressure can be calculated as

 

(2)

 

. (3)

 

 

The function of pressure versus Wigner-Seitz radius can be calculated as

 

, (4)

 

and then for calculation of enthalpy one can use its definition:

 

. (5)

Together equations 4 and 5 define enthalpy function in respect to pressure in parametric form. The minimal possible pressure from equation 4 gives the pressure at which the crystal ruptures (rupture pressure). For the material, for which several phases were calculated, it was possible to consider relative stability of phases and to estimate phase transition pressures at zero temperature. Points of phase transition can be found by solving a system of equations

 

 

(6)

 

in respect to Wigner-Seitz radii aws,1 and aws,2 for both phases. The described fitting and analysis of Equations of States was done with Mathematica code [7].

 

 

Results of testing

 

Results of the described studies can be found through links from the web-site [7], where all described values tabulated and images of calculations with Mathematica code are available.

Lattice constants for calculated metals and oxides are in good agreement with experiment. They are slightly overestimated usually, except Sr. For non-cubic phases of all crystals c/a ratio was usually overestimated as well. The best coincidence with experimental values was found for yttrium metal. This is typical results for PBE functional. Calculated bulk moduli of the crystals are in very well consistent with experimental data.

 

Conclusion

 

In this project I generated pseudopotentials for 46 s-, p- and d- elements, optimized basis sets and occupation numbers of reference states for these elements . The pseudopotentials were tested by computations of basic crystal properties.

 

 

 

 

 

 

 

 

Table 1. Parameters, used to generate pseudopotentials.

 

Atomic No

Element

Core configuration

Reference

Valence

configuration

Matching Radii, Bohrs

Model density

Cut-off radius,

Bohrs

s

p

d

f

3

Li

0s2

1s0.21p0

1.63

0.35

-

-

0.20

17

Cl

0s21s21p6

2s22p5

0.744

0.78

0.84

-

0.235

19

K

0s21s21p62s2

2p63s0.22d0

2.1

0.65

0.7

-

0.40

20

Ca

0s21s21p62s2

2p63s0.22d0

1.77

0.59

0.3

-

0.37

21

Sc

0s21s21p62s2

2p63s12d1

1.72

0.58

0.43

-

0.57

23

V

0s21s21p62s2

2p63s12d3

1.40

0.45

0.30

-

0.40

24

Cr

0s21s21p62s2

2p63s12d4

1.47

0.47

0.33

-

0.46

25

Mn

0s21s21p62s2

2p63s12d5

1.41

0.45

0.30

-

0.40

26

Fe

0s21s21p62s2

2p63s22d5

1.20

0.40

0.25

-

0.40

27

Co

0s21s21p62s2

2p62d63s2

1.20

0.40

0.25

-

0.40

29

Cu

0s21s21p62s2

2p63s12d9

1.20

0.40

0.25

-

0.40

30

Zn

0s21s21p62s2

2p63s12d10

1.15

0.35

0.23

-

0.35

31

Ga

0s21s21p62s22p6

3s12d103p1

1.20

1.10

0.10

-

0.265

33

As

0s21s21p62s22p6

3s22d103p2

1.05

1.05

0.10

-

0.25

34

Se

0s21s21p62s22p6

3s12d103p3

0.98

1.03

0.10

-

0.245

35

Br

0s21s21p62s22p6

3s22d103p5

0.98

0.99

0.10

-

0.235

37

Rb

0s21s21p62s22p63s22d10

3p63d04s0.2

2.20

0.80

1.35

-

0.8

38

Sr

0s21s21p62s22p63s23p62d10

3p63d0.14s0.5

2.00

0.8

1.0

-

2.5

39

Y

0s21s21p62s22p63s23p62d10

3p63d14s1

1.88

0.74

1.04

-

1.20

40

Zr

0s21s21p62s22p63s23p62d10

3p63d24s2

1.14

1.53

0.635

-

2.02

41

Nb

0s21s21p62s22p63s23p62d10

3p63d34s2

1.85

0.60

0.85

-

1.6

43

Tc

0s21s21p62s22p63s23p62d10

3p63d54s2

1.50

0.55

0.70

-

1.

44

Ru

0s21s21p62s22p63s23p62d10

3p63d64s2

1.47

0.54

0.69

-

0.9

45

Rh

0s21s21p62s22p63s23p62d10

3p63d74s2

1.45

0.55

0.65

-

0.9

46

Pd

0s21s21p62s22p63s23p62d10

3p63d84s2

1.40

0.50

0.60

-

0.85

47

Ag

0s21s21p62s22p63s23p62d10

3p63d94s2

1.37

0.50

0.57

-

0.80

48

Cd

0s21s21p62s22p63s23p62d10

3p63d104s2

1.34

0.47

0.54

-

0.75

49

In

0s21s21p62s22p63s22d103p6

3d104s14p1 

1.20

1.50

0.50

-

0.75

50

Sn

0s21s21p62s22p63s22d103p6

3d104s1.84p1.2 

1.20

1.50

0.50

-

0.40

51

Sb

0s21s21p62s22p63s22d103p6

3d104s24p2

1.10

1.34

0.51

-

0.43

53

I

0s21s21p62s22p63s22d103p6

3d104s24p5

1.01

1.20

0.47

-

0.41

55

Cs

0s21s21p62s22p63s22d103p63d104s2

4p65s0.54d0

2.50

1.05

1.35

 

0.8

56

Ba

0s21s21p62s22p63s22d103p63d104s2

4p65s0.54d0

2.20

0.95

1.25

-

0.7

57

La

0s21s21p62s22p63s22d103p63d104s2

4p65s14d1

2.15

0.95

1.30

-

2.05

72

Hf

0s21s21p62s22p63s22d103p64s2 3d104p63f14

5s14d2

1.60

1.90

0.90

1.70

0.69

73

Ta

0s21s21p62s22p63s22d103p64s2 3d104p63f14

5s14d3

1.55

1.85

0.85

1.25

0.67

74

W

0s21s21p62s22p63s22d103p64s2 3d104p63f14

5s14d4

1.50

1.80

0.85

1.25

0.67

75

Re

0s21s21p62s22p63s22d103p64s2 3d104p63f14

5s14d5

1.40

1.70

0.80

1.25

0.66

76

Os

0s21s21p62s22p63s22d103p64s2 3d104p63f14

5s14d6

1.38

1.60

0.75

1.25

0.64

77

Ir

0s21s21p62s22p63s22d103p64s2 3d104p63f14

5s14d7

1.33

1.65

0.70

1.20

0.62

79

Au

0s21s21p62s22p63s22d103p64s2 3d104p63f14

5s14d9

1.30

1.60

0.70

1.10

0.598

80

Hg

0s21s21p62s22p63s22d103p64s2 3d104p63f14

5s14d10

1.25

1.60

0.60

1.10

0.591

81

Tl

0s21s21p62s22p63s22d103p64s2 3d104p63f14

5s14d105p1

1.20

1.50

0.60

1.10

0.58

82

Pb

0s21s21p62s22p63s22d103p64s2 3d104p63f14

5s14d105p3

2.00

0.8

1.0

1.53

3.5

83

Bi

0s21s21p62s22p63s22d103p64s2 3d104p63f14

5s14d105p3

1.10

1.30

0.60

1.00

0.6

 

 

 

References:

 

1.      http://www.wag.caltech.edu/group/quest

2. D. M. Ceperley, B.J. Alder, Phys.Rev.Lett. 45 (1980), 566.

3.J.P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77 (1996), 3865.

4. M. Fuchs, M. Scheffler, Comput. Phys. Commun., Comput. Phys. Commun. 119, 67-98 (1999); http://www.fhi-berlin.mpg.de/th/fhi98md/fhi98PP/

5. D. R. Hamann, Phys. Rev. B 40, 2980 (1989).

6. http://www.wag.caltech.edu/group/pseudopotentials.

7. http://www.wag.caltech.edu/home/heifets/MSCpseudo/Reports/PseudopotentialsProjectHome.htm

8. S.G.Louie, S.Froyen, M.L. Cohen, Phys.Rev.B 26 (1982), 1738.

9. CRC Handbook of Chemistry and Physics, 81th edition, Ed. in chief D. R. Lide, CRC Press, 2000 ( Boca Raton, London, New York, Washington, DC).

10. File /source/seqquest/atoms_pbe/ti_pbe/ti_pbe.atm

11. File /source/seqquest/atoms_pbe/o_tm_pbe/o_tm_pbe.atm

12. P. Vinet, J. Ferrante, J.R. Smith, J.H. Rose, J. Phys. C: Sol. St. Phys. 19 (1986), L467-L473