Coverage for eminus/xc/utils.py: 84.21%
171 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-18 08:43 +0000
1# SPDX-FileCopyrightText: 2021 The eminus developers
2# SPDX-License-Identifier: Apache-2.0
3"""Utility functions for exchange-correlation functionals."""
5import inspect
6import sys
8import numpy as np
10from .. import config
11from ..logger import log
12from ..utils import add_maybe_none
13from .gga_c_chachiyo import gga_c_chachiyo, gga_c_chachiyo_spin
14from .gga_c_pbe import gga_c_pbe, gga_c_pbe_spin
15from .gga_c_pbe_sol import gga_c_pbe_sol, gga_c_pbe_sol_spin
16from .gga_x_chachiyo import gga_x_chachiyo, gga_x_chachiyo_spin
17from .gga_x_pbe import gga_x_pbe, gga_x_pbe_spin
18from .gga_x_pbe_sol import gga_x_pbe_sol, gga_x_pbe_sol_spin
19from .lda_c_chachiyo import lda_c_chachiyo, lda_c_chachiyo_spin
20from .lda_c_chachiyo_mod import lda_c_chachiyo_mod, lda_c_chachiyo_mod_spin
21from .lda_c_pw import lda_c_pw, lda_c_pw_spin
22from .lda_c_pw_mod import lda_c_pw_mod, lda_c_pw_mod_spin
23from .lda_c_vwn import lda_c_vwn, lda_c_vwn_spin
24from .lda_x import lda_x, lda_x_spin
25from .lda_xc_gdsmfb import lda_xc_gdsmfb, lda_xc_gdsmfb_spin
28def get_xc(xc, n_spin, Nspin, dn_spin=None, tau=None, xc_params=None, dens_threshold=0):
29 """Handle and get exchange-correlation functionals.
31 Args:
32 xc: Exchange and correlation identifier.
33 n_spin: Real-space electronic densities per spin channel.
34 Nspin: Number of spin states.
36 Keyword Args:
37 dn_spin: Real-space gradient of densities per spin channel.
38 tau: Real-space kinetic energy densities per spin channel.
39 xc_params: Exchange-correlation functional parameters.
40 dens_threshold: Do not treat densities smaller than the threshold.
42 Returns:
43 Exchange-correlation energy density and potential.
44 """
45 if isinstance(xc, str):
46 xc = parse_functionals(xc)
47 f_exch, f_corr = xc
48 if xc_params is None:
49 xc_params = {}
51 # Only use non-zero values of the density
52 n = np.sum(n_spin, axis=0)
53 nz_mask = np.where(n > dens_threshold)
54 n_nz = n[nz_mask]
55 # Zeta is only needed for non-zero values of the density
56 zeta_nz = get_zeta(n_spin[:, nz_mask])
57 # dn_spin is only needed for non-zero values of the density
58 if dn_spin is not None:
59 dn_spin_nz = dn_spin[:, nz_mask[0], :]
60 else:
61 dn_spin_nz = None
63 def handle_functional(fxc):
64 """Calculate a given functional fxc, same for exchange and correlation."""
65 # Calculate with the libxc extra...
66 if ":" in fxc:
67 from ..extras.libxc import libxc_functional
69 fxc = fxc.split(":")[-1]
70 exc, vxc, vsigma, vtau = libxc_functional(fxc, n_spin, Nspin, dn_spin, tau, xc_params)
71 # ...or use an internal functional
72 else:
73 if Nspin == 2 and fxc != "mock_xc":
74 fxc += "_spin"
75 exc_nz, vxc_nz, vsigma_nz = IMPLEMENTED[fxc](
76 n_nz, zeta=zeta_nz, dn_spin=dn_spin_nz, Nspin=Nspin, **xc_params
77 )
78 # Map the non-zero values back to the right dimension
79 exc = np.zeros_like(n)
80 exc[nz_mask] = exc_nz
81 vxc = np.zeros_like(n_spin)
82 for s in range(Nspin):
83 vxc[s, nz_mask] = vxc_nz[s]
84 if vsigma_nz is not None:
85 vsigma = np.zeros((len(vsigma_nz), len(exc)))
86 for i in range(len(vsigma)):
87 vsigma[i, nz_mask] = vsigma_nz[i]
88 else:
89 vsigma = None
90 # There are no internal meta-GGAs
91 vtau = None
92 return exc, vxc, vsigma, vtau
94 ex, vx, vsigmax, vtaux = handle_functional(f_exch) # Calculate the exchange part
95 ec, vc, vsigmac, vtauc = handle_functional(f_corr) # Calculate the correlation part
96 return ex + ec, vx + vc, add_maybe_none(vsigmax, vsigmac), add_maybe_none(vtaux, vtauc)
99def get_exc(xc, n_spin, Nspin, dn_spin=None, tau=None, xc_params=None, dens_threshold=0):
100 """Get the exchange-correlation energy density.
102 This is a convenience function to interface :func:`~eminus.xc.utils.get_xc`.
104 Args:
105 xc: Exchange and correlation identifier.
106 n_spin: Real-space electronic densities per spin channel.
107 Nspin: Number of spin states.
109 Keyword Args:
110 dn_spin: Real-space gradient of densities per spin channel.
111 tau: Real-space kinetic energy densities per spin channel.
112 xc_params: Exchange-correlation functional parameters.
113 dens_threshold: Do not treat densities smaller than the threshold.
115 Returns:
116 Exchange-correlation energy potential.
117 """
118 exc, _, _, _ = get_xc(xc, n_spin, Nspin, dn_spin, tau, xc_params, dens_threshold)
119 return exc
122def get_vxc(xc, n_spin, Nspin, dn_spin=None, tau=None, xc_params=None, dens_threshold=0):
123 """Get the exchange-correlation potential.
125 This is a convenience function to interface :func:`~eminus.xc.utils.get_xc`.
127 Args:
128 xc: Exchange and correlation identifier.
129 n_spin: Real-space electronic densities per spin channel.
130 Nspin: Number of spin states.
132 Keyword Args:
133 dn_spin: Real-space gradient of densities per spin channel.
134 tau: Real-space kinetic energy densities per spin channel.
135 xc_params: Exchange-correlation functional parameters.
136 dens_threshold: Do not treat densities smaller than the threshold.
138 Returns:
139 Exchange-correlation energy density.
140 """
141 _, vxc, vsigma, vtau = get_xc(xc, n_spin, Nspin, dn_spin, tau, xc_params, dens_threshold)
142 return vxc, vsigma, vtau
145def parse_functionals(xc):
146 """Parse exchange-correlation functional strings to the internal format.
148 Args:
149 xc: Exchange and correlation identifier, separated by a comma.
151 Returns:
152 Exchange and correlation string.
153 """
154 # Check for combined aliases
155 try:
156 # Remove underscores when looking up in the dictionary
157 xc_ = xc.replace("_", "")
158 xc = ALIAS[xc_]
159 except KeyError:
160 pass
162 # Parse functionals
163 functionals = []
164 for f in xc.split(","):
165 if ":" in f or f in IMPLEMENTED:
166 f_xc = f
167 elif not f:
168 f_xc = "mock_xc"
169 else:
170 try:
171 # Remove underscores when looking up in the dictionary
172 f_ = f.replace("_", "")
173 f_xc = XC_MAP[f_]
174 except KeyError:
175 log.exception(f'No functional found for "{f}".')
176 raise
177 functionals.append(f_xc)
179 # If only one or no functional has been parsed append with mock functionals
180 functionals.extend("mock_xc" for _ in range(2 - len(functionals)))
181 return functionals
184def parse_xc_type(xc):
185 """Parse functional strings to identify the corresponding functional type.
187 Args:
188 xc: Exchange and correlation identifier, separated by a comma.
190 Returns:
191 Functional type.
192 """
193 xc_type = []
194 for func in xc:
195 if ":" in func:
196 xc_id = func.split(":")[-1]
197 # Try to parse the functional using pylibxc at first
198 try:
199 family = parse_xc_libxc(xc_id)
200 # Otherwise parse it with PySCF
201 except (ImportError, AssertionError):
202 family = parse_xc_pyscf(xc_id)
204 if family == 1:
205 xc_type.append("lda")
206 elif family == 2:
207 xc_type.append("gga")
208 elif family == 4:
209 xc_type.append("meta-gga")
210 else:
211 msg = "Unsupported functional family."
212 raise NotImplementedError(msg)
213 # Fall back to internal xc functionals
214 elif "gga" in func:
215 xc_type.append("gga")
216 else:
217 xc_type.append("lda")
219 # When mixing functional types use the higher level of theory
220 if xc_type[0] != xc_type[1]:
221 log.warning("Detected mixing of different functional types.")
222 if "meta-gga" in xc_type:
223 return "meta-gga"
224 return "gga"
225 return xc_type[0]
228def parse_xc_libxc(xc_id):
229 """Parse functional type by its ID using pylibxc.
231 Args:
232 xc_id: Functional ID or identifier.
234 Returns:
235 Functional type.
236 """
237 if not config.use_pylibxc:
238 raise AssertionError
239 import pylibxc
241 if not xc_id.isdigit():
242 xc_id = pylibxc.util.xc_functional_get_number(xc_id)
244 func = pylibxc.LibXCFunctional(int(xc_id), 1)
245 if func._needs_laplacian:
246 msg = "meta-GGAs that need a laplacian are not supported."
247 raise NotImplementedError(msg)
248 return func.get_family()
251def parse_xc_pyscf(xc_id):
252 """Parse functional type by its ID using PySCF.
254 Args:
255 xc_id: Functional ID or identifier.
257 Returns:
258 Functional type.
259 """
260 from pyscf.dft.libxc import is_gga, is_lda, is_meta_gga, needs_laplacian, XC_CODES
262 if not xc_id.isdigit():
263 xc_id = XC_CODES[xc_id.upper()]
265 if needs_laplacian(int(xc_id)):
266 msg = "meta-GGAs that need a laplacian are not supported."
267 raise NotImplementedError(msg)
268 # Use the same values as in parse_xc_libxc
269 if is_lda(xc_id):
270 return 1
271 if is_gga(xc_id):
272 return 2
273 if is_meta_gga(xc_id):
274 return 4
275 return -1
278def get_xc_defaults(xc):
279 """Get the default parameters and values for a given set of functionals.
281 Args:
282 xc: Exchange and correlation identifier, separated by a comma.
284 Returns:
285 Default parameters and values.
286 """
287 if isinstance(xc, str):
288 xc = parse_functionals(xc)
290 # Names of special kewyword arguments that should not be used in xc_params
291 SPECIAL_NAMES = ["dn_spin", "Nspin"]
293 params = {}
294 for func in xc:
295 # If pylibxc functionals are used determine the default values from it
296 if ":" in func:
297 # This only works for pylibxc, not with PySCF
298 if not config.use_pylibxc or "pylibxc" not in sys.modules:
299 msg = "ext_params only work with pylibxc as the libxc backend, not with pyscf."
300 raise NotImplementedError(msg)
301 from pylibxc import LibXCFunctional
303 fxc = func.split(":")[-1]
304 try:
305 f_xc = LibXCFunctional(int(fxc), 1)
306 except ValueError:
307 f_xc = LibXCFunctional(fxc, 1)
308 fxc_params = dict(zip(f_xc.get_ext_param_names(), f_xc.get_ext_param_default_values()))
310 # Analyze the signature for implemented functionals
311 if func in IMPLEMENTED:
312 sig = inspect.signature(IMPLEMENTED[func])
313 fxc_params = {
314 param.name: param.default
315 for param in sig.parameters.values()
316 if param.default is not inspect.Parameter.empty
317 }
319 # Remove special names from the parsed parameters
320 for special in SPECIAL_NAMES:
321 if special in fxc_params:
322 del fxc_params[special]
324 # Append all parameters, warn if a parameter has been used before
325 for name in fxc_params:
326 if name in params:
327 log.warning(
328 f'External parameter "{name}" is used in the exchange AND correlation part. '
329 "It will be passed to both functionals if used in xc_params."
330 )
331 params[name] = fxc_params[name]
332 return params
335def get_zeta(n_spin):
336 """Calculate the relative spin polarization.
338 Args:
339 n_spin: Real-space electronic densities per spin channel.
341 Returns:
342 Relative spin polarization.
343 """
344 # If only one spin is given return an array of ones as if the density only is in one channel
345 if len(n_spin) == 1:
346 return np.ones_like(n_spin[0])
347 return (n_spin[0] - n_spin[1]) / (n_spin[0] + n_spin[1])
350def mock_xc(n, Nspin=1, **kwargs):
351 """Mock exchange-correlation functional with no effect (spin-paired).
353 Args:
354 n: Real-space electronic density.
355 Nspin: Number of spin states.
357 Keyword Args:
358 **kwargs: Throwaway arguments.
360 Returns:
361 Mock exchange-correlation energy density and potential.
362 """
363 zeros = np.zeros_like(n)
364 return zeros, np.array([zeros] * Nspin), None
367#: Map functional names with their respective implementation.
368IMPLEMENTED = {
369 i.__name__: i
370 for i in (
371 mock_xc,
372 gga_c_chachiyo,
373 gga_c_chachiyo_spin,
374 gga_c_pbe,
375 gga_c_pbe_spin,
376 gga_c_pbe_sol,
377 gga_c_pbe_sol_spin,
378 gga_x_chachiyo,
379 gga_x_chachiyo_spin,
380 gga_x_pbe,
381 gga_x_pbe_spin,
382 gga_x_pbe_sol,
383 gga_x_pbe_sol_spin,
384 lda_x,
385 lda_x_spin,
386 lda_c_pw,
387 lda_c_pw_spin,
388 lda_c_pw_mod,
389 lda_c_pw_mod_spin,
390 lda_c_vwn,
391 lda_c_vwn_spin,
392 lda_c_chachiyo,
393 lda_c_chachiyo_spin,
394 lda_c_chachiyo_mod,
395 lda_c_chachiyo_mod_spin,
396 lda_xc_gdsmfb,
397 lda_xc_gdsmfb_spin,
398 )
399}
401#: Map functional shorthands to the actual functional names.
402XC_MAP = {
403 # lda_x
404 "1": "lda_x",
405 "s": "lda_x",
406 "lda": "lda_x",
407 "slater": "lda_x",
408 # lda_c_vwn
409 "7": "lda_c_vwn",
410 "vwn": "lda_c_vwn",
411 "vwn5": "lda_c_vwn",
412 # lda_c_pw
413 "12": "lda_c_pw",
414 "pw": "lda_c_pw",
415 "pw92": "lda_c_pw",
416 # lda_c_pw_mod
417 "13": "lda_c_pw_mod",
418 "pwmod": "lda_c_pw_mod",
419 "pw92mod": "lda_c_pw_mod",
420 # gga_x_pbe
421 "101": "gga_x_pbe",
422 "pbex": "gga_x_pbe",
423 # gga_x_pbe_sol
424 "116": "gga_x_pbe_sol",
425 "pbesolx": "gga_x_pbe_sol",
426 # gga_c_pbe
427 "130": "gga_c_pbe",
428 "pbec": "gga_c_pbe",
429 # gga_c_pbe_sol
430 "133": "gga_c_pbe_sol",
431 "pbesolc": "gga_c_pbe_sol",
432 # lda_c_chachiyo
433 "287": "lda_c_chachiyo",
434 "chachiyo": "lda_c_chachiyo",
435 # gga_x_chachiyo
436 "298": "gga_x_chachiyo",
437 "chachiyox": "gga_x_chachiyo",
438 # lda_c_chachiyo_mod
439 "307": "lda_c_chachiyo_mod",
440 "chachiyomod": "lda_c_chachiyo_mod",
441 # gga_c_chachiyo
442 "309": "gga_c_chachiyo",
443 "chachiyoc": "gga_c_chachiyo",
444 # lda_xc_gdsmfb
445 "577": "lda_xc_gdsmfb",
446 "ldaxcgdsmfb": "lda_xc_gdsmfb",
447}
449#: Dictionary of common functional aliases.
450ALIAS = {
451 "svwn": "slater,vwn5",
452 "svwn5": "slater,vwn5",
453 "spw92": "slater,pw92mod",
454 "pbe": "pbex,pbec",
455 "pbesol": "pbesolx,pbesolc",
456 "chachiyo": "chachiyox,chachiyoc",
457}