Coverage for eminus/scf.py: 97.27%
330 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-08 12:59 +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 IMPLEMENTED as ALL_POTENTIALS
29from .potentials import init_pot
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 = pot #: Used potential.
99 self.guess = guess #: Initial wave functions guess.
100 self.etol = etol #: Total energy convergence tolerance.
101 self.gradtol = gradtol #: Gradient norm convergence tolerance.
102 self.sic = sic #: Enables the SIC energy calculation.
103 self.disp = disp #: Enables the dispersion correction calculation.
104 self.opt = opt #: Minimization methods.
105 self.smear_update = 2 #: Steps after the smeared occupations are recalculated.
107 # Initialize other attributes
108 self.energies = Energy() #: Energy object holding energy contributions.
109 self.is_converged = False #: Determines the SCF object convergence.
110 self.W = None #: Unconstrained wave functions.
111 self.xc_params = {} #: Exchange-correlation functional parameters.
112 self.clear()
114 # ### Class properties ###
116 @property
117 def atoms(self):
118 """Atoms object."""
119 return self._atoms
121 @atoms.setter
122 def atoms(self, value):
123 # Build the Atoms object if necessary and make a copy
124 # This way the Atoms objects inside and outside the class are independent but both are build
125 if not value.is_built:
126 self._atoms = copy.deepcopy(value.build())
127 else:
128 self._atoms = copy.deepcopy(value)
130 @property
131 def xc(self):
132 """Exchange-correlation functional."""
133 return self._xc
135 @xc.setter
136 def xc(self, value):
137 self._xc = parse_functionals(value.lower())
138 # Determine the type of functional combinations
139 self._xc_type = parse_xc_type(self._xc)
140 if "mock_xc" in self._xc and "_xc_" not in "".join(self._xc).lower():
141 self._log.warning("Usage of mock functional detected.")
143 @property
144 def xc_params(self):
145 """Exchange-correlation functional parameters."""
146 return self._xc_params
148 @xc_params.setter
149 def xc_params(self, value):
150 # Check if some parameters are unused in the functionals
151 # This also ensures that we print a warning in case of overlapping parameters in x and c
152 if value != {} and value is not None:
153 not_used = value.keys() - self.xc_params_defaults.keys()
154 if len(not_used) > 0:
155 self._log.warning(f"Some xc_params are unused, namely: {', '.join(not_used)}.")
156 self._xc_params = value
158 @property
159 def pot(self):
160 """Used potential."""
161 return self._pot
163 @pot.setter
164 def pot(self, value):
165 if value.lower() in ALL_POTENTIALS:
166 self._pot = value.lower()
167 # Only set the pseudopotential type for GTH pseudopotentials
168 if self._pot == "gth":
169 if "gga" in self._xc_type:
170 self._psp = "pbe"
171 else:
172 self._psp = "pade"
173 # If pot is no supported potential treat it as a path to a directory containing GTH files
174 else:
175 self._log.info(f'Use the path "{value}" to search for GTH pseudopotential files.')
176 self._psp = value
177 self._pot = "gth"
178 # Build the potential
179 if self._pot == "gth":
180 self.gth = GTH(self)
181 self.Vloc = init_pot(self)
183 @property
184 def guess(self):
185 """Initial wave functions guess."""
186 return self._guess
188 @guess.setter
189 def guess(self, value):
190 # Set the guess method
191 value = value.lower()
192 if "rand" in value:
193 self._guess = "random"
194 elif "pseudo" in value:
195 self._guess = "pseudo"
196 else:
197 msg = f"{value} is no valid initial guess."
198 raise ValueError(msg)
199 # Check if a symmetric or unsymmetric guess is selected
200 if "sym" in value and "unsym" not in value:
201 self._symmetric = True
202 else:
203 self._symmetric = False
205 @property
206 def opt(self):
207 """Minimization methods."""
208 return self._opt
210 @opt.setter
211 def opt(self, value):
212 # Set lowercase to all keys
213 value = {k.lower(): v for k, v in value.items()}
214 for opt in value:
215 if opt not in ALL_MINIMIZER:
216 msg = f'No minimizer found for "{opt}".'
217 raise KeyError(msg)
218 self._opt = value
220 @property
221 def verbose(self):
222 """Verbosity level."""
223 return self._verbose
225 @verbose.setter
226 def verbose(self, value):
227 # If no verbosity is given use the one from the Atoms object
228 if value is None:
229 value = self.atoms.verbose
230 self._verbose = get_level(value)
231 self._log.verbose = self._verbose
233 # ### Read-only properties ###
235 @property
236 def kpts(self):
237 """Pass-through to the KPoints object of the Atoms object."""
238 return self.atoms.kpts
240 @property
241 def psp(self):
242 """Pseudopotential path."""
243 return self._psp
245 @property
246 def symmetric(self):
247 """Determines if the initial guess is the same for both spin channels."""
248 return self._symmetric
250 @property
251 def xc_type(self):
252 """Determines the exchange-correlation family."""
253 return self._xc_type
255 @property
256 def xc_params_defaults(self):
257 """Get the default exchange-correlation functional parameters."""
258 return get_xc_defaults(self.xc)
260 # ### Class methods ###
262 def run(self, **kwargs):
263 """Run the self-consistent field (SCF) calculation.
265 Keyword Args:
266 **kwargs: Pass-through keyword arguments.
268 Returns:
269 Total energy.
270 """
271 # Print some information about the calculation
272 if self._log.level <= logging.DEBUG:
273 info()
274 if ":" in "".join(self._xc) and config.use_pylibxc:
275 self._log.info(
276 "This calculation employs Libxc to evaluate density functionals. When using Libxc,"
277 "\nplease cite SoftwareX 7, 1 (2018) (doi:10.1016/j.softx.2017.11.002)."
278 )
279 self._log.debug(
280 f"\n--- Atoms information ---\n{self.atoms}\n"
281 f"\n--- Cell information ---\nCell vectors:\n{self.atoms.a} a0\n"
282 f"Sampling per axis: {self.atoms.s}\n"
283 f"Cut-off energy: {self.atoms.ecut} Eh\n"
284 f"\n--- State information ---\n{self.atoms.occ}\n"
285 f"\n--- Calculation information ---\n{self}\n\n--- SCF data ---"
286 )
288 # Calculate Ewald energy that only depends on the system geometry
289 self.energies.Eewald = get_Eewald(self.atoms)
290 # Build the initial wave function if there is no W to start from
291 if self.W is None:
292 if "random" in self.guess:
293 self.W = guess_random(self, symmetric=self.symmetric)
294 elif "pseudo" in self.guess:
295 self.W = guess_pseudo(self, symmetric=self.symmetric)
297 # Start the minimization procedures
298 self.clear()
299 Etots = []
300 for imin in self.opt:
301 # Call the minimizer
302 self._log.info(f"Start {ALL_MINIMIZER[imin].__name__}...")
303 start = time.perf_counter()
304 Elist = ALL_MINIMIZER[imin](self, self.opt[imin], **kwargs)
305 end = time.perf_counter()
306 # Save the minimizer results
307 self._opt_log[imin] = {}
308 self._opt_log[imin]["iter"] = len(Elist)
309 self._opt_log[imin]["time"] = end - start
310 Etots += Elist
311 # Do not start other minimizations if one converged
312 if self.is_converged:
313 break
314 if self.is_converged:
315 self._log.info(f"SCF converged after {len(Etots)} iterations.")
316 else:
317 self._log.warning("SCF not converged!")
319 # Calculate SIC energy if desired
320 if self.sic:
321 self.energies.Esic = get_Esic(self, self.Y)
322 # Calculate dispersion correction energy if desired
323 if isinstance(self.disp, dict):
324 self.energies.Edisp = get_Edisp(self, **self.disp)
325 elif self.disp:
326 self.energies.Edisp = get_Edisp(self)
328 # Print minimizer timings
329 self._log.debug("\n--- SCF results ---")
330 t_tot = 0
331 for imin in self._opt_log:
332 N = self._opt_log[imin]["iter"]
333 t = self._opt_log[imin]["time"]
334 t_tot += t
335 self._log.debug(
336 f"Minimizer: {imin}"
337 f"\nIterations: {N}"
338 f"\nTime: {t:.5f} s"
339 f"\nTime/Iteration: {t / N:.5f} s"
340 )
341 self._log.info(f"Total SCF time: {t_tot:.5f} s")
342 # Print the S^2 expectation value for unrestricted calculations
343 if self.atoms.unrestricted:
344 self._log.info(f"<S^2> = {get_spin_squared(self):.6e}")
345 # Print energy data
346 if self._log.level <= logging.DEBUG:
347 self._log.debug(
348 "\n--- Energy data ---\n"
349 f"Eigenenergies:\n{get_epsilon(self, self.W)}\n\n{self.energies}"
350 )
351 else:
352 self._log.info(f"Etot = {self.energies.Etot:.9f} Eh")
353 return self.energies.Etot
355 kernel = run
357 def converge_bands(self, **kwargs):
358 """Converge occupied bands after an SCF calculation.
360 Keyword Args:
361 **kwargs: Pass-through keyword arguments.
362 """
363 if not self.is_converged:
364 self._log.warning("The previous calculation has not been converged.")
366 # If new k-points have been set rebuild the atoms object and the potential
367 if not self.atoms.kpts.is_built or (
368 self.W is not None and len(self.W) != self.atoms.kpts.Nk
369 ):
370 self.atoms.build()
371 self.pot = self.pot
372 self.is_converged = False
374 # Build the initial wave function if there is no W to start from
375 if self.W is None or len(self.W) != self.atoms.kpts.Nk:
376 if "random" in self.guess:
377 self.W = guess_random(self, symmetric=self.symmetric)
378 elif "pseudo" in self.guess:
379 self.W = guess_pseudo(self, symmetric=self.symmetric)
381 self._log.info("Minimize occupied band energies...")
382 # Start the minimization procedures
383 Etots = []
384 for imin in self.opt:
385 # Call the minimizer
386 self._log.info(f"Start {BAND_MINIMIZER[imin].__name__}...")
387 start = time.perf_counter()
388 Elist, self.W = BAND_MINIMIZER[imin](self, self.W, self.opt[imin], **kwargs)
389 end = time.perf_counter()
390 # Save the minimizer results
391 self._opt_log[imin] = {}
392 self._opt_log[imin]["iter"] = len(Elist)
393 self._opt_log[imin]["time"] = end - start
394 Etots += Elist
395 # Do not start other minimizations if one converged
396 if self.is_converged:
397 break
398 if self.is_converged:
399 self._log.info(f"Band minimization converged after {len(Etots)} iterations.")
400 else:
401 self._log.warning("Band minimization not converged!")
403 # Print minimizer timings
404 self._log.debug("\n--- Band minimization results ---")
405 t_tot = 0
406 for imin in self._opt_log:
407 N = self._opt_log[imin]["iter"]
408 t = self._opt_log[imin]["time"]
409 t_tot += t
410 self._log.debug(
411 f"Minimizer: {imin}"
412 f"\nIterations: {N}"
413 f"\nTime: {t:.5f} s"
414 f"\nTime/Iteration: {t / N:.5f} s"
415 )
416 self._log.info(f"Total band minimization time: {t_tot:.5f} s")
418 # Converge empty bands automatically if desired
419 if self.atoms.occ.Nempty > 0:
420 self.converge_empty_bands(**kwargs)
421 return self
423 def converge_empty_bands(self, Nempty=None, **kwargs):
424 """Converge unoccupied bands after converging occ. bands.
426 Keyword Args:
427 Nempty: Number of empty states.
428 **kwargs: Pass-through keyword arguments.
429 """
430 if not self.is_converged:
431 self._log.warning("The previous calculation has not been converged.")
432 self.is_converged = False
434 if Nempty is None:
435 Nempty = self.atoms.occ.Nempty
437 # Build the initial wave functions
438 if self.Z is None:
439 if "random" in self.guess:
440 self.Z = guess_random(self, Nempty, symmetric=self.symmetric)
441 elif "pseudo" in self.guess:
442 self.Z = guess_pseudo(self, Nempty, symmetric=self.symmetric)
444 self._log.info("Minimize unoccupied band energies...")
445 # Start the minimization procedures
446 Etots = []
447 for imin in self.opt:
448 # Call the minimizer
449 self._log.info(f"Start {BAND_MINIMIZER[imin].__name__}...")
450 start = time.perf_counter()
451 Elist, self.Z = BAND_MINIMIZER[imin](
452 self, self.Z, self.opt[imin], cost=scf_step_unocc, grad=get_grad_unocc, **kwargs
453 )
454 end = time.perf_counter()
455 # Save the minimizer results
456 self._opt_log[imin] = {}
457 self._opt_log[imin]["iter"] = len(Elist)
458 self._opt_log[imin]["time"] = end - start
459 Etots += Elist
460 # Do not start other minimizations if one converged
461 if self.is_converged:
462 break
463 if self.is_converged:
464 self._log.info(f"Band minimization converged after {len(Etots)} iterations.")
465 else:
466 self._log.warning("Band minimization not converged!")
468 # Print minimizer timings
469 self._log.debug("\n--- Band minimization results ---")
470 t_tot = 0
471 for imin in self._opt_log:
472 N = self._opt_log[imin]["iter"]
473 t = self._opt_log[imin]["time"]
474 t_tot += t
475 self._log.debug(
476 f"Minimizer: {imin}"
477 f"\nIterations: {N}"
478 f"\nTime: {t:.5f} s"
479 f"\nTime/Iteration: {t / N:.5f} s"
480 )
481 self._log.info(f"Total band minimization time: {t_tot:.5f} s")
482 return self
484 def recenter(self, center=None):
485 """Recenter the system inside the cell.
487 Keyword Args:
488 center: Point to center the system around.
489 """
490 atoms = self.atoms
491 # Get the COM before centering the atoms
492 com = center_of_mass(atoms.pos)
493 # Run the recenter method of the atoms object
494 self.atoms.recenter(center=center)
495 if center is None:
496 dr = com - np.sum(atoms.a, axis=0) / 2
497 else:
498 center = np.asarray(center)
499 dr = com - center
501 # Shift orbitals and density
502 if self.W is not None:
503 self.W = atoms.T(self.W, dr=-dr)
504 # Transform the density to the reciprocal space, shift, and transform back
505 n = None
506 if self.n is not None:
507 Jn = atoms.J(self.n)
508 TJn = atoms.T(Jn, dr=-dr)
509 n = np.real(atoms.I(TJn))
511 # Recalculate the potential since it depends on the structure factor
512 self.pot = self.pot
513 # Clear intermediate results to make sure no one uses the unshifted results
514 self.clear()
515 # Set the shifted density after calling the clearing function
516 self.n = n
517 return self
519 def clear(self):
520 """Initialize or clear intermediate results."""
521 self.Y = None # Orthogonal wave functions
522 self.Z = None # Unconstrained wave functions of unoccupied states
523 self.D = None # Orthogonal wave functions of unoccupied states
524 self.n = None #: Electronic density
525 self.n_spin = None # Electronic densities per spin
526 self.dn_spin = None # Gradient of electronic densities per spin
527 self.tau = None # Kinetic energy densities per spin
528 self.phi = None # Hartree field
529 self.exc = None # Exchange-correlation energy density
530 self.vxc = None # Exchange-correlation potential
531 self.vsigma = None # n times d exc/d |dn|^2
532 self.vtau = None # d exc/d tau
533 self._precomputed = {} # Dictionary of pre-computed values not to be saved
534 self._opt_log = {} # Log of the optimization procedure
535 return self
537 def _precompute(self):
538 """Precompute fields stored in the SCF object."""
539 atoms = self.atoms
540 self.Y = orth(atoms, self.W)
541 self.n_spin = get_n_spin(atoms, self.Y)
542 self.n = get_n_total(atoms, self.Y, self.n_spin)
543 if "gga" in self.xc_type:
544 self.dn_spin = get_grad_field(atoms, self.n_spin)
545 if self.xc_type == "meta-gga":
546 self.tau = get_tau(atoms, self.Y)
547 self.phi = get_phi(atoms, self.n)
548 self.exc, self.vxc, self.vsigma, self.vtau = get_xc(
549 self.xc, self.n_spin, atoms.occ.Nspin, self.dn_spin, self.tau, self.xc_params
550 )
551 self._precomputed = {
552 "dn_spin": self.dn_spin,
553 "phi": self.phi,
554 "vxc": self.vxc,
555 "vsigma": self.vsigma,
556 "vtau": self.vtau,
557 }
558 return self
560 def __repr__(self):
561 """Print the most important parameters stored in the SCF object."""
562 # Use chr(10) to create a linebreak since backslashes are not allowed in f-strings
563 return (
564 f"XC functionals: {self.xc}\n"
565 f"Potential: {self.pot}\n"
566 f"{f'GTH files: {self.psp}' + chr(10) if self.pot == 'gth' else ''}"
567 f"Starting guess: {self.guess}\n"
568 f"Symmetric guess: {self.symmetric}\n"
569 f"Energy convergence tolerance: {self.etol} Eh\n"
570 f"Gradient convergence tolerance: {self.gradtol}\n"
571 f"Non-local potential: {self.gth.NbetaNL > 0 if self.pot == 'gth' else 'false'}\n"
572 f"Smearing: {self.atoms.occ.smearing > 0}\n"
573 f"Smearing update cycle: {self.smear_update}"
574 )
577class RSCF(SCF):
578 """SCF class for spin-paired systems.
580 Inherited from :class:`eminus.scf.SCF`.
582 In difference to the SCF class, this class will not build the original Atoms object, only the
583 one attributed to the class. Customized fillings could be overwritten when using this class.
584 """
586 @property
587 def atoms(self):
588 """Atoms object."""
589 return self._atoms
591 @atoms.setter
592 def atoms(self, value):
593 self._atoms = copy.deepcopy(value)
594 self._atoms.unrestricted = False
595 self._atoms.occ.fill()
596 if not self._atoms.is_built:
597 self._atoms = self._atoms.build()
600class USCF(SCF):
601 """SCF class for spin-polarized systems.
603 Inherited from :class:`eminus.scf.SCF`.
605 In difference to the SCF class, this class will not build the original Atoms object, only the
606 one attributed to the class. Customized fillings could be overwritten when using this class.
607 """
609 @property
610 def atoms(self):
611 """Atoms object."""
612 return self._atoms
614 @atoms.setter
615 def atoms(self, value):
616 self._atoms = copy.deepcopy(value)
617 self._atoms.unrestricted = True
618 self._atoms.occ.fill()
619 if not self._atoms.is_built:
620 self._atoms = self._atoms.build()