Coverage for gpaw/ah.py: 97%

1"""Appelbaum-Hamann local pseudo-potential for silicon. 

2 

3See:: 

4 

5 Self-Consistent Pseudopotential for Si 

6 Joel A. Appelbaum and D. R. Hamann 

7 PRB 8, 1777 (1973) 

8 

9""" 

10 

11import numpy as np 

12 

13from gpaw.setup import BaseSetup 

14from gpaw.spline import Spline 

15from gpaw.basis_data import Basis 

16 

17 

18class AppelbaumHamann(BaseSetup): 

19 symbol = 'Si' 

20 

21 def __init__(self, alpha=0.6102, v1=3.042, v2=-1.372): 

22 self.alpha = alpha 

23 self.v1 = v1 

24 self.v2 = v2 

25 

26 self.E = 0.0 

27 self.Z = 14 

28 self.Nc = 10 

29 self.Nv = 4 

30 self.nao = None 

31 nullspline = Spline.from_data(0, 0.5, [0., 0., 0.]) 

32 self.pt_j = [nullspline] 

33 self.ni = 1 

34 self.l_j = [0] 

35 self.f_j = [4] 

36 self.n_j = [1] 

37 self.nct = nullspline 

38 self.tauct = nullspline 

39 self.Nct = 0.0 

40 rc = 4.0 

41 r2_g = np.linspace(0, rc, 100)**2 

42 x_g = np.exp(-alpha * r2_g) 

43 self.ghat_l = [ 

44 Spline.from_data(0, rc, 4 * alpha**1.5 / np.pi**0.5 * x_g), 

45 ] 

46 self.vbar = Spline.from_data( 

47 0, rc, 2 * np.pi**0.5 * (v1 + v2 * r2_g) * x_g, 

48 ) 

49 self.Delta_pL = np.zeros((1, 1)) 

50 self.Delta_iiL = np.zeros((1, 1, 1)) 

51 self.Delta0 = -4 / (4 * np.pi)**0.5 

52 self.ExxC = 0.0 

53 self.lmax = 0 

54 self.K_p = self.M_p = self.MB_p = self.X_p = self.N0_p = np.zeros(1) 

55 self.M_pp = np.zeros((1, 1)) 

56 self.Kc = 0.0 

57 self.MB = 0.0 

58 self.M = 0.0 

59 self.xc_correction = None 

60 self.dO_ii = np.zeros((1, 1)) 

61 self.type = 'ah' 

62 self.fingerprint = None 

63 

64 def build(self, basis): 

65 if basis is None: 

66 basis = Basis.find('Si', 'sz(dzp)') 

67 elif isinstance(basis, str): 

68 basis = Basis.find('Si', basis) 

69 

70 self.basis = basis 

71 self.basis_functions_J = self.basis.tosplines() 

72 self.nao = self.basis.nao 

73 

74 def print_info(self, text): 

75 text('Appelbaum-Hamann pseudo potential') 

76 

77 def calculate_initial_occupation_numbers(self, magmom, hund, charge, 

78 nspins): 

79 assert nspins == 1 

80 return np.array([(2.0, 2.0 / 3, 2.0 / 3, 2.0 / 3)]) 

81 

82 def initialize_density_matrix(self, f_si): 

83 return np.zeros((len(f_si), 1)) 

84 

85 def get_default_nbands(self): 

86 return 3