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