Coverage for gpaw/lcao/generate_ngto_augmented.py: 94%

1from typing import Sequence, Union 

2 

3import numpy as np 

4 

5from gpaw.atom.basis import BasisMaker 

6from gpaw.atom.basis import QuasiGaussian 

7from gpaw.atom.radialgd import EquidistantRadialGridDescriptor 

8from gpaw.atom.configurations import parameters, parameters_extra 

9from gpaw.basis_data import BasisFunction 

10from gpaw.basis_data import parse_basis_name 

11 

12# Module for generating basis sets that compose of usual basis sets 

13# augmented with Gaussian type orbital (GTO). 

14# 

15# GTOs are truncated and represented numerically. 

16 

17 

18def create_GTO_dictionary(l: Union[int, str], exponent: float): 

19 """Dictionary representing Gaussian type orbital. 

20 

21 Parameters 

22 ---------- 

23 l 

24 Angular momentum 

25 exponent 

26 Gaussian exponent 

27 """ 

28 return create_CGTO_dictionary(l, [exponent], [1.0]) 

29 

30 

31def create_CGTO_dictionary(l: Union[int, str], 

32 exponents: Sequence[float], 

33 coefficients: Sequence[float]): 

34 """Dictionary representing contracted Gaussian type orbital. 

35 

36 Parameters 

37 ---------- 

38 l 

39 Angular momentum 

40 exponents 

41 Gaussian exponents 

42 coefficients 

43 Gaussian coefficients 

44 """ 

45 if isinstance(l, str): 

46 l = 'spdfghi'.index(l.lower()) 

47 gto = {'angular_momentum': [l], 

48 'exponents': exponents, 

49 'coefficients': [coefficients]} 

50 return gto 

51 

52 

53def read_gaussian_basis_file(fname): 

54 """Read Gaussian basis set file. 

55 

56 This reads only the first element/atom from the file 

57 as separated with line beginning with '*'. 

58 """ 

59 gtos = [] 

60 description = '' 

61 

62 with open(fname) as fd: 

63 line_i = fd.readlines() 

64 

65 i = 0 

66 Ni = len(line_i) 

67 while True: 

68 line = line_i[i].strip() 

69 if line == '' or line[0] == '*': 

70 pass 

71 elif line[0] == '!': 

72 description += '%s\n' % line[1:].strip() 

73 else: 

74 break 

75 i += 1 

76 description = description.strip() 

77 

78 atom = line_i[i].strip().split()[0] 

79 i += 1 

80 while i < Ni: 

81 line = line_i[i] 

82 if line[0] == '*': 

83 break 

84 i += 1 

85 d = line.split() 

86 l = 'spdfghi'.index(d[0].lower()) 

87 Nj = int(d[1]) 

88 alpha_j = [] 

89 coeff_j = [] 

90 for _ in range(Nj): 

91 line = line_i[i] 

92 d = line.split() 

93 alpha = float(d[0].replace('D', 'E')) 

94 coeff = float(d[1].replace('D', 'E')) 

95 alpha_j.append(alpha) 

96 coeff_j.append(coeff) 

97 i += 1 

98 gto = create_CGTO_dictionary(l, alpha_j, coeff_j) 

99 gtos.append(gto) 

100 

101 return atom, description, gtos 

102 

103 

104def get_ngto(rgd, l, alpha, rcut): 

105 gaussian = QuasiGaussian(alpha, rcut) 

106 psi_g = gaussian(rgd.r_g) * rgd.r_g**l 

107 norm = np.sum(rgd.dr_g * (rgd.r_g * psi_g)**2)**.5 

108 psi_g /= norm 

109 return psi_g 

110 

111 

112def create_ngto(rgd, l, alpha, rmax, tol): 

113 # Get NGTO with the initial (large) rcut=rmax 

114 psiref_g = get_ngto(rgd, l, alpha, rmax) 

115 

116 # Make rcut smaller 

117 

118 # Guess initial rcut where we are close to the tolerance 

119 i = np.where(psiref_g > tol)[0][-1] 

120 rcut = rgd.r_g[i] 

121 psi_g = get_ngto(rgd, l, alpha, rcut) 

122 err = np.max(np.absolute(psi_g - psiref_g)) 

123 

124 # Increase/decrease rcut to find the smallest rcut 

125 # that yields error within the tolerance 

126 if err > tol: 

127 # Increase rcut -> decrease err 

128 for i in range(i, len(rgd.r_g)): 

129 rcut = rgd.r_g[i] 

130 psi_g = get_ngto(rgd, l, alpha, rcut) 

131 err = np.max(np.absolute(psi_g - psiref_g)) 

132 if err < tol: 

133 break 

134 else: 

135 # Decrease rcut -> increase err 

136 for i in range(i, 0, -1): 

137 rcut = rgd.r_g[i] 

138 psi_g = get_ngto(rgd, l, alpha, rcut) 

139 err = np.max(np.absolute(psi_g - psiref_g)) 

140 if err > tol: 

141 i += 1 

142 break 

143 

144 # Construct NGTO with the found rcut 

145 rcut = rgd.r_g[i] 

146 psi_g = get_ngto(rgd, l, alpha, rcut) 

147 

148 # Change norm (maybe unnecessary) 

149 psi_g = psi_g[:(i + 1)] * 0.5 

150 

151 return psi_g 

152 

153 

154def add_ngto(basis, l, coeff_j, alpha_j, tol, label): 

155 rgd = basis.get_grid_descriptor() 

156 rmax = rgd.r_g[-1] 

157 

158 # Create linear combination of NGTO's 

159 psi_g = np.zeros(rgd.r_g.shape) 

160 i_max = 0 

161 for coeff, alpha in zip(coeff_j, alpha_j): 

162 contrib = coeff * create_ngto(rgd, l, alpha, rmax, tol) 

163 i = contrib.size 

164 i_max = max(i, i_max) 

165 psi_g[0:i] += contrib 

166 

167 psi_g = psi_g[0:i_max] 

168 rcut = rgd.r_g[i_max] 

169 

170 # Create associated basis function 

171 bf = BasisFunction(None, l, rcut, psi_g, label) 

172 basis.bf_j.append(bf) 

173 

174 

175def generate_nao_ngto_basis(atom, *, xc, nao, name, 

176 gtos, gto_description=None, 

177 rmax=100.0, tol=0.001): 

178 from dataclasses import replace 

179 # Choose basis sets without semi-core states XXXXXX 

180 if atom == 'Ag': 

181 name = '11.%s' % name 

182 p = parameters_extra 

183 else: 

184 p = parameters 

185 

186 # Generate nao basis 

187 zetacount, polarizationcount = parse_basis_name(nao) 

188 bm = BasisMaker.from_symbol( 

189 atom, name=name, gtxt=None, xc=xc, 

190 generator_run_kwargs=dict(write_xml=False, **p[atom])) 

191 basis = bm.generate(zetacount, polarizationcount, txt=None) 

192 

193 # Increase basis function max radius 

194 # XXX why are we doing this? 

195 assert isinstance(basis.rgd, EquidistantRadialGridDescriptor) 

196 h = basis.rgd.dr_g[0] 

197 assert basis.rgd.r_g[0] == 0.0 

198 N = int(rmax / h) + 1 

199 basis = replace(basis, rgd=EquidistantRadialGridDescriptor(h, N)) 

200 

201 # Add NGTOs 

202 description = [] 

203 msg = 'Augmented with NGTOs' 

204 description.append(msg) 

205 description.append('=' * len(msg)) 

206 description.append('') 

207 if gto_description is not None: 

208 description.append(gto_description) 

209 description.append('') 

210 description.append('NGTO truncation tolerance: %f' % tol) 

211 description.append('Functions: NGTO(l,coeff*alpha + ...)') 

212 

213 for gto in gtos: 

214 assert len(gto['angular_momentum']) == 1 

215 l = gto['angular_momentum'][0] 

216 alpha_j = gto['exponents'] 

217 # Float conversion 

218 alpha_j = [float(a) for a in alpha_j] 

219 for coeff_j in gto['coefficients']: 

220 assert len(alpha_j) == len(coeff_j) 

221 # Float conversion 

222 coeff_j = [float(c) for c in coeff_j] 

223 coeff_alpha_list = [f'{c:+.3f}*{a:.3f}' 

224 for c, a in zip(coeff_j, alpha_j)] 

225 coeff_alpha_label = ''.join(coeff_alpha_list[0:3]) 

226 if len(coeff_alpha_list) > 3: 

227 coeff_alpha_label += '+...' 

228 ngtolabel = 'NGTO({},{})'.format('spdfghi'[l], coeff_alpha_label) 

229 description.append(' ' + ngtolabel) 

230 add_ngto(basis, l, coeff_j, alpha_j, tol, ngtolabel) 

231 

232 basis = replace( 

233 basis, 

234 generatordata=basis.generatordata + '\n\n' + '\n'.join(description)) 

235 

236 basis.write_xml()