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