Coverage for eminus / scf.py: 97.99%
348 statements
« prev ^ index » next coverage.py v7.12.0, created at 2025-11-21 14:20 +0000
1# SPDX-FileCopyrightText: 2021 The eminus developers
2# SPDX-License-Identifier: Apache-2.0
3"""SCF class definition."""
5import copy
6import logging
7import time
9from . import backend as xp
10from . import config
11from .band_minimizer import get_grad_unocc, scf_step_unocc
12from .band_minimizer import IMPLEMENTED as BAND_MINIMIZER
13from .dft import get_epsilon, get_n_spin, get_n_total, get_phi, guess_pseudo, guess_random, orth
14from .energies import Energy, get_Edisp, get_Eewald, get_Esic
15from .gga import get_grad_field, get_tau
16from .gth import GTH
17from .logger import create_logger, get_level
18from .minimizer import IMPLEMENTED as ALL_MINIMIZER
19from .potentials import get_pot_defaults, init_pot
20from .potentials import IMPLEMENTED as ALL_POTENTIALS
21from .tools import center_of_mass, get_spin_squared
22from .utils import BaseObject
23from .version import info
24from .xc import get_xc, get_xc_defaults, parse_functionals, parse_xc_type
27class SCF(BaseObject):
28 """Perform direct minimizations.
30 Args:
31 atoms: Atoms object.
33 Keyword Args:
34 xc: Comma-separated exchange-correlation functional description.
36 Adding "libxc:" before a functional will use the Libxc interface for that functional,
37 e.g., with :code:`libxc:mgga_x_scan,libxc:mgga_c_scan`.
38 pot: Type of potential.
40 Can be one of "GTH" (default), "Coulomb", "Harmonic", or "Ge". Alternatively, a path to
41 a directory containing custom GTH pseudopotential files can be given.
42 guess: Initial guess for the wave functions.
44 Can be one of "random" (default) or "pseudo". Adding "symm" to the string will use the
45 same coefficients for both spin channels, e.g., :code:`symm-rand`.
46 etol: Convergence tolerance of the total energy.
47 gradtol: Convergence tolerance of the gradient norm.
49 This tolerance will only be used in conjugate gradient methods.
50 sic: Calculate the Kohn-Sham Perdew-Zunger SIC energy at the end of the SCF.
51 disp: Calculate a dispersion correction.
53 A dictionary can be used to pass arguments to the respective
54 function, e.g., with :code:`{"version": "d3zero", "atm": False, "xc": "scan"}`.
55 opt: Dictionary to customize the minimization methods.
57 The keys can be chosen out of "sd", "lm", "pclm", "cg", "pccg", and "auto". Defaults to
58 :code:`{"auto": 250}`.
59 verbose: Level of output.
61 Can be one of "critical", "error", "warning", "info" (default), or "debug". An integer
62 value can be used as well, where larger numbers mean more output, starting from 0.
63 Defaults to the verbosity level of the Atoms object.
64 """
66 def __init__(
67 self,
68 atoms,
69 xc="lda,vwn",
70 pot="gth",
71 guess="random",
72 etol=1e-7,
73 gradtol=None,
74 sic=False,
75 disp=False,
76 opt=None,
77 verbose=None,
78 ):
79 """Initialize the SCF object."""
80 # Set opt here, better to not use mutable data types in signatures
81 if opt is None:
82 opt = {"auto": 250}
84 # Set the input parameters (the ordering is important)
85 self.atoms = atoms #: Atoms object.
86 self._log = create_logger(self) #: Logger object.
87 self.verbose = verbose #: Verbosity level.
88 self.xc = xc #: Exchange-correlation functional.
89 self.pot_params = {} #: Potential parameters.
90 self.pot = pot #: Used potential.
91 self.guess = guess #: Initial wave functions guess.
92 self.etol = etol #: Total energy convergence tolerance.
93 self.gradtol = gradtol #: Gradient norm convergence tolerance.
94 self.sic = sic #: Enables the SIC energy calculation.
95 self.disp = disp #: Enables the dispersion correction calculation.
96 self.opt = opt #: Minimization methods.
97 self.smear_update = 2 #: Steps after the smeared occupations are recalculated.
99 # Initialize other attributes
100 self.energies = Energy() #: Energy object holding energy contributions.
101 self.is_converged = False #: Determines the SCF object convergence.
102 self.W = None #: Unconstrained wave functions.
103 self.xc_params = {} #: Exchange-correlation functional parameters.
104 self.clear()
106 # ### Class properties ###
108 @property
109 def atoms(self):
110 """Atoms object."""
111 return self._atoms
113 @atoms.setter
114 def atoms(self, value):
115 # Build the Atoms object if necessary and make a copy
116 # This way the Atoms objects inside and outside the class are independent but both are build
117 if not value.is_built:
118 self._atoms = copy.deepcopy(value.build())
119 else:
120 self._atoms = copy.deepcopy(value)
122 @property
123 def xc(self):
124 """Exchange-correlation functional."""
125 return self._xc
127 @xc.setter
128 def xc(self, value):
129 self._xc = parse_functionals(value.lower())
130 # Determine the type of functional combinations
131 self._xc_type = parse_xc_type(self._xc)
132 if "mock_xc" in self._xc and "_xc_" not in "".join(self._xc).lower():
133 self._log.warning("Usage of mock functional detected.")
135 @property
136 def xc_params(self):
137 """Exchange-correlation functional parameters."""
138 return self._xc_params
140 @xc_params.setter
141 def xc_params(self, value):
142 # Check if some parameters are unused in the functionals
143 # This also ensures that we print a warning in case of overlapping parameters in x and c
144 if value != {} and value is not None:
145 not_used = value.keys() - self.xc_params_defaults.keys()
146 if len(not_used) > 0:
147 self._log.warning(f"Some xc_params are unused, namely: {', '.join(not_used)}.")
148 self._xc_params = value
150 @property
151 def pot(self):
152 """Used potential."""
153 return self._pot
155 @pot.setter
156 def pot(self, value):
157 if value.lower() in ALL_POTENTIALS:
158 self._pot = value.lower()
159 # Only set the pseudopotential type for GTH pseudopotentials
160 if self._pot == "gth":
161 if "gga" in self._xc_type:
162 self._psp = "pbe"
163 else:
164 self._psp = "pade"
165 # If pot is no supported potential treat it as a path to a directory containing GTH files
166 else:
167 self._log.info(f'Use the path "{value}" to search for GTH pseudopotential files.')
168 self._psp = value
169 self._pot = "gth"
170 # Build the potential
171 if self._pot == "gth":
172 self.gth = GTH(self)
173 self.Vloc = init_pot(self, self.pot_params)
175 @property
176 def pot_params(self):
177 """Potential parameters."""
178 return self._pot_params
180 @pot_params.setter
181 def pot_params(self, value):
182 # Check if some parameters are unused in the potential
183 if value != {} and value is not None:
184 not_used = value.keys() - self.pot_params_defaults.keys()
185 if len(not_used) > 0:
186 self._log.warning(f"Some pot_params are unused, namely: {', '.join(not_used)}.")
187 self._pot_params = value
188 # Update the local potential for the new parameters
189 if hasattr(self, "pot"):
190 self.Vloc = init_pot(self, self.pot_params)
192 @property
193 def guess(self):
194 """Initial wave functions guess."""
195 return self._guess
197 @guess.setter
198 def guess(self, value):
199 # Set the guess method
200 value = value.lower()
201 if "rand" in value:
202 self._guess = "random"
203 elif "pseudo" in value:
204 self._guess = "pseudo"
205 else:
206 msg = f"{value} is no valid initial guess."
207 raise ValueError(msg)
208 # Check if a symmetric or unsymmetric guess is selected
209 if "sym" in value and "unsym" not in value:
210 self._symmetric = True
211 else:
212 self._symmetric = False
214 @property
215 def opt(self):
216 """Minimization methods."""
217 return self._opt
219 @opt.setter
220 def opt(self, value):
221 # Set lowercase to all keys
222 value = {k.lower(): v for k, v in value.items()}
223 for opt in value:
224 if opt not in ALL_MINIMIZER:
225 msg = f'No minimizer found for "{opt}".'
226 raise KeyError(msg)
227 self._opt = value
229 @property
230 def verbose(self):
231 """Verbosity level."""
232 return self._verbose
234 @verbose.setter
235 def verbose(self, value):
236 # If no verbosity is given use the one from the Atoms object
237 if value is None:
238 value = self.atoms.verbose
239 self._verbose = get_level(value)
240 self._log.verbose = self._verbose
242 # ### Read-only properties ###
244 @property
245 def kpts(self):
246 """Pass-through to the KPoints object of the Atoms object."""
247 return self.atoms.kpts
249 @property
250 def pot_params_defaults(self):
251 """Get the default potential parameters."""
252 return get_pot_defaults(self.pot)
254 @property
255 def psp(self):
256 """Pseudopotential path."""
257 return self._psp
259 @property
260 def symmetric(self):
261 """Determines if the initial guess is the same for both spin channels."""
262 return self._symmetric
264 @property
265 def xc_type(self):
266 """Determines the exchange-correlation family."""
267 return self._xc_type
269 @property
270 def xc_params_defaults(self):
271 """Get the default exchange-correlation functional parameters."""
272 return get_xc_defaults(self.xc)
274 # ### Class methods ###
276 def run(self, **kwargs):
277 """Run the self-consistent field (SCF) calculation.
279 Keyword Args:
280 **kwargs: Pass-through keyword arguments.
282 Returns:
283 Total energy.
284 """
285 # Print some information about the calculation
286 if self._log.level <= logging.DEBUG:
287 info()
288 if ":" in "".join(self._xc) and config.use_pylibxc:
289 self._log.info(
290 "This calculation employs Libxc to evaluate density functionals. When using Libxc,"
291 "\nplease cite SoftwareX 7, 1 (2018) (doi:10.1016/j.softx.2017.11.002)."
292 )
293 self._log.debug(
294 f"\n--- Atoms information ---\n{self.atoms}\n"
295 f"\n--- Cell information ---\nCell vectors:\n{self.atoms.a} a0\n"
296 f"Sampling per axis: {self.atoms.s}\n"
297 f"Cut-off energy: {self.atoms.ecut} Eh\n"
298 f"\n--- State information ---\n{self.atoms.occ}\n"
299 f"\n--- Calculation information ---\n{self}\n\n--- SCF data ---"
300 )
302 # Calculate Ewald energy that only depends on the system geometry
303 self.energies.Eewald = get_Eewald(self.atoms)
304 # Build the initial wave function if there is no W to start from
305 if self.W is None:
306 if "random" in self.guess:
307 self.W = guess_random(self, symmetric=self.symmetric)
308 elif "pseudo" in self.guess:
309 self.W = guess_pseudo(self, symmetric=self.symmetric)
311 # Start the minimization procedures
312 self.clear()
313 Etots = []
314 for imin in self.opt:
315 # Call the minimizer
316 self._log.info(f"Start {ALL_MINIMIZER[imin].__name__}...")
317 start = time.perf_counter()
318 Elist = ALL_MINIMIZER[imin](self, self.opt[imin], **kwargs)
319 end = time.perf_counter()
320 # Save the minimizer results
321 self._opt_log[imin] = {}
322 self._opt_log[imin]["iter"] = len(Elist)
323 self._opt_log[imin]["time"] = end - start
324 Etots += Elist
325 # Do not start other minimizations if one converged
326 if self.is_converged:
327 break
328 if self.is_converged:
329 self._log.info(f"SCF converged after {len(Etots)} iterations.")
330 else:
331 self._log.warning("SCF not converged!")
333 # Calculate SIC energy if needed
334 if self.sic:
335 self.energies.Esic = get_Esic(self, self.Y)
336 # Calculate dispersion correction energy if needed
337 if isinstance(self.disp, dict):
338 self.energies.Edisp = get_Edisp(self, **self.disp)
339 elif self.disp:
340 self.energies.Edisp = get_Edisp(self)
342 # Print minimizer timings
343 self._log.debug("\n--- SCF results ---")
344 t_tot = 0
345 for imin in self._opt_log:
346 N = self._opt_log[imin]["iter"]
347 t = self._opt_log[imin]["time"]
348 t_tot += t
349 self._log.debug(
350 f"Minimizer: {imin}"
351 f"\nIterations: {N}"
352 f"\nTime: {t:.5f} s"
353 f"\nTime/Iteration: {t / N:.5f} s"
354 )
355 self._log.info(f"Total SCF time: {t_tot:.5f} s")
356 # Print the S^2 expectation value for unrestricted calculations
357 if self.atoms.unrestricted:
358 self._log.info(f"<S^2> = {get_spin_squared(self):.6e}")
359 # Print energy data
360 if self._log.level <= logging.DEBUG:
361 self._log.debug(
362 "\n--- Energy data ---\n"
363 f"Eigenenergies:\n{get_epsilon(self, self.W)}\n\n{self.energies}"
364 )
365 else:
366 self._log.info(f"Etot = {self.energies.Etot:.9f} Eh")
367 return self.energies.Etot
369 kernel = run
371 def converge_bands(self, **kwargs):
372 """Converge occupied bands after an SCF calculation.
374 Keyword Args:
375 **kwargs: Pass-through keyword arguments.
376 """
377 if not self.is_converged:
378 self._log.warning("The previous calculation has not been converged.")
380 # If new k-points have been set rebuild the atoms object and the potential
381 if not self.atoms.kpts.is_built or (
382 self.W is not None and len(self.W) != self.atoms.kpts.Nk
383 ):
384 self.atoms.build()
385 self.pot = self.pot
386 self.is_converged = False
388 # Build the initial wave function if there is no W to start from
389 if self.W is None or len(self.W) != self.atoms.kpts.Nk:
390 if "random" in self.guess:
391 self.W = guess_random(self, symmetric=self.symmetric)
392 elif "pseudo" in self.guess:
393 self.W = guess_pseudo(self, symmetric=self.symmetric)
395 self._log.info("Minimize occupied band energies...")
396 # Start the minimization procedures
397 Etots = []
398 for imin in self.opt:
399 # Call the minimizer
400 self._log.info(f"Start {BAND_MINIMIZER[imin].__name__}...")
401 start = time.perf_counter()
402 Elist, self.W = BAND_MINIMIZER[imin](self, self.W, self.opt[imin], **kwargs)
403 end = time.perf_counter()
404 # Save the minimizer results
405 self._opt_log[imin] = {}
406 self._opt_log[imin]["iter"] = len(Elist)
407 self._opt_log[imin]["time"] = end - start
408 Etots += Elist
409 # Do not start other minimizations if one converged
410 if self.is_converged:
411 break
412 if self.is_converged:
413 self._log.info(f"Band minimization converged after {len(Etots)} iterations.")
414 else:
415 self._log.warning("Band minimization not converged!")
417 # Print minimizer timings
418 self._log.debug("\n--- Band minimization results ---")
419 t_tot = 0
420 for imin in self._opt_log:
421 N = self._opt_log[imin]["iter"]
422 t = self._opt_log[imin]["time"]
423 t_tot += t
424 self._log.debug(
425 f"Minimizer: {imin}"
426 f"\nIterations: {N}"
427 f"\nTime: {t:.5f} s"
428 f"\nTime/Iteration: {t / N:.5f} s"
429 )
430 self._log.info(f"Total band minimization time: {t_tot:.5f} s")
432 # Converge empty bands automatically if needed
433 if self.atoms.occ.Nempty > 0:
434 self.converge_empty_bands(**kwargs)
435 return self
437 def converge_empty_bands(self, Nempty=None, **kwargs):
438 """Converge unoccupied bands after converging occ. bands.
440 Keyword Args:
441 Nempty: Number of empty states.
442 **kwargs: Pass-through keyword arguments.
443 """
444 if not self.is_converged:
445 self._log.warning("The previous calculation has not been converged.")
446 self.is_converged = False
448 if Nempty is None:
449 Nempty = self.atoms.occ.Nempty
451 # Build the initial wave functions
452 if self.Z is None:
453 if "random" in self.guess:
454 self.Z = guess_random(self, Nempty, symmetric=self.symmetric)
455 elif "pseudo" in self.guess:
456 self.Z = guess_pseudo(self, Nempty, symmetric=self.symmetric)
458 self._log.info("Minimize unoccupied band energies...")
459 # Start the minimization procedures
460 Etots = []
461 for imin in self.opt:
462 # Call the minimizer
463 self._log.info(f"Start {BAND_MINIMIZER[imin].__name__}...")
464 start = time.perf_counter()
465 Elist, self.Z = BAND_MINIMIZER[imin](
466 self, self.Z, self.opt[imin], cost=scf_step_unocc, grad=get_grad_unocc, **kwargs
467 )
468 end = time.perf_counter()
469 # Save the minimizer results
470 self._opt_log[imin] = {}
471 self._opt_log[imin]["iter"] = len(Elist)
472 self._opt_log[imin]["time"] = end - start
473 Etots += Elist
474 # Do not start other minimizations if one converged
475 if self.is_converged:
476 break
477 if self.is_converged:
478 self._log.info(f"Band minimization converged after {len(Etots)} iterations.")
479 else:
480 self._log.warning("Band minimization not converged!")
482 # Print minimizer timings
483 self._log.debug("\n--- Band minimization results ---")
484 t_tot = 0
485 for imin in self._opt_log:
486 N = self._opt_log[imin]["iter"]
487 t = self._opt_log[imin]["time"]
488 t_tot += t
489 self._log.debug(
490 f"Minimizer: {imin}"
491 f"\nIterations: {N}"
492 f"\nTime: {t:.5f} s"
493 f"\nTime/Iteration: {t / N:.5f} s"
494 )
495 self._log.info(f"Total band minimization time: {t_tot:.5f} s")
496 return self
498 def recenter(self, center=None):
499 """Recenter the system inside the cell.
501 Keyword Args:
502 center: Point to center the system around.
503 """
504 atoms = self.atoms
505 # Get the COM before centering the atoms
506 com = center_of_mass(atoms.pos)
507 # Run the recenter method of the atoms object
508 self.atoms.recenter(center=center)
509 if center is None:
510 dr = com - xp.sum(atoms.a, axis=0) / 2
511 else:
512 center = xp.asarray(center, dtype=float)
513 dr = com - center
515 # Shift orbitals and density
516 if self.W is not None:
517 self.W = atoms.T(self.W, dr=-dr)
518 # Transform the density to the reciprocal space, shift, and transform back
519 n = None
520 if self.n is not None:
521 Jn = atoms.J(self.n)
522 TJn = atoms.T(Jn, dr=-dr)
523 n = xp.real(atoms.I(TJn))
525 # Recalculate the potential since it depends on the structure factor
526 self.pot = self.pot
527 # Clear intermediate results to make sure no one uses the unshifted results
528 self.clear()
529 # Set the shifted density after calling the clearing function
530 self.n = n
531 return self
533 def clear(self):
534 """Initialize or clear intermediate results."""
535 self.Y = None # Orthogonal wave functions
536 self.Z = None # Unconstrained wave functions of unoccupied states
537 self.D = None # Orthogonal wave functions of unoccupied states
538 self.n = None #: Electronic density
539 self.n_spin = None # Electronic densities per spin
540 self.dn_spin = None # Gradient of electronic densities per spin
541 self.tau = None # Kinetic energy densities per spin
542 self.phi = None # Hartree field
543 self.exc = None # Exchange-correlation energy density
544 self.vxc = None # Exchange-correlation potential
545 self.vsigma = None # n times d exc/d |dn|^2
546 self.vtau = None # d exc/d tau
547 self._precomputed = {} # Dictionary of pre-computed values not to be saved
548 self._opt_log = {} # Log of the optimization procedure
549 return self
551 @staticmethod
552 def callback(scf, step):
553 """Callback function that will get called every SCF iteration.
555 This is just an empty function users can overwrite with their implementation.
556 Remember that the "auto" minimization can call the callback function twice when the pccg
557 step is not preferred.
559 Args:
560 scf: SCF object.
561 step: Optimization step.
562 """
564 def _precompute(self):
565 """Precompute fields stored in the SCF object."""
566 atoms = self.atoms
567 self.Y = orth(atoms, self.W)
568 self.n_spin = get_n_spin(atoms, self.Y)
569 self.n = get_n_total(atoms, self.Y, self.n_spin)
570 if "gga" in self.xc_type:
571 self.dn_spin = get_grad_field(atoms, self.n_spin)
572 if self.xc_type == "meta-gga":
573 self.tau = get_tau(atoms, self.Y)
574 self.phi = get_phi(atoms, self.n)
575 self.exc, self.vxc, self.vsigma, self.vtau = get_xc(
576 self.xc, self.n_spin, atoms.occ.Nspin, self.dn_spin, self.tau, self.xc_params
577 )
578 self._precomputed = {
579 "dn_spin": self.dn_spin,
580 "phi": self.phi,
581 "vxc": self.vxc,
582 "vsigma": self.vsigma,
583 "vtau": self.vtau,
584 }
585 return self
587 def __repr__(self):
588 """Print the most important parameters stored in the SCF object."""
589 # Use chr(10) to create a linebreak since backslashes are not allowed in f-strings
590 return (
591 f"XC functionals: {self.xc}\n"
592 f"Potential: {self.pot}\n"
593 f"{f'GTH files: {self.psp}' + chr(10) if self.pot == 'gth' else ''}"
594 f"Starting guess: {self.guess}\n"
595 f"Symmetric guess: {self.symmetric}\n"
596 f"Energy convergence tolerance: {self.etol} Eh\n"
597 f"Gradient convergence tolerance: {self.gradtol}\n"
598 f"Non-local potential: {self.gth.NbetaNL > 0 if self.pot == 'gth' else 'false'}\n"
599 f"Smearing: {self.atoms.occ.smearing > 0}\n"
600 f"Smearing update cycle: {self.smear_update}"
601 )
604class RSCF(SCF):
605 """SCF class for spin-paired systems.
607 Inherited from :class:`eminus.scf.SCF`.
609 In difference to the SCF class, this class will not build the original Atoms object, only the
610 one attributed to the class. Customized fillings could be overwritten when using this class.
611 """
613 @property
614 def atoms(self):
615 """Atoms object."""
616 return self._atoms
618 @atoms.setter
619 def atoms(self, value):
620 self._atoms = copy.deepcopy(value)
621 self._atoms.unrestricted = False
622 self._atoms.occ.fill()
623 if not self._atoms.is_built:
624 self._atoms = self._atoms.build()
627class USCF(SCF):
628 """SCF class for spin-polarized systems.
630 Inherited from :class:`eminus.scf.SCF`.
632 In difference to the SCF class, this class will not build the original Atoms object, only the
633 one attributed to the class. Customized fillings could be overwritten when using this class.
634 """
636 @property
637 def atoms(self):
638 """Atoms object."""
639 return self._atoms
641 @atoms.setter
642 def atoms(self, value):
643 self._atoms = copy.deepcopy(value)
644 self._atoms.unrestricted = True
645 self._atoms.occ.fill()
646 if not self._atoms.is_built:
647 self._atoms = self._atoms.build()