Coverage for gpaw/pseudopotential.py: 74%
213 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-12 00:18 +0000
1import numpy as np
2from scipy.special import erf
4from gpaw.atom.atompaw import AtomPAW
5from gpaw.atom.radialgd import EquidistantRadialGridDescriptor
6from gpaw.basis_data import Basis, BasisFunction
7from gpaw.setup import BaseSetup, LocalCorrectionVar
8from gpaw.spline import Spline
9from gpaw.utilities import divrl, hartree as hartree_solve
12null_spline = Spline.from_data(0, 1.0, [0., 0., 0.])
15# XXX Not used at the moment; see comment below about rgd splines.
16def projectors_to_splines(rgd, l_j, pt_jg, filter=None):
17 # This function exists because both HGH and SG15 needs to do
18 # exactly the same thing.
19 #
20 # XXX equal-range projectors still required for some reason
21 maxlen = max([len(pt_g) for pt_g in pt_jg])
22 pt_j = []
23 for l, pt1_g in zip(l_j, pt_jg):
24 pt2_g = np.zeros(maxlen)
25 pt2_g[:len(pt1_g)] = pt1_g
26 if filter is not None:
27 filter(rgd, rgd.r_g[maxlen], pt2_g, l=l)
28 pt2_g = divrl(pt2_g, l, rgd.r_g[:maxlen])
29 spline = rgd.spline(pt2_g, rgd.r_g[maxlen - 1], l=l)
30 pt_j.append(spline)
31 return pt_j
34# XXX not used at the moment
35def local_potential_to_spline(rgd, vbar_g, filter=None):
36 vbar_g = vbar_g.copy()
37 rcut = rgd.r_g[len(vbar_g) - 1]
38 if filter is not None:
39 filter(rgd, rcut, vbar_g, l=0)
40 # vbar = Spline.from_data(0, rcut, vbar_g)
41 vbar = rgd.spline(vbar_g, rgd.r_g[len(vbar_g) - 1], l=0)
42 return vbar
45def get_radial_hartree_energy(r_g, rho_g):
46 """Get energy of l=0 compensation charge on equidistant radial grid."""
48 # At least in some cases the zeroth point is moved to 1e-8 or so to
49 # prevent division by zero and the like, so:
50 dr = r_g[2] - r_g[1]
51 rho_r_dr_g = dr * r_g * rho_g
52 vh_r_g = np.zeros(len(r_g)) # "r * vhartree"
53 hartree_solve(0, rho_r_dr_g, r_g, vh_r_g)
54 return 2.0 * np.pi * (rho_r_dr_g * vh_r_g).sum()
57def screen_potential(r, v, charge, rcut=None, a=None):
58 """Split long-range potential into short-ranged contributions.
60 The potential v is a long-ranted potential with the asymptotic form Z/r
61 corresponding to the given charge.
63 Return a potential vscreened and charge distribution rhocomp such that
65 v(r) = vscreened(r) + vHartree[rhocomp](r).
67 The returned quantities are truncated to a reasonable cutoff radius.
68 """
69 vr = v * r + charge
71 if rcut is None:
72 err = 0.0
73 i = len(vr)
74 while err < 1e-4:
75 # Things can be a bit sensitive to the threshold. The O.pz-mt
76 # setup gets 20-30 Bohr long compensation charges if it's 1e-6.
77 i -= 1
78 err = abs(vr[i])
79 i += 1
81 icut = np.searchsorted(r, r[i] * 1.1)
82 else:
83 icut = np.searchsorted(r, rcut)
84 rcut = r[icut]
85 rshort = r[:icut].copy()
86 if rshort[0] < 1e-16:
87 rshort[0] = 1e-10
89 if a is None:
90 a = rcut / 5.0 # XXX why is this so important?
91 vcomp = np.zeros_like(rshort)
92 vcomp = charge * erf(rshort / (np.sqrt(2.0) * a)) / rshort
93 # XXX divide by r
94 rhocomp = charge * (np.sqrt(2.0 * np.pi) * a)**(-3) * \
95 np.exp(-0.5 * (rshort / a)**2)
96 vscreened = v[:icut] + vcomp
97 return vscreened, rhocomp
100def figure_out_valence_states(ppdata):
101 from gpaw.atom.configurations import configurations
102 from ase.data import chemical_symbols
103 # ppdata.symbol may not be a chemical symbol so use Z
104 chemical_symbol = chemical_symbols[ppdata.Z]
105 Z, config = configurations[chemical_symbol]
106 assert Z == ppdata.Z
108 # Okay, we need to figure out occupations f_ln when we don't know
109 # any info about existing states on the pseudopotential.
110 #
111 # The plan is to loop over all states and count until only the correct
112 # number of valence electrons "remain".
113 nelectrons = 0
114 ncore = ppdata.Z - ppdata.Nv
116 energies = [c[3] for c in config]
117 args = np.argsort(energies)
118 config = list(np.array(config, dtype=object)[args])
120 nelectrons = 0
121 ncore = ppdata.Z - ppdata.Nv
122 assert ppdata.Nv > 0
123 iterconfig = iter(config)
124 if ncore > 0:
125 for n, l, occ, eps in iterconfig:
126 nelectrons += occ
127 if nelectrons == ncore:
128 break
129 elif nelectrons >= ncore:
130 raise ValueError('Cannot figure out what states should exist '
131 'on this pseudopotential.')
133 f_ln = {}
134 l_j = []
135 f_j = []
136 n_j = []
137 for n, l, occ, eps in iterconfig:
138 f_ln.setdefault(l, []).append(occ)
139 l_j.append(l)
140 f_j.append(occ)
141 n_j.append(n)
142 lmax = max(f_ln.keys())
143 f_ln = [f_ln.get(l, []) for l in range(lmax + 1)]
144 return n_j, l_j, f_j, f_ln
147def generate_basis_functions(ppdata):
148 class SimpleBasis(Basis):
149 def __init__(self, symbol, l_j, n_j):
150 rgd = EquidistantRadialGridDescriptor(0.02, 160)
151 super().__init__(symbol, 'simple', rgd=rgd, generatordata='simple')
152 bf_j = self.bf_j
153 rcgauss = rgd.r_g[-1] / 3.0
154 gauss_g = np.exp(-(rgd.r_g / rcgauss)**2.0)
155 for l, n in zip(l_j, n_j):
156 phit_g = rgd.r_g**l * gauss_g
157 norm = (rgd.integrate(phit_g**2) / (4 * np.pi))**0.5
158 phit_g /= norm
159 bf = BasisFunction(n, l, rgd.r_g[-1], phit_g, 'gaussian')
160 bf_j.append(bf)
161 # l_orb_J = [state.l for state in self.data['states']]
162 b1 = SimpleBasis(ppdata.symbol, ppdata.l_orb_J, ppdata.n_j)
163 apaw = AtomPAW(ppdata.symbol, [ppdata.f_ln], h=0.05, rcut=9.0,
164 basis={ppdata.symbol: b1},
165 setups={ppdata.symbol: ppdata},
166 maxiter=60,
167 txt=None)
168 basis = apaw.extract_basis_functions(ppdata.n_j, ppdata.l_j)
169 return basis
172def pseudoplot(pp, show=True):
173 import matplotlib.pyplot as plt
175 fig = plt.figure()
176 wfsax = fig.add_subplot(221)
177 ptax = fig.add_subplot(222)
178 vax = fig.add_subplot(223)
179 rhoax = fig.add_subplot(224)
181 def spline2grid(spline):
182 rcut = spline.get_cutoff()
183 r = np.linspace(0.0, rcut, 2000)
184 return r, spline.map(r)
186 for phit in pp.phit_j:
187 r, y = spline2grid(phit)
188 wfsax.plot(r, y, label='wf l=%d' % phit.get_angular_momentum_number())
190 for pt in pp.pt_j:
191 r, y = spline2grid(pt)
192 ptax.plot(r, y, label='pr l=%d' % pt.get_angular_momentum_number())
194 for ghat in pp.ghat_l:
195 r, y = spline2grid(ghat)
196 rhoax.plot(r, y, label='cc l=%d' % ghat.get_angular_momentum_number())
198 r, y = spline2grid(pp.vbar)
199 vax.plot(r, y, label='vbar')
201 vax.set_ylabel('potential')
202 rhoax.set_ylabel('density')
203 wfsax.set_ylabel('wfs')
204 ptax.set_ylabel('projectors')
206 for ax in [vax, rhoax, wfsax, ptax]:
207 ax.legend()
209 if show:
210 plt.show()
213class PseudoPotential(BaseSetup):
214 is_pseudo = True
216 def __init__(self, data, basis=None, filter=None):
217 self.data = data
219 self.N0_q = None
221 self.filename = None
222 self.fingerprint = None
223 self.symbol = data.symbol
224 self.type = data.name
226 self.Z = data.Z
227 self.Nv = data.Nv
228 self.Nc = data.Nc
230 self.f_j = data.f_j
231 self.n_j = data.n_j
232 self.l_j = data.l_j
233 self.l_orb_J = data.l_orb_J
234 self.nj = len(data.l_j)
235 self.nq = self.nj * (self.nj + 1) // 2
237 self.ni = sum([2 * l + 1 for l in data.l_j])
238 # self.pt_j = projectors_to_splines(data.rgd, data.l_j, data.pt_jg,
239 # filter=filter)
240 self.pt_j = data.get_projectors()
242 if len(self.pt_j) == 0:
243 assert False # not sure yet about the consequences of
244 # cleaning this up in the other classes
245 self.l_j = [0]
246 self.pt_j = [null_spline]
248 if basis is None:
249 basis = data.create_basis_functions()
251 self.basis_functions_J = basis.tosplines()
253 # We declare (for the benefit of the wavefunctions reuse method)
254 # that we have no PAW projectors as such. This makes the
255 # 'paw' wfs reuse method a no-op.
256 self.pseudo_partial_waves_j = []
258 self.basis = basis
259 self.nao = sum([2 * phit.get_angular_momentum_number() + 1
260 for phit in self.basis_functions_J])
262 self.Nct = 0.0
263 self.nct = null_spline
265 self.lmax = 0
267 self.xc_correction = None
269 r, l_comp, g_comp = data.get_compensation_charge_functions()
270 assert l_comp == [0] # Presumably only spherical charges
271 self.ghat_l = [
272 Spline.from_data(l, r[-1], g) for l, g in zip(l_comp, g_comp)
273 ]
274 self.rcgauss = data.rcgauss
276 # accuracy is rather sensitive to this
277 # self.vbar = local_potential_to_spline(data.rgd, data.vbar_g,
278 # filter=filter)
279 self.vbar = data.get_local_potential()
280 # XXX HGH and UPF use different radial grids, and this for
281 # some reason makes it difficult to use the exact same code to
282 # construct vbar and projectors. This should be fixed since
283 # either type of rgd should be able to always produce a valid
284 # and equivalent spline transparently.
286 _np = self.ni * (self.ni + 1) // 2
287 self.Delta0 = data.Delta0
288 self.Delta_pL = np.zeros((_np, 1))
290 self.E = 0.0
291 self.Kc = 0.0
292 self.M = -data.Eh_compcharge
293 self.M_p = np.zeros(_np)
294 self.M_pp = np.zeros((_np, _np))
295 self.M_wpp = {}
297 self.K_p = data.expand_hamiltonian_matrix()
298 self.MB = 0.0
299 self.MB_p = np.zeros(_np)
300 self.dO_ii = np.zeros((self.ni, self.ni))
302 # We don't really care about these variables
303 self.rcutfilter = None
304 self.rcore = None
306 self.N0_p = np.zeros(_np) # not really implemented
307 if hasattr(data, 'nabla_iiv'):
308 self.nabla_iiv = data.nabla_iiv
309 else:
310 self.nabla_iiv = None
311 self.rxnabla_iiv = None
312 self.phicorehole_g = None
313 self.rgd = data.rgd
314 self.rcut_j = data.rcut_j
315 self.tauct = None
316 self.Delta_iiL = np.zeros((self.ni, self.ni, 1))
317 self.B_ii = None
318 self.dC_ii = None
319 self.X_p = None
320 self.X_wp = {}
321 self.X_pg = None
322 self.ExxC = None
323 self.ExxC_w = {}
324 self.dEH0 = 0.0
325 self.dEH_p = np.zeros(_np)
326 self.extra_xc_data = {}
328 self.wg_lg = None
329 self.g_lg = None
330 self.local_corr = LocalCorrectionVar(None)