Coverage for gpaw/setup.py: 97%
893 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-08 00:17 +0000
1from __future__ import annotations
2import functools
3from io import StringIO
4from math import pi, sqrt
5import ase.units as units
6import numpy as np
7from ase.data import chemical_symbols
9from gpaw import debug
10from gpaw.basis_data import Basis, BasisFunction
11from gpaw.gaunt import gaunt, nabla
12from gpaw.overlap import OverlapCorrections
13from gpaw.setup_data import SetupData, search_for_file
14from gpaw.spline import Spline
15from gpaw.utilities import pack_density, unpack_hermitian
16from gpaw.xc import XC
17from gpaw.new import zips
18from gpaw.xc.ri.spherical_hse_kernel import RadialHSE
19from gpaw.core.atom_arrays import AtomArraysLayout
22class WrongMagmomForHundsRuleError(ValueError):
23 """
24 Custom error for catching bad magnetic moments in Hund's rule calculation
25 """
28def create_setup(symbol, xc='LDA', lmax=0,
29 type='paw', basis=None, setupdata=None,
30 filter=None, world=None, backwards_compatible=True):
31 if isinstance(xc, str):
32 xc = XC(xc)
34 if isinstance(type, str) and ':' in type:
35 from gpaw.hubbard import parse_hubbard_string
36 type, hubbard_u = parse_hubbard_string(type)
37 else:
38 hubbard_u = None
40 if setupdata is None:
41 if type == 'hgh' or type == 'hgh.sc':
42 lmax = 0
43 from gpaw.hgh import HGHSetupData, sc_setups, setups
44 if type == 'hgh.sc':
45 table = sc_setups
46 else:
47 table = setups
48 parameters = table[symbol]
49 setupdata = HGHSetupData(parameters)
50 elif type == 'ah':
51 from gpaw.ah import AppelbaumHamann
52 ah = AppelbaumHamann()
53 ah.build(basis)
54 return ah
55 elif type == 'ae':
56 from gpaw.ae import HydrogenAllElectronSetup
57 assert symbol == 'H'
58 ae = HydrogenAllElectronSetup()
59 ae.build(basis)
60 return ae
61 elif type == 'ghost':
62 from gpaw.lcao.bsse import GhostSetupData
63 setupdata = GhostSetupData(symbol)
64 elif type == 'sg15':
65 from gpaw.upf import read_sg15
66 upfname = f'{symbol}_ONCV_PBE-*.upf'
67 try:
68 upfpath, source = search_for_file(upfname, world=world)
69 except RuntimeError:
70 raise OSError('Could not find pseudopotential file %s '
71 'in any GPAW search path. '
72 'Please install the SG15 setups using, '
73 'e.g., "gpaw install-data".' % upfname)
74 setupdata = read_sg15(upfpath)
75 if xc.get_setup_name() != 'PBE':
76 raise ValueError('SG15 pseudopotentials support only the PBE '
77 'functional. This calculation would use '
78 'the %s functional.' % xc.get_setup_name())
79 else:
80 setupdata = SetupData.find_and_read_path(symbol,
81 xc.get_setup_name(),
82 setuptype=type,
83 world=world)
85 if hasattr(setupdata, 'build'):
86 # It is not so nice that we have hubbard_u floating around here.
87 # For example, none of the other setup types are aware
88 # of hubbard u, so they silently ignore it!
89 if isinstance(setupdata, SetupData):
90 kwargs = dict(backwards_compatible=backwards_compatible)
91 else:
92 kwargs = {}
93 setup = LeanSetup(
94 setupdata.build(xc, lmax, basis, filter, **kwargs),
95 hubbard_u=hubbard_u)
96 return setup
97 else:
98 return setupdata
101def correct_occ_numbers(f_j,
102 degeneracy_j,
103 jsorted,
104 correction: float,
105 eps=1e-12) -> None:
106 """Correct f_j ndarray in-place."""
108 if correction > 0:
109 # Add electrons to the lowest eigenstates:
110 for j in jsorted:
111 c = min(correction, degeneracy_j[j] - f_j[j])
112 f_j[j] += c
113 correction -= c
114 if correction < eps:
115 break
116 elif correction < 0:
117 # Add electrons to the highest eigenstates:
118 for j in jsorted[::-1]:
119 c = min(-correction, f_j[j])
120 f_j[j] -= c
121 correction += c
122 if correction > -eps:
123 break
126class LocalCorrectionVar:
127 """Class holding data for local the calculation of local corr."""
128 def __init__(self, s=None):
129 """Initialize our data."""
130 for work_key in ('nq', 'lcut', 'n_qg', 'nt_qg', 'nc_g', 'nct_g',
131 'rgd2', 'Delta_lq', 'T_Lqp'):
132 if s is None or not hasattr(s, work_key):
133 setattr(self, work_key, None)
134 else:
135 setattr(self, work_key, getattr(s, work_key))
138class BaseSetup:
139 """Mixin-class for setups.
141 This makes it possible to inherit the most important methods without
142 the cumbersome constructor of the ordinary Setup class.
144 Maybe this class will be removed in the future, or it could be
145 made a proper base class with attributes and so on."""
147 orbital_free = False
148 hubbard_u = None # XXX remove me
149 is_pseudo = False
151 def print_info(self, text):
152 self.data.print_info(text, self)
154 def get_basis_description(self):
155 return self.basis.get_description()
157 def get_partial_waves_for_atomic_orbitals(self):
158 """Get those states phit that represent a real atomic state.
160 This typically corresponds to the (truncated) partial waves (PAW) or
161 a single-zeta basis."""
163 # XXX ugly hack for pseudopotentials:
164 if not hasattr(self, 'pseudo_partial_waves_j'):
165 return []
167 # The zip may cut off part of phit_j if there are more states than
168 # projectors. This should be the correct behaviour for all the
169 # currently supported PAW/pseudopotentials.
170 partial_waves_j = []
171 for n, phit in zip(self.n_j, self.pseudo_partial_waves_j):
172 if n > 0:
173 partial_waves_j.append(phit)
175 assert all(n > 0 for n in self.n_j[:len(partial_waves_j)])
177 return partial_waves_j
179 def calculate_initial_occupation_numbers(self, magmom, hund, charge,
180 nspins, f_j=None):
181 """If f_j is specified, custom occupation numbers will be used.
183 Hund rules disabled if so."""
184 nao = self.nao
185 f_sI = np.zeros((nspins, nao))
187 assert (not hund) or f_j is None
188 if f_j is None:
189 f_j = self.f_j
190 f_j = np.array(f_j, float)
191 l_j = np.array(self.l_j)
193 if hasattr(self, 'data') and hasattr(self.data, 'eps_j'):
194 eps_j = np.array(self.data.eps_j)
195 else:
196 eps_j = np.ones(len(self.n_j))
197 # Bound states:
198 for j, n in enumerate(self.n_j):
199 if n > 0:
200 eps_j[j] = -1.0
202 deg_j = 2 * (2 * l_j + 1)
204 # Sort after:
205 #
206 # 1) empty state (f == 0)
207 # 2) open shells (d - f)
208 # 3) eigenvalues (e)
210 states = []
211 for j, (f, d, e) in enumerate(zips(f_j, deg_j, eps_j, strict=False)):
212 if e < 0.0:
213 states.append((f == 0, d - f, e, j))
214 states.sort()
215 jsorted = [j for _, _, _, j in states]
217 # if len(l_j) == 0:
218 # l_j = np.ones(1)
220 # distribute the charge to the radial orbitals
221 if nspins == 1:
222 assert magmom == 0.0
223 f_sj = np.array([f_j])
224 if not self.orbital_free:
225 correct_occ_numbers(f_sj[0], deg_j, jsorted, -charge)
226 else:
227 # ofdft degeneracy of one orbital is infinite
228 f_sj[0, 0] -= charge
229 else:
230 nval = f_j.sum() - charge
231 if np.abs(magmom) > nval:
232 raise RuntimeError(
233 'Magnetic moment larger than number ' +
234 f'of valence electrons (|{magmom:g}| > {nval:g})')
235 f_sj = 0.5 * np.array([f_j, f_j])
236 nup = 0.5 * (nval + magmom)
237 ndn = 0.5 * (nval - magmom)
238 deg_j //= 2
239 correct_occ_numbers(f_sj[0], deg_j, jsorted, nup - f_sj[0].sum())
240 correct_occ_numbers(f_sj[1], deg_j, jsorted, ndn - f_sj[1].sum())
242 for j, (n1, l1, f, f_s) in enumerate(zip(self.n_j, self.l_j,
243 f_j, f_sj.T)):
244 if not f_s.any() or n1 < 0:
245 continue
246 I = 0
247 for bf in self.basis.bf_j:
248 l = bf.l
249 n = bf.n
250 if (n, l) == (n1, l1):
251 break
252 I += 2 * l + 1
253 else: # no break
254 raise ValueError(f'Bad basis for {self.symbol}')
256 degeneracy = 2 * l + 1
258 if hund:
259 # Use Hunds rules:
260 # assert f == int(f)
261 f = int(f)
262 f_sI[0, I:I + min(f, degeneracy)] = 1.0 # spin up
263 f_sI[1, I:I + max(f - degeneracy, 0)] = 1.0 # spin down
264 if f < degeneracy:
265 magmom -= f
266 else:
267 magmom -= 2 * degeneracy - f
268 else:
269 for s in range(nspins):
270 f_sI[s, I:I + degeneracy] = f_s[s] / degeneracy
272 if hund and magmom != 0:
273 raise WrongMagmomForHundsRuleError(
274 f'Bad magnetic moment {magmom:g} for {self.symbol} atom!')
276 # print('fsi=', f_si)
277 return f_sI
279 def get_hunds_rule_moment(self, charge=0):
280 for M in range(10):
281 try:
282 self.calculate_initial_occupation_numbers(M, True, charge, 2)
283 except WrongMagmomForHundsRuleError:
284 pass
285 else:
286 return M
287 raise RuntimeError
289 def initialize_density_matrix(self, f_sI):
290 nspins, nao = f_sI.shape
291 ni = self.ni
293 D_sii = np.zeros((nspins, ni, ni))
294 D_sp = np.zeros((nspins, ni * (ni + 1) // 2))
296 I = 0
297 for bf in self.basis.bf_j:
298 f_sm = f_sI[:, I:I + 2 * bf.l + 1]
299 if f_sm.any():
300 i = 0
301 for n, l in zip(self.n_j, self.l_j):
302 if (n, l) == (bf.n, bf.l):
303 break
304 i += 2 * l + 1
305 else: # no break
306 raise ValueError(f'Bad basis for {self.symbol}')
308 if i < ni:
309 for m in range(2 * l + 1):
310 D_sii[:, i + m, i + m] = f_sm[:, m]
312 I += 2 * bf.l + 1
314 for s in range(nspins):
315 D_sp[s] = pack_density(D_sii[s])
316 return D_sp
318 def get_partial_waves(self):
319 """Return spline representation of partial waves and densities."""
321 l_j = self.l_j
323 # cutoffs
324 rcut2 = 2 * max(self.rcut_j)
325 gcut2 = self.rgd.ceil(rcut2)
327 data = self.data
329 # Construct splines:
330 nc_g = data.nc_g.copy()
331 nct_g = data.nct_g.copy()
332 tauc_g = data.tauc_g
333 tauct_g = data.tauct_g
334 nc = self.rgd.spline(nc_g, rcut2, points=1000)
335 nct = self.rgd.spline(nct_g, rcut2, points=1000)
336 if tauc_g is None:
337 tauc_g = np.zeros(nct_g.shape)
338 tauct_g = tauc_g
339 tauc = self.rgd.spline(tauc_g, rcut2, points=1000)
340 tauct = self.rgd.spline(tauct_g, rcut2, points=1000)
341 phi_j = []
342 phit_j = []
343 for j, (phi_g, phit_g) in enumerate(zips(data.phi_jg, data.phit_jg)):
344 l = l_j[j]
345 phi_g = phi_g.copy()
346 phit_g = phit_g.copy()
347 phi_g[gcut2:] = phit_g[gcut2:] = 0.0
348 phi_j.append(self.rgd.spline(phi_g, rcut2, l, points=100))
349 phit_j.append(self.rgd.spline(phit_g, rcut2, l, points=100))
350 return phi_j, phit_j, nc, nct, tauc, tauct
352 def four_phi_integrals(self):
353 """Calculate four-phi integral.
355 Calculate the integral over the product of four all electron
356 functions in the augmentation sphere, i.e.::
358 /
359 | d vr ( phi_i1 phi_i2 phi_i3 phi_i4
360 / - phit_i1 phit_i2 phit_i3 phit_i4 ),
362 where phi_i1 is an all electron function and phit_i1 is its
363 smooth partner.
364 """
365 if hasattr(self, 'I4_pp'):
366 return self.I4_pp
368 # radial grid
369 r2dr_g = self.rgd.r_g**2 * self.rgd.dr_g
371 phi_jg = self.data.phi_jg
372 phit_jg = self.data.phit_jg
374 # compute radial parts
375 nj = len(self.l_j)
376 R_jjjj = np.empty((nj, nj, nj, nj))
377 for j1 in range(nj):
378 for j2 in range(nj):
379 for j3 in range(nj):
380 for j4 in range(nj):
381 R_jjjj[j1, j2, j3, j4] = np.dot(
382 r2dr_g,
383 phi_jg[j1] * phi_jg[j2] *
384 phi_jg[j3] * phi_jg[j4] -
385 phit_jg[j1] * phit_jg[j2] *
386 phit_jg[j3] * phit_jg[j4])
388 # prepare for angular parts
389 L_i = []
390 j_i = []
391 for j, l in enumerate(self.l_j):
392 for m in range(2 * l + 1):
393 L_i.append(l**2 + m)
394 j_i.append(j)
395 ni = len(L_i)
396 # j_i is the list of j values
397 # L_i is the list of L (=l**2+m for 0<=m<2*l+1) values
398 # https://gpaw.readthedocs.io/devel/overview.html
400 G_LLL = gaunt(max(self.l_j))
402 # calculate the integrals
403 _np = ni * (ni + 1) // 2 # length for packing
404 self.I4_pp = np.empty((_np, _np))
405 p1 = 0
406 for i1 in range(ni):
407 L1 = L_i[i1]
408 j1 = j_i[i1]
409 for i2 in range(i1, ni):
410 L2 = L_i[i2]
411 j2 = j_i[i2]
412 p2 = 0
413 for i3 in range(ni):
414 L3 = L_i[i3]
415 j3 = j_i[i3]
416 for i4 in range(i3, ni):
417 L4 = L_i[i4]
418 j4 = j_i[i4]
419 self.I4_pp[p1, p2] = (np.dot(G_LLL[L1, L2],
420 G_LLL[L3, L4]) *
421 R_jjjj[j1, j2, j3, j4])
422 p2 += 1
423 p1 += 1
425 # To unpack into I4_iip do:
426 # from gpaw.utilities import unpack_hermitian
427 # I4_iip = np.empty((ni, ni, _np)):
428 # for p in range(_np):
429 # I4_iip[..., p] = unpack_hermitian(I4_pp[:, p])
431 return self.I4_pp
433 def get_default_nbands(self):
434 assert len(self.l_orb_J) == len(self.n_j), (self.l_orb_J, self.n_j)
435 return sum([2 * l + 1 for (l, n) in zips(self.l_orb_J, self.n_j)
436 if n > 0])
438 def calculate_coulomb_corrections(self, wn_lqg, wnt_lqg, wg_lg, wnc_g,
439 wmct_g):
440 """Calculate "Coulomb" energies."""
441 # Can we reduce the excessive parameter passing?
442 # Seems so ....
443 # Added instance variables
444 # T_Lqp = self.local_corr.T_Lqp
445 # n_qg = self.local_corr.n_qg
446 # Delta_lq = self.local_corr.Delta_lq
447 # nt_qg = self.local_corr.nt_qg
448 # Local variables derived from instance variables
449 _np = self.ni * (self.ni + 1) // 2 # change to inst. att.?
450 mct_g = self.local_corr.nct_g + self.Delta0 * self.g_lg[0] # s.a.
451 rdr_g = self.local_corr.rgd2.r_g * \
452 self.local_corr.rgd2.dr_g # change to inst. att.?
454 A_q = 0.5 * (np.dot(wn_lqg[0], self.local_corr.nc_g) + np.dot(
455 self.local_corr.n_qg, wnc_g))
456 A_q -= sqrt(4 * pi) * self.Z * np.dot(self.local_corr.n_qg, rdr_g)
457 A_q -= 0.5 * (np.dot(wnt_lqg[0], mct_g) +
458 np.dot(self.local_corr.nt_qg, wmct_g))
459 A_q -= 0.5 * (np.dot(mct_g, wg_lg[0]) +
460 np.dot(self.g_lg[0], wmct_g)) * \
461 self.local_corr.Delta_lq[0]
462 M_p = np.dot(A_q, self.local_corr.T_Lqp[0])
464 A_lqq = []
465 for l in range(2 * self.local_corr.lcut + 1):
466 A_qq = 0.5 * np.dot(self.local_corr.n_qg, np.transpose(wn_lqg[l]))
467 A_qq -= 0.5 * np.dot(self.local_corr.nt_qg,
468 np.transpose(wnt_lqg[l]))
469 if l <= self.lmax:
470 A_qq -= 0.5 * np.outer(self.local_corr.Delta_lq[l],
471 np.dot(wnt_lqg[l], self.g_lg[l]))
472 A_qq -= 0.5 * np.outer(np.dot(self.local_corr.nt_qg,
473 wg_lg[l]),
474 self.local_corr.Delta_lq[l])
475 A_qq -= 0.5 * np.dot(self.g_lg[l], wg_lg[l]) * \
476 np.outer(self.local_corr.Delta_lq[l],
477 self.local_corr.Delta_lq[l])
478 A_lqq.append(A_qq)
480 M_pp = np.zeros((_np, _np))
481 L = 0
482 for l in range(2 * self.local_corr.lcut + 1):
483 for m in range(2 * l + 1): # m?
484 M_pp += np.dot(np.transpose(self.local_corr.T_Lqp[L]),
485 np.dot(A_lqq[l], self.local_corr.T_Lqp[L]))
486 L += 1
488 return M_p, M_pp # , V_p
490 def calculate_integral_potentials(self, func):
491 """Calculates a set of potentials using func."""
492 wg_lg = [func(self.g_lg[l], l)
493 for l in range(self.lmax + 1)]
494 wn_lqg = [np.array([func(self.local_corr.n_qg[q], l)
495 for q in range(self.nq)])
496 for l in range(2 * self.local_corr.lcut + 1)]
497 wnt_lqg = [np.array([func(self.local_corr.nt_qg[q], l)
498 for q in range(self.nq)])
499 for l in range(2 * self.local_corr.lcut + 1)]
500 wnc_g = func(self.local_corr.nc_g, l=0)
501 wnct_g = func(self.local_corr.nct_g, l=0)
502 wmct_g = wnct_g + self.Delta0 * wg_lg[0]
503 return wg_lg, wn_lqg, wnt_lqg, wnc_g, wnct_g, wmct_g
505 def calculate_yukawa_interaction(self, gamma):
506 """Calculate and return the Yukawa based interaction."""
508 # Solves the radial screened poisson equation for density n_g
509 def Yuk(n_g, l):
510 """Solve radial screened poisson for density n_g."""
511 return self.local_corr.rgd2.yukawa(n_g, l, gamma) * \
512 self.local_corr.rgd2.r_g * self.local_corr.rgd2.dr_g
514 return self.calculate_vvx_interactions(Yuk)
516 def calculate_erfc_interaction(self, omega):
517 """Calculate and return erfc based valence valence
518 exchange interactions."""
519 hse = RadialHSE(self.local_corr.rgd2, omega).screened_coulomb_dv
521 def erfc_interaction(n_g, l):
522 return hse(n_g, l)
523 return self.calculate_vvx_interactions(erfc_interaction)
525 def calculate_vvx_interactions(self, interaction):
526 """Calculate valence valence interactions for generic
527 interaction."""
528 (wg_lg, wn_lqg, wnt_lqg, wnc_g, wnct_g, wmct_g) = \
529 self.calculate_integral_potentials(interaction)
530 return self.calculate_coulomb_corrections(
531 wn_lqg, wnt_lqg, wg_lg, wnc_g, wmct_g)[1]
534class LeanSetup(BaseSetup):
535 """Setup class with minimal attribute set.
537 A setup-like class must define at least the attributes of this
538 class in order to function in a calculation."""
539 def __init__(self, s, hubbard_u=None):
540 """Copies precisely the necessary attributes of the Setup s."""
541 # Hubbard U is poked onto the setup in hacky ways.
542 # This needs cleaning.
543 self.hubbard_u = hubbard_u
545 self.N0_q = s.N0_q # Required for LDA+U.
546 self.type = s.type # required for writing to file
547 self.fingerprint = s.fingerprint # also req. for writing
548 self.filename = s.filename
550 self.symbol = s.symbol
551 self.Z = s.Z
552 self.Nv = s.Nv
553 self.Nc = s.Nc
555 self.ni = s.ni
556 self.nao = s.nao
558 self.pt_j = s.pt_j
560 self.pseudo_partial_waves_j = s.pseudo_partial_waves_j
561 self.basis_functions_J = s.basis_functions_J
563 self.Nct = s.Nct
564 self.nct = s.nct
566 self.lmax = s.lmax
567 self.ghat_l = s.ghat_l
568 self.vbar = s.vbar
570 self.Delta_pL = s.Delta_pL
571 self.Delta0 = s.Delta0
573 self.E = s.E
574 self.Kc = s.Kc
576 self.M = s.M
577 self.M_p = s.M_p
578 self.M_pp = s.M_pp
579 self.M_wpp = s.M_wpp
580 self.K_p = s.K_p
581 self.MB = s.MB
582 self.MB_p = s.MB_p
584 self.dO_ii = s.dO_ii
586 self.xc_correction = s.xc_correction
588 # Required to calculate initial occupations
589 self.f_j = s.f_j
590 self.n_j = s.n_j
591 self.l_j = s.l_j
592 self.l_orb_J = s.l_orb_J
593 self.nj = len(s.l_j)
594 self.nq = s.nq
596 self.data = s.data
598 # Below are things which are not really used all that much,
599 # i.e. shouldn't generally be necessary. Maybe we can make a system
600 # involving dictionaries for these "optional" parameters
602 # Required by print_info
603 self.rcutfilter = s.rcutfilter
604 self.rcore = s.rcore
605 self.basis = s.basis # we don't need nao if we use this instead
607 # XXX figure out better way to store these.
608 # Refactoring: We should delete this and use psit_j. However
609 # the code depends on psit_j being the *basis* functions sometimes.
610 if hasattr(s, 'pseudo_partial_waves_j'):
611 self.pseudo_partial_waves_j = s.pseudo_partial_waves_j
612 # Can also get rid of the phit_j splines if need be
614 self.N0_p = s.N0_p # req. by estimate_magnetic_moments
615 self.nabla_iiv = s.nabla_iiv # req. by lrtddft
616 self.rxnabla_iiv = s.rxnabla_iiv # req. by lrtddft and lrtddft2
618 # XAS stuff
619 self.phicorehole_g = s.phicorehole_g # should be optional
621 # Required to get all electron density
622 self.rgd = s.rgd
623 self.rcut_j = s.rcut_j
625 self.tauct = s.tauct # required by TPSS, MGGA
627 self.Delta_iiL = s.Delta_iiL # required with external potential
629 self.B_ii = s.B_ii # required for exact inverse overlap operator
630 self.dC_ii = s.dC_ii # required by time-prop tddft with apply_inverse
632 # Required by exx
633 self.X_p = s.X_p
634 self.X_wp = s.X_wp
635 self.ExxC = s.ExxC
636 self.ExxC_w = s.ExxC_w
638 # Required by yukawa rsf
639 self.X_pg = s.X_pg
641 # Required by electrostatic correction
642 self.dEH0 = s.dEH0
643 self.dEH_p = s.dEH_p
645 # Required by utilities/kspot.py (AllElectronPotential)
646 self.g_lg = s.g_lg
648 # Probably empty dictionary, required by GLLB
649 self.extra_xc_data = s.extra_xc_data
651 self.orbital_free = s.orbital_free
653 # Stuff required by Yukawa RSF to calculate Mg_pp at runtime
654 # the calcualtion of Mg_pp at rt is needed for dscf
655 if hasattr(s, 'local_corr'):
656 self.local_corr = s.local_corr
657 else:
658 self.local_corr = LocalCorrectionVar(s)
661class Setup(BaseSetup):
662 """Attributes:
664 ========== =====================================================
665 Name Description
666 ========== =====================================================
667 ``Z`` Charge
668 ``type`` Type-name of setup (eg. 'paw')
669 ``symbol`` Chemical element label (eg. 'Mg')
670 ``xcname`` Name of xc
671 ``data`` Container class for information on the the atom, eg.
672 Nc, Nv, n_j, l_j, f_j, eps_j, rcut_j.
673 It defines the radial grid by ng and beta, from which
674 r_g = beta * arange(ng) / (ng - arange(ng)).
675 It stores pt_jg, phit_jg, phi_jg, vbar_g
676 ========== =====================================================
679 Attributes for making PAW corrections
681 ============= ==========================================================
682 Name Description
683 ============= ==========================================================
684 ``Delta0`` Constant in compensation charge expansion coeff.
685 ``Delta_iiL`` Linear term in compensation charge expansion coeff.
686 ``Delta_pL`` Packed version of ``Delta_iiL``.
687 ``dO_ii`` Overlap coefficients
688 ``B_ii`` Projector function overlaps B_ii = <pt_i | pt_i>
689 ``dC_ii`` Inverse overlap coefficients
690 ``E`` Reference total energy of atom
691 ``M`` Constant correction to Coulomb energy
692 ``M_p`` Linear correction to Coulomb energy
693 ``M_pp`` 2nd order correction to Coulomb energy and Exx energy
694 ``M_wpp`` 2nd order correction to erfc screened Coulomb energy and Exx
695 energy for given w.
696 ``Kc`` Core kinetic energy
697 ``K_p`` Linear correction to kinetic energy
698 ``ExxC`` Core Exx energy
699 ``X_p`` Linear correction to Exx energy
700 ``MB`` Constant correction due to vbar potential
701 ``MB_p`` Linear correction due to vbar potential
702 ``dEH0`` Constant correction due to average electrostatic potential
703 ``dEH_p`` Linear correction due to average electrostatic potential
704 ``I4_iip`` Correction to integrals over 4 all electron wave functions
705 ``Nct`` Analytical integral of the pseudo core density ``nct``
706 ============= ==========================================================
708 It also has the attribute ``xc_correction`` which is an XCCorrection class
709 instance capable of calculating the corrections due to the xc functional.
712 Splines:
714 ========== ============================================
715 Name Description
716 ========== ============================================
717 ``pt_j`` Projector functions
718 ``phit_j`` Pseudo partial waves
719 ``vbar`` vbar potential
720 ``nct`` Pseudo core density
721 ``ghat_l`` Compensation charge expansion functions
722 ``tauct`` Pseudo core kinetic energy density
723 ========== ============================================
724 """
725 def __init__(self, data, xc, lmax=0, basis=None, filter=None,
726 backwards_compatible=True):
727 self.type = data.name
729 if not data.is_compatible(xc):
730 raise ValueError('Cannot use %s setup with %s functional' %
731 (data.setupname, xc.get_setup_name()))
733 self.symbol = data.symbol
734 self.data = data
736 self.Nc = data.Nc
737 self.Nv = data.Nv
738 self.Z = data.Z
739 l_j = self.l_j = data.l_j
740 self.l_orb_J = data.l_orb_J
741 n_j = self.n_j = data.n_j
742 self.f_j = data.f_j
743 self.eps_j = data.eps_j
744 nj = self.nj = len(l_j)
745 self.nq = nj * (nj + 1) // 2
746 rcut_j = self.rcut_j = data.rcut_j
748 self.ExxC = data.ExxC
749 self.ExxC_w = data.ExxC_w
750 self.X_p = data.X_p
751 self.X_wp = data.X_wp
753 self.X_pg = data.X_pg
755 self.orbital_free = data.orbital_free
757 pt_jg = data.pt_jg
758 phit_jg = data.phit_jg
759 phi_jg = data.phi_jg
761 self.fingerprint = data.fingerprint
762 self.filename = data.filename
764 rgd = self.rgd = data.rgd
765 r_g = rgd.r_g
766 dr_g = rgd.dr_g
768 self.lmax = lmax
770 # Attributes for run-time calculation of _Mg_pp
771 self.local_corr = LocalCorrectionVar(data)
773 rcutmax = max(rcut_j)
774 rcut2 = 2 * rcutmax
775 gcut2 = self.gcut2 = rgd.ceil(rcut2)
777 vbar_g = data.vbar_g
779 # if float(data.version) < 0.7 and data.generator_version < 2:
780 if data.generator_version < 2:
781 # Old-style Fourier-filtered datasets.
782 # Find Fourier-filter cutoff radius:
783 gcutfilter = rgd.get_cutoff(pt_jg[0])
785 elif filter:
786 rc = rcutmax
787 vbar_g = vbar_g.copy()
788 filter(rgd, rc, vbar_g)
790 pt_jg = [pt_g.copy() for pt_g in pt_jg]
791 for l, pt_g in zips(l_j, pt_jg):
792 filter(rgd, rc, pt_g, l)
794 for l in range(max(l_j) + 1):
795 J = [j for j, lj in enumerate(l_j) if lj == l]
796 A_nn = [[rgd.integrate(phit_jg[j1] * pt_jg[j2]) / 4 / pi
797 for j1 in J] for j2 in J]
798 B_nn = np.linalg.inv(A_nn)
799 pt_ng = np.dot(B_nn, [pt_jg[j] for j in J])
800 for n, j in enumerate(J):
801 pt_jg[j] = pt_ng[n]
802 gcutfilter = rgd.get_cutoff(pt_jg[0])
803 else:
804 gcutfilter = rgd.ceil(max(rcut_j))
806 if (vbar_g[gcutfilter:] != 0.0).any():
807 gcutfilter = rgd.get_cutoff(vbar_g)
808 assert r_g[gcutfilter] < 2.0 * max(rcut_j)
810 self.rcutfilter = rcutfilter = r_g[gcutfilter]
812 ni = 0
813 i = 0
814 j = 0
815 jlL_i = []
816 for l, n in zips(l_j, n_j):
817 for m in range(2 * l + 1):
818 jlL_i.append((j, l, l**2 + m))
819 i += 1
820 j += 1
821 ni = i
822 self.ni = ni
824 _np = ni * (ni + 1) // 2
826 lcut = max(l_j)
827 if 2 * lcut < lmax:
828 lcut = (lmax + 1) // 2
829 self.local_corr.lcut = lcut
831 self.B_ii = self.calculate_projector_overlaps(pt_jg)
833 self.fcorehole = data.fcorehole
834 self.lcorehole = data.lcorehole
836 # Construct splines:
837 self.vbar = rgd.spline(vbar_g, rcutfilter)
839 rcore, nc_g, nct_g, nct = self.construct_core_densities(
840 data, backwards_compatible)
841 self.rcore = rcore
842 self.nct = nct
844 # Construct splines for core kinetic energy density:
845 tauct_g = data.tauct_g
846 if tauct_g is not None:
847 self.tauct = rgd.spline(tauct_g, self.rcore)
848 else:
849 self.tauct = None
851 self.pt_j = self.create_projectors(pt_jg, rcutfilter)
853 partial_waves = self.create_basis_functions(phit_jg, rcut2, gcut2)
854 self.pseudo_partial_waves_j = partial_waves.tosplines()
856 if basis is None:
857 basis = partial_waves
858 basis_functions_J = self.pseudo_partial_waves_j
859 else:
860 basis_functions_J = basis.tosplines()
862 self.basis_functions_J = basis_functions_J
863 self.basis = basis
865 self.nao = 0
866 for bf in self.basis_functions_J:
867 l = bf.get_angular_momentum_number()
868 self.nao += 2 * l + 1
870 rgd2 = self.local_corr.rgd2 = rgd.new(gcut2)
871 r_g = rgd2.r_g
872 dr_g = rgd2.dr_g
873 phi_jg = np.array([phi_g[:gcut2].copy() for phi_g in phi_jg])
874 phit_jg = np.array([phit_g[:gcut2].copy() for phit_g in phit_jg])
875 self.local_corr.nc_g = nc_g = nc_g[:gcut2].copy()
876 self.local_corr.nct_g = nct_g = nct_g[:gcut2].copy()
877 vbar_g = vbar_g[:gcut2].copy()
879 extra_xc_data = dict(data.extra_xc_data)
880 # Cut down the GLLB related extra data
881 for key, item in extra_xc_data.items():
882 if len(item) == rgd.N:
883 extra_xc_data[key] = item[:gcut2].copy()
884 self.extra_xc_data = extra_xc_data
886 self.phicorehole_g = data.phicorehole_g
887 if self.phicorehole_g is not None:
888 self.phicorehole_g = self.phicorehole_g[:gcut2].copy()
890 self.local_corr.T_Lqp = self.calculate_T_Lqp(lcut, _np, nj, jlL_i)
891 # set the attributes directly?
892 (self.g_lg, self.local_corr.n_qg, self.local_corr.nt_qg,
893 self.local_corr.Delta_lq, self.Lmax, self.Delta_pL, self.Delta0,
894 self.N0_p) = self.get_compensation_charges(phi_jg, phit_jg, _np,
895 self.local_corr.T_Lqp)
897 # Solves the radial poisson equation for density n_g
898 def H(n_g, l):
899 return rgd2.poisson(n_g, l) * r_g * dr_g
901 (wg_lg, wn_lqg, wnt_lqg, wnc_g, wnct_g, wmct_g) = \
902 self.calculate_integral_potentials(H)
903 self.wg_lg = wg_lg
905 rdr_g = r_g * dr_g
906 dv_g = r_g * rdr_g
907 A = 0.5 * np.dot(nc_g, wnc_g)
908 A -= sqrt(4 * pi) * self.Z * np.dot(rdr_g, nc_g)
909 mct_g = nct_g + self.Delta0 * self.g_lg[0]
910 # wmct_g = wnct_g + self.Delta0 * wg_lg[0]
911 A -= 0.5 * np.dot(mct_g, wmct_g)
912 self.M = A
913 self.MB = -np.dot(dv_g * nct_g, vbar_g)
915 AB_q = -np.dot(self.local_corr.nt_qg, dv_g * vbar_g)
916 self.MB_p = np.dot(AB_q, self.local_corr.T_Lqp[0])
918 # Correction for average electrostatic potential:
919 #
920 # dEH = dEH0 + dot(D_p, dEH_p)
921 #
922 self.dEH0 = sqrt(4 * pi) * (wnc_g - wmct_g -
923 sqrt(4 * pi) * self.Z * r_g * dr_g).sum()
924 dEh_q = (wn_lqg[0].sum(1) - wnt_lqg[0].sum(1) -
925 self.local_corr.Delta_lq[0] * wg_lg[0].sum())
926 self.dEH_p = np.dot(dEh_q, self.local_corr.T_Lqp[0]) * sqrt(4 * pi)
928 M_p, M_pp = self.calculate_coulomb_corrections(wn_lqg, wnt_lqg,
929 wg_lg, wnc_g, wmct_g)
930 self.M_p = M_p
931 self.M_pp = M_pp
933 self.M_wpp = {}
934 for omega in self.ExxC_w:
935 self.M_wpp[omega] = self.calculate_erfc_interaction(omega)
937 self.Kc = data.e_kinetic_core - data.e_kinetic
938 self.M -= data.e_electrostatic
939 self.E = data.e_total
941 Delta0_ii = unpack_hermitian(self.Delta_pL[:, 0].copy())
942 self.dO_ii = data.get_overlap_correction(Delta0_ii)
943 self.dC_ii = self.get_inverse_overlap_coefficients(self.B_ii,
944 self.dO_ii)
946 self.Delta_iiL = np.zeros((ni, ni, self.Lmax))
947 for L in range(self.Lmax):
948 self.Delta_iiL[:, :, L] = unpack_hermitian(
949 self.Delta_pL[:, L].copy())
951 self.Nct = data.get_smooth_core_density_integral(self.Delta0)
952 self.K_p = data.get_linear_kinetic_correction(self.local_corr.T_Lqp[0])
954 self.ghat_l = [rgd2.spline(g_g, rcut2, l, 50)
955 for l, g_g in enumerate(self.g_lg)]
957 self.xc_correction = data.get_xc_correction(rgd2, xc, gcut2, lcut)
958 self.nabla_iiv = self.get_derivative_integrals(rgd2, phi_jg, phit_jg)
959 self.rxnabla_iiv = self.get_magnetic_integrals(rgd2, phi_jg, phit_jg)
961 def create_projectors(self, pt_jg, rcut):
962 pt_j = []
963 for j, pt_g in enumerate(pt_jg):
964 l = self.l_j[j]
965 pt_j.append(self.rgd.spline(pt_g, rcut, l))
966 return pt_j
968 def get_inverse_overlap_coefficients(self, B_ii, dO_ii):
969 ni = len(B_ii)
970 xO_ii = np.dot(B_ii, dO_ii)
971 return -np.dot(dO_ii, np.linalg.inv(np.identity(ni) + xO_ii))
973 def calculate_T_Lqp(self, lcut, _np, nj, jlL_i):
974 Lcut = (2 * lcut + 1)**2
975 G_LLL = gaunt(max(self.l_j))[:, :, :Lcut]
976 LGcut = G_LLL.shape[2]
977 T_Lqp = np.zeros((Lcut, self.nq, _np))
978 p = 0
979 i1 = 0
980 for j1, l1, L1 in jlL_i:
981 for j2, l2, L2 in jlL_i[i1:]:
982 if j1 < j2:
983 q = j2 + j1 * nj - j1 * (j1 + 1) // 2
984 else:
985 q = j1 + j2 * nj - j2 * (j2 + 1) // 2
986 T_Lqp[:LGcut, q, p] = G_LLL[L1, L2]
987 p += 1
988 i1 += 1
989 return T_Lqp
991 def calculate_projector_overlaps(self, pt_jg):
992 """Compute projector function overlaps B_ii = <pt_i | pt_i>."""
993 nj = len(pt_jg)
994 B_jj = np.zeros((nj, nj))
995 for j1, pt1_g in enumerate(pt_jg):
996 for j2, pt2_g in enumerate(pt_jg):
997 B_jj[j1, j2] = self.rgd.integrate(pt1_g * pt2_g) / (4 * pi)
998 B_ii = np.zeros((self.ni, self.ni))
999 i1 = 0
1000 for j1, l1 in enumerate(self.l_j):
1001 for m1 in range(2 * l1 + 1):
1002 i2 = 0
1003 for j2, l2 in enumerate(self.l_j):
1004 for m2 in range(2 * l2 + 1):
1005 if l1 == l2 and m1 == m2:
1006 B_ii[i1, i2] = B_jj[j1, j2]
1007 i2 += 1
1008 i1 += 1
1009 return B_ii
1011 def get_compensation_charges(self, phi_jg, phit_jg, _np, T_Lqp):
1012 lmax = self.lmax
1013 gcut2 = self.gcut2
1014 rcut_j = self.rcut_j
1015 nq = self.nq
1017 rgd = self.local_corr.rgd2
1018 r_g = rgd.r_g
1019 dr_g = rgd.dr_g
1021 g_lg = self.data.create_compensation_charge_functions(lmax)
1023 n_qg = np.zeros((nq, gcut2))
1024 nt_qg = np.zeros((nq, gcut2))
1025 gcut_q = np.zeros(nq, dtype=int)
1026 N0_q = np.zeros(nq)
1027 q = 0 # q: common index for j1, j2
1028 for j1 in range(self.nj):
1029 for j2 in range(j1, self.nj):
1030 n_qg[q] = phi_jg[j1] * phi_jg[j2]
1031 nt_qg[q] = phit_jg[j1] * phit_jg[j2]
1033 gcut = rgd.ceil(min(rcut_j[j1], rcut_j[j2]))
1034 N0_q[q] = sum(n_qg[q, :gcut] * r_g[:gcut]**2 * dr_g[:gcut])
1035 gcut_q[q] = gcut
1037 q += 1
1039 self.gcut_q = gcut_q
1040 self.N0_q = N0_q
1042 Delta_lq = np.zeros((lmax + 1, nq))
1043 for l in range(lmax + 1):
1044 Delta_lq[l] = np.dot(n_qg - nt_qg, r_g**(2 + l) * dr_g)
1046 Lmax = (lmax + 1)**2
1047 Delta_pL = np.zeros((_np, Lmax))
1048 for l in range(lmax + 1):
1049 L = l**2
1050 for m in range(2 * l + 1):
1051 delta_p = np.dot(Delta_lq[l], T_Lqp[L + m])
1052 Delta_pL[:, L + m] = delta_p
1054 Delta0 = np.dot(self.local_corr.nc_g - self.local_corr.nct_g,
1055 r_g**2 * dr_g) - self.Z / sqrt(4 * pi)
1057 # Electron density inside augmentation sphere. Used for estimating
1058 # atomic magnetic moment:
1059 N0_p = N0_q @ T_Lqp[0] * sqrt(4 * pi)
1061 return (g_lg[:, :gcut2].copy(), n_qg, nt_qg,
1062 Delta_lq, Lmax, Delta_pL, Delta0, N0_p)
1064 def get_derivative_integrals(self, rgd, phi_jg, phit_jg):
1065 """Calculate PAW-correction matrix elements of nabla.
1067 ::
1069 / _ _ d _ ~ _ d ~ _
1070 | dr [phi (r) -- phi (r) - phi (r) -- phi (r)]
1071 / 1 dx 2 1 dx 2
1073 and similar for y and z."""
1074 # lmax needs to be at least 1 for evaluating
1075 # the Gaunt coefficients from derivatives
1076 lmax = max(1, max(self.l_j))
1077 G_LLL = gaunt(lmax)
1078 Y_LLv = nabla(lmax)
1080 r_g = rgd.r_g
1081 dr_g = rgd.dr_g
1082 nabla_iiv = np.empty((self.ni, self.ni, 3))
1083 if debug:
1084 nabla_iiv[:] = np.nan
1085 i1 = 0
1086 for j1 in range(self.nj):
1087 l1 = self.l_j[j1]
1088 nm1 = 2 * l1 + 1
1089 i2 = 0
1090 for j2 in range(self.nj):
1091 l2 = self.l_j[j2]
1092 nm2 = 2 * l2 + 1
1093 f1f2or = np.dot(phi_jg[j1] * phi_jg[j2] -
1094 phit_jg[j1] * phit_jg[j2], r_g * dr_g)
1095 dphidr_g = np.empty_like(phi_jg[j2])
1096 rgd.derivative_spline(phi_jg[j2], dphidr_g)
1097 dphitdr_g = np.empty_like(phit_jg[j2])
1098 rgd.derivative_spline(phit_jg[j2], dphitdr_g)
1099 f1df2dr = np.dot(phi_jg[j1] * dphidr_g -
1100 phit_jg[j1] * dphitdr_g, r_g**2 * dr_g)
1101 for v in range(3):
1102 Lv = 1 + (v + 2) % 3
1103 G_12 = G_LLL[Lv, l1**2:l1**2 + nm1, l2**2:l2**2 + nm2]
1104 Y_12 = Y_LLv[l1**2:l1**2 + nm1, l2**2:l2**2 + nm2, v]
1105 nabla_iiv[i1:i1 + nm1, i2:i2 + nm2, v] = (
1106 sqrt(4 * pi / 3) * (f1df2dr - l2 * f1f2or) * G_12
1107 + f1f2or * Y_12)
1108 i2 += nm2
1109 i1 += nm1
1110 if debug:
1111 assert not np.any(np.isnan(nabla_iiv))
1112 return nabla_iiv
1114 def get_magnetic_integrals(self, rgd, phi_jg, phit_jg):
1115 """Calculate PAW-correction matrix elements of r x nabla.
1117 ::
1119 / _ _ _ ~ _ ~ _
1120 | dr [phi (r) O phi (r) - phi (r) O phi (r)]
1121 / 1 x 2 1 x 2
1123 d d
1124 where O = y -- - z --
1125 x dz dy
1127 and similar for y and z."""
1128 # lmax needs to be at least 1 for evaluating
1129 # the Gaunt coefficients from derivatives
1130 lmax = max(1, max(self.l_j))
1131 G_LLL = gaunt(lmax)
1132 Y_LLv = nabla(2 * lmax)
1134 r_g = rgd.r_g
1135 dr_g = rgd.dr_g
1136 rxnabla_iiv = np.empty((self.ni, self.ni, 3))
1137 if debug:
1138 rxnabla_iiv[:] = np.nan
1139 i1 = 0
1140 for j1 in range(self.nj):
1141 l1 = self.l_j[j1]
1142 nm1 = 2 * l1 + 1
1143 i2 = 0
1144 for j2 in range(self.nj):
1145 l2 = self.l_j[j2]
1146 nm2 = 2 * l2 + 1
1147 f1f2or = np.dot(phi_jg[j1] * phi_jg[j2] -
1148 phit_jg[j1] * phit_jg[j2], r_g**2 * dr_g)
1149 for v in range(3):
1150 v1 = (v + 1) % 3
1151 v2 = (v + 2) % 3
1152 Lv1 = 1 + (v1 + 2) % 3
1153 Lv2 = 1 + (v2 + 2) % 3
1154 # term from radial wfs does not contribute
1155 # term from spherical harmonics derivatives
1156 G_12 = np.zeros((nm1, nm2))
1157 G_12 += np.dot(G_LLL[Lv1, l1**2:l1**2 + nm1, :],
1158 Y_LLv[:, l2**2:l2**2 + nm2, v2])
1159 G_12 -= np.dot(G_LLL[Lv2, l1**2:l1**2 + nm1, :],
1160 Y_LLv[:, l2**2:l2**2 + nm2, v1])
1161 rxnabla_iiv[i1:i1 + nm1, i2:i2 + nm2, v] = (
1162 sqrt(4 * pi / 3) * f1f2or * G_12)
1163 i2 += nm2
1164 i1 += nm1
1165 if debug:
1166 assert not np.any(np.isnan(rxnabla_iiv))
1167 return rxnabla_iiv
1169 def construct_core_densities(self, setupdata, backwards_compatible=True):
1170 rcore = self.data.find_core_density_cutoff(setupdata.nc_g)
1171 nct = self.rgd.spline(
1172 setupdata.nct_g, rcore,
1173 points=None if backwards_compatible else 256,
1174 backwards_compatible=backwards_compatible)
1175 return rcore, setupdata.nc_g, setupdata.nct_g, nct
1177 def create_basis_functions(self, phit_jg, rcut2, gcut2):
1178 # Cutoff for atomic orbitals used for initial guess:
1179 rcut3 = 8.0 # XXXXX Should depend on the size of the atom!
1180 gcut3 = self.rgd.ceil(rcut3)
1182 # We cut off the wave functions smoothly at rcut3 by the
1183 # following replacement:
1184 #
1185 # /
1186 # | f(r), r < rcut2
1187 # f(r) <- < f(r) - a(r) f(rcut3) - b(r) f'(rcut3), rcut2 < r < rcut3
1188 # | 0, r > rcut3
1189 # \
1190 #
1191 # where a(r) and b(r) are 4. order polynomials:
1192 #
1193 # a(rcut2) = 0, a'(rcut2) = 0, a''(rcut2) = 0,
1194 # a(rcut3) = 1, a'(rcut3) = 0
1195 # b(rcut2) = 0, b'(rcut2) = 0, b''(rcut2) = 0,
1196 # b(rcut3) = 0, b'(rcut3) = 1
1197 #
1198 r_g = self.rgd.r_g
1199 x = (r_g[gcut2:gcut3] - rcut2) / (rcut3 - rcut2)
1200 a_g = 4 * x**3 * (1 - 0.75 * x)
1201 b_g = x**3 * (x - 1) * (rcut3 - rcut2)
1203 basis_functions_J = []
1204 n_J = []
1205 for j, phit_g in enumerate(phit_jg):
1206 n = self.n_j[j]
1207 if n > 0:
1208 n_J.append(n)
1209 l = self.l_j[j]
1210 phit_g = phit_g.copy()
1211 phit = phit_g[gcut3]
1212 dphitdr = ((phit - phit_g[gcut3 - 1]) /
1213 (r_g[gcut3] - r_g[gcut3 - 1]))
1214 phit_g[gcut2:gcut3] -= phit * a_g + dphitdr * b_g
1215 phit_g[gcut3:] = 0.0
1216 basis_function = self.rgd.spline(phit_g, rcut3, l, points=100)
1217 basis_functions_J.append(basis_function)
1218 basis = PartialWaveBasis(self.symbol, basis_functions_J, n_J)
1219 return basis
1222class PartialWaveBasis(Basis): # yuckkk
1223 def __init__(self, symbol, phit_J, n_J):
1224 self._basis_functions_J = phit_J
1225 super().__init__(
1226 symbol, 'partial-waves',
1227 bf_j=[BasisFunction(n, phit.get_angular_momentum_number())
1228 for n, phit in zip(n_J, phit_J)])
1230 def tosplines(self):
1231 return self._basis_functions_J
1233 def get_description(self):
1234 template = 'Using partial waves for %s as LCAO basis'
1235 string = template % self.symbol
1236 return string
1239class Setups(list):
1240 """Collection of Setup objects. One for each distinct atom.
1242 Non-distinct atoms are those with the same atomic number, setup, and basis.
1244 Class attributes:
1246 ``nvalence`` Number of valence electrons.
1247 ``nao`` Number of atomic orbitals.
1248 ``Eref`` Reference energy.
1249 ``core_charge`` Core hole charge.
1250 """
1252 def __init__(self, Z_a, setup_types, basis_sets, xc, *,
1253 filter=None,
1254 world=None,
1255 backwards_compatible=True):
1256 list.__init__(self)
1257 symbols = [chemical_symbols[Z] for Z in Z_a]
1258 type_a = types2atomtypes(symbols, setup_types, default='paw')
1259 basis_a = types2atomtypes(symbols, basis_sets, default=None)
1261 for a, _type in enumerate(type_a):
1262 # Make basis files correspond to setup files.
1263 #
1264 # If the setup has a name (i.e. non-default _type), then
1265 # prepend that name to the basis name.
1266 #
1267 # Typically people might specify '11' as the setup but just
1268 # 'dzp' for the basis set. Here we adjust to
1269 # obtain, say, '11.dzp' which loads the correct basis set.
1270 #
1271 # There will be no way to obtain the original 'dzp' with
1272 # a custom-named setup except by loading directly from
1273 # BasisData.
1274 #
1275 # Due to the "szp(dzp)" syntax this is complicated!
1276 # The name has to go as "szp(name.dzp)".
1277 basis = basis_a[a]
1278 if isinstance(basis, str):
1279 if isinstance(_type, str):
1280 setupname = _type
1281 else:
1282 setupname = _type.name # _type is an object like SetupData
1283 # Drop DFT+U specification from type string if it is there:
1284 if hasattr(setupname, 'swapcase'):
1285 setupname = setupname.split(':')[0]
1287 # Basis names inherit setup names except default setups
1288 # and ghost atoms.
1289 if setupname != 'paw' and setupname != 'ghost':
1290 if setupname:
1291 if '(' in basis:
1292 reduced, name = basis.split('(')
1293 assert name.endswith(')')
1294 name = name[:-1]
1295 fullname = f'{reduced}({setupname}.{name})'
1296 else:
1297 fullname = f'{setupname}.{basis_a[a]}'
1298 basis_a[a] = fullname
1300 # Construct necessary PAW-setup objects:
1301 self.setups = {}
1302 natoms = {}
1303 Mcumulative = 0
1304 self.M_a = []
1305 self.id_a = list(zips(Z_a, type_a, basis_a))
1306 for id in self.id_a:
1307 setup = self.setups.get(id)
1308 if setup is None:
1309 Z, type, basis = id
1310 symbol = chemical_symbols[Z]
1311 setupdata = None
1312 if not isinstance(type, str):
1313 setupdata = type
1314 # Basis may be None (meaning that the setup decides), a string
1315 # (meaning we load the basis set now from a file) or an actual
1316 # pre-created Basis object (meaning we just pass it along)
1317 if isinstance(basis, str):
1318 basis = Basis.find(symbol, basis, world=world)
1319 setup = create_setup(symbol, xc, 2, type,
1320 basis, setupdata=setupdata,
1321 filter=filter, world=world,
1322 backwards_compatible=backwards_compatible)
1323 self.setups[id] = setup
1324 natoms[id] = 0
1325 natoms[id] += 1
1326 self.append(setup)
1327 self.M_a.append(Mcumulative)
1328 Mcumulative += setup.nao
1330 # Sum up ...
1331 self.nvalence = 0 # number of valence electrons
1332 self.nao = 0 # number of atomic orbitals
1333 self.Eref = 0.0 # reference energy
1334 self.core_charge = 0.0 # core hole charge
1335 for id, setup in self.setups.items():
1336 n = natoms[id]
1337 self.Eref += n * setup.E
1338 self.core_charge += n * (setup.Z - setup.Nv - setup.Nc)
1339 self.nvalence += n * setup.Nv
1340 self.nao += n * setup.nao
1342 self.dS = OverlapCorrections(self)
1344 self.backwards_compatible = backwards_compatible
1346 def __str__(self):
1347 # Write PAW setup information in order of appearance:
1348 ids = set()
1349 s = 'species:\n'
1350 for id in self.id_a:
1351 if id in ids:
1352 continue
1353 ids.add(id)
1354 setup = self.setups[id]
1355 output = StringIO()
1356 setup.print_info(functools.partial(print, file=output))
1357 txt = output.getvalue()
1358 txt += ' # ' + setup.get_basis_description().replace('\n',
1359 '\n # ')
1360 txt = txt.replace('\n', '\n ')
1361 s += ' ' + txt + '\n\n'
1363 s += f'Reference energy: {self.Eref * units.Hartree:.6f} # eV\n'
1364 return s
1366 def set_symmetry(self, symmetry):
1367 """Find rotation matrices for spherical harmonics."""
1368 # XXX It is ugly that we set self.atomrotations from here;
1369 # it would be better to return it to the caller.
1370 from gpaw.atomrotations import AtomRotations
1371 self.atomrotations = AtomRotations(self.setups, self.id_a, symmetry)
1373 def empty_atomic_matrix(self, ns, atom_partition, dtype=float):
1374 Dshapes_a = [(ns, setup.ni * (setup.ni + 1) // 2)
1375 for setup in self]
1376 return atom_partition.arraydict(Dshapes_a, dtype)
1378 def estimate_dedecut(self, ecut):
1379 from gpaw.utilities.ekin import dekindecut, ekin
1380 dedecut = 0.0
1381 e = {}
1382 for id in self.id_a:
1383 if id not in e:
1384 G, de, e0 = ekin(self.setups[id])
1385 e[id] = -dekindecut(G, de, ecut)
1386 dedecut += e[id]
1387 return dedecut
1389 def basis_indices(self):
1390 return FunctionIndices([setup.basis_functions_J for setup in self])
1392 def projector_indices(self):
1393 return FunctionIndices([setup.pt_j for setup in self])
1395 def create_pseudo_core_densities(self, domain, positions, atomdist,
1396 xp=np):
1397 spline_aj = []
1398 for setup in self:
1399 if setup.nct is None:
1400 spline_aj.append([])
1401 else:
1402 spline_aj.append([setup.nct])
1403 if self.backwards_compatible and hasattr(domain, 'ecut'):
1404 integrals = None
1405 else:
1406 # 0.0 will skip normalization:
1407 integrals = [setup.Nct if abs(setup.Nct) > 1e-12 else 0.0
1408 for setup in self]
1409 return domain.atom_centered_functions(
1410 spline_aj, positions,
1411 atomdist=atomdist,
1412 integrals=integrals,
1413 cut=True, xp=xp)
1415 def create_pseudo_core_ked(self,
1416 domain,
1417 positions,
1418 atomdist):
1419 return domain.atom_centered_functions(
1420 [[setup.tauct] for setup in self],
1421 positions,
1422 atomdist=atomdist,
1423 cut=True)
1425 def create_local_potentials(self, domain, positions, atomdist, xp=np):
1426 return domain.atom_centered_functions(
1427 [[setup.vbar] for setup in self],
1428 positions,
1429 atomdist=atomdist,
1430 xp=xp)
1432 def create_compensation_charges(self, domain, positions, atomdist=None,
1433 xp=np):
1434 if self.backwards_compatible and hasattr(domain, 'ecut'):
1435 integral = None
1436 else:
1437 integral = sqrt(4 * pi)
1438 return domain.atom_centered_functions(
1439 [setup.ghat_l for setup in self], positions,
1440 atomdist=atomdist,
1441 integrals=integral,
1442 xp=xp)
1444 def get_overlap_corrections(self, atomdist, xp, dtype=np.float64):
1445 if atomdist is getattr(self, '_atomdist', None):
1446 if self.dS_aii.data.dtype == dtype:
1447 return self.dS_aii
1448 self._atomdist = atomdist
1449 dS_aii = AtomArraysLayout([setup.dO_ii.shape for setup in self],
1450 atomdist=atomdist, dtype=dtype).empty()
1451 for a, dS_ii in dS_aii.items():
1452 dS_ii[:] = self[a].dO_ii
1453 self.dS_aii = dS_aii.to_xp(xp)
1454 return self.dS_aii
1456 def partial_wave_corrections(self) -> list[list[Spline]]:
1457 splines: dict[Setup, list[Spline]] = {}
1458 dphi_aj = []
1459 for setup in self:
1460 dphi_j = splines.get(setup)
1461 if dphi_j is None:
1462 rcut = max(setup.rcut_j) * 1.1
1463 gcut = setup.rgd.ceil(rcut)
1464 dphi_j = []
1465 for l, phi_g, phit_g in zips(setup.l_j,
1466 setup.data.phi_jg,
1467 setup.data.phit_jg):
1468 dphi_g = (phi_g - phit_g)[:gcut]
1469 dphi_j.append(setup.rgd.spline(dphi_g, rcut, l,
1470 points=200))
1471 splines[setup] = dphi_j
1472 dphi_aj.append(dphi_j)
1474 return dphi_aj
1477class FunctionIndices:
1478 def __init__(self, f_aj):
1479 nm_a = [0]
1480 for f_j in f_aj:
1481 nm = sum([2 * f.get_angular_momentum_number() + 1 for f in f_j])
1482 nm_a.append(nm)
1483 self.M_a = np.cumsum(nm_a)
1484 self.nm_a = np.array(nm_a[1:])
1485 self.max = self.M_a[-1]
1487 def __getitem__(self, a):
1488 return self.M_a[a], self.M_a[a + 1]
1491class CachedYukawaInteractions:
1492 def __init__(self, omega):
1493 self.omega = omega
1494 self._cache = {}
1496 def get_Mg_pp(self, setup):
1497 if setup not in self._cache:
1498 Mg_pp = setup.calculate_yukawa_interaction(self.omega)
1499 self._cache[setup] = Mg_pp
1500 return self._cache[setup]
1503def types2atomtypes(symbols, types, default):
1504 """Map a types identifier to a list with a type id for each atom.
1506 types can be a single str, or a dictionary mapping chemical
1507 symbols and/or atom numbers to a type identifier.
1508 If both a symbol key and atomnumber key relates to the same atom, then
1509 the atomnumber key is dominant.
1511 If types is a dictionary and contains the string 'default', this will
1512 be used as default type, otherwize input arg ``default`` is used as
1513 default.
1514 """
1515 natoms = len(symbols)
1516 if isinstance(types, str):
1517 return [types] * natoms
1519 # If present, 'default' will map to the default type,
1520 # else use the input default
1521 type_a = [types.get('default', default)] * natoms
1523 # First symbols ...
1524 for symbol, type in types.items():
1525 # Types are given either by strings or they are objects that
1526 # have a 'symbol' attribute (SetupData, Pseudopotential, Basis, etc.).
1527 assert isinstance(type, str) or hasattr(type, 'symbol')
1528 if isinstance(symbol, str):
1529 for a, symbol2 in enumerate(symbols):
1530 if symbol == symbol2:
1531 type_a[a] = type
1533 # and then atom indices
1534 for a, type in types.items():
1535 if isinstance(a, int):
1536 type_a[a] = type
1538 return type_a
1541if __name__ == '__main__':
1542 print("""\
1543This is not the setup.py you are looking for! This setup.py defines a
1544Setup class used to hold the atomic data needed for a specific atom.
1545For building the GPAW code you must use the setup.py setuptools script
1546at the root of the code tree. Just do "cd .." and you will be at the
1547right place.""")
1548 raise SystemExit