Coverage for gpaw/test/radial/test_yukawa_radial.py: 100%

1"""Small test for local parts using yukawa potential""" 

2import pytest 

3from math import factorial as fac 

4from numpy import exp, sqrt, pi, absolute 

5from gpaw.atom.radialgd import AERadialGridDescriptor 

6 

7# Values from Rico, Lopez, Ramirez, Ema, Theor Chem Acc (2013) 132:1304 

8# Table 1 and 2. 

9# 

10# Normalized STO: N = n + L : lambda (gamma) = 0.34 

11# n n' L Zeta Zeta' I(n,n',L,Zeta,Zeta') 

12# 8 6 4 8.4 3.5 0.793087136533672 

13# 2 1 5 8.4 3.5 0.271337849250231 

14# 3 4 6 8.4 3.5 0.170346973782779 

15# 3 2 8 8.4 3.5 0.141826311951003 

16# 2 1 10 8.4 3.5 9.005003835363871(-2) 

17 

18# 

19# Zeta = lambda (gamma) 

20# n = 3, L=5, n'=2, L'=5, zeta'=3.5 

21# 0.10 9.503545229396430(-6) 

22# 0.25 5.869262722600324(-4) 

23# 0.50 1.194846721531269(-2) 

24# 0.75 5.829196062489732(-2) 

25# 1.00 0.153003818405446 

26 

27 

28@pytest.fixture 

29def rgd(): 

30 N = 1200 

31 beta = 0.4 * 600 / N 

32 return AERadialGridDescriptor(beta / N, 1.0 / N, N) 

33 

34 

35basic_sto_list = [ # See table 1 of Rico, Lopez, Ramirez, Ema, 

36 # Theor Chem Acc (2013) 132:1304 

37 { 

38 'n': [8, 6], # n, n' 

39 'Zeta': [8.4, 3.5], # Z, Z' 

40 'L': 4, # L 

41 'I': 0.793087136533672 # [X_LM^n|X_LM^n'] 

42 }, 

43 # { 

44 # 'n': [2, 1], 

45 # 'Zeta': [8.4, 3.5], 

46 # 'L': 5, 

47 # 'I': 0.271337849250231 

48 # }, 

49 # { # values are completely of 

50 # 'n': [3, 4], 

51 # 'Zeta': [8.4, 3.5], 

52 # 'L': 6, 

53 # 'I': 0.170346973782779 

54 # }, 

55 # { 

56 # 'n': [3, 2], # n, n' 

57 # 'Zeta': [8.4, 3.5], # Z, Z' 

58 # 'L': 8, # L 

59 # 'I': 0.141826311951003 # [X_LM^n|X_LM^n'] 

60 # }, 

61 # { 

62 # 'n': [2, 1], 

63 # 'Zeta': [8.4, 3.5], 

64 # 'L': 10, 

65 # 'I': 9.005003835363871e-2 

66 # } 

67] 

68 

69gamma_sto_list = [ # See table 2 of Rico, Lopez, Ramirez, Ema, 

70 # Theor Chem Acc (2013) 132:1304 

71 { 

72 'zeta': 0.10, 

73 'I': 9.503545229396430e-6 

74 }, 

75 # { 

76 # 'zeta': 0.25, 

77 # 'I': 5.869262722600324e-4 

78 # }, 

79 # { 

80 # 'zeta': 0.50, 

81 # 'I': 1.194846721531269e-2 

82 # }, 

83 # { 

84 # 'zeta': 0.75, 

85 # 'I': 5.829196062489732e-2 

86 # }, 

87 # { 

88 # 'zeta': 1.00, 

89 # 'I': 0.153003818405446 

90 # }, 

91] 

92 

93 

94def radial_sto(n, zeta, l, r): 

95 """Build radial part of slater type orbital""" 

96# Stolen from gpaw all_electron.py (intialize_wave_functions) 

97# 

98# STOs are defined as 

99# 

100# N Y_l^m r^{n-1} exp(-zeta * r) 

101# 

102# N = (2 zeta)^n sqrt(2 zeta / (2n)!) 

103 assert n > 0 

104 radial = r**(n + l - 1.0) 

105 radial *= exp(-zeta * r) 

106 N = (2 * zeta)**(n + l) * sqrt(2 * zeta / fac(2 * (n + l))) 

107 radial *= N # perhaps also do Y_l^m normalization 

108 return radial 

109 

110 

111def test_different_sto(rgd): 

112 """Check integration of two different STO.""" 

113 gamma = 0.34 

114 for test_sto in basic_sto_list: 

115 wave_functions = [] 

116 l_sto = test_sto['L'] 

117 for n, zeta in zip(test_sto['n'], test_sto['Zeta']): 

118 radial = radial_sto(n, zeta, l_sto, rgd.r_g) 

119 norm = sqrt((2 * l_sto + 1) / (4 * pi)) # m = 0 

120 radial *= norm 

121 wave_functions.append(radial) 

122 yukawa_wf = rgd.yukawa(wave_functions[1], l_sto, gamma) 

123 test_int = rgd.integrate(wave_functions[0] * yukawa_wf, -1) / \ 

124 (2 * l_sto + 1) 

125# print(u"{:7.5e}||{:7.5e}||{:7.5e}".format(test_sto['I'], I, 

126# absolute(I - test_sto['I']) * 100.0 / test_sto['I'])) 

127 assert absolute(test_int - test_sto['I']) / test_sto['I'] < 0.001 

128 

129 

130def test_poisson(rgd): 

131 """Check integration against rgd.poisson.""" 

132 gamma = 1e-5 

133 for test_sto in basic_sto_list: 

134 wave_functions = [] 

135 l_sto = test_sto['L'] 

136 for n, zeta in zip(test_sto['n'], test_sto['Zeta']): 

137 radial = radial_sto(n, zeta, l_sto, rgd.r_g) 

138 norm = sqrt((2 * l_sto + 1) / (4 * pi)) # m = 0 

139 radial *= norm 

140 wave_functions.append(radial) 

141 yukawa_wf = rgd.yukawa(wave_functions[1], l_sto, gamma) 

142 poisson_wf = rgd.poisson(wave_functions[1], l_sto) 

143 I_p = rgd.integrate(wave_functions[0] * poisson_wf, -1) 

144 I_y = rgd.integrate(wave_functions[0] * yukawa_wf, -1) 

145 assert absolute(I_y - I_p) / I_p < 5e-4 

146 

147 

148def test_same_sto(rgd): 

149 """Check integration of nearly identical STO.""" 

150 n = 3 

151 l_sto = 5 

152 np = 2 

153 zetap = 3.5 

154 for test_sto in gamma_sto_list: 

155 zeta = test_sto['zeta'] 

156 radial = radial_sto(n, zeta, l_sto, rgd.r_g) 

157 norm = sqrt((2 * l_sto + 1) / (4 * pi)) # m = 0 

158 radial *= norm 

159 radial2 = radial_sto(np, zetap, l_sto, rgd.r_g) 

160 radial2 *= norm 

161 yukawa_wf = rgd.yukawa(radial2, l_sto, zeta) 

162 test_int = rgd.integrate(radial * yukawa_wf, -1) / (2 * l_sto + 1) 

163 assert absolute(test_int - test_sto['I']) / test_sto['I'] < 0.001