Coverage for gpaw/atom/generator2.py: 85%
915 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1from __future__ import annotations
3import sys
4from functools import partial
5from math import exp, log, pi, sqrt
6from typing import Any, Dict, List, Optional, Tuple, Union, TYPE_CHECKING
8import numpy as np
9from ase.data import atomic_numbers, chemical_symbols
10from ase.units import Ha
11from gpaw import __version__ as version
12from gpaw.atom.aeatom import (AllElectronAtom, Channel, GaussianBasis, colors,
13 parse_ld_str)
14from gpaw.basis_data import Basis, BasisFunction
15from gpaw.gaunt import gaunt
16from gpaw.setup_data import SetupData
17from gpaw.typing import Array2D
18from gpaw.utilities import pack_hermitian
19from gpaw.xc.ri.ribasis import generate_ri_basis
20from gpaw.xc.ri.spherical_hse_kernel import RadialHSE
21from scipy.linalg import eigh
22from scipy.optimize import fsolve
23from scipy.special import erf
25if TYPE_CHECKING:
26 from matplotlib import pyplot as plt
29class DatasetGenerationError(Exception):
30 pass
33parameters: dict[str,
34 Union[tuple[str, float | list[float]],
35 tuple[str, float | list[float], dict[str, Any]]]] = {
36 # 1-2:
37 'H1': ('1s,s,p', 0.9),
38 'He2': ('1s,s,p', 1.5),
39 # 3-10:
40 'Li1': ('2s,s,2p', 2.1),
41 'Li3': ('1s,2s,2p,p,d', 1.5),
42 'Be2': ('2s,s,2p', 1.5),
43 'Be4': ('1s,2s,2p,p,d', 1.4),
44 'B3': ('2s,s,2p,p,d', 1.2),
45 'C4': ('2s,s,2p,p,d', 1.2),
46 'N5': ('2s,s,2p,p,d', [1.2, 1.3], {'r0': 1.1}),
47 'O6': ('2s,s,2p,p,d,F', 1.2),
48 'F7': ('2s,s,2p,p,d', [1.2, 1.4]),
49 'Ne8': ('2s,s,2p,p,d', 1.8),
50 # 11-18:
51 'Na1': ('3s,s,3p', 2.6),
52 'Na9': ('2s,3s,2p,3p,d,F', 2.3),
53 'Mg2': ('3s,s,3p,D', 2.6),
54 'Mg10': ('2s,3s,2p,3p,d,F', [2.0, 1.8]),
55 'Al3': ('3s,s,3p,p,d,F', 2.1),
56 'Si4': ('3s,s,3p,p,d,F', 1.9),
57 'P5': ('3s,s,3p,p,d,F', 1.7),
58 'S6': ('3s,s,3p,p,d,F', 1.6),
59 'Cl7': ('3s,s,3p,p,d,F', 1.5),
60 'Ar8': ('3s,s,3p,p,d,F', 1.5),
61 # 19-36:
62 'K1': ('4s,s,4p,D', 3.5),
63 'K9': ('3s,4s,3p,4p,d,d,F', 2.1),
64 'Ca2': ('4s,s,4p', 3.1),
65 'Ca10': ('3s,4s,3p,4p,3d,d,F', 2.1),
66 'Sc3': ('4s,s,4p,p,3d,d', 2.7),
67 'Sc11': ('3s,4s,3p,4p,3d,d,F', 2.3),
68 'Ti4': ('4s,s,4p,p,3d,d', 2.7),
69 'Ti12': ('3s,4s,3p,4p,3d,d,F', [2.2, 2.2, 2.3]),
70 'V5': ('4s,s,4p,p,3d,d', 2.6),
71 'V13': ('3s,4s,3p,4p,3d,d,F', [2.1, 2.1, 2.3]),
72 'Cr6': ('4s,s,4p,p,3d,d', 2.5),
73 'Cr14': ('3s,4s,3p,4p,3d,d,F', [2.1, 2.1, 2.3]),
74 'Mn7': ('4s,s,4p,p,3d,d', 2.4),
75 'Mn15': ('3s,4s,3p,4p,3d,d,F', [2.0, 2.0, 2.2]),
76 'Fe8': ('4s,s,4p,p,3d,d', 2.2),
77 'Fe16': ('3s,4s,3p,4p,3d,d,F', 2.1),
78 'Co9': ('4s,s,4p,p,3d,d', 2.2),
79 'Co17': ('3s,4s,3p,4p,3d,d,F', 2.1),
80 'Ni10': ('4s,s,4p,p,3d,d', 2.1),
81 'Ni18': ('3s,4s,3p,4p,3d,d,F', 2.0),
82 'Cu11': ('4s,s,4p,p,3d,d', 2.1),
83 'Cu19': ('3s,4s,3p,4p,3d,d,F', 1.9),
84 'Zn12': ('4s,s,4p,p,3d', 2.1),
85 'Zn20': ('3s,4s,3p,4p,3d,d,F', 1.9),
86 'Ga3': ('4s,s,4p,p,d,F', 2.2),
87 'Ga13': ('4s,s,4p,p,3d,d,F', 2.2),
88 'Ge4': ('4s,s,4p,p,d,F', 2.1),
89 'Ge14': ('4s,s,4p,p,3d,d,F', 2.1),
90 'As5': ('4s,s,4p,p,d,F', 2.0),
91 'Se6': ('4s,s,4p,p,d,F', 2.1),
92 'Br7': ('4s,s,4p,p,d,F', 2.1),
93 'Kr8': ('4s,s,4p,p,d,F', 2.1),
94 # 37-54:
95 'Rb1': ('5s,s,5p', 3.6),
96 'Rb9': ('4s,5s,4p,5p,d,d,F', 2.5),
97 'Sr2': ('5s,s,5p', 3.3),
98 'Sr10': ('4s,5s,4p,5p,4d,d,F', 2.5),
99 'Y3': ('5s,s,5p,p,4d,d', 3.1),
100 'Y11': ('4s,5s,4p,5p,4d,d,F', 2.5),
101 'Zr4': ('5s,s,5p,p,4d,d', 3.0),
102 'Zr12': ('4s,5s,4p,5p,4d,d,F', 2.5),
103 'Nb5': ('5s,s,5p,p,4d,d', 2.9),
104 'Nb13': ('4s,5s,4p,5p,4d,d,F', [2.4, 2.4, 2.5]),
105 'Mo6': ('5s,s,5p,p,4d,d', 2.8),
106 'Mo14': ('4s,5s,4p,5p,4d,d,F', 2.3),
107 'Tc7': ('5s,s,5p,p,4d,d', 2.7),
108 'Tc15': ('4s,5s,4p,5p,4d,d,F', 2.3),
109 'Ru8': ('5s,s,5p,p,4d,d', 2.6),
110 'Ru16': ('4s,5s,4p,5p,4d,d,F', 2.3),
111 'Rh9': ('5s,s,5p,p,4d,d', 2.5),
112 'Rh17': ('4s,5s,4p,5p,4d,d,F', 2.3),
113 'Pd10': ('5s,s,5p,p,4d,d', 2.4),
114 'Pd18': ('4s,5s,4p,5p,4d,d,F', 2.3),
115 'Ag11': ('5s,s,5p,p,4d,d', 2.4),
116 'Ag19': ('4s,5s,4p,5p,4d,d,F', 2.3),
117 'Cd12': ('5s,s,5p,p,4d,d', 2.4),
118 'Cd20': ('4s,5s,4p,5p,4d,d,F', 2.3),
119 'In13': ('5s,s,5p,p,4d,d,F', 2.6),
120 'Sn14': ('5s,s,5p,p,4d,d,F', 2.5),
121 'Sb15': ('5s,s,5p,p,4d,d,F', 2.5),
122 'Te6': ('5s,6s,5p,p,d,d,F', 2.5),
123 'I7': ('5s,s,5p,p,d,F', 2.4),
124 'Xe8': ('5s,s,5p,p,d,F', 2.3),
125 # 55-56:
126 'Cs1': ('6s,s,6p,5d', [4.3, 4.6, 4.0]),
127 'Cs9': ('5s,6s,5p,6p,5d,0.5d,F', 3.2),
128 'Ba2': ('6s,s,6p,5d', 3.9),
129 'Ba10': ('5s,6s,5p,6p,5d,d,F', 2.2),
130 # 57-71:
131 'La11': ('5s,6s,5p,6p,5d,d,4f,f,G', 2.2),
132 'Ce12': ('5s,6s,5p,6p,5d,d,4f,f,G', 2.2),
133 'Pr13': ('5s,6s,5p,6p,5d,d,4f,f,G', 2.2),
134 'Nd14': ('5s,6s,5p,6p,5d,d,4f,f,G', 2.2),
135 'Pm15': ('5s,6s,5p,6p,5d,d,4f,f,G', 2.2),
136 'Sm16': ('5s,6s,5p,6p,5d,d,4f,f,G', 2.2),
137 'Eu17': ('5s,6s,5p,6p,5d,d,4f,f,G', 2.2),
138 'Gd18': ('5s,6s,5p,6p,5d,d,4f,f,G', 2.2),
139 'Tb19': ('5s,6s,5p,6p,5d,d,4f,f,G', 2.2),
140 'Dy20': ('5s,6s,5p,6p,5d,d,4f,f,G', 2.2),
141 'Ho21': ('5s,6s,5p,6p,5d,d,4f,f,G', 2.2),
142 'Er22': ('5s,6s,5p,6p,5d,d,4f,f,G', 2.2),
143 'Tm23': ('5s,6s,5p,6p,5d,d,4f,f,G', 2.2),
144 'Yb24': ('5s,6s,5p,6p,5d,d,4f,f,G', 2.2),
145 'Lu25': ('5s,6s,5p,6p,5d,d,4f,f,G', 2.2),
146 # 72-86:
147 'Hf4': ('6s,s,6p,p,5d,d', 2.9),
148 'Hf12': ('5s,6s,5p,6p,5d,d,F', 2.4),
149 'Ta5': ('6s,s,6p,p,5d,d', 2.8),
150 'Ta13': ('5s,6s,5p,6p,5d,d,F', 2.4),
151 'W6': ('6s,s,6p,p,5d,d', 2.7),
152 'W14': ('5s,6s,5p,6p,5d,d,F', 2.4),
153 'Re7': ('6s,s,6p,p,5d,d', 2.6),
154 'Re15': ('5s,6s,5p,6p,5d,d,F', 2.4),
155 'Os8': ('6s,s,6p,p,5d,d', 2.6),
156 'Os16': ('5s,6s,5p,6p,5d,d,F', 2.4),
157 'Ir9': ('6s,s,6p,p,5d,d', 2.6),
158 'Ir17': ('5s,6s,5p,6p,5d,d,F', 2.4),
159 'Pt10': ('6s,s,6p,p,5d,d', 2.5),
160 'Pt18': ('5s,6s,5p,6p,5d,d,F', 2.3),
161 'Au11': ('6s,s,6p,p,5d,d', 2.5),
162 'Au19': ('5s,6s,5p,6p,5d,d,F', 2.3),
163 'Hg12': ('6s,s,6p,p,5d,d', 2.5),
164 'Hg20': ('5s,6s,5p,6p,5d,d,F', 2.3),
165 'Tl13': ('6s,s,6p,p,5d,d,F', 2.8),
166 'Pb14': ('6s,s,6p,p,5d,d,F', 2.6),
167 'Bi5': ('6s,s,6p,p,d,F', 2.8),
168 'Bi15': ('6s,s,6p,p,5d,d,F', 2.6),
169 'Po6': ('6s,s,6p,p,d,F', 2.7),
170 'At7': ('6s,s,6p,p,d,F', 2.6),
171 'Rn8': ('6s,s,6p,p,d,F', 2.6),
172 # 87-88:
173 'Fr1': ('6s,s,6p,5d', 4.5),
174 'Fr9': ('6s,7s,6p,7p,6d,d,F', [2.7, 2.5]),
175 'Ra2': ('6s,s,6p,5d', 4.5),
176 'Ra10': ('6s,7s,6p,7p,6d,d,F', [2.7, 2.5]),
177 # 89-102:
178 'Ac11': ('6s,7s,6p,7p,6d,d,5f,f,G', 2.5),
179 'Th12': ('6s,7s,6p,7p,6d,d,5f,f,G', 2.4),
180 'Pa13': ('6s,7s,6p,7p,6d,d,5f,f,G', 2.5),
181 'U14': ('6s,7s,6p,7p,6d,d,5f,f,G', 2.5),
182 'Np15': ('6s,7s,6p,7p,6d,d,5f,f,G', 2.5),
183 'Pu16': ('6s,7s,6p,7p,6d,d,5f,f,G', 2.5),
184 'Am17': ('6s,7s,6p,7p,6d,d,5f,f,G', 2.5),
185 'Cm18': ('6s,7s,6p,7p,6d,d,5f,f,G', 2.5),
186 'Bk19': ('6s,7s,6p,7p,6d,d,5f,f,G', 2.5),
187 'Cf20': ('6s,7s,6p,7p,6d,d,5f,f,G', 2.5),
188 'Es21': ('6s,7s,6p,7p,6d,d,5f,f,G', 2.5),
189 'Fm22': ('6s,7s,6p,7p,6d,d,5f,f,G', 2.5),
190 'Md23': ('6s,7s,6p,7p,6d,d,5f,f,G', 2.5),
191 'No24': ('6s,7s,6p,7p,6d,d,5f,f,G', 2.5)}
194default = [0,
195 1, 2,
196 1, 2, 3, 4, 5, 6, 7, 8,
197 1, 2, 3, 4, 5, 6, 7, 8,
198 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 3, 4, 5, 6, 7, 8,
199 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 6, 7, 8,
200 1, 2, 11,
201 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25,
202 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 6, 7, 8,
203 9, 10,
204 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24]
207semicore = [0,
208 1, 2,
209 3, 4, 3, 4, 5, 6, 7, 8,
210 9, 10, 3, 4, 5, 6, 7, 8,
211 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 13, 14, 5, 6, 7, 8,
212 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 13, 14, 15, 6, 7, 8,
213 9, 10, 11,
214 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25,
215 12, 13, 14, 15, 16, 17, 18, 19, 20, 13, 14, 15, 6, 7, 8]
218def get_number_of_electrons(symbol, name):
219 Z = atomic_numbers[symbol]
220 if name == 'default':
221 return default[Z]
222 return semicore[Z]
225class PAWWaves:
226 def __init__(self, rgd, l, rcut):
227 self.rgd = rgd
228 self.l = l
229 self.rcut = rcut
231 self.n_n = []
232 self.e_n = []
233 self.f_n = []
234 self.phi_ng = []
235 self.phit_ng = None
236 self.pt_ng = None
238 def __len__(self):
239 return len(self.n_n)
241 def add(self, phi_g, n, e, f):
242 self.phi_ng.append(phi_g)
243 self.n_n.append(n)
244 self.e_n.append(e)
245 self.f_n.append(f)
247 def pseudize(self, pseudizer, vtr_g, vr_g, rcmax, ecut):
248 rgd = self.rgd
249 r_g = rgd.r_g
250 phi_ng = self.phi_ng = np.array(self.phi_ng)
251 N = len(phi_ng)
252 phit_ng = self.phit_ng = rgd.empty(N)
253 pt_ng = self.pt_ng = rgd.empty(N)
254 gc = rgd.ceil(self.rcut)
255 gcmax = rgd.ceil(rcmax)
257 l = self.l
259 dgdr_g = 1 / rgd.dr_g
260 d2gdr2_g = rgd.d2gdr2()
262 self.nt_g = rgd.zeros()
263 for n, phi_g in enumerate(phi_ng):
264 phit_ng[n], c0 = pseudizer(phi_g, gc, self.l,
265 divergent=self.n_n[n] <= 0,
266 ecut=ecut)
267 a_g, dadg_g, d2adg2_g = rgd.zeros(3)
268 a_g[1:] = self.phit_ng[n, 1:] / r_g[1:]**l
269 a_g[0] = c0
270 dadg_g[1:-1] = 0.5 * (a_g[2:] - a_g[:-2])
271 d2adg2_g[1:-1] = a_g[2:] - 2 * a_g[1:-1] + a_g[:-2]
272 q_g = (vtr_g - self.e_n[n] * r_g) * self.phit_ng[n]
273 q_g -= 0.5 * r_g**l * (
274 (2 * (l + 1) * dgdr_g + r_g * d2gdr2_g) * dadg_g +
275 r_g * d2adg2_g * dgdr_g**2)
276 q_g[gcmax:] = 0
277 rgd.cut(q_g, self.rcut)
278 q_g[1:] /= r_g[1:]
279 if l == 0:
280 q_g[0] = q_g[1]
281 pt_ng[n] = q_g
283 self.nt_g += self.f_n[n] / 4 / pi * phit_ng[n]**2
285 self.dS_nn = (rgd.integrate(phi_ng[:, None] * phi_ng) -
286 rgd.integrate(phit_ng[:, None] * phit_ng)) / (4 * pi)
287 self.Q = np.dot(self.f_n, self.dS_nn.diagonal())
288 A_nn = rgd.integrate(phit_ng[:, None] * pt_ng) / (4 * pi)
289 self.dH_nn = np.dot(self.dS_nn, np.diag(self.e_n)) - A_nn
290 self.dH_nn *= 0.5
291 self.dH_nn += self.dH_nn.T.copy()
292 pt_ng[:] = np.dot(np.linalg.inv(A_nn.T), pt_ng)
294 def construct_projectors(self, vtr_g, rcmax):
295 pass
297 def calculate_kinetic_energy_correction(self, vr_g, vtr_g):
298 if len(self) == 0:
299 return
300 self.dekin_nn = (self.rgd.integrate(self.phit_ng[:, None] *
301 self.phit_ng *
302 vtr_g, -1) / (4 * pi) -
303 self.rgd.integrate(self.phi_ng[:, None] *
304 self.phi_ng *
305 vr_g, -1) / (4 * pi) +
306 self.dH_nn)
309class PAWSetupGenerator:
310 def __init__(self, aea, projectors, *,
311 core_hole=None,
312 fd=None,
313 yukawa_gamma=0.0,
314 ecut: float = None,
315 omega=None):
316 """fd: stream
317 Text output.
319 yukawa_gamma: separation parameter for RSF"""
321 self.aea = aea
323 self.fd = fd or sys.stdout
324 self.yukawa_gamma = yukawa_gamma
325 self.exxcc_w: Dict[float, float] = {}
326 self.exxcv_wii: Dict[float, Array2D] = {}
327 self.omega = omega
328 self.ecut = ecut
330 self.core_hole: Optional[Tuple[int, int, float]]
331 if core_hole:
332 state, occ = core_hole.split(',')
333 n0 = int(state[0])
334 l0 = 'spdf'.find(state[1])
335 occ0 = float(occ)
336 aea.add(n0, l0, -occ0)
337 self.core_hole = (n0, l0, occ0)
338 else:
339 self.core_hole = None
341 self.l0: Optional[int] = None
342 if projectors[-1].isupper():
343 assert projectors[-2] == ',', projectors
344 self.l0 = 'SPDFG'.find(projectors[-1])
345 projectors = projectors[:-2]
347 self.lmax = -1
348 self.states: Dict[int, List[Union[None, int, float]]] = {}
349 for s in projectors.split(','):
350 l = 'spdf'.find(s[-1])
351 n: Union[None, int, float]
352 if len(s) == 1:
353 n = None
354 elif '.' in s:
355 n = float(s[:-1])
356 else:
357 n = int(s[:-1])
358 if l in self.states:
359 self.states[l].append(n)
360 else:
361 self.states[l] = [n]
362 if l > self.lmax:
363 self.lmax = l
365 # Add empty bound states:
366 for l, nn in self.states.items():
367 for n in nn:
368 if (isinstance(n, int) and
369 (l not in aea.f_lsn or n - l > len(aea.f_lsn[l][0]))):
370 aea.add(n, l, 0)
372 aea.initialize()
373 aea.run()
374 aea.refine()
376 self.rgd = aea.rgd
378 def pseudize(a_g, gc, l=0, points=4, ecut=None, divergent=False):
379 if ecut is None or divergent:
380 return self.rgd.pseudize(a_g, gc, l, points)
381 Gcut = (2 * ecut)**0.5
382 return self.rgd.pseudize_smooth(a_g, gc, l, points,
383 Gcut=Gcut)
385 self.pseudizer = pseudize
387 self.vtr_g = None
389 self.basis = None
391 self.log('\nGenerating PAW', aea.xc.name, 'setup for', aea.symbol)
393 def construct_shape_function(self, alpha=None, rc=None, eps=1e-10):
394 """Build shape-function for compensation charge."""
396 self.alpha = alpha
398 if self.alpha is None:
399 if not isinstance(rc, (float, int)):
400 rc = min(rc)
401 rc = 1.5 * rc
403 def spillage(alpha):
404 """Fraction of gaussian charge outside rc."""
405 x = alpha * rc**2
406 return 1 - erf(sqrt(x)) + 2 * sqrt(x / pi) * exp(-x)
408 def f(alpha):
409 return log(spillage(alpha[0])) - log(eps)
411 self.alpha = fsolve(f, 7.0)[0]
413 self.alpha = round(self.alpha, 1)
414 elif alpha < 0:
415 rc = min(rc)
416 self.log('Shape function: (sin(pi*r/rc)/r)^2, rc={rc:.2f} Bohr'
417 .format(rc=rc))
418 self.ghat_g = np.sinc(self.rgd.r_g / rc)**2 * (pi / 2 / rc**3)
419 self.ghat_g[self.rgd.ceil(rc):] = 0.0
420 self.alpha = -rc
421 return
423 self.log('Shape function: exp(-alpha*r^2), alpha=%.1f Bohr^-2' %
424 self.alpha)
426 self.ghat_g = (np.exp(-self.alpha * self.rgd.r_g**2) *
427 (self.alpha / pi)**1.5)
429 def log(self, *args, **kwargs):
430 print(file=self.fd, *args, **kwargs)
432 def calculate_core_density(self):
433 self.nc_g = self.rgd.zeros()
434 self.tauc_g = self.rgd.zeros()
435 self.ncore = 0
436 self.nvalence = 0
437 self.ekincore = 0.0
438 for l, ch in enumerate(self.aea.channels):
439 for n, f in enumerate(ch.f_n):
440 if (l <= self.lmax and
441 any(n + l + 1 == nn
442 for nn in self.states[l]
443 if isinstance(nn, int))):
444 self.nvalence += f
445 else:
446 self.nc_g += f * ch.calculate_density(n)
447 self.tauc_g += f * ch.calculate_kinetic_energy_density(n)
448 self.ncore += f
449 self.ekincore += f * ch.e_n[n]
451 self.ekincore -= self.rgd.integrate(self.nc_g * self.aea.vr_sg[0], -1)
453 self.log('Core electrons:', self.ncore)
454 self.log('Valence electrons:', self.nvalence)
456 def find_local_potential(self, r0, P, dv0=None):
457 self.r0 = r0
458 self.nderiv0 = P
459 if self.l0 is None:
460 self.find_polynomial_potential(r0, P, dv0)
461 else:
462 self.match_local_potential(r0, P)
464 def find_polynomial_potential(self, r0, P, dv0):
465 self.log('Constructing smooth local potential for r < %.3f' % r0)
466 g0 = self.rgd.ceil(r0)
467 self.vtr_g = self.pseudizer(self.aea.vr_sg[0], g0, 1, P)[0]
468 if dv0:
469 x = self.rgd.r_g[:g0] / r0
470 dv_g = dv0 * (1 - 2 * x**2 + x**4)
471 self.vtr_g[:g0] += x * r0 * dv_g
473 def match_local_potential(self, r0, P):
474 l0 = self.l0
475 self.log('Local potential matching %s-scattering at e=0.0 eV' %
476 'spdfg'[l0] + ' and r=%.2f Bohr' % r0)
478 g0 = self.rgd.ceil(r0)
479 gc = g0 + 20
481 e0 = 0.0
483 ch = Channel(l0)
484 phi_g = self.rgd.zeros()
485 a = ch.integrate_outwards(phi_g, self.rgd, self.aea.vr_sg[0], gc, e0,
486 self.aea.scalar_relativistic, self.aea.Z)[1]
487 phi_g[1:gc] /= self.rgd.r_g[1:gc]
488 phi_g[0] = a
489 phit_g, c = self.pseudizer(phi_g, g0, l=l0, points=P, divergent=True)
491 dgdr_g = 1 / self.rgd.dr_g
492 d2gdr2_g = self.rgd.d2gdr2()
493 a_g = phit_g.copy()
494 a_g[1:] /= self.rgd.r_g[1:]**l0
495 a_g[0] = c
496 dadg_g = self.rgd.zeros()
497 d2adg2_g = self.rgd.zeros()
498 dadg_g[1:-1] = 0.5 * (a_g[2:] - a_g[:-2])
499 d2adg2_g[1:-1] = a_g[2:] - 2 * a_g[1:-1] + a_g[:-2]
500 q_g = (((l0 + 1) * dgdr_g + 0.5 * self.rgd.r_g * d2gdr2_g) * dadg_g +
501 0.5 * self.rgd.r_g * d2adg2_g * dgdr_g**2)
502 q_g[:g0] /= a_g[:g0]
503 q_g += e0 * self.rgd.r_g
504 q_g[0] = 0.0
506 self.vtr_g = self.aea.vr_sg[0].copy()
507 self.vtr_g[0] = 0.0
508 self.vtr_g[1:g0] = q_g[1:g0]
510 def add_waves(self, rc):
511 if isinstance(rc, float):
512 radii = [rc]
513 else:
514 radii = list(rc)
516 if self.lmax >= 0:
517 radii += [radii[-1]] * (self.lmax + 1 - len(radii))
518 del radii[self.lmax + 1:] # remove unused radii
520 self.rcmax = max(radii)
522 self.waves_l = []
523 for l in range(self.lmax + 1):
524 rcut = radii[l]
525 waves = PAWWaves(self.rgd, l, rcut)
526 e = -1.0
527 for n in self.states[l]:
528 if isinstance(n, int):
529 # Bound state:
530 ch = self.aea.channels[l]
531 e = ch.e_n[n - l - 1]
532 f = ch.f_n[n - l - 1]
533 phi_g = ch.phi_ng[n - l - 1]
534 else:
535 if n is None:
536 e += 1.0
537 else:
538 e = n
539 n = -1
540 f = 0.0
541 phi_g = self.rgd.zeros()
542 gc = self.rgd.round(2.5 * self.rcmax)
543 ch = Channel(l)
544 a = ch.integrate_outwards(phi_g, self.rgd,
545 self.aea.vr_sg[0], gc, e,
546 self.aea.scalar_relativistic,
547 self.aea.Z)[1]
548 phi_g[1:gc + 1] /= self.rgd.r_g[1:gc + 1]
549 phi_g[0] = a
550 phi_g /= (self.rgd.integrate(phi_g**2) / (4 * pi))**0.5
552 waves.add(phi_g, n, e, f)
553 self.waves_l.append(waves)
555 def pseudize(self, type='poly', nderiv=6, rcore=None):
556 self.Q = -self.aea.Z + self.ncore
558 self.nt_g = self.rgd.zeros()
559 if type == 'poly':
560 pseudizer = partial(self.pseudizer, points=nderiv)
561 elif type == 'nc':
562 def pseudizer(a_g, gc, l=0, ecut=None, divergent=False):
563 return self.rgd.pseudize_normalized(a_g, gc, l, points=nderiv)
564 for waves in self.waves_l:
565 waves.pseudize(pseudizer, self.vtr_g, self.aea.vr_sg[0],
566 2.0 * self.rcmax, ecut=self.ecut)
567 self.nt_g += waves.nt_g
568 self.Q += waves.Q
570 self.construct_pseudo_core_density(rcore)
571 self.calculate_potentials()
572 self.summarize()
574 def construct_pseudo_core_density(self, rcore):
575 if rcore is None:
576 rcore = self.rcmax * 0.8
577 else:
578 assert abs(rcore) <= self.rcmax
580 if self.ncore == 0:
581 self.nct_g = self.rgd.zeros()
582 self.tauct_g = self.rgd.zeros()
583 elif rcore > 0.0:
584 # Make sure pseudo density is monotonically decreasing:
585 while True:
586 gcore = self.rgd.round(rcore)
587 self.nct_g = self.rgd.pseudize(self.nc_g, gcore)[0]
588 nt_g = self.nt_g + self.nct_g
589 dntdr_g = self.rgd.derivative(nt_g)[:gcore]
590 if dntdr_g.max() < 0.0:
591 break
592 rcore -= 0.01
594 rcore *= 1.2
595 gcore = self.rgd.round(rcore)
596 self.nct_g = self.pseudizer(self.nc_g, gcore)[0]
597 nt_g = self.nt_g + self.nct_g
599 self.log('Constructing smooth pseudo core density for r < %.3f' %
600 rcore)
601 self.nt_g = nt_g
603 self.tauct_g = self.pseudizer(self.tauc_g, gcore)[0]
604 else:
605 rcore *= -1
606 gcore = self.rgd.round(rcore)
607 nt_g = self.pseudizer(self.aea.n_sg[0], gcore)[0]
608 self.nct_g = nt_g - self.nt_g
609 self.nt_g = nt_g
611 self.log('Constructing NLCC-style smooth pseudo core density for '
612 'r < %.3f' % rcore)
614 self.tauct_g = self.pseudizer(self.tauc_g, gcore)[0]
616 self.npseudocore = self.rgd.integrate(self.nct_g)
617 self.log('Pseudo core electrons: %.6f' % self.npseudocore)
618 self.Q -= self.npseudocore
620 def calculate_potentials(self):
621 self.rhot_g = self.nt_g + self.Q * self.ghat_g
622 self.vHtr_g = self.rgd.poisson(self.rhot_g)
624 self.vxct_g = self.rgd.zeros()
625 nt_sg = self.nt_g.reshape((1, -1))
626 self.exct = self.aea.xc.calculate_spherical(
627 self.rgd, nt_sg, self.vxct_g.reshape((1, -1)))
629 self.v0r_g = self.vtr_g - self.vHtr_g - self.vxct_g * self.rgd.r_g
630 self.v0r_g[self.rgd.round(self.rcmax):] = 0.0
632 def summarize(self):
633 self.log('\nProjectors:')
634 self.log(' state occ energy norm rcut')
635 self.log(' nl [Hartree] [eV] [electrons] [Bohr]')
636 self.log('----------------------------------------------------------')
637 for l, waves in enumerate(self.waves_l):
638 for n, e, f, ds in zip(waves.n_n, waves.e_n, waves.f_n,
639 waves.dS_nn.diagonal()):
640 if n == -1:
641 self.log(' %s %10.6f %10.5f %19.2f' %
642 ('spdf'[l], e, e * Ha, waves.rcut))
643 else:
644 self.log(
645 ' %d%s %5.2f %10.6f %10.5f %5.3f %9.2f' %
646 (n, 'spdf'[l], f, e, e * Ha, 1 - ds,
647 waves.rcut))
648 self.log()
650 def construct_projectors(self, rcore):
651 for waves in self.waves_l:
652 waves.construct_projectors(self.vtr_g, 2.45 * self.rcmax)
653 waves.calculate_kinetic_energy_correction(self.aea.vr_sg[0],
654 self.vtr_g)
656 def check_all(self,
657 tol: float = 0.05 # eV
658 ) -> bool:
659 tol /= Ha
660 self.log(('Checking eigenvalues of %s pseudo atom using ' +
661 'a Gaussian basis set:') % self.aea.symbol)
662 self.log(' AE [eV] PS [eV] error [eV]')
664 ok = True
666 for l in range(4):
667 try:
668 e_b = self.check(l)
669 except RuntimeError:
670 self.log('Singular overlap matrix!')
671 ok = False
672 continue
674 n0 = self.number_of_core_states(l)
676 if l < len(self.aea.channels):
677 e0_b = self.aea.channels[l].e_n
678 extra = 6
679 nae = len(self.aea.channels[l].f_n)
680 for n in range(1 + l, nae + 1 + l + extra):
681 if n - 1 - l < nae:
682 f = self.aea.channels[l].f_n[n - 1 - l]
683 self.log('%2d%s %2d' % (n, 'spdf'[l], f), end='')
684 else:
685 self.log(' ', end='')
686 self.log(' %15.3f' % (e0_b[n - 1 - l] * Ha), end='')
687 if n - 1 - l - n0 >= 0:
688 self.log('%15.3f' * 2 %
689 (e_b[n - 1 - l - n0] * Ha,
690 (e_b[n - 1 - l - n0] - e0_b[n - 1 - l]) *
691 Ha))
692 else:
693 self.log()
695 errors = abs(e_b[:nae - n0] - e0_b[n0:nae])
696 if (errors > tol).any():
697 self.log('Error in bound %s-states!' % 'spdf'[l])
698 ok = False
699 errors = abs(e_b[nae - n0:nae - n0 + extra] -
700 e0_b[nae:nae + extra])
701 if (not self.aea.scalar_relativistic and
702 (errors > 10 * tol).any()):
703 self.log('Error in %s-states!' % 'spdf'[l])
704 ok = False
706 return ok
708 def number_of_core_states(self, l):
709 n0 = 0
710 if l < len(self.waves_l):
711 waves = self.waves_l[l]
712 if len(waves) > 0:
713 n0 = waves.n_n[0] - l - 1
714 if n0 < 0 and l < len(self.aea.channels):
715 n0 = (self.aea.channels[l].f_n > 0).sum()
716 elif l < len(self.aea.channels):
717 n0 = (self.aea.channels[l].f_n > 0).sum()
718 return n0
720 def check(self, l):
721 basis = self.aea.channels[0].basis
722 eps = basis.eps
723 alpha_B = basis.alpha_B
725 basis = GaussianBasis(l, alpha_B, self.rgd, eps)
726 H_bb = basis.calculate_potential_matrix(self.vtr_g)
727 H_bb += basis.T_bb
728 S_bb = np.eye(len(basis))
730 if l < len(self.waves_l):
731 waves = self.waves_l[l]
732 if len(waves) > 0:
733 P_bn = self.rgd.integrate(basis.basis_bg[:, None] *
734 waves.pt_ng) / (4 * pi)
735 H_bb += np.dot(np.dot(P_bn, waves.dH_nn), P_bn.T)
736 S_bb += np.dot(np.dot(P_bn, waves.dS_nn), P_bn.T)
738 e_b, _ = eigh(H_bb, S_bb)
739 return e_b
741 def test_convergence(self,
742 ax: 'plt.Axes',
743 show: bool = True) -> None:
744 rgd = self.rgd
745 r_g = rgd.r_g
746 G_k, nt_k = self.rgd.fft(self.nt_g * r_g)
747 rhot_k = self.rgd.fft(self.rhot_g * r_g)[1]
748 ghat_k = self.rgd.fft(self.ghat_g * r_g)[1]
749 vt_k = self.rgd.fft(self.vtr_g)[1]
750 phi_k = self.rgd.fft(self.waves_l[0].phit_ng[0] * r_g)[1]
751 eee_k = 0.5 * nt_k**2 * (4 * pi)**2 / (2 * pi)**3
752 ecc_k = 0.5 * rhot_k**2 * (4 * pi)**2 / (2 * pi)**3
753 egg_k = 0.5 * ghat_k**2 * (4 * pi)**2 / (2 * pi)**3
754 ekin_k = 0.5 * phi_k**2 * G_k**4 / (2 * pi)**3
755 evt_k = nt_k * vt_k * G_k**2 * 4 * pi / (2 * pi)**3
757 eee = 0.5 * rgd.integrate(self.nt_g * rgd.poisson(self.nt_g), -1)
758 ecc = 0.5 * rgd.integrate(self.rhot_g * self.vHtr_g, -1)
759 egg = 0.5 * rgd.integrate(self.ghat_g * rgd.poisson(self.ghat_g), -1)
760 ekin = self.aea.ekin - self.ekincore - self.waves_l[0].dekin_nn[0, 0]
761 evt = rgd.integrate(self.nt_g * self.vtr_g, -1)
763 errors = 10.0**np.arange(-4, 0) / Ha
764 self.log('\nConvergence of energy:')
765 self.log('plane-wave cutoff (wave-length) [ev (Bohr)]\n ', end='')
766 for de in errors:
767 self.log('%14.4f' % (de * Ha), end='')
768 for label, e_k, e0 in [
769 ('e-e', eee_k, eee),
770 ('c-c', ecc_k, ecc),
771 ('g-g', egg_k, egg),
772 ('kin', ekin_k, ekin),
773 ('vt', evt_k, evt)]:
774 self.log('\n%3s: ' % label, end='')
775 e_k = (np.add.accumulate(e_k) - 0.5 * e_k[0] - 0.5 * e_k) * G_k[1]
776 k = len(e_k) - 1
777 for de in errors:
778 while abs(e_k[k] - e_k[-1]) < de:
779 k -= 1
780 G = k * G_k[1]
781 ecut = 0.5 * G**2
782 h = pi / G
783 self.log(f' {ecut * Ha:6.1f} ({h:4.2f})', end='')
784 ax.semilogy(G_k, abs(e_k - e_k[-1]) * Ha, label=label)
785 self.log()
786 ax.axis(xmax=20)
787 ax.set_xlabel('G')
788 ax.set_ylabel('[eV]')
789 ax.legend()
790 if show:
791 plt.show()
793 def plot(
794 self,
795 *,
796 potential_components: 'plt.Axes' | None = None,
797 partial_waves: 'plt.Axes' | None = None,
798 projectors: 'plt.Axes' | None = None,
799 ) -> None:
800 if potential_components is not None:
801 from .plot_dataset import (
802 plot_potential_components,
803 get_plot_pot_comps_params_from_generator as get_pc_args)
804 plot_potential_components(potential_components, *get_pc_args(self))
805 if partial_waves is not None:
806 from .plot_dataset import (
807 plot_partial_waves,
808 get_plot_pwaves_params_from_generator as get_ppw_args)
810 plot_partial_waves(partial_waves, *get_ppw_args(self))
811 if projectors is not None:
812 from .plot_dataset import (
813 plot_projectors,
814 get_plot_projs_params_from_generator as get_pp_args)
816 plot_projectors(projectors, *get_pp_args(self))
818 def create_basis_set(self, tailnorm=0.0005, scale=200.0, splitnorm=0.16,
819 tag=None, ri=None):
820 rgd = self.rgd
821 name = 'dzp' if not tag else f'{tag}.dzp'
823 # We print text to stdout and put it in the basis-set file
824 txt = 'Basis functions:\n'
826 # Bound states:
827 bf_j = []
828 for l, waves in enumerate(self.waves_l):
829 for i, n in enumerate(waves.n_n):
830 if n > 0:
831 tn = tailnorm
832 if waves.f_n[i] == 0:
833 tn = min(0.05, tn * 20) # no need for long tail
834 if l == 3:
835 tn = 0.01
836 phit_g, ronset, rc, de = self.create_basis_function(
837 l, i, tn, scale)
838 bf = BasisFunction(n, l, rc, phit_g, 'bound state')
839 bf_j.append(bf)
841 txt += '%d%s bound state:\n' % (n, 'spdf'[l])
842 txt += (' cutoff: %.3f to %.3f Bohr (tail-norm=%f)\n' %
843 (ronset, rc, tn))
844 txt += ' eigenvalue shift: %.3f eV\n' % (de * Ha)
846 # Split valence:
847 for l, waves in enumerate(self.waves_l):
848 # Find the largest n that is occupied:
849 n0 = None
850 for f, n in zip(waves.f_n, waves.n_n):
851 if n > 0 and f > 0:
852 n0 = n
853 if n0 is None:
854 continue
856 for bf in bf_j:
857 if bf.l == l and bf.n == n0:
858 break
860 # Radius and l-value used for polarization function below:
861 rcpol = bf.rc
862 lpol = l + 1
864 phit_g = bf.phit_g
866 # Find cutoff radius:
867 n_g = np.add.accumulate(phit_g**2 * rgd.r_g**2 * rgd.dr_g)
868 norm = n_g[-1]
869 gc = (norm - n_g > splitnorm * norm).sum()
870 rc = rgd.r_g[gc]
872 phit2_g = self.pseudizer(phit_g, gc, l, 2)[0] # "split valence"
873 bf = BasisFunction(n, l, rc, phit_g - phit2_g, 'split valence')
874 bf_j.append(bf)
876 txt += '%d%s split valence:\n' % (n0, 'spdf'[l])
877 txt += f' cutoff: {rc:.3f} Bohr (tail-norm={splitnorm:f})\n'
879 # Polarization:
880 if lpol < 4:
881 gcpol = rgd.round(rcpol)
882 alpha = 1 / (0.25 * rcpol)**2
884 # Gaussian that is continuous and has a continuous
885 # derivative at rcpol:
886 phit_g = np.exp(-alpha * rgd.r_g**2) * rgd.r_g**lpol
887 phit_g -= self.pseudizer(phit_g, gcpol, lpol, 2)[0]
888 phit_g[gcpol:] = 0.0
890 bf = BasisFunction(None, lpol, rcpol, phit_g, 'polarization')
891 bf_j.append(bf)
892 txt += f'l={lpol} polarization functions:\n'
893 txt += f' cutoff: {rcpol:.3f} Bohr '
894 txt += f'(r^{lpol} exp(-{alpha:.3f}*r^2))\n'
896 self.log(txt)
898 basis = Basis(self.aea.symbol, name, rgd=rgd, bf_j=bf_j,
899 generatordata=txt,
900 generatorattrs=dict(tailnorm=tailnorm,
901 scale=scale,
902 splitnorm=splitnorm))
904 if ri:
905 basis = generate_ri_basis(basis, ri)
907 self.basis = basis
908 return basis
910 def create_basis_function(self, l, n, tailnorm, scale):
911 rgd = self.rgd
912 waves = self.waves_l[l]
914 # Find cutoff radii:
915 n_g = np.add.accumulate(waves.phit_ng[n]**2 * rgd.r_g**2 * rgd.dr_g)
916 norm = n_g[-1]
917 g2 = (norm - n_g > tailnorm * norm).sum()
918 r2 = rgd.r_g[g2]
919 r1 = max(0.6 * r2, waves.rcut)
920 g1 = rgd.ceil(r1)
921 # Set up confining potential:
922 r = rgd.r_g[g1:g2]
923 vtr_g = self.vtr_g.copy()
924 vtr_g[g1:g2] += scale * np.exp((r2 - r1) / (r1 - r)) / (r - r2)**2
925 vtr_g[g2:] = np.inf
927 # Nonlocal PAW stuff:
928 pt_ng = waves.pt_ng
929 dH_nn = waves.dH_nn
930 dS_nn = waves.dS_nn
931 N = len(pt_ng)
933 u_g = rgd.zeros()
934 u_ng = rgd.zeros(N)
935 duodr_n = np.empty(N)
936 a_n = np.empty(N)
938 e = waves.e_n[n]
939 e0 = e
940 ch = Channel(l)
941 while True:
942 duodr, a = ch.integrate_outwards(u_g, rgd, vtr_g, g1, e,
943 scalar_relativistic=False)
945 for n in range(N):
946 duodr_n[n], a_n[n] = ch.integrate_outwards(
947 u_ng[n], rgd, vtr_g, g1, e,
948 scalar_relativistic=False, pt_g=pt_ng[n])
950 A_nn = (dH_nn - e * dS_nn) / (4 * pi)
951 B_nn = rgd.integrate(pt_ng[:, None] * u_ng, -1)
952 c_n = rgd.integrate(pt_ng * u_g, -1)
953 d_n = np.linalg.solve(np.dot(A_nn, B_nn) + np.eye(N),
954 np.dot(A_nn, c_n))
955 u_g[:g1 + 1] -= np.dot(d_n, u_ng[:, :g1 + 1])
956 a -= np.dot(d_n, a_n)
957 duodr -= np.dot(duodr_n, d_n)
958 uo = u_g[g1]
960 duidr = ch.integrate_inwards(u_g, rgd, vtr_g, g1, e, gmax=g2,
961 scalar_relativistic=False)
962 ui = u_g[g1]
963 A = duodr / uo - duidr / ui
964 u_g[g1:] *= uo / ui
965 x = (norm / rgd.integrate(u_g**2, -2) * (4 * pi))**0.5
966 u_g *= x
967 a *= x
969 if abs(A) < 1e-5:
970 break
972 e += 0.5 * A * u_g[g1]**2
974 u_g[1:] /= rgd.r_g[1:]
975 u_g[0] = a * 0.0**l
976 return u_g, r1, r2, e - e0
978 def logarithmic_derivative(self, l, energies, rcut):
979 rgd = self.rgd
980 ch = Channel(l)
981 gcut = rgd.round(rcut)
983 N = 0
984 if l < len(self.waves_l):
985 # Nonlocal PAW stuff:
986 waves = self.waves_l[l]
987 if len(waves) > 0:
988 pt_ng = waves.pt_ng
989 dH_nn = waves.dH_nn
990 dS_nn = waves.dS_nn
991 N = len(pt_ng)
993 u_g = rgd.zeros()
994 u_ng = rgd.zeros(N)
995 dudr_n = np.empty(N)
997 logderivs = []
998 d0 = 42.0
999 offset = 0
1000 for e in energies:
1001 dudr = ch.integrate_outwards(u_g, rgd, self.vtr_g, gcut, e,
1002 scalar_relativistic=False)[0]
1003 u = u_g[gcut]
1005 if N:
1006 for n in range(N):
1007 dudr_n[n] = ch.integrate_outwards(
1008 u_ng[n], rgd, self.vtr_g, gcut, e,
1009 scalar_relativistic=False, pt_g=pt_ng[n])[0]
1011 A_nn = (dH_nn - e * dS_nn) / (4 * pi)
1012 B_nn = rgd.integrate(pt_ng[:, None] * u_ng, -1)
1013 c_n = rgd.integrate(pt_ng * u_g, -1)
1014 d_n = np.linalg.solve(np.dot(A_nn, B_nn) + np.eye(N),
1015 np.dot(A_nn, c_n))
1016 u -= np.dot(u_ng[:, gcut], d_n)
1017 dudr -= np.dot(dudr_n, d_n)
1019 d1 = np.arctan(dudr / u) / pi + offset
1020 if d1 > d0:
1021 offset -= 1
1022 d1 -= 1
1023 logderivs.append(d1)
1024 d0 = d1
1026 return np.array(logderivs)
1028 def make_paw_setup(self, tag=None):
1029 aea = self.aea
1030 setup = SetupData(aea.symbol, aea.xc.name, tag, readxml=False,
1031 generator_version=3)
1033 setup.id_j = []
1035 J = [] # new reordered j and i indices
1036 I = []
1038 # Bound states:
1039 j = 0
1040 i = 0
1041 for l, waves in enumerate(self.waves_l):
1042 for n, f, e, phi_g, phit_g, pt_g in zip(waves.n_n, waves.f_n,
1043 waves.e_n, waves.phi_ng,
1044 waves.phit_ng,
1045 waves.pt_ng):
1046 if n != -1:
1047 setup.append(n, l, f, e, waves.rcut, phi_g, phit_g, pt_g)
1048 id = '%d%s' % (n, 'spdf'[l])
1049 setup.id_j.append(id)
1050 J.append(j)
1051 I.extend(range(i, i + 2 * l + 1))
1052 j += 1
1053 i += 2 * l + 1
1055 # Excited states:
1056 j = 0
1057 i = 0
1058 for l, waves in enumerate(self.waves_l):
1059 ne = 0
1060 for n, f, e, phi_g, phit_g, pt_g in zip(waves.n_n, waves.f_n,
1061 waves.e_n, waves.phi_ng,
1062 waves.phit_ng,
1063 waves.pt_ng):
1064 if n == -1:
1065 setup.append(n, l, f, e, waves.rcut, phi_g, phit_g, pt_g)
1066 ne += 1
1067 id = '%s%d' % ('spdf'[l], ne)
1068 setup.id_j.append(id)
1069 J.append(j)
1070 I.extend(range(i, i + 2 * l + 1))
1071 j += 1
1072 i += 2 * l + 1
1074 nj = sum(len(waves) for waves in self.waves_l)
1075 e_kin_jj = np.zeros((nj, nj))
1076 j1 = 0
1077 for waves in self.waves_l:
1078 j2 = j1 + len(waves)
1079 e_kin_jj[j1:j2, j1:j2] = waves.dekin_nn
1080 j1 = j2
1081 setup.e_kin_jj = e_kin_jj[J][:, J].copy()
1083 setup.nc_g = self.nc_g * sqrt(4 * pi)
1084 setup.nct_g = self.nct_g * sqrt(4 * pi)
1085 setup.e_kinetic_core = self.ekincore
1086 setup.vbar_g = self.v0r_g * sqrt(4 * pi)
1087 setup.vbar_g[1:] /= self.rgd.r_g[1:]
1088 setup.vbar_g[0] = setup.vbar_g[1]
1089 setup.vt_g = self.vtr_g * sqrt(4 * pi)
1090 setup.vt_g[1:] /= self.rgd.r_g[1:]
1091 setup.vt_g[0] = setup.vt_g[1]
1092 setup.Z = aea.Z
1093 setup.Nc = self.ncore
1094 setup.Nv = self.nvalence
1095 setup.e_kinetic = aea.ekin
1096 setup.e_xc = aea.exc
1097 setup.e_electrostatic = aea.eH + aea.eZ
1098 setup.e_total = aea.exc + aea.ekin + aea.eH + aea.eZ
1099 setup.rgd = self.rgd
1100 if self.alpha > 0:
1101 setup.shape_function = {'type': 'gauss',
1102 'rc': 1 / sqrt(self.alpha)}
1103 else:
1104 setup.shape_function = {'type': 'sinc',
1105 'rc': -self.alpha}
1107 self.calculate_exx_integrals()
1108 setup.ExxC = self.exxcc
1109 setup.ExxC_w = self.exxcc_w # erfc screened core contributions
1110 setup.X_p = pack_hermitian(self.exxcv_ii[I][:, I])
1111 setup.X_wp = {omega: pack_hermitian(self.exxcv_wii[omega][I][:, I])
1112 for omega in self.exxcv_wii}
1114 if self.yukawa_gamma > 0.0:
1115 self.calculate_yukawa_integrals()
1116 setup.X_pg = pack_hermitian(self.exxgcv_ii[I][:, I])
1118 setup.tauc_g = self.tauc_g * (4 * pi)**0.5
1119 setup.tauct_g = self.tauct_g * (4 * pi)**0.5
1121 if self.aea.scalar_relativistic:
1122 reltype = 'scalar-relativistic'
1123 else:
1124 reltype = 'non-relativistic'
1125 attrs = [('type', reltype),
1126 ('version', 3),
1127 ('name', 'gpaw-%s' % version)]
1128 setup.generatorattrs = attrs
1129 setup.version = '0.7'
1131 setup.l0 = self.l0
1132 setup.e0 = 0.0
1133 setup.r0 = self.r0
1134 setup.nderiv0 = self.nderiv0
1136 setup.basis = self.basis
1138 if self.core_hole:
1139 n, l, occ = self.core_hole
1140 phi_g = self.aea.channels[l].phi_ng[n - l - 1]
1141 setup.ncorehole = n
1142 setup.lcorehole = l
1143 setup.fcorehole = occ
1144 setup.phicorehole_g = phi_g
1145 setup.has_corehole = True
1146 ch = self.aea.channels[l]
1147 eig = ch.e_n[n - 1 - l]
1148 setup.core_hole_e = eig
1149 setup.core_hole_e_kin = eig - self.rgd.integrate(
1150 ch.calculate_density(n - 1 - l) * self.aea.vr_sg[0], -1)
1151 return setup
1153 def find_core_states(self):
1154 # Find core states:
1155 core = []
1156 lmax = 0
1157 for l, ch in enumerate(self.aea.channels):
1158 for n, phi_g in enumerate(ch.phi_ng):
1159 if (l >= len(self.waves_l) or
1160 (l < len(self.waves_l) and
1161 n + l + 1 not in self.waves_l[l].n_n)):
1162 core.append((l, phi_g))
1163 if l > lmax:
1164 lmax = l
1166 lmax = max(lmax, len(self.waves_l) - 1)
1167 G_LLL = gaunt(lmax)
1169 return lmax, core, G_LLL
1171 def core_core_exchange(self, interaction):
1172 # Calculate core contribution to EXX energy per interaction kernel
1173 (lmax, core, G_LLL) = self.find_core_states()
1174 E = 0.0
1175 j1 = 0
1176 for l1, phi1_g in core:
1177 f = 1.0
1178 for l2, phi2_g in core[j1:]:
1179 n_g = phi1_g * phi2_g
1180 for l in range((l1 + l2) % 2, l1 + l2 + 1, 2):
1181 G = (G_LLL[l1**2:(l1 + 1)**2,
1182 l2**2:(l2 + 1)**2,
1183 l**2:(l + 1)**2]**2).sum()
1184 vr_g = interaction(n_g, l)
1185 e = f * self.rgd.integrate(vr_g * n_g, -1) / 4 / pi
1186 E -= e * G
1187 f = 2.0
1188 j1 += 1
1189 return E
1191 def calculate_exx_integrals(self):
1193 # Calculate core-valence contribution to EXX energy:
1194 ni = sum(len(waves) * (2 * l + 1)
1195 for l, waves in enumerate(self.waves_l))
1197 self.exxcc = self.core_core_exchange(self.rgd.poisson)
1198 self.log('EXX (core-core):', self.exxcc, 'Hartree')
1200 if self.omega is not None:
1201 hse = RadialHSE(self.rgd, self.omega)
1203 self.exxcc_w = {self.omega:
1204 self.core_core_exchange(hse.screened_coulomb)}
1205 self.log(f'EXX omega={self.omega} (core-core):',
1206 self.exxcc_w[self.omega], 'Hartree')
1208 X_ii = self.calculate_exx_cv_integrals(ni, hse.screened_coulomb)
1209 self.exxcv_wii = {self.omega: X_ii}
1211 self.exxcv_ii = self.calculate_exx_cv_integrals(ni, self.rgd.poisson)
1213 def calculate_yukawa_integrals(self):
1214 """Wrapper to calculate the rsf core-valence contribution."""
1215 ni = sum(len(waves) * (2 * l + 1)
1216 for l, waves in enumerate(self.waves_l))
1218 def interaction(n_g, l):
1219 return self.rgd.yukawa(n_g, l, self.yukawa_gamma)
1221 self.exxgcv_ii = self.calculate_exx_cv_integrals(ni, interaction)
1223 def calculate_exx_cv_integrals(self, ni, interaction):
1224 (lmax, core, G_LLL) = self.find_core_states()
1225 cv_ii = np.zeros((ni, ni))
1227 i1 = 0
1228 for l1, waves1 in enumerate(self.waves_l):
1229 for phi1_g in waves1.phi_ng:
1230 i2 = 0
1231 for l2, waves2 in enumerate(self.waves_l):
1232 for phi2_g in waves2.phi_ng:
1233 X_mm = cv_ii[i1:i1 + 2 * l1 + 1,
1234 i2:i2 + 2 * l2 + 1]
1235 if (l1 + l2) % 2 == 0:
1236 for lc, phi_g in core:
1237 n_g = phi1_g * phi_g
1238 for l in range((l1 + lc) % 2,
1239 max(l1, l2) + lc + 1, 2):
1240 n2c = phi2_g * phi_g
1241 vr_g = interaction(n2c, l)
1242 e = (self.rgd.integrate(vr_g * n_g, -1) /
1243 (4 * pi))
1244 for mc in range(2 * lc + 1):
1245 for m in range(2 * l + 1):
1246 G_L = G_LLL[:,
1247 lc**2 + mc,
1248 l**2 + m]
1249 X_mm += np.outer(
1250 G_L[l1**2:(l1 + 1)**2],
1251 G_L[l2**2:(l2 + 1)**2]) * e
1252 i2 += 2 * l2 + 1
1253 i1 += 2 * l1 + 1
1254 return cv_ii
1257def get_parameters(symbol, args):
1258 Z = atomic_numbers.get(symbol)
1259 if Z is None:
1260 Z = float(symbol)
1261 symbol = chemical_symbols[int(round(Z))]
1263 if args.electrons:
1264 par = parameters[symbol + str(args.electrons)]
1265 else:
1266 par = parameters[symbol + str(default[int(round(Z))])]
1268 projectors, radii = par[:2]
1269 if len(par) == 3:
1270 extra = par[2]
1271 else:
1272 extra = {}
1274 if args.configuration:
1275 configuration = args.configuration
1276 else:
1277 configuration = None
1279 if args.projectors:
1280 projectors = args.projectors
1282 if args.radius:
1283 radii = [float(r) for r in args.radius.split(',')]
1285 if isinstance(radii, float):
1286 radii = [radii]
1288 if args.pseudize:
1289 type, nderiv = args.pseudize.split(',')
1290 pseudize = (type, int(nderiv))
1291 else:
1292 pseudize = ('poly', 4)
1294 if args.zero_potential:
1295 x = args.zero_potential.split(',')
1296 nderiv0 = int(x[0])
1297 r0 = float(x[1])
1298 else:
1299 if projectors[-1].isupper():
1300 nderiv0 = 5
1301 r0 = extra.get('r0', min(radii) * 0.9)
1302 else:
1303 nderiv0 = 2
1304 r0 = extra.get('r0', min(radii))
1306 if args.pseudo_core_density_radius:
1307 rcore = args.pseudo_core_density_radius
1308 else:
1309 rcore = extra.get('rcore')
1311 if args.nlcc:
1312 rcore *= -1
1314 return dict(symbol=symbol,
1315 Z=Z,
1316 xc=args.xc_functional,
1317 configuration=configuration,
1318 projectors=projectors,
1319 radii=radii,
1320 scalar_relativistic=args.scalar_relativistic, alpha=args.alpha,
1321 r0=r0, v0=None, nderiv0=nderiv0,
1322 pseudize=pseudize, rcore=rcore,
1323 core_hole=args.core_hole,
1324 output=None, # args.output,
1325 yukawa_gamma=args.gamma,
1326 ecut=args.ecut,
1327 omega=args.omega)
1330def generate(symbol,
1331 Z,
1332 projectors,
1333 radii,
1334 r0, v0,
1335 nderiv0,
1336 xc='LDA',
1337 scalar_relativistic=True,
1338 pseudize=('poly', 4),
1339 configuration=None,
1340 alpha=None,
1341 rcore=None,
1342 core_hole=None,
1343 *,
1344 output=None,
1345 yukawa_gamma=0.0,
1346 ecut=None,
1347 omega=None):
1348 aea = AllElectronAtom(symbol, xc, Z=Z,
1349 configuration=configuration,
1350 scalar_relativistic=scalar_relativistic)
1351 gen = PAWSetupGenerator(aea, projectors,
1352 core_hole=core_hole,
1353 fd=output,
1354 yukawa_gamma=yukawa_gamma,
1355 ecut=ecut,
1356 omega=omega)
1358 gen.construct_shape_function(alpha, radii, eps=1e-10)
1359 gen.calculate_core_density()
1360 gen.find_local_potential(r0, nderiv0, v0)
1361 gen.add_waves(radii)
1362 gen.pseudize(pseudize[0], pseudize[1], rcore=rcore)
1363 gen.construct_projectors(rcore)
1364 return gen
1367def generate_all():
1368 if len(sys.argv) > 1:
1369 atoms = sys.argv[1:]
1370 else:
1371 atoms = None
1372 functionals = ['LDA', 'PBE', 'revPBE', 'RPBE']
1373 for name in parameters:
1374 n = sum(c.isalpha() for c in name)
1375 symbol = name[:n]
1376 electrons = name[n:]
1377 if atoms and symbol not in atoms:
1378 continue
1379 print(name, symbol, electrons)
1380 for xc in functionals:
1381 argv = [symbol, '-swf', xc, '-e', electrons, '-t', electrons + 'e']
1382 if xc == 'PBE':
1383 argv.append('-b')
1384 main(argv)
1387class CLICommand:
1388 """Create PAW dataset."""
1390 @staticmethod
1391 def add_arguments(parser):
1392 add = parser.add_argument
1393 add('symbol')
1394 add('-f', '--xc-functional', type=str, default='LDA',
1395 help='Exchange-Correlation functional (default value LDA)',
1396 metavar='<XC>')
1397 add('-C', '--configuration',
1398 help='Example for Nd: "[Xe]6s1,5d1,4f4".')
1399 add('-P', '--projectors',
1400 help='Projector functions - use comma-separated - ' +
1401 'nl values, where n can be principal quantum number ' +
1402 '(integer) or energy (floating point number). ' +
1403 'Example: 2s,0.5s,2p,0.5p,0.0d.')
1404 add('-r', '--radius',
1405 help='Cutoff radius for all projectors or comma-separated list '
1406 'of radii for each angular momentum channel. '
1407 'Example: 1.2 or 1.2,1.1,1.1. Units: Bohr.')
1408 add('-0', '--zero-potential',
1409 metavar='nderivs,radius',
1410 help='Parameters for zero potential.')
1411 add('-c', '--pseudo-core-density-radius', type=float,
1412 metavar='radius',
1413 help='Radius for pseudizing core density.')
1414 add('-z', '--pseudize',
1415 metavar='type,nderivs',
1416 help='Parameters for pseudizing wave functions.')
1417 add('-p', '--plot',
1418 const=True,
1419 default=False,
1420 metavar='FILE',
1421 nargs='?',
1422 help='Show plots of the setup; '
1423 'if a filename is supplied, write the plots thereto instead of '
1424 '`plt.show()`-ing them')
1425 add('-S', '--separate-figures',
1426 action='store_true',
1427 help='If not plotting to a file, '
1428 'plot the plots in separate figure windows/tabs, '
1429 'instead of as subplots/panels in the same figure')
1430 add('-l', '--logarithmic-derivatives',
1431 metavar='spdfg,e1:e2:de,radius',
1432 help='Plot logarithmic derivatives. ' +
1433 'Example: -l spdf,-1:1:0.05,1.3. ' +
1434 'Energy range and/or radius can be left out.')
1435 add('-w', '--write', action='store_true',
1436 help='Write setup to file <symbol>.<XC> '
1437 'or, with --tag, <symbol>.<TAG>.<XC>.')
1438 add('-s', '--scalar-relativistic', action='store_true',
1439 help='Perform a scalar-relativistic calculation. '
1440 'Default is a non-scalar-relativistic.')
1441 add('-n', '--no-check', action='store_true',
1442 help='Disable error checks. This allows saving files that would '
1443 'normally be considered invalid.')
1444 add('-t', '--tag', type=str, help='Include TAG in output filename.')
1445 add('-a', '--alpha', type=float,
1446 help='Decay of shape function exp(-alpha r²) used for '
1447 'compensation charges. Default is based on --radius.')
1448 add('-g', '--gamma', type=float, default=0.0,
1449 help='Yukawa gamma parameter for range-separated functionals. '
1450 'Default is zero (no range separation).')
1451 add('-b', '--create-basis-set', action='store_true',
1452 help='Create a rudimentary basis set.')
1453 add('--nlcc', action='store_true',
1454 help='Use NLCC-style pseudo core density '
1455 '(for vdW-DF functionals).')
1456 add('--core-hole', metavar='STATE,OCC',
1457 help='State and occupation of core hole. Default is no '
1458 'core hole. Example: "3s,0.5."')
1459 add('-e', '--electrons', type=int,
1460 help='Number of valence electrons. Requires a corresponding '
1461 'entry in the global parameters dictionary (read source code).')
1462 # --output does not currently work since it never closes or flushes.
1463 # add('-o', '--output', metavar='FILE',
1464 # help='Write log to FILE.')
1465 add('--ecut', type=float, help='Minimize Fourier components above '
1466 'ECUT for pseudo wave functions.')
1467 add('--ri', type=str,
1468 help='Calculate also resolution of identity basis.')
1469 add('--omega', type=float, default=None,
1470 help='Calculate core-core and core-valence contributions'
1471 ' to erfc screen HSE-potential')
1473 @staticmethod
1474 def run(args):
1475 main(args)
1478def main(args):
1479 kwargs = get_parameters(args.symbol, args)
1480 gen = generate(**kwargs)
1482 should_plot_dataset = args.logarithmic_derivatives or args.plot
1484 if not args.no_check:
1485 if not gen.check_all():
1486 raise DatasetGenerationError
1488 if args.create_basis_set:
1489 basis = gen.create_basis_set(tag=args.tag, ri=args.ri)
1490 basis.write_xml() # XXX: should this only happen if `.write`?
1491 else:
1492 basis = None
1494 if args.write or should_plot_dataset:
1495 setup = gen.make_paw_setup(args.tag)
1496 parameters = []
1497 for key, value in kwargs.items():
1498 if value is not None:
1499 parameters.append(f'{key}={value!r}')
1500 setup.generatordata = ',\n '.join(parameters)
1501 if args.write:
1502 setup.write_xml()
1503 else:
1504 setup = None
1506 if not args.create_basis_set and args.ri:
1507 raise ValueError('Basis set must be created in order to create the '
1508 'RI-basis set as well')
1510 if should_plot_dataset:
1511 from matplotlib import pyplot as plt
1512 from .plot_dataset import plot_dataset
1514 assert setup is not None
1515 ax_objs, plot_fname = plot_dataset(
1516 setup,
1517 basis=basis,
1518 gen=gen,
1519 plot_potential_components=bool(args.plot),
1520 plot_partial_waves=bool(args.plot),
1521 plot_projectors=bool(args.plot),
1522 plot_logarithmic_derivatives=args.logarithmic_derivatives,
1523 separate_figures=args.separate_figures,
1524 savefig=(None if args.plot in (True, False) else args.plot),
1525 )
1526 if ax_objs and plot_fname is None:
1527 plt.show()
1530def plot_log_derivs(gen: PAWSetupGenerator,
1531 ld_str: str,
1532 plot: bool,
1533 ax: 'plt.Axes') -> None:
1534 """Make nice log-derivs plot."""
1536 r = 1.1 * gen.rcmax
1537 emin = min(min(wave.e_n) for wave in gen.waves_l) - 0.8
1538 emax = max(max(wave.e_n) for wave in gen.waves_l) + 0.8
1539 lvalues, energies, r = parse_ld_str(ld_str, (emin, emax, 0.05), r)
1540 emin = energies[0]
1541 de = energies[1] - emin
1543 error = 0.0
1544 aea = gen.aea
1545 for l in lvalues:
1546 efix = []
1547 # Fixed points:
1548 if l < len(gen.waves_l):
1549 efix.extend(gen.waves_l[l].e_n)
1550 if l == gen.l0:
1551 efix.append(0.0)
1553 ld1 = aea.logarithmic_derivative(l, energies, r)
1554 ld2 = gen.logarithmic_derivative(l, energies, r)
1555 for e in efix:
1556 i = int((e - emin) / de)
1557 if 0 <= i < len(energies):
1558 ld1 -= round(ld1[i] - ld2[i])
1559 if plot:
1560 ldfix = ld1[i]
1561 ax.plot([energies[i]], [ldfix], 'x' + colors[l])
1563 if plot:
1564 ax.plot(energies, ld1, colors[l], label='spdfg'[l])
1565 ax.plot(energies, ld2, '--' + colors[l])
1567 error = abs(ld1 - ld2).sum() * de
1568 print('Logarithmic derivative error:', l, error)
1570 if plot:
1571 ax.set_title(f'Logarithmic derivatives: {aea.symbol} {aea.xc.name}')
1572 ax.set_xlabel('energy [Ha]')
1573 ax.set_ylabel(r'$\arctan(d\log\phi_{\ell\epsilon}(r)/dr)/\pi'
1574 r'|_{r=r_c}$')
1575 ax.legend(loc='best')