Coverage for eminus/scf.py: 97.99%
348 statements
« prev ^ index » next coverage.py v7.7.0, created at 2025-03-18 09: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
9import numpy as np
11from . import config
12from .band_minimizer import get_grad_unocc, scf_step_unocc
13from .band_minimizer import IMPLEMENTED as BAND_MINIMIZER
14from .dft import (
15 get_epsilon,
16 get_n_spin,
17 get_n_total,
18 get_phi,
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) or "pseudo". Adding "symm" to the string will use the
54 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 "pseudo" in value:
213 self._guess = "pseudo"
214 else:
215 msg = f"{value} is no valid initial guess."
216 raise ValueError(msg)
217 # Check if a symmetric or unsymmetric guess is selected
218 if "sym" in value and "unsym" not in value:
219 self._symmetric = True
220 else:
221 self._symmetric = False
223 @property
224 def opt(self):
225 """Minimization methods."""
226 return self._opt
228 @opt.setter
229 def opt(self, value):
230 # Set lowercase to all keys
231 value = {k.lower(): v for k, v in value.items()}
232 for opt in value:
233 if opt not in ALL_MINIMIZER:
234 msg = f'No minimizer found for "{opt}".'
235 raise KeyError(msg)
236 self._opt = value
238 @property
239 def verbose(self):
240 """Verbosity level."""
241 return self._verbose
243 @verbose.setter
244 def verbose(self, value):
245 # If no verbosity is given use the one from the Atoms object
246 if value is None:
247 value = self.atoms.verbose
248 self._verbose = get_level(value)
249 self._log.verbose = self._verbose
251 # ### Read-only properties ###
253 @property
254 def kpts(self):
255 """Pass-through to the KPoints object of the Atoms object."""
256 return self.atoms.kpts
258 @property
259 def pot_params_defaults(self):
260 """Get the default potential parameters."""
261 return get_pot_defaults(self.pot)
263 @property
264 def psp(self):
265 """Pseudopotential path."""
266 return self._psp
268 @property
269 def symmetric(self):
270 """Determines if the initial guess is the same for both spin channels."""
271 return self._symmetric
273 @property
274 def xc_type(self):
275 """Determines the exchange-correlation family."""
276 return self._xc_type
278 @property
279 def xc_params_defaults(self):
280 """Get the default exchange-correlation functional parameters."""
281 return get_xc_defaults(self.xc)
283 # ### Class methods ###
285 def run(self, **kwargs):
286 """Run the self-consistent field (SCF) calculation.
288 Keyword Args:
289 **kwargs: Pass-through keyword arguments.
291 Returns:
292 Total energy.
293 """
294 # Print some information about the calculation
295 if self._log.level <= logging.DEBUG:
296 info()
297 if ":" in "".join(self._xc) and config.use_pylibxc:
298 self._log.info(
299 "This calculation employs Libxc to evaluate density functionals. When using Libxc,"
300 "\nplease cite SoftwareX 7, 1 (2018) (doi:10.1016/j.softx.2017.11.002)."
301 )
302 self._log.debug(
303 f"\n--- Atoms information ---\n{self.atoms}\n"
304 f"\n--- Cell information ---\nCell vectors:\n{self.atoms.a} a0\n"
305 f"Sampling per axis: {self.atoms.s}\n"
306 f"Cut-off energy: {self.atoms.ecut} Eh\n"
307 f"\n--- State information ---\n{self.atoms.occ}\n"
308 f"\n--- Calculation information ---\n{self}\n\n--- SCF data ---"
309 )
311 # Calculate Ewald energy that only depends on the system geometry
312 self.energies.Eewald = get_Eewald(self.atoms)
313 # Build the initial wave function if there is no W to start from
314 if self.W is None:
315 if "random" in self.guess:
316 self.W = guess_random(self, symmetric=self.symmetric)
317 elif "pseudo" in self.guess:
318 self.W = guess_pseudo(self, symmetric=self.symmetric)
320 # Start the minimization procedures
321 self.clear()
322 Etots = []
323 for imin in self.opt:
324 # Call the minimizer
325 self._log.info(f"Start {ALL_MINIMIZER[imin].__name__}...")
326 start = time.perf_counter()
327 Elist = ALL_MINIMIZER[imin](self, self.opt[imin], **kwargs)
328 end = time.perf_counter()
329 # Save the minimizer results
330 self._opt_log[imin] = {}
331 self._opt_log[imin]["iter"] = len(Elist)
332 self._opt_log[imin]["time"] = end - start
333 Etots += Elist
334 # Do not start other minimizations if one converged
335 if self.is_converged:
336 break
337 if self.is_converged:
338 self._log.info(f"SCF converged after {len(Etots)} iterations.")
339 else:
340 self._log.warning("SCF not converged!")
342 # Calculate SIC energy if needed
343 if self.sic:
344 self.energies.Esic = get_Esic(self, self.Y)
345 # Calculate dispersion correction energy if needed
346 if isinstance(self.disp, dict):
347 self.energies.Edisp = get_Edisp(self, **self.disp)
348 elif self.disp:
349 self.energies.Edisp = get_Edisp(self)
351 # Print minimizer timings
352 self._log.debug("\n--- SCF results ---")
353 t_tot = 0
354 for imin in self._opt_log:
355 N = self._opt_log[imin]["iter"]
356 t = self._opt_log[imin]["time"]
357 t_tot += t
358 self._log.debug(
359 f"Minimizer: {imin}"
360 f"\nIterations: {N}"
361 f"\nTime: {t:.5f} s"
362 f"\nTime/Iteration: {t / N:.5f} s"
363 )
364 self._log.info(f"Total SCF time: {t_tot:.5f} s")
365 # Print the S^2 expectation value for unrestricted calculations
366 if self.atoms.unrestricted:
367 self._log.info(f"<S^2> = {get_spin_squared(self):.6e}")
368 # Print energy data
369 if self._log.level <= logging.DEBUG:
370 self._log.debug(
371 "\n--- Energy data ---\n"
372 f"Eigenenergies:\n{get_epsilon(self, self.W)}\n\n{self.energies}"
373 )
374 else:
375 self._log.info(f"Etot = {self.energies.Etot:.9f} Eh")
376 return self.energies.Etot
378 kernel = run
380 def converge_bands(self, **kwargs):
381 """Converge occupied bands after an SCF calculation.
383 Keyword Args:
384 **kwargs: Pass-through keyword arguments.
385 """
386 if not self.is_converged:
387 self._log.warning("The previous calculation has not been converged.")
389 # If new k-points have been set rebuild the atoms object and the potential
390 if not self.atoms.kpts.is_built or (
391 self.W is not None and len(self.W) != self.atoms.kpts.Nk
392 ):
393 self.atoms.build()
394 self.pot = self.pot
395 self.is_converged = False
397 # Build the initial wave function if there is no W to start from
398 if self.W is None or len(self.W) != self.atoms.kpts.Nk:
399 if "random" in self.guess:
400 self.W = guess_random(self, symmetric=self.symmetric)
401 elif "pseudo" in self.guess:
402 self.W = guess_pseudo(self, symmetric=self.symmetric)
404 self._log.info("Minimize occupied band energies...")
405 # Start the minimization procedures
406 Etots = []
407 for imin in self.opt:
408 # Call the minimizer
409 self._log.info(f"Start {BAND_MINIMIZER[imin].__name__}...")
410 start = time.perf_counter()
411 Elist, self.W = BAND_MINIMIZER[imin](self, self.W, self.opt[imin], **kwargs)
412 end = time.perf_counter()
413 # Save the minimizer results
414 self._opt_log[imin] = {}
415 self._opt_log[imin]["iter"] = len(Elist)
416 self._opt_log[imin]["time"] = end - start
417 Etots += Elist
418 # Do not start other minimizations if one converged
419 if self.is_converged:
420 break
421 if self.is_converged:
422 self._log.info(f"Band minimization converged after {len(Etots)} iterations.")
423 else:
424 self._log.warning("Band minimization not converged!")
426 # Print minimizer timings
427 self._log.debug("\n--- Band minimization results ---")
428 t_tot = 0
429 for imin in self._opt_log:
430 N = self._opt_log[imin]["iter"]
431 t = self._opt_log[imin]["time"]
432 t_tot += t
433 self._log.debug(
434 f"Minimizer: {imin}"
435 f"\nIterations: {N}"
436 f"\nTime: {t:.5f} s"
437 f"\nTime/Iteration: {t / N:.5f} s"
438 )
439 self._log.info(f"Total band minimization time: {t_tot:.5f} s")
441 # Converge empty bands automatically if needed
442 if self.atoms.occ.Nempty > 0:
443 self.converge_empty_bands(**kwargs)
444 return self
446 def converge_empty_bands(self, Nempty=None, **kwargs):
447 """Converge unoccupied bands after converging occ. bands.
449 Keyword Args:
450 Nempty: Number of empty states.
451 **kwargs: Pass-through keyword arguments.
452 """
453 if not self.is_converged:
454 self._log.warning("The previous calculation has not been converged.")
455 self.is_converged = False
457 if Nempty is None:
458 Nempty = self.atoms.occ.Nempty
460 # Build the initial wave functions
461 if self.Z is None:
462 if "random" in self.guess:
463 self.Z = guess_random(self, Nempty, symmetric=self.symmetric)
464 elif "pseudo" in self.guess:
465 self.Z = guess_pseudo(self, Nempty, symmetric=self.symmetric)
467 self._log.info("Minimize unoccupied band energies...")
468 # Start the minimization procedures
469 Etots = []
470 for imin in self.opt:
471 # Call the minimizer
472 self._log.info(f"Start {BAND_MINIMIZER[imin].__name__}...")
473 start = time.perf_counter()
474 Elist, self.Z = BAND_MINIMIZER[imin](
475 self, self.Z, self.opt[imin], cost=scf_step_unocc, grad=get_grad_unocc, **kwargs
476 )
477 end = time.perf_counter()
478 # Save the minimizer results
479 self._opt_log[imin] = {}
480 self._opt_log[imin]["iter"] = len(Elist)
481 self._opt_log[imin]["time"] = end - start
482 Etots += Elist
483 # Do not start other minimizations if one converged
484 if self.is_converged:
485 break
486 if self.is_converged:
487 self._log.info(f"Band minimization converged after {len(Etots)} iterations.")
488 else:
489 self._log.warning("Band minimization not converged!")
491 # Print minimizer timings
492 self._log.debug("\n--- Band minimization results ---")
493 t_tot = 0
494 for imin in self._opt_log:
495 N = self._opt_log[imin]["iter"]
496 t = self._opt_log[imin]["time"]
497 t_tot += t
498 self._log.debug(
499 f"Minimizer: {imin}"
500 f"\nIterations: {N}"
501 f"\nTime: {t:.5f} s"
502 f"\nTime/Iteration: {t / N:.5f} s"
503 )
504 self._log.info(f"Total band minimization time: {t_tot:.5f} s")
505 return self
507 def recenter(self, center=None):
508 """Recenter the system inside the cell.
510 Keyword Args:
511 center: Point to center the system around.
512 """
513 atoms = self.atoms
514 # Get the COM before centering the atoms
515 com = center_of_mass(atoms.pos)
516 # Run the recenter method of the atoms object
517 self.atoms.recenter(center=center)
518 if center is None:
519 dr = com - np.sum(atoms.a, axis=0) / 2
520 else:
521 center = np.asarray(center)
522 dr = com - center
524 # Shift orbitals and density
525 if self.W is not None:
526 self.W = atoms.T(self.W, dr=-dr)
527 # Transform the density to the reciprocal space, shift, and transform back
528 n = None
529 if self.n is not None:
530 Jn = atoms.J(self.n)
531 TJn = atoms.T(Jn, dr=-dr)
532 n = np.real(atoms.I(TJn))
534 # Recalculate the potential since it depends on the structure factor
535 self.pot = self.pot
536 # Clear intermediate results to make sure no one uses the unshifted results
537 self.clear()
538 # Set the shifted density after calling the clearing function
539 self.n = n
540 return self
542 def clear(self):
543 """Initialize or clear intermediate results."""
544 self.Y = None # Orthogonal wave functions
545 self.Z = None # Unconstrained wave functions of unoccupied states
546 self.D = None # Orthogonal wave functions of unoccupied states
547 self.n = None #: Electronic density
548 self.n_spin = None # Electronic densities per spin
549 self.dn_spin = None # Gradient of electronic densities per spin
550 self.tau = None # Kinetic energy densities per spin
551 self.phi = None # Hartree field
552 self.exc = None # Exchange-correlation energy density
553 self.vxc = None # Exchange-correlation potential
554 self.vsigma = None # n times d exc/d |dn|^2
555 self.vtau = None # d exc/d tau
556 self._precomputed = {} # Dictionary of pre-computed values not to be saved
557 self._opt_log = {} # Log of the optimization procedure
558 return self
560 @staticmethod
561 def callback(scf, step):
562 """Callback function that will get called every SCF iteration.
564 This is just an empty function users can overwrite with their implementation.
565 Remember that the "auto" minimization can call the callback function twice when the pccg
566 step is not preferred.
568 Args:
569 scf: SCF object.
570 step: Optimization step.
571 """
573 def _precompute(self):
574 """Precompute fields stored in the SCF object."""
575 atoms = self.atoms
576 self.Y = orth(atoms, self.W)
577 self.n_spin = get_n_spin(atoms, self.Y)
578 self.n = get_n_total(atoms, self.Y, self.n_spin)
579 if "gga" in self.xc_type:
580 self.dn_spin = get_grad_field(atoms, self.n_spin)
581 if self.xc_type == "meta-gga":
582 self.tau = get_tau(atoms, self.Y)
583 self.phi = get_phi(atoms, self.n)
584 self.exc, self.vxc, self.vsigma, self.vtau = get_xc(
585 self.xc, self.n_spin, atoms.occ.Nspin, self.dn_spin, self.tau, self.xc_params
586 )
587 self._precomputed = {
588 "dn_spin": self.dn_spin,
589 "phi": self.phi,
590 "vxc": self.vxc,
591 "vsigma": self.vsigma,
592 "vtau": self.vtau,
593 }
594 return self
596 def __repr__(self):
597 """Print the most important parameters stored in the SCF object."""
598 # Use chr(10) to create a linebreak since backslashes are not allowed in f-strings
599 return (
600 f"XC functionals: {self.xc}\n"
601 f"Potential: {self.pot}\n"
602 f"{f'GTH files: {self.psp}' + chr(10) if self.pot == 'gth' else ''}"
603 f"Starting guess: {self.guess}\n"
604 f"Symmetric guess: {self.symmetric}\n"
605 f"Energy convergence tolerance: {self.etol} Eh\n"
606 f"Gradient convergence tolerance: {self.gradtol}\n"
607 f"Non-local potential: {self.gth.NbetaNL > 0 if self.pot == 'gth' else 'false'}\n"
608 f"Smearing: {self.atoms.occ.smearing > 0}\n"
609 f"Smearing update cycle: {self.smear_update}"
610 )
613class RSCF(SCF):
614 """SCF class for spin-paired systems.
616 Inherited from :class:`eminus.scf.SCF`.
618 In difference to the SCF class, this class will not build the original Atoms object, only the
619 one attributed to the class. Customized fillings could be overwritten when using this class.
620 """
622 @property
623 def atoms(self):
624 """Atoms object."""
625 return self._atoms
627 @atoms.setter
628 def atoms(self, value):
629 self._atoms = copy.deepcopy(value)
630 self._atoms.unrestricted = False
631 self._atoms.occ.fill()
632 if not self._atoms.is_built:
633 self._atoms = self._atoms.build()
636class USCF(SCF):
637 """SCF class for spin-polarized systems.
639 Inherited from :class:`eminus.scf.SCF`.
641 In difference to the SCF class, this class will not build the original Atoms object, only the
642 one attributed to the class. Customized fillings could be overwritten when using this class.
643 """
645 @property
646 def atoms(self):
647 """Atoms object."""
648 return self._atoms
650 @atoms.setter
651 def atoms(self, value):
652 self._atoms = copy.deepcopy(value)
653 self._atoms.unrestricted = True
654 self._atoms.occ.fill()
655 if not self._atoms.is_built:
656 self._atoms = self._atoms.build()