1"""This module defines an ASE-calculator interface to GPAW. 2 3The central object that glues everything together. 4""" 5 6import warnings 7from typing import Any, Dict 8 9import numpy as np 10from ase import Atoms 11from ase.calculators.calculator import Calculator, kpts2ndarray 12from ase.dft.bandgap import bandgap 13from ase.units import Bohr, Ha 14from ase.utils import plural 15from ase.utils.timing import Timer 16 17import gpaw 18import gpaw.mpi as mpi 19import gpaw.wavefunctions.pw as pw 20from gpaw.band_descriptor import BandDescriptor 21from gpaw.density import RealSpaceDensity 22from gpaw.dos import DOSCalculator 23from gpaw.eigensolvers import get_eigensolver 24from gpaw.external import PointChargePotential 25from gpaw.forces import calculate_forces 26from gpaw.grid_descriptor import GridDescriptor 27from gpaw.hamiltonian import RealSpaceHamiltonian 28from gpaw.hybrids import HybridXC 29from gpaw.io import Reader, Writer 30from gpaw.io.logger import GPAWLogger 31from gpaw.jellium import create_background_charge 32from gpaw.kohnsham_layouts import get_KohnSham_layouts 33from gpaw.kpt_descriptor import KPointDescriptor 34from gpaw.kpt_refine import create_kpoint_descriptor_with_refinement 35from gpaw.matrix import suggest_blocking 36from gpaw.occupations import ParallelLayout, create_occ_calc 37from gpaw.output import (print_cell, print_parallelization_details, 38 print_positions) 39from gpaw.scf import SCFLoop 40from gpaw.setup import Setups 41from gpaw.stress import calculate_stress 42from gpaw.symmetry import Symmetry 43from gpaw.utilities import check_atoms_too_close 44from gpaw.utilities.gpts import get_number_of_grid_points 45from gpaw.utilities.grid import GridRedistributor 46from gpaw.utilities.memory import MemNode, maxrss 47from gpaw.utilities.partition import AtomPartition 48from gpaw.wavefunctions.mode import create_wave_function_mode 49from gpaw.xc import XC 50from gpaw.xc.kernel import XCKernel 51from gpaw.xc.sic import SIC 52from gpaw.typing import Array1D 53 54 55class GPAW(Calculator): 56 """This is the ASE-calculator frontend for doing a GPAW calculation.""" 57 58 implemented_properties = ['energy', 'forces', 'stress', 'dipole', 59 'magmom', 'magmoms'] 60 61 default_parameters: Dict[str, Any] = { 62 'mode': 'fd', 63 'xc': 'LDA', 64 'occupations': None, 65 'poissonsolver': None, 66 'h': None, # Angstrom 67 'gpts': None, 68 'kpts': [(0.0, 0.0, 0.0)], 69 'nbands': None, 70 'charge': 0, 71 'setups': {}, 72 'basis': {}, 73 'spinpol': None, 74 'filter': None, 75 'mixer': None, 76 'eigensolver': None, 77 'background_charge': None, 78 'experimental': {'reuse_wfs_method': 'paw', 79 'niter_fixdensity': 0, 80 'magmoms': None, 81 'soc': None, 82 'kpt_refine': None}, 83 'external': None, 84 'random': False, 85 'hund': False, 86 'maxiter': 333, 87 'idiotproof': True, 88 'symmetry': {'point_group': True, 89 'time_reversal': True, 90 'symmorphic': True, 91 'tolerance': 1e-7, 92 'do_not_symmetrize_the_density': None}, # deprecated 93 'convergence': {'energy': 0.0005, # eV / electron 94 'density': 1.0e-4, 95 'eigenstates': 4.0e-8, # eV^2 96 'bands': 'occupied', 97 'forces': np.inf}, # eV / Ang 98 'verbose': 0, 99 'fixdensity': False, # deprecated 100 'dtype': None} # deprecated 101 102 default_parallel: Dict[str, Any] = { 103 'kpt': None, 104 'domain': None, 105 'band': None, 106 'order': 'kdb', 107 'stridebands': False, 108 'augment_grids': False, 109 'sl_auto': False, 110 'sl_default': None, 111 'sl_diagonalize': None, 112 'sl_inverse_cholesky': None, 113 'sl_lcao': None, 114 'sl_lrtddft': None, 115 'use_elpa': False, 116 'elpasolver': '2stage', 117 'buffer_size': None} 118 119 def __init__(self, 120 restart=None, 121 *, 122 label=None, 123 timer=None, 124 communicator=None, 125 txt='?', 126 parallel=None, 127 **kwargs): 128 129 if txt == '?': 130 txt = '-' if restart is None else None 131 132 self.parallel = dict(self.default_parallel) 133 if parallel: 134 for key in parallel: 135 if key not in self.default_parallel: 136 allowed = ', '.join(list(self.default_parallel.keys())) 137 raise TypeError('Unexpected keyword "{}" in "parallel" ' 138 'dictionary. Must be one of: {}' 139 .format(key, allowed)) 140 self.parallel.update(parallel) 141 142 if timer is None: 143 self.timer = Timer() 144 else: 145 self.timer = timer 146 147 self.scf = None 148 self.wfs = None 149 self.density = None 150 self.hamiltonian = None 151 self.spos_ac = None # XXX store this in some better way. 152 153 self.observers = [] # XXX move to self.scf 154 self.initialized = False 155 156 self.world = communicator 157 if self.world is None: 158 self.world = mpi.world 159 elif not hasattr(self.world, 'new_communicator'): 160 self.world = mpi.world.new_communicator(np.asarray(self.world)) 161 162 self.log = GPAWLogger(world=self.world) 163 self.log.fd = txt 164 165 self.reader = None 166 167 Calculator.__init__(self, restart, label=label, **kwargs) 168 169 def fixed_density(self, *, 170 update_fermi_level: bool = False, 171 communicator=None, 172 txt='-', 173 parallel: Dict[str, Any] = None, 174 **kwargs) -> 'GPAW': 175 """Create new calculator and do SCF calculation with fixed density. 176 177 Returns a new GPAW object fully converged. 178 179 Useful for band-structure calculations. Given a ground-state 180 calculation, ``gs_calc``, one can do:: 181 182 bs_calc = gs_calc.fixed_density(kpts=<path>, 183 symmetry='off') 184 bs = bs_calc.get_band_structure() 185 """ 186 assert not update_fermi_level # for now ... 187 188 for key in kwargs: 189 if key not in {'nbands', 'occupations', 'poissonsolver', 'kpts', 190 'eigensolver', 'random', 'maxiter', 'basis', 191 'symmetry', 'convergence', 'verbose'}: 192 raise TypeError(f'{key:!r} is an invalid keyword argument') 193 194 params = self.parameters.copy() 195 params.update(kwargs) 196 197 if params['h'] is None: 198 # Backwards compatibility 199 params['gpts'] = self.density.gd.N_c 200 201 calc = GPAW(communicator=communicator, 202 txt=txt, 203 parallel=parallel, 204 **params) 205 calc.initialize(self.atoms) 206 calc.density.initialize_from_other_density(self.density, 207 calc.wfs.kptband_comm) 208 calc.density.fixed = True 209 calc.wfs.fermi_levels = self.wfs.fermi_levels 210 if calc.hamiltonian.xc.type == 'GLLB': 211 new_response = calc.hamiltonian.xc.response 212 old_response = self.hamiltonian.xc.response 213 new_response.initialize_from_other_response(old_response) 214 new_response.fix_potential = True 215 calc.calculate(system_changes=[]) 216 return calc 217 218 def __del__(self): 219 # Write timings and close reader if necessary. 220 # If we crashed in the constructor (e.g. a bad keyword), we may not 221 # have the normally expected attributes: 222 if hasattr(self, 'timer') and not self.log.fd.closed: 223 self.timer.write(self.log.fd) 224 225 if hasattr(self, 'reader') and self.reader is not None: 226 self.reader.close() 227 228 def write(self, filename, mode=''): 229 self.log(f'Writing to {filename} (mode={mode!r})\n') 230 writer = Writer(filename, self.world) 231 self._write(writer, mode) 232 writer.close() 233 self.world.barrier() 234 235 def _write(self, writer, mode): 236 from ase.io.trajectory import write_atoms 237 writer.write(version=3, gpaw_version=gpaw.__version__, 238 ha=Ha, bohr=Bohr) 239 240 write_atoms(writer.child('atoms'), self.atoms) 241 writer.child('results').write(**self.results) 242 writer.child('parameters').write(**self.todict()) 243 244 self.density.write(writer.child('density')) 245 self.hamiltonian.write(writer.child('hamiltonian')) 246 # self.occupations.write(writer.child('occupations')) 247 self.scf.write(writer.child('scf')) 248 self.wfs.write(writer.child('wave_functions'), mode == 'all') 249 250 return writer 251 252 def _set_atoms(self, atoms): 253 check_atoms_too_close(atoms) 254 self.atoms = atoms 255 # GPAW works in terms of the scaled positions. We want to 256 # extract the scaled positions in only one place, and that is 257 # here. No other place may recalculate them, or we might end up 258 # with rounding errors and inconsistencies. 259 self.spos_ac = atoms.get_scaled_positions() % 1.0 260 261 def read(self, filename): 262 from ase.io.trajectory import read_atoms 263 self.log('Reading from {}'.format(filename)) 264 265 self.reader = reader = Reader(filename) 266 267 atoms = read_atoms(reader.atoms) 268 self._set_atoms(atoms) 269 270 res = reader.results 271 self.results = dict((key, res.get(key)) for key in res.keys()) 272 if self.results: 273 self.log('Read {}'.format(', '.join(sorted(self.results)))) 274 275 self.log('Reading input parameters:') 276 # XXX param 277 self.parameters = self.get_default_parameters() 278 dct = {} 279 for key, value in reader.parameters.asdict().items(): 280 if key in {'txt', 'fixdensity'}: 281 continue # old gpw-files may have these 282 if (isinstance(value, dict) and 283 isinstance(self.parameters[key], dict)): 284 self.parameters[key].update(value) 285 else: 286 self.parameters[key] = value 287 dct[key] = self.parameters[key] 288 289 self.log.print_dict(dct) 290 self.log() 291 292 self.initialize(reading=True) 293 294 self.density.read(reader) 295 self.hamiltonian.read(reader) 296 self.scf.read(reader) 297 self.wfs.read(reader) 298 299 # We need to do this in a better way: XXX 300 from gpaw.utilities.partition import AtomPartition 301 atom_partition = AtomPartition(self.wfs.gd.comm, 302 np.zeros(len(self.atoms), dtype=int)) 303 self.density.atom_partition = atom_partition 304 self.hamiltonian.atom_partition = atom_partition 305 rank_a = self.density.gd.get_ranks_from_positions(self.spos_ac) 306 new_atom_partition = AtomPartition(self.density.gd.comm, rank_a) 307 for obj in [self.density, self.hamiltonian]: 308 obj.set_positions_without_ruining_everything(self.spos_ac, 309 new_atom_partition) 310 if new_atom_partition != atom_partition: 311 for kpt in self.wfs.kpt_u: 312 kpt.projections = kpt.projections.redist(new_atom_partition) 313 self.wfs.atom_partition = new_atom_partition 314 315 self.hamiltonian.xc.read(reader) 316 317 return reader 318 319 def check_state(self, atoms, tol=1e-12): 320 system_changes = Calculator.check_state(self, atoms, tol) 321 if 'positions' not in system_changes: 322 if self.hamiltonian: 323 if self.hamiltonian.vext: 324 if self.hamiltonian.vext.vext_g is None: 325 # QMMM atoms have moved: 326 system_changes.append('positions') 327 return system_changes 328 329 def calculate(self, atoms=None, properties=['energy'], 330 system_changes=['cell']): 331 for _ in self.icalculate(atoms, properties, system_changes): 332 pass 333 334 def icalculate(self, atoms=None, properties=['energy'], 335 system_changes=['cell']): 336 """Calculate things.""" 337 338 Calculator.calculate(self, atoms) 339 atoms = self.atoms 340 341 if system_changes: 342 self.log('System changes:', ', '.join(system_changes), '\n') 343 if system_changes == ['positions']: 344 # Only positions have changed: 345 self.density.reset() 346 else: 347 # Drastic changes: 348 self.wfs = None 349 self.density = None 350 self.hamiltonian = None 351 self.scf = None 352 self.initialize(atoms) 353 354 self.set_positions(atoms) 355 356 if not self.initialized: 357 self.initialize(atoms) 358 self.set_positions(atoms) 359 360 if not (self.wfs.positions_set and self.hamiltonian.positions_set): 361 self.set_positions(atoms) 362 363 yield 364 365 if not self.scf.converged: 366 print_cell(self.wfs.gd, self.atoms.pbc, self.log) 367 368 with self.timer('SCF-cycle'): 369 yield from self.scf.irun( 370 self.wfs, self.hamiltonian, 371 self.density, 372 self.log, self.call_observers) 373 374 self.log('\nConverged after {} iterations.\n' 375 .format(self.scf.niter)) 376 377 e_free = self.hamiltonian.e_total_free 378 e_extrapolated = self.hamiltonian.e_total_extrapolated 379 self.results['energy'] = e_extrapolated * Ha 380 self.results['free_energy'] = e_free * Ha 381 382 dipole_v = self.density.calculate_dipole_moment() * Bohr 383 self.log('Dipole moment: ({:.6f}, {:.6f}, {:.6f}) |e|*Ang\n' 384 .format(*dipole_v)) 385 self.results['dipole'] = dipole_v 386 387 if self.wfs.nspins == 2 or not self.density.collinear: 388 totmom_v, magmom_av = self.density.estimate_magnetic_moments() 389 self.log('Total magnetic moment: ({:.6f}, {:.6f}, {:.6f})' 390 .format(*totmom_v)) 391 self.log('Local magnetic moments:') 392 symbols = self.atoms.get_chemical_symbols() 393 for a, mom_v in enumerate(magmom_av): 394 self.log('{:4} {:2} ({:9.6f}, {:9.6f}, {:9.6f})' 395 .format(a, symbols[a], *mom_v)) 396 self.log() 397 self.results['magmom'] = totmom_v[2] 398 self.results['magmoms'] = magmom_av[:, 2].copy() 399 else: 400 self.results['magmom'] = 0.0 401 self.results['magmoms'] = np.zeros(len(self.atoms)) 402 403 self.summary() 404 405 self.call_observers(self.scf.niter, final=True) 406 407 if 'forces' in properties: 408 with self.timer('Forces'): 409 F_av = calculate_forces(self.wfs, self.density, 410 self.hamiltonian, self.log) 411 self.results['forces'] = F_av * (Ha / Bohr) 412 413 if 'stress' in properties: 414 with self.timer('Stress'): 415 try: 416 stress = calculate_stress(self).flat[[0, 4, 8, 5, 2, 1]] 417 except NotImplementedError: 418 # Our ASE Calculator base class will raise 419 # PropertyNotImplementedError for us. 420 pass 421 else: 422 self.results['stress'] = stress * (Ha / Bohr**3) 423 424 def summary(self): 425 self.hamiltonian.summary(self.wfs, self.log) 426 self.density.summary(self.atoms, self.results.get('magmom', 0.0), 427 self.log) 428 self.wfs.summary(self.log) 429 if len(self.wfs.fermi_levels) == 1: 430 try: 431 bandgap(self, 432 output=self.log.fd, 433 efermi=self.wfs.fermi_level * Ha) 434 except ValueError: 435 pass 436 self.log.fd.flush() 437 438 def set(self, **kwargs): 439 """Change parameters for calculator. 440 441 Examples:: 442 443 calc.set(xc='PBE') 444 calc.set(nbands=20, kpts=(4, 1, 1)) 445 """ 446 447 # Verify that keys are consistent with default ones. 448 for key in kwargs: 449 if key != 'txt' and key not in self.default_parameters: 450 raise TypeError('Unknown GPAW parameter: {}'.format(key)) 451 452 if key in ['convergence', 'symmetry', 453 'experimental'] and isinstance(kwargs[key], dict): 454 # For values that are dictionaries, verify subkeys, too. 455 default_dict = self.default_parameters[key] 456 for subkey in kwargs[key]: 457 if subkey not in default_dict: 458 allowed = ', '.join(list(default_dict.keys())) 459 raise TypeError('Unknown subkeyword "{}" of keyword ' 460 '"{}". Must be one of: {}' 461 .format(subkey, key, allowed)) 462 463 changed_parameters = Calculator.set(self, **kwargs) 464 465 for key in ['setups', 'basis']: 466 if key in changed_parameters: 467 dct = changed_parameters[key] 468 if isinstance(dct, dict) and None in dct: 469 dct['default'] = dct.pop(None) 470 warnings.warn('Please use {key}={dct}' 471 .format(key=key, dct=dct)) 472 473 # We need to handle txt early in order to get logging up and running: 474 if 'txt' in changed_parameters: 475 self.log.fd = changed_parameters.pop('txt') 476 477 if not changed_parameters: 478 return {} 479 480 self.initialized = False 481 self.scf = None 482 self.results = {} 483 484 self.log('Input parameters:') 485 self.log.print_dict(changed_parameters) 486 self.log() 487 488 for key in changed_parameters: 489 if key in ['eigensolver', 'convergence'] and self.wfs: 490 self.wfs.set_eigensolver(None) 491 492 if key in ['mixer', 'verbose', 'txt', 'hund', 'random', 493 'eigensolver', 'idiotproof']: 494 continue 495 496 if key in ['convergence', 'fixdensity', 'maxiter']: 497 continue 498 499 # Check nested arguments 500 if key in ['experimental']: 501 changed_parameters2 = changed_parameters[key] 502 for key2 in changed_parameters2: 503 if key2 in ['kpt_refine', 'magmoms', 'soc']: 504 self.wfs = None 505 elif key2 in ['reuse_wfs_method', 'niter_fixdensity']: 506 continue 507 else: 508 raise TypeError('Unknown keyword argument:', key2) 509 continue 510 511 # More drastic changes: 512 if self.wfs: 513 self.wfs.set_orthonormalized(False) 514 if key in ['xc', 'poissonsolver']: 515 self.hamiltonian = None 516 elif key in ['occupations', 'width']: 517 pass 518 elif key in ['external', 'charge', 'background_charge']: 519 self.hamiltonian = None 520 self.density = None 521 self.wfs = None 522 elif key in ['kpts', 'nbands', 'symmetry']: 523 self.wfs = None 524 elif key in ['h', 'gpts', 'setups', 'spinpol', 'dtype', 'mode']: 525 self.density = None 526 self.hamiltonian = None 527 self.wfs = None 528 elif key in ['basis']: 529 self.wfs = None 530 else: 531 raise TypeError('Unknown keyword argument: "%s"' % key) 532 533 def initialize_positions(self, atoms=None): 534 """Update the positions of the atoms.""" 535 self.log('Initializing position-dependent things.\n') 536 if atoms is None: 537 atoms = self.atoms 538 else: 539 atoms = atoms.copy() 540 self._set_atoms(atoms) 541 542 mpi.synchronize_atoms(atoms, self.world) 543 544 rank_a = self.wfs.gd.get_ranks_from_positions(self.spos_ac) 545 atom_partition = AtomPartition(self.wfs.gd.comm, rank_a, name='gd') 546 self.wfs.set_positions(self.spos_ac, atom_partition) 547 self.density.set_positions(self.spos_ac, atom_partition) 548 self.hamiltonian.set_positions(self.spos_ac, atom_partition) 549 550 def set_positions(self, atoms=None): 551 """Update the positions of the atoms and initialize wave functions.""" 552 self.initialize_positions(atoms) 553 554 nlcao, nrand = self.wfs.initialize(self.density, self.hamiltonian, 555 self.spos_ac) 556 if nlcao + nrand: 557 self.log('Creating initial wave functions:') 558 if nlcao: 559 self.log(' ', plural(nlcao, 'band'), 'from LCAO basis set') 560 if nrand: 561 self.log(' ', plural(nrand, 'band'), 'from random numbers') 562 self.log() 563 564 self.wfs.eigensolver.reset() 565 self.scf.reset() 566 occ_name = getattr(self.wfs.occupations, "name", None) 567 if occ_name == 'mom': 568 # Initialize MOM reference orbitals 569 self.wfs.occupations.initialize_reference_orbitals() 570 print_positions(self.atoms, self.log, self.density.magmom_av) 571 572 def initialize(self, atoms=None, reading=False): 573 """Inexpensive initialization.""" 574 575 self.log('Initialize ...\n') 576 577 if atoms is None: 578 atoms = self.atoms 579 else: 580 atoms = atoms.copy() 581 self._set_atoms(atoms) 582 583 par = self.parameters 584 585 natoms = len(atoms) 586 587 cell_cv = atoms.get_cell() / Bohr 588 number_of_lattice_vectors = cell_cv.any(axis=1).sum() 589 if number_of_lattice_vectors < 3: 590 raise ValueError( 591 'GPAW requires 3 lattice vectors. Your system has {}.' 592 .format(number_of_lattice_vectors)) 593 594 pbc_c = atoms.get_pbc() 595 assert len(pbc_c) == 3 596 magmom_a = atoms.get_initial_magnetic_moments() 597 598 if par.experimental.get('magmoms') is not None: 599 magmom_av = np.array(par.experimental['magmoms'], float) 600 collinear = False 601 else: 602 magmom_av = np.zeros((natoms, 3)) 603 magmom_av[:, 2] = magmom_a 604 collinear = True 605 606 mpi.synchronize_atoms(atoms, self.world) 607 608 # Generate new xc functional only when it is reset by set 609 # XXX sounds like this should use the _changed_keywords dictionary. 610 if self.hamiltonian is None or self.hamiltonian.xc is None: 611 if isinstance(par.xc, (str, dict, XCKernel)): 612 xc = XC(par.xc, collinear=collinear, atoms=atoms) 613 else: 614 xc = par.xc 615 else: 616 xc = self.hamiltonian.xc 617 618 if par.fixdensity: 619 warnings.warn( 620 'The fixdensity keyword has been deprecated. ' 621 'Please use the GPAW.fixed_density() method instead.', 622 DeprecationWarning) 623 624 mode = par.mode 625 if isinstance(mode, str): 626 mode = {'name': mode} 627 if isinstance(mode, dict): 628 mode = create_wave_function_mode(**mode) 629 630 if par.dtype == complex: 631 warnings.warn('Use mode={}(..., force_complex_dtype=True) ' 632 'instead of dtype=complex'.format(mode.name.upper())) 633 mode.force_complex_dtype = True 634 del par['dtype'] 635 par.mode = mode 636 637 if xc.orbital_dependent and mode.name == 'lcao': 638 raise ValueError('LCAO mode does not support ' 639 'orbital-dependent XC functionals.') 640 641 realspace = (mode.name != 'pw' and mode.interpolation != 'fft') 642 643 self.create_setups(mode, xc) 644 645 if not realspace: 646 pbc_c = np.ones(3, bool) 647 648 if par.hund: 649 if natoms != 1: 650 raise ValueError('hund=True arg only valid for single atoms!') 651 spinpol = True 652 magmom_av[0, 2] = self.setups[0].get_hunds_rule_moment(par.charge) 653 654 if collinear: 655 magnetic = magmom_av.any() 656 657 spinpol = par.spinpol 658 if spinpol is None: 659 spinpol = magnetic 660 elif magnetic and not spinpol: 661 raise ValueError('Non-zero initial magnetic moment for a ' + 662 'spin-paired calculation!') 663 nspins = 1 + int(spinpol) 664 665 if spinpol: 666 self.log('Spin-polarized calculation.') 667 self.log('Magnetic moment: {:.6f}\n'.format(magmom_av.sum())) 668 else: 669 self.log('Spin-paired calculation\n') 670 else: 671 nspins = 1 672 self.log('Non-collinear calculation.') 673 self.log('Magnetic moment: ({:.6f}, {:.6f}, {:.6f})\n' 674 .format(*magmom_av.sum(0))) 675 676 self.create_symmetry(magmom_av, cell_cv, reading) 677 678 if par.gpts is not None: 679 if par.h is not None: 680 raise ValueError("""You can't use both "gpts" and "h"!""") 681 N_c = np.array(par.gpts) 682 h = None 683 else: 684 h = par.h 685 if h is not None: 686 h /= Bohr 687 if h is None and reading: 688 shape = self.reader.density.proxy('density').shape[-3:] 689 N_c = 1 - pbc_c + shape 690 elif h is None and self.density is not None: 691 N_c = self.density.gd.N_c 692 else: 693 N_c = get_number_of_grid_points(cell_cv, h, mode, realspace, 694 self.symmetry, self.log) 695 696 self.setups.set_symmetry(self.symmetry) 697 698 if isinstance(par.background_charge, dict): 699 background = create_background_charge(**par.background_charge) 700 else: 701 background = par.background_charge 702 703 nao = self.setups.nao 704 nvalence = self.setups.nvalence - par.charge 705 if par.background_charge is not None: 706 nvalence += background.charge 707 708 M = np.linalg.norm(magmom_av.sum(0)) 709 710 nbands = par.nbands 711 712 orbital_free = any(setup.orbital_free for setup in self.setups) 713 if orbital_free: 714 nbands = 1 715 716 if isinstance(nbands, str): 717 if nbands == 'nao': 718 nbands = nao 719 elif nbands[-1] == '%': 720 basebands = (nvalence + M) / 2 721 nbands = int(np.ceil(float(nbands[:-1]) / 100 * basebands)) 722 else: 723 raise ValueError('Integer expected: Only use a string ' 724 'if giving a percentage of occupied bands') 725 726 if nbands is None: 727 # Number of bound partial waves: 728 nbandsmax = sum(setup.get_default_nbands() 729 for setup in self.setups) 730 nbands = int(np.ceil((1.2 * (nvalence + M) / 2))) + 4 731 if nbands > nbandsmax: 732 nbands = nbandsmax 733 if mode.name == 'lcao' and nbands > nao: 734 nbands = nao 735 elif nbands <= 0: 736 nbands = max(1, int(nvalence + M + 0.5) // 2 + (-nbands)) 737 738 if nbands > nao and mode.name == 'lcao': 739 raise ValueError('Too many bands for LCAO calculation: ' 740 '%d bands and only %d atomic orbitals!' % 741 (nbands, nao)) 742 743 if nvalence < 0: 744 raise ValueError( 745 'Charge %f is not possible - not enough valence electrons' % 746 par.charge) 747 748 if nvalence > 2 * nbands and not orbital_free: 749 raise ValueError('Too few bands! Electrons: %f, bands: %d' 750 % (nvalence, nbands)) 751 752 if self.scf is None: 753 self.create_scf(nvalence, mode) 754 755 if not collinear: 756 nbands *= 2 757 758 if not self.wfs: 759 self.create_wave_functions(mode, realspace, 760 nspins, collinear, nbands, nao, 761 nvalence, self.setups, 762 cell_cv, pbc_c, N_c, 763 xc) 764 else: 765 self.wfs.set_setups(self.setups) 766 767 occ = self.create_occupations(cell_cv, magmom_av[:, 2].sum(), 768 orbital_free) 769 self.wfs.occupations = occ 770 771 if not self.wfs.eigensolver: 772 self.create_eigensolver(xc, nbands, mode) 773 774 if self.density is None and not reading: 775 assert not par.fixdensity, 'No density to fix!' 776 777 olddens = None 778 if (self.density is not None and 779 (self.density.gd.parsize_c != self.wfs.gd.parsize_c).any()): 780 # Domain decomposition has changed, so we need to 781 # reinitialize density and hamiltonian: 782 if par.fixdensity: 783 olddens = self.density 784 785 self.density = None 786 self.hamiltonian = None 787 788 if self.density is None: 789 self.create_density(realspace, mode, background, h) 790 791 # XXXXXXXXXX if setups change, then setups.core_charge may change. 792 # But that parameter was supplied in Density constructor! 793 # This surely is a bug! 794 self.density.initialize(self.setups, self.timer, 795 magmom_av, par.hund) 796 self.density.set_mixer(par.mixer) 797 if self.density.mixer.driver.name == 'dummy' or par.fixdensity: 798 self.log('No density mixing\n') 799 else: 800 self.log(self.density.mixer, '\n') 801 self.density.fixed = par.fixdensity 802 self.density.log = self.log 803 804 if olddens is not None: 805 self.density.initialize_from_other_density(olddens, 806 self.wfs.kptband_comm) 807 808 if self.hamiltonian is None: 809 self.create_hamiltonian(realspace, mode, xc) 810 811 xc.initialize(self.density, self.hamiltonian, self.wfs) 812 813 description = xc.get_description() 814 if description is not None: 815 self.log('XC parameters: {}\n' 816 .format('\n '.join(description.splitlines()))) 817 818 if xc.type == 'GLLB' and olddens is not None: 819 xc.heeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeelp(olddens) 820 821 self.print_memory_estimate(maxdepth=3) 822 823 print_parallelization_details(self.wfs, self.hamiltonian, self.log) 824 825 self.log('Number of atoms:', natoms) 826 self.log('Number of atomic orbitals:', self.wfs.setups.nao) 827 self.log('Number of bands in calculation:', self.wfs.bd.nbands) 828 self.log('Number of valence electrons:', self.wfs.nvalence) 829 830 n = par.convergence.get('bands', 'occupied') 831 if isinstance(n, int) and n < 0: 832 n += self.wfs.bd.nbands 833 self.log('Bands to converge:', n) 834 835 self.log(flush=True) 836 837 self.timer.print_info(self) 838 839 if gpaw.dry_run: 840 self.dry_run() 841 842 if (realspace and 843 self.hamiltonian.poisson.get_description() == 'FDTD+TDDFT'): 844 self.hamiltonian.poisson.set_density(self.density) 845 self.hamiltonian.poisson.print_messages(self.log) 846 self.log.fd.flush() 847 848 self.initialized = True 849 self.log('... initialized\n') 850 851 def create_setups(self, mode, xc): 852 if self.parameters.filter is None and mode.name != 'pw': 853 gamma = 1.6 854 N_c = self.parameters.get('gpts') 855 if N_c is None: 856 h = (self.parameters.h or 0.2) / Bohr 857 else: 858 icell_vc = np.linalg.inv(self.atoms.cell) 859 h = ((icell_vc**2).sum(0)**-0.5 / N_c).max() / Bohr 860 861 def filter(rgd, rcut, f_r, l=0): 862 gcut = np.pi / h - 2 / rcut / gamma 863 ftmp = rgd.filter(f_r, rcut * gamma, gcut, l) 864 f_r[:] = ftmp[:len(f_r)] 865 else: 866 filter = self.parameters.filter 867 868 Z_a = self.atoms.get_atomic_numbers() 869 self.setups = Setups(Z_a, 870 self.parameters.setups, self.parameters.basis, 871 xc, filter, self.world) 872 self.log(self.setups) 873 874 def create_grid_descriptor(self, N_c, cell_cv, pbc_c, 875 domain_comm, parsize_domain): 876 return GridDescriptor(N_c, cell_cv, pbc_c, domain_comm, 877 parsize_domain) 878 879 def create_occupations(self, cell_cv, magmom, orbital_free): 880 dct = self.parameters.occupations 881 882 if dct is None: 883 if orbital_free: 884 dct = {'name': 'orbital-free'} 885 else: 886 if self.atoms.pbc.any(): 887 dct = {'name': 'fermi-dirac', 888 'width': 0.1} # eV 889 else: 890 dct = {'width': 0.0} 891 elif not isinstance(dct, dict): 892 return dct 893 894 if self.wfs.nspins == 1: 895 dct.pop('fixmagmom', None) 896 897 kwargs = dct.copy() 898 name = kwargs.pop('name', '') 899 if name == 'mom': 900 from gpaw.mom import OccupationsMOM 901 occ = OccupationsMOM(self.wfs, **kwargs) 902 903 self.log(occ) 904 return occ 905 906 occ = create_occ_calc( 907 dct, 908 parallel_layout=ParallelLayout(self.wfs.bd, 909 self.wfs.kd.comm, 910 self.wfs.gd.comm), 911 fixed_magmom_value=magmom, 912 rcell=np.linalg.inv(cell_cv).T, 913 monkhorst_pack_size=self.wfs.kd.N_c, 914 bz2ibzmap=self.wfs.kd.bz2ibz_k) 915 916 self.log('Occupation numbers:', occ, '\n') 917 return occ 918 919 def create_scf(self, nvalence, mode): 920 # if mode.name == 'lcao': 921 # niter_fixdensity = 0 922 # else: 923 # niter_fixdensity = 2 924 925 nv = max(nvalence, 1) 926 cc = self.parameters.convergence 927 self.scf = SCFLoop( 928 cc.get('eigenstates', 4.0e-8) / Ha**2 * nv, 929 cc.get('energy', 0.0005) / Ha * nv, 930 cc.get('density', 1.0e-4) * nv, 931 cc.get('forces', np.inf) / (Ha / Bohr), 932 self.parameters.maxiter, 933 # XXX make sure niter_fixdensity value is *always* set from default 934 # Subdictionary defaults seem to not be set when user provides 935 # e.g. {}. We should change that so it works like the ordinary 936 # parameters. 937 self.parameters.experimental.get('niter_fixdensity', 0), 938 nv) 939 self.log(self.scf) 940 941 def create_symmetry(self, magmom_av, cell_cv, reading): 942 symm = self.parameters.symmetry 943 if symm == 'off': 944 symm = {'point_group': False, 'time_reversal': False} 945 946 if 'do_not_symmetrize_the_density' in symm: 947 symm = symm.copy() 948 dnstd = symm.pop('do_not_symmetrize_the_density') 949 if dnstd is not None: 950 info = 'Ignoring your "do_not_symmetrize_the_density" keyword!' 951 if reading: 952 self.log(info) 953 else: 954 warnings.warn(info + ' Please remove it.') 955 956 if self.parameters.external is not None: 957 symm = symm.copy() 958 symm['point_group'] = False 959 960 if reading and self.reader.version <= 1: 961 symm['allow_invert_aperiodic_axes'] = False 962 963 m_av = magmom_av.round(decimals=3) # round off 964 id_a = [id + tuple(m_v) for id, m_v in zip(self.setups.id_a, m_av)] 965 self.symmetry = Symmetry(id_a, cell_cv, self.atoms.pbc, **symm) 966 self.symmetry.analyze(self.spos_ac) 967 968 def create_eigensolver(self, xc, nbands, mode): 969 # Number of bands to converge: 970 nbands_converge = self.parameters.convergence.get('bands', 'occupied') 971 if nbands_converge == 'all': 972 nbands_converge = nbands 973 elif isinstance(nbands_converge, int): 974 if nbands_converge < 0: 975 nbands_converge += nbands 976 eigensolver = get_eigensolver(self.parameters.eigensolver, mode, 977 self.parameters.convergence) 978 eigensolver.nbands_converge = nbands_converge 979 # XXX Eigensolver class doesn't define an nbands_converge property 980 981 if isinstance(xc, SIC): 982 eigensolver.blocksize = 1 983 984 self.wfs.set_eigensolver(eigensolver) 985 986 self.log('Eigensolver\n ', self.wfs.eigensolver, '\n') 987 988 def create_density(self, realspace, mode, background, h): 989 gd = self.wfs.gd 990 991 big_gd = gd.new_descriptor(comm=self.world) 992 # Check whether grid is too small. 8 is smallest admissible. 993 # (we decide this by how difficult it is to make the tests pass) 994 # (Actually it depends on stencils! But let the user deal with it) 995 N_c = big_gd.get_size_of_global_array(pad=True) 996 too_small = np.any(N_c / big_gd.parsize_c < 8) 997 if (self.parallel['augment_grids'] and not too_small and 998 mode.name != 'pw'): 999 aux_gd = big_gd 1000 else: 1001 aux_gd = gd 1002 1003 redistributor = GridRedistributor(self.world, 1004 self.wfs.kptband_comm, 1005 gd, aux_gd) 1006 1007 # Construct grid descriptor for fine grids for densities 1008 # and potentials: 1009 finegd = aux_gd.refine() 1010 1011 kwargs = dict( 1012 gd=gd, finegd=finegd, 1013 nspins=self.wfs.nspins, 1014 collinear=self.wfs.collinear, 1015 charge=self.parameters.charge + self.wfs.setups.core_charge, 1016 redistributor=redistributor, 1017 background_charge=background) 1018 1019 if realspace: 1020 self.density = RealSpaceDensity(stencil=mode.interpolation, 1021 **kwargs) 1022 else: 1023 if h is None: 1024 ecut = 2 * self.wfs.pd.ecut 1025 else: 1026 ecut = 0.5 * (np.pi / h)**2 1027 self.density = pw.ReciprocalSpaceDensity(ecut=ecut, 1028 **kwargs) 1029 1030 self.log(self.density, '\n') 1031 1032 def create_hamiltonian(self, realspace, mode, xc): 1033 dens = self.density 1034 kwargs = dict( 1035 gd=dens.gd, finegd=dens.finegd, 1036 nspins=dens.nspins, 1037 collinear=dens.collinear, 1038 setups=dens.setups, 1039 timer=self.timer, 1040 xc=xc, 1041 world=self.world, 1042 redistributor=dens.redistributor, 1043 vext=self.parameters.external, 1044 psolver=self.parameters.poissonsolver) 1045 if realspace: 1046 self.hamiltonian = RealSpaceHamiltonian(stencil=mode.interpolation, 1047 **kwargs) 1048 xc.set_grid_descriptor(self.hamiltonian.finegd) 1049 else: 1050 # This code will work if dens.redistributor uses 1051 # ordinary density.gd as aux_gd 1052 gd = dens.finegd 1053 1054 xc_redist = None 1055 if self.parallel['augment_grids']: 1056 from gpaw.grid_descriptor import BadGridError 1057 try: 1058 aux_gd = gd.new_descriptor(comm=self.world) 1059 except BadGridError as err: 1060 import warnings 1061 warnings.warn('Ignoring augment_grids: {}' 1062 .format(err)) 1063 else: 1064 bcast_comm = dens.redistributor.broadcast_comm 1065 xc_redist = GridRedistributor(self.world, bcast_comm, 1066 gd, aux_gd) 1067 1068 self.hamiltonian = pw.ReciprocalSpaceHamiltonian( 1069 pd2=dens.pd2, pd3=dens.pd3, realpbc_c=self.atoms.pbc, 1070 xc_redistributor=xc_redist, 1071 **kwargs) 1072 xc.set_grid_descriptor(self.hamiltonian.xc_gd) 1073 1074 self.hamiltonian.soc = self.parameters.experimental.get('soc') 1075 self.log(self.hamiltonian, '\n') 1076 1077 def create_kpoint_descriptor(self, nspins): 1078 par = self.parameters 1079 1080 # Zero cell vectors that are not periodic so that ASE's 1081 # kpts2ndarray can handle 1-d and 2-d correctly: 1082 atoms = Atoms(cell=self.atoms.cell * self.atoms.pbc[:, np.newaxis], 1083 pbc=self.atoms.pbc) 1084 bzkpts_kc = kpts2ndarray(par.kpts, atoms) 1085 1086 kpt_refine = par.experimental.get('kpt_refine') 1087 if kpt_refine is None: 1088 kd = KPointDescriptor(bzkpts_kc, nspins) 1089 1090 self.timer.start('Set symmetry') 1091 kd.set_symmetry(self.atoms, self.symmetry, comm=self.world) 1092 self.timer.stop('Set symmetry') 1093 1094 else: 1095 self.timer.start('Set k-point refinement') 1096 kd = create_kpoint_descriptor_with_refinement( 1097 kpt_refine, 1098 bzkpts_kc, nspins, self.atoms, 1099 self.symmetry, comm=self.world, 1100 timer=self.timer) 1101 self.timer.stop('Set k-point refinement') 1102 # Update quantities which might have changed, if symmetry 1103 # was changed 1104 self.symmetry = kd.symmetry 1105 self.setups.set_symmetry(kd.symmetry) 1106 1107 self.log(kd) 1108 1109 return kd 1110 1111 def create_wave_functions(self, mode, realspace, 1112 nspins, collinear, nbands, nao, nvalence, 1113 setups, cell_cv, pbc_c, N_c, xc): 1114 par = self.parameters 1115 1116 kd = self.create_kpoint_descriptor(nspins) 1117 1118 parallelization = mpi.Parallelization(self.world, 1119 kd.nibzkpts) 1120 1121 parsize_kpt = self.parallel['kpt'] 1122 parsize_domain = self.parallel['domain'] 1123 parsize_bands = self.parallel['band'] 1124 1125 if isinstance(xc, HybridXC): 1126 parsize_kpt = 1 1127 parsize_domain = self.world.size 1128 parsize_bands = 1 1129 1130 ndomains = None 1131 if parsize_domain is not None: 1132 ndomains = np.prod(parsize_domain) 1133 parallelization.set(kpt=parsize_kpt, 1134 domain=ndomains, 1135 band=parsize_bands) 1136 comms = parallelization.build_communicators() 1137 domain_comm = comms['d'] 1138 kpt_comm = comms['k'] 1139 band_comm = comms['b'] 1140 kptband_comm = comms['D'] 1141 domainband_comm = comms['K'] 1142 1143 self.comms = comms 1144 1145 kd.set_communicator(kpt_comm) 1146 1147 parstride_bands = self.parallel['stridebands'] 1148 if parstride_bands: 1149 raise RuntimeError('stridebands is unreliable') 1150 1151 bd = BandDescriptor(nbands, band_comm, parstride_bands) 1152 1153 # Construct grid descriptor for coarse grids for wave functions: 1154 gd = self.create_grid_descriptor(N_c, cell_cv, pbc_c, 1155 domain_comm, parsize_domain) 1156 1157 if hasattr(self, 'time') or mode.force_complex_dtype or not collinear: 1158 dtype = complex 1159 else: 1160 if kd.gamma: 1161 dtype = float 1162 else: 1163 dtype = complex 1164 1165 wfs_kwargs = dict(gd=gd, nvalence=nvalence, setups=setups, 1166 bd=bd, dtype=dtype, world=self.world, kd=kd, 1167 kptband_comm=kptband_comm, timer=self.timer) 1168 1169 if self.parallel['sl_auto']: 1170 # Choose scalapack parallelization automatically 1171 1172 for key, val in self.parallel.items(): 1173 if (key.startswith('sl_') and key != 'sl_auto' and 1174 val is not None): 1175 raise ValueError("Cannot use 'sl_auto' together " 1176 "with '%s'" % key) 1177 1178 max_scalapack_cpus = bd.comm.size * gd.comm.size 1179 sl_default = suggest_blocking(nbands, max_scalapack_cpus) 1180 else: 1181 sl_default = self.parallel['sl_default'] 1182 1183 if mode.name == 'lcao': 1184 assert collinear 1185 # Layouts used for general diagonalizer 1186 sl_lcao = self.parallel['sl_lcao'] 1187 if sl_lcao is None: 1188 sl_lcao = sl_default 1189 1190 elpasolver = None 1191 if self.parallel['use_elpa']: 1192 elpasolver = self.parallel['elpasolver'] 1193 lcaoksl = get_KohnSham_layouts(sl_lcao, 'lcao', 1194 gd, bd, domainband_comm, dtype, 1195 nao=nao, timer=self.timer, 1196 elpasolver=elpasolver) 1197 1198 self.wfs = mode(lcaoksl, **wfs_kwargs) 1199 1200 elif mode.name == 'fd' or mode.name == 'pw': 1201 # Use (at most) all available LCAO for initialization 1202 lcaonbands = min(nbands, nao) 1203 1204 try: 1205 lcaobd = BandDescriptor(lcaonbands, band_comm, 1206 parstride_bands) 1207 except RuntimeError: 1208 initksl = None 1209 else: 1210 # Layouts used for general diagonalizer 1211 # (LCAO initialization) 1212 sl_lcao = self.parallel['sl_lcao'] 1213 if sl_lcao is None: 1214 sl_lcao = sl_default 1215 initksl = get_KohnSham_layouts(sl_lcao, 'lcao', 1216 gd, lcaobd, domainband_comm, 1217 dtype, nao=nao, 1218 timer=self.timer) 1219 1220 reuse_wfs_method = par.experimental.get('reuse_wfs_method', 'paw') 1221 sl = (domainband_comm,) + (self.parallel['sl_diagonalize'] or 1222 sl_default or 1223 (1, 1, None)) 1224 self.wfs = mode(sl, initksl, 1225 reuse_wfs_method=reuse_wfs_method, 1226 collinear=collinear, 1227 **wfs_kwargs) 1228 else: 1229 self.wfs = mode(self, collinear=collinear, **wfs_kwargs) 1230 1231 self.log(self.wfs, '\n') 1232 1233 def dry_run(self): 1234 # Can be overridden like in gpaw.atom.atompaw 1235 print_cell(self.wfs.gd, self.atoms.pbc, self.log) 1236 print_positions(self.atoms, self.log, self.density.magmom_av) 1237 self.log.fd.flush() 1238 1239 # Write timing info now before the interpreter shuts down: 1240 self.__del__() 1241 1242 # Disable timing output during shut-down: 1243 del self.timer 1244 1245 raise SystemExit 1246 1247 def get_atomic_electrostatic_potentials(self) -> Array1D: 1248 r"""Return the electrostatic potential at the atomic sites. 1249 1250 Return list of energies in eV, one for each atom: 1251 1252 .. math:: 1253 1254 Y_{00} 1255 \int d\mathbf{r} 1256 \tilde{v}_H(\mathbf{r}) 1257 \hat{g}_{00}^a(\mathbf{r} - \mathbf{R}^a) 1258 1259 """ 1260 ham = self.hamiltonian 1261 dens = self.density 1262 self.initialize_positions() 1263 dens.interpolate_pseudo_density() 1264 dens.calculate_pseudo_charge() 1265 ham.update(dens) 1266 W_aL = ham.calculate_atomic_hamiltonians(dens) 1267 return np.array([W_L[0] / (4 * np.pi)**0.5 * Ha 1268 for W_L in W_aL.values()]) 1269 1270 def linearize_to_xc(self, newxc): 1271 """Linearize Hamiltonian to difference XC functional. 1272 1273 Used in real time TDDFT to perform calculations with various kernels. 1274 """ 1275 if isinstance(newxc, str): 1276 newxc = XC(newxc) 1277 self.log('Linearizing xc-hamiltonian to ' + str(newxc)) 1278 newxc.initialize(self.density, self.hamiltonian, self.wfs) 1279 self.hamiltonian.linearize_to_xc(newxc, self.density) 1280 1281 def attach(self, function, n=1, *args, **kwargs): 1282 """Register observer function. 1283 1284 Call *function* using *args* and 1285 *kwargs* as arguments. 1286 1287 If *n* is positive, then 1288 *function* will be called every *n* iterations + the 1289 final iteration if it would not be otherwise 1290 1291 If *n* is negative, then *function* will only be 1292 called on iteration *abs(n)*. 1293 1294 If *n* is 0, then *function* will only be called 1295 on convergence""" 1296 1297 try: 1298 slf = function.__self__ 1299 except AttributeError: 1300 pass 1301 else: 1302 if slf is self: 1303 # function is a bound method of self. Store the name 1304 # of the method and avoid circular reference: 1305 function = function.__func__.__name__ 1306 1307 # Replace self in args with another unique reference 1308 # to avoid circular reference 1309 if not hasattr(self, 'self_ref'): 1310 self.self_ref = object() 1311 self_ = self.self_ref 1312 args = tuple([self_ if arg is self else arg for arg in args]) 1313 1314 self.observers.append((function, n, args, kwargs)) 1315 1316 def call_observers(self, iter, final=False): 1317 """Call all registered callback functions.""" 1318 for function, n, args, kwargs in self.observers: 1319 call = False 1320 # Call every n iterations, including the last 1321 if n > 0: 1322 if ((iter % n) == 0) != final: 1323 call = True 1324 # Call only on iteration n 1325 elif n < 0 and not final: 1326 if iter == abs(n): 1327 call = True 1328 # Call only on convergence 1329 elif n == 0 and final: 1330 call = True 1331 if call: 1332 if isinstance(function, str): 1333 function = getattr(self, function) 1334 # Replace self reference with self 1335 self_ = self.self_ref 1336 args = tuple([self if arg is self_ else arg for arg in args]) 1337 function(*args, **kwargs) 1338 1339 def get_reference_energy(self): 1340 return self.wfs.setups.Eref * Ha 1341 1342 def get_homo_lumo(self, spin=None): 1343 """Return HOMO and LUMO eigenvalues. 1344 1345 By default, return the true HOMO-LUMO eigenvalues (spin=None). 1346 1347 If spin is 0 or 1, return HOMO-LUMO eigenvalues taken among 1348 only those states with the given spin.""" 1349 return self.wfs.get_homo_lumo(spin) * Ha 1350 1351 def estimate_memory(self, mem): 1352 """Estimate memory use of this object.""" 1353 for name, obj in [('Density', self.density), 1354 ('Hamiltonian', self.hamiltonian), 1355 ('Wavefunctions', self.wfs)]: 1356 obj.estimate_memory(mem.subnode(name)) 1357 1358 def print_memory_estimate(self, log=None, maxdepth=-1): 1359 """Print estimated memory usage for PAW object and components. 1360 1361 maxdepth is the maximum nesting level of displayed components. 1362 1363 The PAW object must be initialize()'d, but needs not have large 1364 arrays allocated.""" 1365 # NOTE. This should work with "--dry-run=N" 1366 # 1367 # However, the initial overhead estimate is wrong if this method 1368 # is called within a real mpirun/gpaw-python context. 1369 if log is None: 1370 log = self.log 1371 log('Memory estimate:') 1372 1373 mem_init = maxrss() # initial overhead includes part of Hamiltonian! 1374 log(' Process memory now: %.2f MiB' % (mem_init / 1024.0**2)) 1375 1376 mem = MemNode('Calculator', 0) 1377 mem.indent = ' ' 1378 try: 1379 self.estimate_memory(mem) 1380 except AttributeError as m: 1381 log('Attribute error: %r' % m) 1382 log('Some object probably lacks estimate_memory() method') 1383 log('Memory breakdown may be incomplete') 1384 mem.calculate_size() 1385 mem.write(log.fd, maxdepth=maxdepth, depth=1) 1386 log() 1387 1388 def converge_wave_functions(self): 1389 """Converge the wave-functions if not present.""" 1390 1391 if self.scf and self.scf.converged: 1392 if isinstance(self.wfs.kpt_u[0].psit_nG, np.ndarray): 1393 return 1394 if self.wfs.kpt_u[0].psit_nG is not None: 1395 self.wfs.initialize_wave_functions_from_restart_file() 1396 return 1397 1398 if not self.initialized: 1399 self.initialize() 1400 1401 self.set_positions() 1402 1403 self.scf.converged = False 1404 fixed = self.density.fixed 1405 self.density.fixed = True 1406 self.calculate(system_changes=[]) 1407 self.density.fixed = fixed 1408 1409 def diagonalize_full_hamiltonian(self, nbands=None, ecut=None, 1410 scalapack=None, 1411 expert=False): 1412 if not self.initialized: 1413 self.initialize() 1414 nbands = self.wfs.diagonalize_full_hamiltonian( 1415 self.hamiltonian, self.atoms, self.log, 1416 nbands, ecut, scalapack, expert) 1417 self.parameters.nbands = nbands 1418 1419 def get_number_of_bands(self) -> int: 1420 """Return the number of bands.""" 1421 return self.wfs.bd.nbands 1422 1423 def get_xc_functional(self) -> str: 1424 """Return the XC-functional identifier. 1425 1426 'LDA', 'PBE', ...""" 1427 1428 xc = self.parameters.get('xc', 'LDA') 1429 if isinstance(xc, dict): 1430 xc = xc['name'] 1431 return xc 1432 1433 def get_number_of_spins(self): 1434 return self.wfs.nspins 1435 1436 def get_spin_polarized(self): 1437 """Is it a spin-polarized calculation?""" 1438 return self.wfs.nspins == 2 1439 1440 def get_bz_k_points(self): 1441 """Return the k-points.""" 1442 return self.wfs.kd.bzk_kc.copy() 1443 1444 def get_ibz_k_points(self): 1445 """Return k-points in the irreducible part of the Brillouin zone.""" 1446 return self.wfs.kd.ibzk_kc.copy() 1447 1448 def get_bz_to_ibz_map(self): 1449 """Return indices from BZ to IBZ.""" 1450 return self.wfs.kd.bz2ibz_k.copy() 1451 1452 def get_k_point_weights(self): 1453 """Weights of the k-points. 1454 1455 The sum of all weights is one.""" 1456 1457 return self.wfs.kd.weight_k 1458 1459 def get_pseudo_density(self, spin=None, gridrefinement=1, 1460 pad=True, broadcast=True): 1461 """Return pseudo-density array. 1462 1463 If *spin* is not given, then the total density is returned. 1464 Otherwise, the spin up or down density is returned (spin=0 or 1465 1).""" 1466 1467 if gridrefinement == 1: 1468 nt_sG = self.density.nt_sG 1469 gd = self.density.gd 1470 elif gridrefinement == 2: 1471 if self.density.nt_sg is None: 1472 self.density.interpolate_pseudo_density() 1473 nt_sG = self.density.nt_sg 1474 gd = self.density.finegd 1475 else: 1476 raise NotImplementedError 1477 1478 if spin is None: 1479 if self.density.nspins == 1: 1480 nt_G = nt_sG[0] 1481 else: 1482 nt_G = nt_sG.sum(axis=0) 1483 else: 1484 if self.density.nspins == 1: 1485 nt_G = 0.5 * nt_sG[0] 1486 else: 1487 nt_G = nt_sG[spin] 1488 1489 nt_G = gd.collect(nt_G, broadcast=broadcast) 1490 1491 if nt_G is None: 1492 return None 1493 1494 if pad: 1495 nt_G = gd.zero_pad(nt_G) 1496 1497 return nt_G / Bohr**3 1498 1499 get_pseudo_valence_density = get_pseudo_density # Don't use this one! 1500 1501 def get_effective_potential(self, spin=0, pad=True, broadcast=True): 1502 """Return pseudo effective-potential.""" 1503 vt_G = self.hamiltonian.gd.collect(self.hamiltonian.vt_sG[spin], 1504 broadcast=broadcast) 1505 if vt_G is None: 1506 return None 1507 1508 if pad: 1509 vt_G = self.hamiltonian.gd.zero_pad(vt_G) 1510 return vt_G * Ha 1511 1512 def get_electrostatic_potential(self): 1513 """Return the electrostatic potential. 1514 1515 This is the potential from the pseudo electron density and the 1516 PAW-compensation charges. So, the electrostatic potential will 1517 only be correct outside the PAW augmentation spheres. 1518 """ 1519 1520 ham = self.hamiltonian 1521 dens = self.density 1522 self.initialize_positions() 1523 dens.interpolate_pseudo_density() 1524 dens.calculate_pseudo_charge() 1525 return ham.get_electrostatic_potential(dens) * Ha 1526 1527 def get_pseudo_density_corrections(self): 1528 """Integrated density corrections. 1529 1530 Returns the integrated value of the difference between the pseudo- 1531 and the all-electron densities at each atom. These are the numbers 1532 you should add to the result of doing e.g. Bader analysis on the 1533 pseudo density.""" 1534 if self.wfs.nspins == 1: 1535 return np.array([self.density.get_correction(a, 0) 1536 for a in range(len(self.atoms))]) 1537 else: 1538 return np.array([[self.density.get_correction(a, spin) 1539 for a in range(len(self.atoms))] 1540 for spin in range(2)]) 1541 1542 def get_all_electron_density(self, spin=None, gridrefinement=2, 1543 pad=True, broadcast=True, collect=True, 1544 skip_core=False): 1545 """Return reconstructed all-electron density array.""" 1546 n_sG, gd = self.density.get_all_electron_density( 1547 self.atoms, gridrefinement=gridrefinement, skip_core=skip_core) 1548 if spin is None: 1549 if self.density.nspins == 1: 1550 n_G = n_sG[0] 1551 else: 1552 n_G = n_sG.sum(axis=0) 1553 else: 1554 if self.density.nspins == 1: 1555 n_G = 0.5 * n_sG[0] 1556 else: 1557 n_G = n_sG[spin] 1558 1559 if collect: 1560 n_G = gd.collect(n_G, broadcast=broadcast) 1561 1562 if n_G is None: 1563 return None 1564 1565 if pad: 1566 n_G = gd.zero_pad(n_G) 1567 1568 return n_G / Bohr**3 1569 1570 def get_fermi_level(self): 1571 """Return the Fermi-level.""" 1572 assert self.wfs.fermi_levels is not None 1573 if len(self.wfs.fermi_levels) != 1: 1574 raise ValueError('There are two Fermi-levels!') 1575 return self.wfs.fermi_levels[0] * Ha 1576 1577 def get_fermi_levels(self): 1578 """Return the Fermi-levels in case of fixed-magmom.""" 1579 assert self.wfs.fermi_levels is not None 1580 if len(self.wfs.fermi_levels) != 2: 1581 raise ValueError('There is only one Fermi-level!') 1582 return self.wfs.fermi_levels * Ha 1583 1584 def get_wigner_seitz_densities(self, spin): 1585 """Get the weight of the spin-density in Wigner-Seitz cells 1586 around each atom. 1587 1588 The density assigned to each atom is relative to the neutral atom, 1589 i.e. the density sums to zero. 1590 """ 1591 from gpaw.analyse.wignerseitz import wignerseitz 1592 atom_index = wignerseitz(self.wfs.gd, self.atoms) 1593 1594 nt_G = self.density.nt_sG[spin] 1595 weight_a = np.empty(len(self.atoms)) 1596 for a in range(len(self.atoms)): 1597 # XXX Optimize! No need to integrate in zero-region 1598 smooth = self.wfs.gd.integrate(np.where(atom_index == a, 1599 nt_G, 0.0)) 1600 correction = self.density.get_correction(a, spin) 1601 weight_a[a] = smooth + correction 1602 1603 return weight_a 1604 1605 def get_dos(self, spin=0, npts=201, width=None): 1606 """The total DOS. 1607 1608 Fold eigenvalues with Gaussians, and put on an energy grid. 1609 1610 returns an (energies, dos) tuple, where energies are relative to the 1611 vacuum level for non-periodic systems, and the average potential for 1612 periodic systems. 1613 """ 1614 if width is None: 1615 width = 0.1 1616 1617 w_k = self.wfs.kd.weight_k 1618 Nb = self.wfs.bd.nbands 1619 energies = np.empty(len(w_k) * Nb) 1620 weights = np.empty(len(w_k) * Nb) 1621 x = 0 1622 for k, w in enumerate(w_k): 1623 energies[x:x + Nb] = self.get_eigenvalues(k, spin) 1624 weights[x:x + Nb] = w 1625 x += Nb 1626 1627 from gpaw.utilities.dos import fold 1628 return fold(energies, weights, npts, width) 1629 1630 def get_wigner_seitz_ldos(self, a, spin=0, npts=201, width=None): 1631 """The Local Density of States, using a Wigner-Seitz basis function. 1632 1633 Project wave functions onto a Wigner-Seitz box at atom ``a``, and 1634 use this as weight when summing the eigenvalues.""" 1635 if width is None: 1636 width = 0.1 1637 1638 from gpaw.utilities.dos import fold, raw_wignerseitz_LDOS 1639 energies, weights = raw_wignerseitz_LDOS(self, a, spin) 1640 return fold(energies * Ha, weights, npts, width) 1641 1642 def get_orbital_ldos(self, a, 1643 spin=0, angular='spdf', npts=201, width=None, 1644 nbands=None, spinorbit=False): 1645 """The Local Density of States, using atomic orbital basis functions. 1646 1647 Project wave functions onto an atom orbital at atom ``a``, and 1648 use this as weight when summing the eigenvalues. 1649 1650 The atomic orbital has angular momentum ``angular``, which can be 1651 's', 'p', 'd', 'f', or any combination (e.g. 'sdf'). 1652 1653 An integer value for ``angular`` can also be used to specify a specific 1654 projector function to project onto. 1655 1656 Setting nbands limits the number of bands included. This speeds up the 1657 calculation if one has many bands in the calculator but is only 1658 interested in the DOS at low energies. 1659 """ 1660 from gpaw.utilities.dos import fold, raw_orbital_LDOS 1661 if width is None: 1662 width = 0.1 1663 1664 if not spinorbit: 1665 energies, weights = raw_orbital_LDOS(self, a, spin, angular, 1666 nbands) 1667 else: 1668 raise DeprecationWarning( 1669 'Please use GPAW.dos(soc=True, ...).raw_pdos(...)') 1670 1671 return fold(energies * Ha, weights, npts, width) 1672 1673 def get_lcao_dos(self, atom_indices=None, basis_indices=None, 1674 npts=201, width=None): 1675 """Get density of states projected onto orbitals in LCAO mode. 1676 1677 basis_indices is a list of indices of basis functions on which 1678 to project. To specify all basis functions on a set of atoms, 1679 you can supply atom_indices instead. Both cannot be given 1680 simultaneously.""" 1681 1682 both_none = atom_indices is None and basis_indices is None 1683 neither_none = atom_indices is not None and basis_indices is not None 1684 if both_none or neither_none: 1685 raise ValueError('Please give either atom_indices or ' 1686 'basis_indices but not both') 1687 1688 if width is None: 1689 width = 0.1 1690 1691 if self.wfs.S_qMM is None: 1692 from gpaw.utilities.dos import RestartLCAODOS 1693 lcaodos = RestartLCAODOS(self) 1694 else: 1695 from gpaw.utilities.dos import LCAODOS 1696 lcaodos = LCAODOS(self) 1697 1698 if atom_indices is not None: 1699 basis_indices = lcaodos.get_atom_indices(atom_indices) 1700 1701 eps_n, w_n = lcaodos.get_subspace_pdos(basis_indices) 1702 from gpaw.utilities.dos import fold 1703 return fold(eps_n * Ha, w_n, npts, width) 1704 1705 def get_all_electron_ldos(self, mol, spin=0, npts=201, width=None, 1706 wf_k=None, P_aui=None, lc=None, raw=False): 1707 """The Projected Density of States, using all-electron wavefunctions. 1708 1709 Projects onto a pseudo_wavefunctions (wf_k) corresponding to some band 1710 n and uses P_aui ([paw.nuclei[a].P_uni[:,n,:] for a in atoms]) to 1711 obtain the all-electron overlaps. 1712 Instead of projecting onto a wavefunction, a molecular orbital can 1713 be specified by a linear combination of weights (lc) 1714 """ 1715 from gpaw.utilities.dos import all_electron_LDOS, fold 1716 1717 if raw: 1718 return all_electron_LDOS(self, mol, spin, lc=lc, 1719 wf_k=wf_k, P_aui=P_aui) 1720 if width is None: 1721 width = 0.1 1722 1723 energies, weights = all_electron_LDOS(self, mol, spin, 1724 lc=lc, wf_k=wf_k, P_aui=P_aui) 1725 return fold(energies * Ha, weights, npts, width) 1726 1727 def get_pseudo_wave_function(self, band=0, kpt=0, spin=0, broadcast=True, 1728 pad=True, periodic=False): 1729 """Return pseudo-wave-function array. 1730 1731 Units: 1/Angstrom^(3/2) 1732 """ 1733 if self.wfs.mode == 'lcao' and not self.wfs.positions_set: 1734 self.initialize_positions() 1735 1736 if pad: 1737 psit_G = self.get_pseudo_wave_function(band, kpt, spin, broadcast, 1738 pad=False, 1739 periodic=periodic) 1740 if psit_G is None: 1741 return 1742 else: 1743 return self.wfs.gd.zero_pad(psit_G) 1744 1745 psit_G = self.wfs.get_wave_function_array(band, kpt, spin, 1746 periodic=periodic) 1747 if broadcast: 1748 if self.wfs.world.rank != 0: 1749 psit_G = self.wfs.gd.empty(dtype=self.wfs.dtype, 1750 global_array=True) 1751 psit_G = np.ascontiguousarray(psit_G) 1752 self.wfs.world.broadcast(psit_G, 0) 1753 return psit_G / Bohr**1.5 1754 elif self.wfs.world.rank == 0: 1755 return psit_G / Bohr**1.5 1756 1757 def get_eigenvalues(self, kpt=0, spin=0, broadcast=True): 1758 """Return eigenvalue array.""" 1759 assert 0 <= kpt < self.wfs.kd.nibzkpts, kpt 1760 eps_n = self.wfs.collect_eigenvalues(kpt, spin) 1761 if broadcast: 1762 if self.wfs.world.rank != 0: 1763 eps_n = np.empty(self.wfs.bd.nbands) 1764 self.wfs.world.broadcast(eps_n, 0) 1765 return eps_n * Ha 1766 1767 def get_occupation_numbers(self, kpt=0, spin=0, broadcast=True): 1768 """Return occupation array.""" 1769 f_n = self.wfs.collect_occupations(kpt, spin) 1770 if broadcast: 1771 if self.wfs.world.rank != 0: 1772 f_n = np.empty(self.wfs.bd.nbands) 1773 self.wfs.world.broadcast(f_n, 0) 1774 return f_n 1775 1776 def get_xc_difference(self, xc): 1777 if isinstance(xc, (str, dict)): 1778 xc = XC(xc) 1779 xc.set_grid_descriptor(self.density.finegd) 1780 xc.initialize(self.density, self.hamiltonian, self.wfs) 1781 xc.set_positions(self.spos_ac) 1782 if xc.orbital_dependent: 1783 self.converge_wave_functions() 1784 return self.hamiltonian.get_xc_difference(xc, self.density) * Ha 1785 1786 def initial_wannier(self, initialwannier, kpointgrid, fixedstates, 1787 edf, spin, nbands): 1788 """Initial guess for the shape of wannier functions. 1789 1790 Use initial guess for wannier orbitals to determine rotation 1791 matrices U and C. 1792 """ 1793 from ase.dft.wannier import rotation_from_projection 1794 proj_knw = self.get_projections(initialwannier, spin) 1795 U_kww = [] 1796 C_kul = [] 1797 for fixed, proj_nw in zip(fixedstates, proj_knw): 1798 U_ww, C_ul = rotation_from_projection(proj_nw[:nbands], 1799 fixed, 1800 ortho=True) 1801 U_kww.append(U_ww) 1802 C_kul.append(C_ul) 1803 1804 U_kww = np.asarray(U_kww) 1805 return C_kul, U_kww 1806 1807 def get_wannier_localization_matrix(self, nbands, dirG, kpoint, 1808 nextkpoint, G_I, spin): 1809 """Calculate integrals for maximally localized Wannier functions.""" 1810 1811 # Due to orthorhombic cells, only one component of dirG is non-zero. 1812 k_kc = self.wfs.kd.bzk_kc 1813 G_c = k_kc[nextkpoint] - k_kc[kpoint] - G_I 1814 1815 return self.get_wannier_integrals(spin, kpoint, 1816 nextkpoint, G_c, nbands) 1817 1818 def get_wannier_integrals(self, s, k, k1, G_c, nbands=None): 1819 """Calculate integrals for maximally localized Wannier functions.""" 1820 1821 assert s <= self.wfs.nspins 1822 kpt_rank, u = divmod(k + len(self.wfs.kd.ibzk_kc) * s, 1823 len(self.wfs.kpt_u)) 1824 kpt_rank1, u1 = divmod(k1 + len(self.wfs.kd.ibzk_kc) * s, 1825 len(self.wfs.kpt_u)) 1826 kpt_u = self.wfs.kpt_u 1827 1828 # XXX not for the kpoint/spin parallel case 1829 assert self.wfs.kd.comm.size == 1 1830 1831 # If calc is a save file, read in tar references to memory 1832 # For lcao mode just initialize the wavefunctions from the 1833 # calculated lcao coefficients 1834 if self.wfs.mode == 'lcao': 1835 self.wfs.initialize_wave_functions_from_lcao() 1836 else: 1837 self.wfs.initialize_wave_functions_from_restart_file() 1838 1839 # Get pseudo part 1840 Z_nn = self.wfs.gd.wannier_matrix(kpt_u[u].psit_nG, 1841 kpt_u[u1].psit_nG, G_c, nbands) 1842 1843 # Add corrections 1844 self.add_wannier_correction(Z_nn, G_c, u, u1, nbands) 1845 1846 self.wfs.gd.comm.sum(Z_nn) 1847 1848 return Z_nn 1849 1850 def add_wannier_correction(self, Z_nn, G_c, u, u1, nbands=None): 1851 r"""Calculate the correction to the wannier integrals. 1852 1853 See: (Eq. 27 ref1):: 1854 1855 -i G.r 1856 Z = <psi | e |psi > 1857 nm n m 1858 1859 __ __ 1860 ~ \ a \ a* a a 1861 Z = Z + ) exp[-i G . R ] ) P dO P 1862 nmx nmx /__ x /__ ni ii' mi' 1863 1864 a ii' 1865 1866 Note that this correction is an approximation that assumes the 1867 exponential varies slowly over the extent of the augmentation sphere. 1868 1869 ref1: Thygesen et al, Phys. Rev. B 72, 125119 (2005) 1870 """ 1871 1872 if nbands is None: 1873 nbands = self.wfs.bd.nbands 1874 1875 P_ani = self.wfs.kpt_u[u].P_ani 1876 P1_ani = self.wfs.kpt_u[u1].P_ani 1877 for a, P_ni in P_ani.items(): 1878 P_ni = P_ani[a][:nbands] 1879 P1_ni = P1_ani[a][:nbands] 1880 dO_ii = self.wfs.setups[a].dO_ii 1881 e = np.exp(-2.j * np.pi * np.dot(G_c, self.spos_ac[a])) 1882 Z_nn += e * np.dot(np.dot(P_ni.conj(), dO_ii), P1_ni.T) 1883 1884 def get_projections(self, locfun, spin=0): 1885 """Project wave functions onto localized functions 1886 1887 Determine the projections of the Kohn-Sham eigenstates 1888 onto specified localized functions of the format:: 1889 1890 locfun = [[spos_c, l, sigma], [...]] 1891 1892 spos_c can be an atom index, or a scaled position vector. l is 1893 the angular momentum, and sigma is the (half-) width of the 1894 radial gaussian. 1895 1896 Return format is:: 1897 1898 f_kni = <psi_kn | f_i> 1899 1900 where psi_kn are the wave functions, and f_i are the specified 1901 localized functions. 1902 1903 As a special case, locfun can be the string 'projectors', in which 1904 case the bound state projectors are used as localized functions. 1905 """ 1906 1907 wfs = self.wfs 1908 1909 if locfun == 'projectors': 1910 f_kin = [] 1911 for kpt in wfs.kpt_u: 1912 if kpt.s == spin: 1913 f_in = [] 1914 for a, P_ni in kpt.P_ani.items(): 1915 i = 0 1916 setup = wfs.setups[a] 1917 for l, n in zip(setup.l_j, setup.n_j): 1918 if n >= 0: 1919 for j in range(i, i + 2 * l + 1): 1920 f_in.append(P_ni[:, j]) 1921 i += 2 * l + 1 1922 f_kin.append(f_in) 1923 f_kni = np.array(f_kin).transpose(0, 2, 1) 1924 return f_kni.conj() 1925 1926 from math import factorial as fac 1927 1928 from gpaw.lfc import LocalizedFunctionsCollection as LFC 1929 from gpaw.spline import Spline 1930 1931 nkpts = len(wfs.kd.ibzk_kc) 1932 nbf = np.sum([2 * l + 1 for pos, l, a in locfun]) 1933 f_kni = np.zeros((nkpts, wfs.bd.nbands, nbf), wfs.dtype) 1934 1935 spos_xc = [] 1936 splines_x = [] 1937 for spos_c, l, sigma in locfun: 1938 if isinstance(spos_c, int): 1939 spos_c = self.spos_ac[spos_c] 1940 spos_xc.append(spos_c) 1941 alpha = .5 * Bohr**2 / sigma**2 1942 r = np.linspace(0, 10. * sigma, 500) 1943 f_g = (fac(l) * (4 * alpha)**(l + 3 / 2.) * 1944 np.exp(-alpha * r**2) / 1945 (np.sqrt(4 * np.pi) * fac(2 * l + 1))) 1946 splines_x.append([Spline(l, rmax=r[-1], f_g=f_g)]) 1947 1948 lf = LFC(wfs.gd, splines_x, wfs.kd, dtype=wfs.dtype) 1949 lf.set_positions(spos_xc) 1950 1951 assert wfs.gd.comm.size == 1 1952 k = 0 1953 f_ani = lf.dict(wfs.bd.nbands) 1954 for kpt in wfs.kpt_u: 1955 if kpt.s != spin: 1956 continue 1957 lf.integrate(kpt.psit_nG[:], f_ani, kpt.q) 1958 i1 = 0 1959 for x, f_ni in f_ani.items(): 1960 i2 = i1 + f_ni.shape[1] 1961 f_kni[k, :, i1:i2] = f_ni 1962 i1 = i2 1963 k += 1 1964 1965 return f_kni.conj() 1966 1967 def get_number_of_grid_points(self): 1968 return self.wfs.gd.N_c 1969 1970 def get_number_of_iterations(self): 1971 return self.scf.niter 1972 1973 def get_number_of_electrons(self): 1974 return self.wfs.setups.nvalence - self.density.charge 1975 1976 def get_electrostatic_corrections(self): 1977 """Calculate PAW correction to average electrostatic potential.""" 1978 dEH_a = np.zeros(len(self.atoms)) 1979 for a, D_sp in self.density.D_asp.items(): 1980 setup = self.wfs.setups[a] 1981 dEH_a[a] = setup.dEH0 + np.dot(setup.dEH_p, D_sp.sum(0)) 1982 self.wfs.gd.comm.sum(dEH_a) 1983 return dEH_a * Ha * Bohr**3 1984 1985 def get_nonselfconsistent_energies(self, type='beefvdw'): 1986 from gpaw.xc.bee import BEEFEnsemble 1987 if type not in ['beefvdw', 'mbeef', 'mbeefvdw']: 1988 raise NotImplementedError('Not implemented for type = %s' % type) 1989 assert self.scf.converged 1990 bee = BEEFEnsemble(self) 1991 x = bee.create_xc_contributions('exch') 1992 c = bee.create_xc_contributions('corr') 1993 if type == 'beefvdw': 1994 return np.append(x, c) 1995 elif type == 'mbeef': 1996 return x.flatten() 1997 elif type == 'mbeefvdw': 1998 return np.append(x.flatten(), c) 1999 2000 def embed(self, q_p, rc=0.2, rc2=np.inf, width=1.0): 2001 """Embed QM region in point-charges.""" 2002 pc = PointChargePotential(q_p, rc=rc, rc2=rc2, width=width) 2003 self.set(external=pc) 2004 return pc 2005 2006 def dos(self, 2007 soc: bool = False, 2008 theta: float = 0.0, 2009 phi: float = 0.0, 2010 shift_fermi_level: bool = True) -> DOSCalculator: 2011 """Create DOS-calculator. 2012 2013 Default is to shift_fermi_level to 0.0 eV. For soc=True, angles 2014 can be given in degrees. 2015 """ 2016 return DOSCalculator.from_calculator( 2017 self, soc=soc, 2018 theta=theta, phi=phi, 2019 shift_fermi_level=shift_fermi_level) 2020