Coverage for eminus/scf.py: 97.99%

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 get_epsilon, get_n_spin, get_n_total, get_phi, guess_pseudo, guess_random, orth

15from .energies import Energy, get_Edisp, get_Eewald, get_Esic

16from .gga import get_grad_field, get_tau

17from .gth import GTH

18from .logger import create_logger, get_level

19from .minimizer import IMPLEMENTED as ALL_MINIMIZER

20from .potentials import get_pot_defaults, init_pot

21from .potentials import IMPLEMENTED as ALL_POTENTIALS

22from .tools import center_of_mass, get_spin_squared

23from .utils import BaseObject

24from .version import info

25from .xc import get_xc, get_xc_defaults, parse_functionals, parse_xc_type

28class SCF(BaseObject):

29 """Perform direct minimizations.

31 Args:

32 atoms: Atoms object.

34 Keyword Args:

35 xc: Comma-separated exchange-correlation functional description.

37 Adding "libxc:" before a functional will use the Libxc interface for that functional,

38 e.g., with :code:`libxc:mgga_x_scan,libxc:mgga_c_scan`.

39 pot: Type of potential.

41 Can be one of "GTH" (default), "Coulomb", "Harmonic", or "Ge". Alternatively, a path to

42 a directory containing custom GTH pseudopotential files can be given.

43 guess: Initial guess for the wave functions.

45 Can be one of "random" (default) or "pseudo". Adding "symm" to the string will use the

46 same coefficients for both spin channels, e.g., :code:`symm-rand`.

47 etol: Convergence tolerance of the total energy.

48 gradtol: Convergence tolerance of the gradient norm.

50 This tolerance will only be used in conjugate gradient methods.

51 sic: Calculate the Kohn-Sham Perdew-Zunger SIC energy at the end of the SCF.

52 disp: Calculate a dispersion correction.

54 A dictionary can be used to pass arguments to the respective

55 function, e.g., with :code:`{"version": "d3zero", "atm": False, "xc": "scan"}`.

56 opt: Dictionary to customize the minimization methods.

58 The keys can be chosen out of "sd", "lm", "pclm", "cg", "pccg", and "auto". Defaults to

59 :code:`{"auto": 250}`.

60 verbose: Level of output.

62 Can be one of "critical", "error", "warning", "info" (default), or "debug". An integer

63 value can be used as well, where larger numbers mean more output, starting from 0.

64 Defaults to the verbosity level of the Atoms object.

65 """

67 def __init__(

68 self,

69 atoms,

70 xc="lda,vwn",

71 pot="gth",

72 guess="random",

73 etol=1e-7,

74 gradtol=None,

75 sic=False,

76 disp=False,

77 opt=None,

78 verbose=None,

79 ):

80 """Initialize the SCF object."""

81 # Set opt here, better to not use mutable data types in signatures

82 if opt is None:

83 opt = {"auto": 250}

85 # Set the input parameters (the ordering is important)

86 self.atoms = atoms #: Atoms object.

87 self._log = create_logger(self) #: Logger object.

88 self.verbose = verbose #: Verbosity level.

89 self.xc = xc #: Exchange-correlation functional.

90 self.pot_params = {} #: Potential parameters.

91 self.pot = pot #: Used potential.

92 self.guess = guess #: Initial wave functions guess.

93 self.etol = etol #: Total energy convergence tolerance.

94 self.gradtol = gradtol #: Gradient norm convergence tolerance.

95 self.sic = sic #: Enables the SIC energy calculation.

96 self.disp = disp #: Enables the dispersion correction calculation.

97 self.opt = opt #: Minimization methods.

98 self.smear_update = 2 #: Steps after the smeared occupations are recalculated.

100 # Initialize other attributes

101 self.energies = Energy() #: Energy object holding energy contributions.

102 self.is_converged = False #: Determines the SCF object convergence.

103 self.W = None #: Unconstrained wave functions.

104 self.xc_params = {} #: Exchange-correlation functional parameters.

105 self.clear()

106

107 # ### Class properties ###

108

109 @property

110 def atoms(self):

111 """Atoms object."""

112 return self._atoms

113

114 @atoms.setter

115 def atoms(self, value):

116 # Build the Atoms object if necessary and make a copy

117 # This way the Atoms objects inside and outside the class are independent but both are build

118 if not value.is_built:

119 self._atoms = copy.deepcopy(value.build())

120 else:

121 self._atoms = copy.deepcopy(value)

122

123 @property

124 def xc(self):

125 """Exchange-correlation functional."""

126 return self._xc

127

128 @xc.setter

129 def xc(self, value):

130 self._xc = parse_functionals(value.lower())

131 # Determine the type of functional combinations

132 self._xc_type = parse_xc_type(self._xc)

133 if "mock_xc" in self._xc and "_xc_" not in "".join(self._xc).lower():

134 self._log.warning("Usage of mock functional detected.")

135

136 @property

137 def xc_params(self):

138 """Exchange-correlation functional parameters."""

139 return self._xc_params

140

141 @xc_params.setter

142 def xc_params(self, value):

143 # Check if some parameters are unused in the functionals

144 # This also ensures that we print a warning in case of overlapping parameters in x and c

145 if value != {} and value is not None:

146 not_used = value.keys() - self.xc_params_defaults.keys()

147 if len(not_used) > 0:

148 self._log.warning(f"Some xc_params are unused, namely: {', '.join(not_used)}.")

149 self._xc_params = value

150

151 @property

152 def pot(self):

153 """Used potential."""

154 return self._pot

155

156 @pot.setter

157 def pot(self, value):

158 if value.lower() in ALL_POTENTIALS:

159 self._pot = value.lower()

160 # Only set the pseudopotential type for GTH pseudopotentials

161 if self._pot == "gth":

162 if "gga" in self._xc_type:

163 self._psp = "pbe"

164 else:

165 self._psp = "pade"

166 # If pot is no supported potential treat it as a path to a directory containing GTH files

167 else:

168 self._log.info(f'Use the path "{value}" to search for GTH pseudopotential files.')

169 self._psp = value

170 self._pot = "gth"

171 # Build the potential

172 if self._pot == "gth":

173 self.gth = GTH(self)

174 self.Vloc = init_pot(self, self.pot_params)

175

176 @property

177 def pot_params(self):

178 """Potential parameters."""

179 return self._pot_params

180

181 @pot_params.setter

182 def pot_params(self, value):

183 # Check if some parameters are unused in the potential

184 if value != {} and value is not None:

185 not_used = value.keys() - self.pot_params_defaults.keys()

186 if len(not_used) > 0:

187 self._log.warning(f"Some pot_params are unused, namely: {', '.join(not_used)}.")

188 self._pot_params = value

189 # Update the local potential for the new parameters

190 if hasattr(self, "pot"):

191 self.Vloc = init_pot(self, self.pot_params)

192

193 @property

194 def guess(self):

195 """Initial wave functions guess."""

196 return self._guess

197

198 @guess.setter

199 def guess(self, value):

200 # Set the guess method

201 value = value.lower()

202 if "rand" in value:

203 self._guess = "random"

204 elif "pseudo" in value:

205 self._guess = "pseudo"

206 else:

207 msg = f"{value} is no valid initial guess."

208 raise ValueError(msg)

209 # Check if a symmetric or unsymmetric guess is selected

210 if "sym" in value and "unsym" not in value:

211 self._symmetric = True

212 else:

213 self._symmetric = False

214

215 @property

216 def opt(self):

217 """Minimization methods."""

218 return self._opt

219

220 @opt.setter

221 def opt(self, value):

222 # Set lowercase to all keys

223 value = {k.lower(): v for k, v in value.items()}

224 for opt in value:

225 if opt not in ALL_MINIMIZER:

226 msg = f'No minimizer found for "{opt}".'

227 raise KeyError(msg)

228 self._opt = value

229

230 @property

231 def verbose(self):

232 """Verbosity level."""

233 return self._verbose

234

235 @verbose.setter

236 def verbose(self, value):

237 # If no verbosity is given use the one from the Atoms object

238 if value is None:

239 value = self.atoms.verbose

240 self._verbose = get_level(value)

241 self._log.verbose = self._verbose

242

243 # ### Read-only properties ###

244

245 @property

246 def kpts(self):

247 """Pass-through to the KPoints object of the Atoms object."""

248 return self.atoms.kpts

249

250 @property

251 def pot_params_defaults(self):

252 """Get the default potential parameters."""

253 return get_pot_defaults(self.pot)

254

255 @property

256 def psp(self):

257 """Pseudopotential path."""

258 return self._psp

259

260 @property

261 def symmetric(self):

262 """Determines if the initial guess is the same for both spin channels."""

263 return self._symmetric

264

265 @property

266 def xc_type(self):

267 """Determines the exchange-correlation family."""

268 return self._xc_type

269

270 @property

271 def xc_params_defaults(self):

272 """Get the default exchange-correlation functional parameters."""

273 return get_xc_defaults(self.xc)

274

275 # ### Class methods ###

276

277 def run(self, **kwargs):

278 """Run the self-consistent field (SCF) calculation.

279

280 Keyword Args:

281 **kwargs: Pass-through keyword arguments.

282

283 Returns:

284 Total energy.

285 """

286 # Print some information about the calculation

287 if self._log.level <= logging.DEBUG:

288 info()

289 if ":" in "".join(self._xc) and config.use_pylibxc:

290 self._log.info(

291 "This calculation employs Libxc to evaluate density functionals. When using Libxc,"

292 "\nplease cite SoftwareX 7, 1 (2018) (doi:10.1016/j.softx.2017.11.002)."

293 )

294 self._log.debug(

295 f"\n--- Atoms information ---\n{self.atoms}\n"

296 f"\n--- Cell information ---\nCell vectors:\n{self.atoms.a} a0\n"

297 f"Sampling per axis: {self.atoms.s}\n"

298 f"Cut-off energy: {self.atoms.ecut} Eh\n"

299 f"\n--- State information ---\n{self.atoms.occ}\n"

300 f"\n--- Calculation information ---\n{self}\n\n--- SCF data ---"

301 )

302

303 # Calculate Ewald energy that only depends on the system geometry

304 self.energies.Eewald = get_Eewald(self.atoms)

305 # Build the initial wave function if there is no W to start from

306 if self.W is None:

307 if "random" in self.guess:

308 self.W = guess_random(self, symmetric=self.symmetric)

309 elif "pseudo" in self.guess:

310 self.W = guess_pseudo(self, symmetric=self.symmetric)

311

312 # Start the minimization procedures

313 self.clear()

314 Etots = []

315 for imin in self.opt:

316 # Call the minimizer

317 self._log.info(f"Start {ALL_MINIMIZER[imin].__name__}...")

318 start = time.perf_counter()

319 Elist = ALL_MINIMIZER[imin](self, self.opt[imin], **kwargs)

320 end = time.perf_counter()

321 # Save the minimizer results

322 self._opt_log[imin] = {}

323 self._opt_log[imin]["iter"] = len(Elist)

324 self._opt_log[imin]["time"] = end - start

325 Etots += Elist

326 # Do not start other minimizations if one converged

327 if self.is_converged:

328 break

329 if self.is_converged:

330 self._log.info(f"SCF converged after {len(Etots)} iterations.")

331 else:

332 self._log.warning("SCF not converged!")

333

334 # Calculate SIC energy if needed

335 if self.sic:

336 self.energies.Esic = get_Esic(self, self.Y)

337 # Calculate dispersion correction energy if needed

338 if isinstance(self.disp, dict):

339 self.energies.Edisp = get_Edisp(self, **self.disp)

340 elif self.disp:

341 self.energies.Edisp = get_Edisp(self)

342

343 # Print minimizer timings

344 self._log.debug("\n--- SCF results ---")

345 t_tot = 0

346 for imin in self._opt_log:

347 N = self._opt_log[imin]["iter"]

348 t = self._opt_log[imin]["time"]

349 t_tot += t

350 self._log.debug(

351 f"Minimizer: {imin}"

352 f"\nIterations: {N}"

353 f"\nTime: {t:.5f} s"

354 f"\nTime/Iteration: {t / N:.5f} s"

355 )

356 self._log.info(f"Total SCF time: {t_tot:.5f} s")

357 # Print the S^2 expectation value for unrestricted calculations

358 if self.atoms.unrestricted:

359 self._log.info(f"<S^2> = {get_spin_squared(self):.6e}")

360 # Print energy data

361 if self._log.level <= logging.DEBUG:

362 self._log.debug(

363 "\n--- Energy data ---\n"

364 f"Eigenenergies:\n{get_epsilon(self, self.W)}\n\n{self.energies}"

365 )

366 else:

367 self._log.info(f"Etot = {self.energies.Etot:.9f} Eh")

368 return self.energies.Etot

369

370 kernel = run

371

372 def converge_bands(self, **kwargs):

373 """Converge occupied bands after an SCF calculation.

374

375 Keyword Args:

376 **kwargs: Pass-through keyword arguments.

377 """

378 if not self.is_converged:

379 self._log.warning("The previous calculation has not been converged.")

380

381 # If new k-points have been set rebuild the atoms object and the potential

382 if not self.atoms.kpts.is_built or (

383 self.W is not None and len(self.W) != self.atoms.kpts.Nk

384 ):

385 self.atoms.build()

386 self.pot = self.pot

387 self.is_converged = False

388

389 # Build the initial wave function if there is no W to start from

390 if self.W is None or len(self.W) != self.atoms.kpts.Nk:

391 if "random" in self.guess:

392 self.W = guess_random(self, symmetric=self.symmetric)

393 elif "pseudo" in self.guess:

394 self.W = guess_pseudo(self, symmetric=self.symmetric)

395

396 self._log.info("Minimize occupied band energies...")

397 # Start the minimization procedures

398 Etots = []

399 for imin in self.opt:

400 # Call the minimizer

401 self._log.info(f"Start {BAND_MINIMIZER[imin].__name__}...")

402 start = time.perf_counter()

403 Elist, self.W = BAND_MINIMIZER[imin](self, self.W, self.opt[imin], **kwargs)

404 end = time.perf_counter()

405 # Save the minimizer results

406 self._opt_log[imin] = {}

407 self._opt_log[imin]["iter"] = len(Elist)

408 self._opt_log[imin]["time"] = end - start

409 Etots += Elist

410 # Do not start other minimizations if one converged

411 if self.is_converged:

412 break

413 if self.is_converged:

414 self._log.info(f"Band minimization converged after {len(Etots)} iterations.")

415 else:

416 self._log.warning("Band minimization not converged!")

417

418 # Print minimizer timings

419 self._log.debug("\n--- Band minimization results ---")

420 t_tot = 0

421 for imin in self._opt_log:

422 N = self._opt_log[imin]["iter"]

423 t = self._opt_log[imin]["time"]

424 t_tot += t

425 self._log.debug(

426 f"Minimizer: {imin}"

427 f"\nIterations: {N}"

428 f"\nTime: {t:.5f} s"

429 f"\nTime/Iteration: {t / N:.5f} s"

430 )

431 self._log.info(f"Total band minimization time: {t_tot:.5f} s")

432

433 # Converge empty bands automatically if needed

434 if self.atoms.occ.Nempty > 0:

435 self.converge_empty_bands(**kwargs)

436 return self

437

438 def converge_empty_bands(self, Nempty=None, **kwargs):

439 """Converge unoccupied bands after converging occ. bands.

440

441 Keyword Args:

442 Nempty: Number of empty states.

443 **kwargs: Pass-through keyword arguments.

444 """

445 if not self.is_converged:

446 self._log.warning("The previous calculation has not been converged.")

447 self.is_converged = False

448

449 if Nempty is None:

450 Nempty = self.atoms.occ.Nempty

451

452 # Build the initial wave functions

453 if self.Z is None:

454 if "random" in self.guess:

455 self.Z = guess_random(self, Nempty, symmetric=self.symmetric)

456 elif "pseudo" in self.guess:

457 self.Z = guess_pseudo(self, Nempty, symmetric=self.symmetric)

458

459 self._log.info("Minimize unoccupied band energies...")

460 # Start the minimization procedures

461 Etots = []

462 for imin in self.opt:

463 # Call the minimizer

464 self._log.info(f"Start {BAND_MINIMIZER[imin].__name__}...")

465 start = time.perf_counter()

466 Elist, self.Z = BAND_MINIMIZER[imin](

467 self, self.Z, self.opt[imin], cost=scf_step_unocc, grad=get_grad_unocc, **kwargs

468 )

469 end = time.perf_counter()

470 # Save the minimizer results

471 self._opt_log[imin] = {}

472 self._opt_log[imin]["iter"] = len(Elist)

473 self._opt_log[imin]["time"] = end - start

474 Etots += Elist

475 # Do not start other minimizations if one converged

476 if self.is_converged:

477 break

478 if self.is_converged:

479 self._log.info(f"Band minimization converged after {len(Etots)} iterations.")

480 else:

481 self._log.warning("Band minimization not converged!")

482

483 # Print minimizer timings

484 self._log.debug("\n--- Band minimization results ---")

485 t_tot = 0

486 for imin in self._opt_log:

487 N = self._opt_log[imin]["iter"]

488 t = self._opt_log[imin]["time"]

489 t_tot += t

490 self._log.debug(

491 f"Minimizer: {imin}"

492 f"\nIterations: {N}"

493 f"\nTime: {t:.5f} s"

494 f"\nTime/Iteration: {t / N:.5f} s"

495 )

496 self._log.info(f"Total band minimization time: {t_tot:.5f} s")

497 return self

498

499 def recenter(self, center=None):

500 """Recenter the system inside the cell.

501

502 Keyword Args:

503 center: Point to center the system around.

504 """

505 atoms = self.atoms

506 # Get the COM before centering the atoms

507 com = center_of_mass(atoms.pos)

508 # Run the recenter method of the atoms object

509 self.atoms.recenter(center=center)

510 if center is None:

511 dr = com - np.sum(atoms.a, axis=0) / 2

512 else:

513 center = np.asarray(center)

514 dr = com - center

515

516 # Shift orbitals and density

517 if self.W is not None:

518 self.W = atoms.T(self.W, dr=-dr)

519 # Transform the density to the reciprocal space, shift, and transform back

520 n = None

521 if self.n is not None:

522 Jn = atoms.J(self.n)

523 TJn = atoms.T(Jn, dr=-dr)

524 n = np.real(atoms.I(TJn))

525

526 # Recalculate the potential since it depends on the structure factor

527 self.pot = self.pot

528 # Clear intermediate results to make sure no one uses the unshifted results

529 self.clear()

530 # Set the shifted density after calling the clearing function

531 self.n = n

532 return self

533

534 def clear(self):

535 """Initialize or clear intermediate results."""

536 self.Y = None # Orthogonal wave functions

537 self.Z = None # Unconstrained wave functions of unoccupied states

538 self.D = None # Orthogonal wave functions of unoccupied states

539 self.n = None #: Electronic density

540 self.n_spin = None # Electronic densities per spin

541 self.dn_spin = None # Gradient of electronic densities per spin

542 self.tau = None # Kinetic energy densities per spin

543 self.phi = None # Hartree field

544 self.exc = None # Exchange-correlation energy density

545 self.vxc = None # Exchange-correlation potential

546 self.vsigma = None # n times d exc/d |dn|^2

547 self.vtau = None # d exc/d tau

548 self._precomputed = {} # Dictionary of pre-computed values not to be saved

549 self._opt_log = {} # Log of the optimization procedure

550 return self

551

552 @staticmethod

553 def callback(scf, step):

554 """Callback function that will get called every SCF iteration.

555

556 This is just an empty function users can overwrite with their implementation.

557 Remember that the "auto" minimization can call the callback function twice when the pccg

558 step is not preferred.

559

560 Args:

561 scf: SCF object.

562 step: Optimization step.

563 """

564

565 def _precompute(self):

566 """Precompute fields stored in the SCF object."""

567 atoms = self.atoms

568 self.Y = orth(atoms, self.W)

569 self.n_spin = get_n_spin(atoms, self.Y)

570 self.n = get_n_total(atoms, self.Y, self.n_spin)

571 if "gga" in self.xc_type:

572 self.dn_spin = get_grad_field(atoms, self.n_spin)

573 if self.xc_type == "meta-gga":

574 self.tau = get_tau(atoms, self.Y)

575 self.phi = get_phi(atoms, self.n)

576 self.exc, self.vxc, self.vsigma, self.vtau = get_xc(

577 self.xc, self.n_spin, atoms.occ.Nspin, self.dn_spin, self.tau, self.xc_params

578 )

579 self._precomputed = {

580 "dn_spin": self.dn_spin,

581 "phi": self.phi,

582 "vxc": self.vxc,

583 "vsigma": self.vsigma,

584 "vtau": self.vtau,

585 }

586 return self

587

588 def __repr__(self):

589 """Print the most important parameters stored in the SCF object."""

590 # Use chr(10) to create a linebreak since backslashes are not allowed in f-strings

591 return (

592 f"XC functionals: {self.xc}\n"

593 f"Potential: {self.pot}\n"

594 f"{f'GTH files: {self.psp}' + chr(10) if self.pot == 'gth' else ''}"

595 f"Starting guess: {self.guess}\n"

596 f"Symmetric guess: {self.symmetric}\n"

597 f"Energy convergence tolerance: {self.etol} Eh\n"

598 f"Gradient convergence tolerance: {self.gradtol}\n"

599 f"Non-local potential: {self.gth.NbetaNL > 0 if self.pot == 'gth' else 'false'}\n"

600 f"Smearing: {self.atoms.occ.smearing > 0}\n"

601 f"Smearing update cycle: {self.smear_update}"

602 )

603

604

605class RSCF(SCF):

606 """SCF class for spin-paired systems.

607

608 Inherited from :class:`eminus.scf.SCF`.

609

610 In difference to the SCF class, this class will not build the original Atoms object, only the

611 one attributed to the class. Customized fillings could be overwritten when using this class.

612 """

613

614 @property

615 def atoms(self):

616 """Atoms object."""

617 return self._atoms

618

619 @atoms.setter

620 def atoms(self, value):

621 self._atoms = copy.deepcopy(value)

622 self._atoms.unrestricted = False

623 self._atoms.occ.fill()

624 if not self._atoms.is_built:

625 self._atoms = self._atoms.build()

626

627

628class USCF(SCF):

629 """SCF class for spin-polarized systems.

630

631 Inherited from :class:`eminus.scf.SCF`.

632

633 In difference to the SCF class, this class will not build the original Atoms object, only the

634 one attributed to the class. Customized fillings could be overwritten when using this class.

635 """

636

637 @property

638 def atoms(self):

639 """Atoms object."""

640 return self._atoms

641

642 @atoms.setter

643 def atoms(self, value):

644 self._atoms = copy.deepcopy(value)

645 self._atoms.unrestricted = True

646 self._atoms.occ.fill()

647 if not self._atoms.is_built:

648 self._atoms = self._atoms.build()