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