Coverage for gpaw/atom/check.py: 12%

1import optparse 

2import traceback 

3from pathlib import Path 

4from typing import Any 

5 

6import ase.db 

7import numpy as np 

8from ase import Atoms 

9from ase.data import atomic_numbers, covalent_radii 

10from ase.optimize import BFGS 

11from gpaw import GPAW, PW, KohnShamConvergenceError 

12 

13cutoffs = [200, 250, 300, 400, 500, 600, 700, 800, 1500] 

14 

15 

16def check(con, name: str, lcao=True): 

17 params: dict[str, Any] = dict(xc='PBE', 

18 symmetry='off') 

19 

20 if '.' in name: 

21 symbol, _, setup = name.partition('.') 

22 params['setups'] = setup 

23 else: 

24 symbol = name 

25 

26 for h in [0.16, 0.18, 0.2]: 

27 a = 16 * h 

28 atoms = Atoms(symbol, cell=(a, a, 2 * a), pbc=True) 

29 atoms.calc = GPAW(mode='fd', 

30 h=h, 

31 txt=f'{name}-eggbox-{h:.2f}.txt', 

32 **params) 

33 energies = [] 

34 try: 

35 for i in range(4): 

36 energies.append(atoms.get_potential_energy()) 

37 atoms.positions += h / 6 

38 except KohnShamConvergenceError: 

39 continue 

40 atoms.positions -= h / 6 

41 eegg = np.ptp(energies) 

42 con.write(atoms, name=name, test='eggbox', eegg=eegg, h=h) 

43 

44 a = 4.0 

45 atoms = Atoms(symbol, cell=(a, a, a), pbc=True) 

46 for ecut in cutoffs: 

47 atoms.calc = GPAW(mode=PW(ecut), 

48 txt=f'{name}-pw-{ecut:04}.txt', 

49 **params) 

50 try: 

51 atoms.get_potential_energy() 

52 except KohnShamConvergenceError: 

53 continue 

54 con.write(atoms, name=name, test='pw1', ecut=ecut) 

55 

56 for g in [20, 24, 28]: 

57 atoms.calc = GPAW(mode='fd', 

58 gpts=(g, g, g), 

59 txt=f'{name}-fd-{g}.txt', 

60 **params) 

61 try: 

62 atoms.get_potential_energy() 

63 except KohnShamConvergenceError: 

64 continue 

65 con.write(atoms, name=name, test='fd1', gpts=g) 

66 

67 if lcao: 

68 for g in [20, 24, 28]: 

69 id = con.reserve(name=name, test='lcao1', gpts=g) 

70 if id is None: 

71 continue 

72 atoms.calc = GPAW(gpts=(g, g, g), 

73 mode='lcao', basis='dzp', 

74 txt=f'{name}-lcao-{g}.txt', 

75 **params) 

76 atoms.get_potential_energy() 

77 con.write(atoms, name=name, test='lcao1', gpts=g) 

78 del con[id] 

79 

80 Z = atomic_numbers[symbol] 

81 d = 2 * covalent_radii[Z] 

82 atoms = Atoms(symbol * 2, cell=(a, a, 2 * a), pbc=True, 

83 positions=[(0, 0, 0), (0, 0, d)]) 

84 for ecut in cutoffs: 

85 atoms.calc = GPAW(mode=PW(ecut), 

86 txt=f'{name}2-pw-{ecut:04}.txt', 

87 **params) 

88 try: 

89 atoms.get_potential_energy() 

90 except KohnShamConvergenceError: 

91 continue 

92 con.write(atoms, name=name, test='pw2', ecut=ecut) 

93 

94 if 0: 

95 id = con.reserve(name=name, test='relax') 

96 if id is not None: 

97 atoms.calc = GPAW(mode=PW(1500), 

98 txt=f'{name}2-relax.txt', 

99 **params) 

100 BFGS(atoms).run(fmax=0.02) 

101 con.write(atoms, name=name, test='relax') 

102 del con[id] 

103 

104 for g in [20, 24, 28]: 

105 atoms.calc = GPAW(mode='fd', 

106 gpts=(g, g, 2 * g), 

107 txt=f'{name}2-fd-{g}.txt', 

108 **params) 

109 try: 

110 atoms.get_potential_energy() 

111 except KohnShamConvergenceError: 

112 continue 

113 con.write(atoms, name=name, test='fd2', gpts=g) 

114 

115 if lcao: 

116 for g in [20, 24, 28]: 

117 id = con.reserve(name=name, test='lcao2', gpts=g) 

118 if id is None: 

119 continue 

120 atoms.calc = GPAW(gpts=(g, g, 2 * g), 

121 mode='lcao', basis='dzp', 

122 txt=f'{name}2-lcao-{g}.txt', 

123 **params) 

124 atoms.get_potential_energy() 

125 con.write(atoms, name=name, test='lcao2', gpts=g) 

126 del con[id] 

127 

128 

129def solve(energies, de): 

130 for i1 in range(len(energies) - 3, -1, -1): 

131 if energies[i1] > de: 

132 break 

133 c1 = cutoffs[i1] 

134 c2 = cutoffs[i1 + 1] 

135 # a * exp(-b * i) 

136 e1 = energies[i1] 

137 e2 = energies[i1 + 1] 

138 b = np.log(e1 / e2) / (c2 - c1) 

139 a = e1 * np.exp(b * c1) 

140 return np.log(a / de) / b 

141 

142 

143def summary(con, name): 

144 eegg = [row.get('eegg', np.nan) 

145 for row in con.select(name=name, test='eggbox', sort='h')] 

146 ecut = np.array([row.energy for row in con.select(name=name, 

147 test='pw1', 

148 sort='ecut')]) 

149 ecut2 = np.array([row.get('energy', np.nan) 

150 for row in con.select(name=name, test='pw2', 

151 sort='ecut')]) 

152 eg = np.array([row.energy for row in con.select(name=name, 

153 test='fd1', sort='gpts')]) 

154 eg2 = np.array([row.get('energy', np.nan) 

155 for row in con.select(name=name, test='fd2', sort='gpts')]) 

156 eL = np.array([row.energy for row in con.select(name=name, 

157 test='lcao1', 

158 sort='gpts')]) 

159 eL2 = np.array([row.get('energy', np.nan) 

160 for row in con.select(name=name, 

161 test='lcao2', 

162 sort='gpts')]) 

163 

164 decut = ecut - 0.5 * ecut2 

165 energies = abs(ecut - ecut[-1]) 

166 denergies = abs(decut - decut[-1]) 

167 assert len(energies) == len(cutoffs) 

168 

169 eg -= ecut[-1] 

170 eg2 -= ecut2[-1] 

171 deg = eg - 0.5 * eg2 

172 

173 eL -= ecut[-1] 

174 eL2 -= ecut2[-1] 

175 deL = eL - 0.5 * eL2 

176 

177 return (energies, denergies, 

178 abs(eg), abs(deg), abs(eL), abs(deL), 

179 eegg) 

180 

181 

182all_names = [ 

183 'H', 'He', 'Li', 'Be', 'B', 'C', 'N', 'O', 'F', 'Ne', 'Na', 'Na.1', 

184 'Mg', 'Mg.2', 'Al', 'Si', 'P', 'S', 'Cl', 'Ar', 'K', 'Ca', 'Sc', 'Ti', 

185 'V', 'V.5', 'Cr', 'Mn', 'Mn.7', 'Fe', 'Co', 'Ni', 'Ni.10', 'Cu', 'Zn', 

186 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', 'Rb', 'Sr', 'Y', 'Zr', 'Nb', 'Nb.5', 

187 'Mo', 'Mo.6', 'Ru', 'Ru.8', 'Rh', 'Rh.9', 'Pd', 'Pd.10', 'Ag', 'Ag.11', 

188 'Cd', 'In', 'Sn', 'Sb', 'Te', 'Te.16', 'I', 'Xe', 'Cs', 'Ba', 'Hf', 

189 'Ta', 'Ta.5', 'W', 'W.6', 'Re', 'Os', 'Os.8', 'Ir', 'Ir.9', 'Pt', 

190 'Pt.10', 'Au', 'Hg', 'Tl', 'Pb', 'Bi', 'Rn'] 

191 

192new_names = [ 

193 'Cr.14', 

194 'La', 'Ce', 'Pr', 'Nd', 'Pm', 'Sm', 'Eu', 'Gd', 

195 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu'] 

196 

197all_names += new_names 

198 

199 

200def main(): 

201 parser = optparse.OptionParser(usage='python -m gpaw.atom.check ' 

202 '[options] symbol ...', 

203 description='Check dataset.') 

204 parser.add_option('-s', '--summary', action='store_true') 

205 parser.add_option('-v', '--verbose', action='store_true') 

206 parser.add_option('-p', '--plot', action='store_true') 

207 parser.add_option('-l', '--lcao', action='store_true') 

208 parser.add_option('-d', '--database', default='check.db') 

209 parser.add_option('--datasets', default='.') 

210 parser.add_option('-e', '--energy-difference', type=float, default=0.01) 

211 opts, names = parser.parse_args() 

212 if not names: 

213 names = [Path.cwd().name] 

214 con = ase.db.connect(opts.database) 

215 if opts.datasets: 

216 from gpaw import setup_paths 

217 setup_paths[:0] = opts.datasets.split(',') 

218 if opts.summary: 

219 for name in names: 

220 try: 

221 E, dE, eegg, ecut, decut, eg, deg, eL, deL = summary( 

222 con, name, opts.energy_difference) 

223 except Exception as ex: 

224 if opts.verbose: 

225 print(name) 

226 traceback.print_exc() 

227 else: 

228 print('{} {}: {}'.format(name, 

229 ex.__class__.__name__, ex)) 

230 else: 

231 print('{0:5} {1:6.1f} {2:6.1f} {2:6.1f}' 

232 .format(name, ecut, decut), 

233 ''.join(f'{e:7.4f}' for e in eegg), 

234 ''.join(f'{e:7.3f}' for e in eg), 

235 ''.join(f'{e:7.3f}' for e in deg), 

236 ''.join(f'{e:7.3f}' for e in eL), 

237 ''.join(f'{e:7.3f}' for e in deL)) 

238 if opts.plot: 

239 for name in names: 

240 E, dE, eegg, ecut, decut, eg, deg, eL, deL = summary( 

241 con, name, opts.energy_difference) 

242 import matplotlib.pyplot as plt 

243 plt.semilogy(cutoffs[:-1], E[:-1]) 

244 plt.semilogy(cutoffs[:-1], dE[:-1]) 

245 plt.semilogy([solve(E, de) for de in eg], eg, 's') 

246 plt.semilogy([solve(dE, de) for de in deg], deg, 'o') 

247 plt.semilogy([solve(E, de) for de in eL], eL, '+') 

248 plt.semilogy([solve(dE, de) for de in deL], deL, 'x') 

249 plt.show() 

250 else: 

251 for name in names: 

252 check(con, name, opts.lcao) 

253 

254 

255if __name__ == '__main__': 

256 main()