from __future__ import annotations
import warnings
from collections.abc import Callable
from pathlib import Path
from typing import IO, TYPE_CHECKING, Any, Sequence
import numpy as np
from ase import Atoms
from ase.geometry import cell_to_cellpar
from ase.units import Bohr, Ha
from gpaw import GPAW_NO_C_EXTENSION
from gpaw.core import UGArray, UGDesc
from gpaw.core.atom_arrays import AtomDistribution
from gpaw.densities import Densities
from gpaw.electrostatic_potential import ElectrostaticPotential
from gpaw.gpu import as_np
from gpaw.gpu import cupy as cp
from gpaw.mpi import MPIComm
from gpaw.mpi import broadcast as bcast
from gpaw.mpi import broadcast_float, receive, send, serial_comm
from gpaw.new import trace
from gpaw.new.density import Density
from gpaw.new.energies import DFTEnergies
from gpaw.new.gpw import read_gpw
from gpaw.new.ibzwfs import IBZWaveFunctions
from gpaw.new.logger import Logger, indent
from gpaw.new.potential import Potential
from gpaw.new.scf import SCFLoop
from gpaw.setup import Setups
from gpaw.typing import Array1D, Array2D
from gpaw.utilities import (check_atoms_too_close,
check_atoms_too_close_to_boundary)
from gpaw.utilities.timing import simpletimer
if TYPE_CHECKING:
from gpaw.dft import (XC, Eigensolver, ExtensionInput, KptsType, Mixer,
Mode, MonkhorstPack, Occupations, Parameters,
PoissonSolver, Symmetry)
class ReuseWaveFunctionsError(Exception):
"""Reusing the old wave functions after cell-change failed.
Most likekly, the number of k-points changed.
"""
class NonsenseError(Exception):
"""Operation doesn't make sense."""
class CalculationModeError(Exception):
"""Calculation mode does not match what is expected from a given method.
For example, if a method only works in collinear mode and receives a
calculator in non-collinear mode, this exception should be raised.
"""
units = {'energy': Ha,
'free_energy': Ha,
'forces': Ha / Bohr,
'stress': Ha / Bohr**3,
'dipole': Bohr,
'magmom': 1.0,
'magmoms': 1.0,
'non_collinear_magmom': 1.0,
'non_collinear_magmoms': 1.0}
[docs]
class DFT:
def __init__(
self,
atoms: Atoms,
*,
mode: str | dict | Mode,
basis: str | dict[str | int | None, str] | None = None,
charge: float | None = None,
convergence: dict | None = None,
eigensolver: str | dict | Eigensolver | None = None,
experimental: dict | None = None,
extensions: Sequence[ExtensionInput] | None = None,
gpts: Sequence[int] | None = None,
h: float | None = None,
hund: bool | None = None,
interpolation: int | None = None,
kpts: KptsType | MonkhorstPack | None = None,
magmoms: Sequence[float] | Sequence[Sequence[float]] | None = None,
maxiter: int | None = None,
mixer: dict | Mixer | None = None,
nbands: int | str | None = None,
occupations: dict | Occupations | None = None,
parallel: dict | None = None,
poissonsolver: dict | PoissonSolver | None = None,
random: bool | None = None,
setups: str | dict | None = None,
soc: bool | None = None,
spinpol: bool | None = None,
symmetry: str | dict | Symmetry | None = None,
xc: str | dict | XC | None = None,
txt: str | Path | IO[str] | Logger | None = '-',
communicator: MPIComm | None = None,
converge=True,
_components=None,
_parameters=None):
"""Create a DFT object.
See :class:`gpaw.dft.Parameters` for the complete list of parameters.
Parameters
==========
atoms:
ASE-Atoms object.
txt:
Text log-file. Use ``None`` for no logging and ``'-'`` for using
standard out.
communicator:
MPI-communicator. Default is to use ``gpaw.mpi.world``.
"""
from gpaw.dft import PARAMETER_NAMES, Parameters
self.atoms = atoms.copy()
if _components is None:
if _parameters is None:
kwargs = {k: v for k, v in locals().items()
if k in PARAMETER_NAMES}
params = Parameters(**kwargs)
else:
params = _parameters
_components = params.create_dft_calculation_components(
self.atoms, communicator, txt)
(self.ibzwfs, self.density, self.potential,
self.setups, self.scf_loop, self.pot_calc,
self.log, self.params, energies) = _components
self.comm = self.ibzwfs.comm
self.results: dict[str, Any] = {}
self.relpos_ac = self.pot_calc.relpos_ac
self.energies = energies or DFTEnergies()
self.forces_have_been_printed = False
if converge:
self.converge()
[docs]
@classmethod
def from_parameters(cls,
atoms: Atoms,
params: Parameters,
comm,
log: Logger | str | None = None,
converge=True):
return cls(
atoms,
mode='',
communicator=comm,
txt=log,
converge=converge,
_parameters=params)
[docs]
@classmethod
def from_components(cls,
atoms: Atoms,
ibzwfs: IBZWaveFunctions,
density: Density,
potential: Potential,
setups: Setups,
scf_loop: SCFLoop,
pot_calc,
log: Logger,
params: Parameters,
energies: DFTEnergies | None = None,
converge=True):
return cls(
atoms,
mode='',
_components=(ibzwfs, density, potential,
setups, scf_loop, pot_calc,
log, params, energies),
converge=converge)
def __getattr__(self, name):
matches = [ext
for ext in self.pot_calc.extensions
if ext.name == name]
if len(matches) != 1:
raise AttributeError
return matches[0]
[docs]
def ase_calculator(self):
"""Create ASE-compatible GPAW calculator.
"""
from gpaw.new.ase_interface import ASECalculator
return ASECalculator(params=self.params,
log=self.log,
dft=self,
atoms=self.atoms)
[docs]
def move_atoms(self, atoms) -> DFT:
check_atoms_too_close(atoms)
self.atoms = atoms
self.relpos_ac = np.ascontiguousarray(atoms.get_scaled_positions())
self.comm.broadcast(self.relpos_ac, 0)
atomdist = self.density.D_asii.layout.atomdist
grid = self.density.nt_sR.desc
rank_a = grid.ranks_from_fractional_positions(self.relpos_ac)
atomdist = AtomDistribution(rank_a, atomdist.comm)
self.pot_calc.move(self.relpos_ac, atomdist)
self.ibzwfs.move(self.relpos_ac, atomdist)
self.density.move(self.relpos_ac, atomdist)
if self.ibzwfs.has_wave_functions():
self.density.update(self.ibzwfs)
self.potential.move(atomdist)
self.scf_loop.hamiltonian.move(self.relpos_ac)
if self.params.experimental.get('paw_corr_mixer', False):
for basemixer in self.scf_loop.mixer.basemixers:
basemixer.dotprod.atomdist = atomdist
self.potential, self.energies, _ = self.pot_calc.calculate(
self.density, self.ibzwfs, self.potential.vHt_x)
mm_av = self.results.get('non_collinear_magmoms')
if mm_av is None:
_, mm_av = self.density.calculate_magnetic_moments()
write_atoms(atoms, mm_av, self.density.nt_sR.desc, self.log)
self.results = {}
self.forces_have_been_printed = False
return self
[docs]
@trace
def iconverge(self, maxiter=None, calculate_forces=None):
if calculate_forces is None:
calculate_forces = self._calculate_forces
self.ibzwfs.make_sure_wfs_are_read_from_gpw_file()
if self.params.mode.name == 'pw':
if self.params.mode.dtype == 'single':
if not (GPAW_NO_C_EXTENSION or self.ibzwfs.xp is cp):
raise NotImplementedError(
'Single precision is GPU or pure python only.')
for ctx in self.scf_loop.iterate(self.ibzwfs,
self.density,
self.potential,
self.energies,
self.pot_calc,
maxiter=maxiter,
calculate_forces=calculate_forces,
log=self.log):
self.energies = ctx.energies
self.potential = ctx.potential
yield ctx
[docs]
@trace
def converge(self,
maxiter=None,
steps=99999999999999999):
"""Converge to self-consistent solution of Kohn-Sham equation."""
for step, _ in enumerate(self.iconverge(maxiter),
start=1):
if step == steps:
break
else: # no break
self.log('SCF steps:', step)
[docs]
def calculate_energy(self) -> float:
"""Calculate total energy in eV."""
self.results['free_energy'] = broadcast_float(
self.energies.total_free, self.comm)
self.results['energy'] = broadcast_float(
self.energies.total_extrapolated, self.comm)
self.log('Energy contributions relative to reference atoms:',
f'(reference = {self.setups.Eref * Ha:.6f})\n')
self.energies.summary(self.log)
return self.results['energy'] * units['energy']
[docs]
def calculate_dipole(self):
"""Calculate dipole moment in Angstrom."""
if 'dipole' in self.results:
return self.results['dipole'] * units['dipole']
dipole_v = self.density.calculate_dipole_moment(self.relpos_ac)
x, y, z = dipole_v * Bohr
self.log(f'dipole moment: ({x:.6f}, {y:.6f}, {z:.6f}) # |e|*Ang\n')
self.results['dipole'] = dipole_v
return self.results['dipole'] * units['dipole']
[docs]
def magmoms(self) -> tuple[Array1D, Array2D]:
mm_v, mm_av = self.density.calculate_magnetic_moments()
self.results['magmom'] = mm_v[2]
self.results['magmoms'] = mm_av[:, 2].copy()
self.results['non_collinear_magmoms'] = mm_av
self.results['non_collinear_magmom'] = mm_v
if self.density.ncomponents > 1:
x, y, z = mm_v
self.log(f'Total magnetic moment: ({x:.6f}, {y:.6f}, {z:.6f})\n')
self.log.table(
'Local magnetic moments',
header=['x', 'y', 'z'],
rows=[[f'{mm:.6f}' for mm in mm_v] for mm_v in mm_av],
comments=[f'{setup.symbol:2} {a:3}'
for a, setup in enumerate(self.setups)])
self.log()
return mm_v, mm_av
[docs]
def calculate_forces(self):
"""Calculate atomic forces in eV / Angstrom."""
forcetimer = simpletimer()
if 'forces' not in self.results:
self._calculate_forces()
if self.forces_have_been_printed:
return self.results['forces'] * units['forces']
self.forces_have_been_printed = True
F_av = self.results['forces'] * units['forces']
self.log.table(
'\nForces',
comment='eV/Ang',
header=['x', 'y', 'z'],
rows=[[f'{f:.5f}' for f in F_v] for F_v in F_av],
comments=[f'{setup.symbol:2} {a:3}'
for a, setup in enumerate(self.setups)])
self.log(f'\nForce computation time: {forcetimer():.3f} # s')
self.log.fd.flush()
return F_av
def _calculate_forces(self):
xc = self.pot_calc.xc
assert not hasattr(xc.xc, 'setup_force_corrections')
# Force from projector functions (and basis set, hybrids):
F_av = self.ibzwfs.forces(self.potential, self.scf_loop.hamiltonian,
self.density.D_asii)
if xc.name == 'QNA':
getattr(xc.xc, 'add_forces', lambda F_av: None)(F_av) # QNA
pot_calc = self.pot_calc
Q_aL = self.density.calculate_compensation_charge_coefficients()
Fcc_av, Fnct_av, Ftauct_av, Fvbar_av, Fext_av = \
pot_calc.force_contributions(Q_aL,
self.density,
self.potential)
# Force from compensation charges:
F_av += Fcc_av
# Force from smooth core charge:
for a, dF_v in Fnct_av.items():
F_av[a] += dF_v[:, 0]
if Ftauct_av is not None:
# Force from smooth core ked:
for a, dF_v in Ftauct_av.items():
F_av[a] += dF_v[:, 0]
# Force from zero potential:
for a, dF_v in Fvbar_av.items():
F_av[a] += dF_v[:, 0]
F_av = as_np(F_av)
domain_comm = Q_aL.layout.atomdist.comm
domain_comm.sum(F_av)
# Force from extensions (only from rank 0)
F_av += Fext_av
F_av = self.ibzwfs.ibz.symmetries.symmetrize_forces(F_av)
self.comm.broadcast(F_av, 0)
self.results['forces'] = F_av
return F_av
[docs]
def calculate_stress(self) -> None:
"""Calculate stress in eV / Angstrom^3."""
stresstimer = simpletimer()
if 'stress' in self.results:
return self.results['stress'] * units['stress']
stress_vv = self.pot_calc.stress(
self.ibzwfs, self.density, self.potential)
self.log.table(
'\nStress tensor',
comment='eV/Ang^3',
header=list('xyz'),
rows=[[f'{s:.6f}' for s in s_v]
for s_v in stress_vv * units['stress']])
self.results['stress'] = stress_vv.flat[[0, 4, 8, 5, 2, 1]]
self.log(f'\nStress computation time: {stresstimer():.3f} # s\n')
return self.results['stress'] * units['stress']
[docs]
def write_converged(self) -> None:
self.ibzwfs.write_summary(self.log)
vacuum_level = self.potential.get_vacuum_level()
if not np.isnan(vacuum_level):
self.log(f'vacuum-level: {vacuum_level:.3f} # V')
try:
wf1, wf2 = self.workfunctions(_vacuum_level=vacuum_level)
except NonsenseError:
pass
else:
self.log(
f'Workfunctions: {wf1:.3f}, {wf2:.3f} eV (top, bottom)')
self.log.fd.flush()
[docs]
def workfunctions(self,
*, _vacuum_level=None) -> tuple[float, float]:
if _vacuum_level is None:
_vacuum_level = self.potential.get_vacuum_level()
if np.isnan(_vacuum_level):
raise NonsenseError('No vacuum')
try:
correction = self.pot_calc.poisson_solver.dipole_layer_correction()
except NotImplementedError:
raise NonsenseError('No dipole layer')
correction *= Ha
fermi_level = self.ibzwfs.fermi_level * Ha
wf = _vacuum_level - fermi_level
return wf - correction, wf + correction
[docs]
def vacuum_level(self) -> float:
return self.potential.get_vacuum_level()
[docs]
def electrostatic_potential(self) -> ElectrostaticPotential:
return ElectrostaticPotential.from_calculation(self)
[docs]
def densities(self) -> Densities:
return Densities.from_calculation(self)
[docs]
def wave_function(self, band: int, kpt=0, spin=None,
periodic=False,
broadcast=True) -> UGArray | None:
psit_nR = self.wave_functions(n1=band, n2=band + 1, kpt=kpt, spin=spin,
periodic=periodic, broadcast=broadcast)
if psit_nR is not None:
return psit_nR[0]
return None
[docs]
def wave_functions(self, n1=0, n2=None, kpt=0, spin=None,
periodic=False,
broadcast=True,
_pad=True) -> UGArray | None:
ibzwfs = self.ibzwfs
collinear = ibzwfs.collinear
if collinear:
if spin is None:
spin = 0
else:
assert spin is None or spin == 0
spin = 0
kpt_comm = ibzwfs.kpt_comm
krank = ibzwfs.rank_ks[kpt][spin]
if krank == kpt_comm.rank:
wfs = ibzwfs._get_wfs(kpt, spin)
wfs = wfs.collect_bands(n1, n2)
if wfs is not None:
basis = getattr(self.scf_loop.hamiltonian, 'basis', None)
grid = self.density.nt_sR.desc
if collinear:
wfs = wfs.to_uniform_grid_wave_functions(grid, basis)
psit_nR = wfs.psit_nX
else:
psit_nsG = wfs.psit_nX
grid = grid.new(kpt=psit_nsG.desc.kpt_c,
dtype=psit_nsG.desc.dtype)
psit_nR = psit_nsG.ifft(grid=grid)
if not psit_nR.desc.pbc.all() and _pad:
psit_nR = psit_nR.to_pbc_grid()
if periodic:
psit_nR.multiply_by_eikr(-psit_nR.desc.kpt_c)
psit_nR = psit_nR.gather()
if krank != 0 and psit_nR is not None:
send(psit_nR, 0, kpt_comm)
psit_nR = None
else:
psit_nR = None
elif self.comm.rank == 0:
psit_nR = receive(krank, kpt_comm)
else:
psit_nR = None
if broadcast:
psit_nR = bcast(psit_nR, 0, comm=self.comm)
if psit_nR is None:
return None
return psit_nR.scaled(cell=Bohr, values=Bohr**-1.5)
[docs]
def gather(self, txt='-') -> DFT | None:
"""Gather calculation data from DFT object
on master and return new DFT object
(only on master, None everywhere else)."""
atoms = self.atoms
params = self.params
mode = params.mode
comm = self.comm
# gather data on master
ibz, ncomponents, wfs_u, fermi_levels = self.ibzwfs.gather()
dH_asp, vt_sR, dedtaut_sR, vHt_x = self.potential.gather()
D_asp, nt_sR, taut_sR = self.density.gather()
# only create new dft object on master
if comm.rank != 0:
return None
params.parallel = {'kpt': 1, 'band': 1, 'domain': 1}
builder = params.dft_component_builder(atoms, log=txt,
comm=serial_comm)
# make new wfs on master
if mode == 'lcao':
from gpaw.new.lcao.ibzwfs import LCAOIBZWaveFunctions
IBZWFs = LCAOIBZWaveFunctions
else:
from gpaw.new.pwfd.ibzwfs import PWFDIBZWaveFunctions
IBZWFs = PWFDIBZWaveFunctions # type: ignore
ibzwfs = IBZWFs(
ibz=ibz,
ncomponents=ncomponents,
wfs_u=wfs_u,
kpt_comm=serial_comm,
kpt_band_comm=serial_comm,
comm=serial_comm)
ibzwfs.fermi_levels = fermi_levels
potential = Potential(vt_sR=vt_sR, dH_asii=dH_asp.to_full(),
dedtaut_sR=dedtaut_sR, vHt_x=vHt_x,
e_stress=self.potential.e_stress)
density = Density.from_data_and_setups(
nt_sR, taut_sR, D_asp.to_full(),
builder.params.charge,
builder.setups,
builder.get_pseudo_core_densities(),
builder.get_pseudo_core_ked())
dft = DFT.from_components(
atoms, ibzwfs, density, potential,
builder.setups,
builder.create_scf_loop(),
builder.create_potential_calculator(),
builder.log,
params=params,
energies=self.energies,
converge=False)
dft.results = self.results.copy()
return dft
[docs]
def change(self,
*,
xc=None,
eigensolver=None,
mixer=None,
occupations=None,
convergence=None) -> None:
from gpaw.dft import XC, Eigensolver, Mixer, Occupations
# build kwargs
allargs = {'xc': xc, 'eigensolver': eigensolver,
'mixer': mixer, 'occupations': occupations,
'convergence': convergence}
kwargs = {key: val for key, val in allargs.items() if val is not None}
atoms = self.atoms
params = self.params
log = self.log
prop_update = {'xc': XC.from_param,
'eigensolver': Eigensolver.from_param,
'mixer': Mixer.from_param,
'occupations': Occupations.from_param}
# update params
for prop in kwargs:
val = kwargs[prop]
update_func = prop_update.get(prop, None)
# actually change property
if update_func is None:
# update dictionary
getattr(params, prop).update(val)
else:
setattr(params, prop, update_func(val))
# rebuild
builder = params.dft_component_builder(atoms, log=log,
comm=self.comm)
scf_loop = builder.create_scf_loop()
pot_calc = builder.create_potential_calculator()
xcfunc_n = pot_calc.xc
# checks compatibility
if 'xc' in kwargs:
# check that 'base' functional is the same
xcfunc_o = self.pot_calc.xc
assert xcfunc_n.get_setup_name() == xcfunc_o.get_setup_name()
self.scf_loop = scf_loop
self.pot_calc = pot_calc
if 'occupations' in kwargs:
dens = self.density
wfs = self.ibzwfs
# update occupations numbers
nelectrons = dens.nvalence - dens.charge + self.pot_calc.charge
wfs.calculate_occs(self.scf_loop.occ_calc, nelectrons,
fix_fermi_level=self.scf_loop.fix_fermi_level)
prop_log = {'xc': xcfunc_n.name,
'eigensolver': self.scf_loop.eigensolver,
'mixer': self.scf_loop.mixer,
'occupations': self.scf_loop.occ_calc,
'convergence': allargs['convergence']}
# assemble log
for prop in kwargs:
log(f'Changed {prop}:')
log(indent(prop_log[prop]))
log('Reusing wavefunctions.')
self.results = {}
[docs]
def new(self,
atoms: Atoms,
params: Parameters,
log=None) -> DFT:
"""Create new DFT object."""
if params.mode.name != 'pw':
raise ReuseWaveFunctionsError
ibzwfs = self.ibzwfs
# if ibzwfs.domain_comm.size != 1:
# raise ReuseWaveFunctionsError
check_atoms_too_close(atoms)
check_atoms_too_close_to_boundary(atoms)
builder = params.dft_component_builder(atoms, log=log,
comm=self.comm)
kpt_kc = builder.ibz.kpt_kc
old_kpt_kc = ibzwfs.ibz.kpt_kc
if len(kpt_kc) != len(old_kpt_kc):
raise ReuseWaveFunctionsError
if abs(kpt_kc - old_kpt_kc).max() > 1e-9:
raise ReuseWaveFunctionsError
log('Interpolating wave functions to new cell')
scf_loop = builder.create_scf_loop()
pot_calc = builder.create_potential_calculator()
density = self.density.new(builder.grid,
builder.interpolation_desc,
builder.relpos_ac,
builder.atomdist)
density.normalize(pot_calc.charge)
# Make sure all have exactly the same density.
# Not quite sure it is needed???
# At the moment we skip it on GPU's because it doesn't
# work!
if density.nt_sR.xp is np:
ibzwfs.kpt_band_comm.broadcast(density.nt_sR.data, 0)
potential, energies, _ = pot_calc.calculate(density)
ibzwfs.make_sure_wfs_are_read_from_gpw_file()
old_ibzwfs = ibzwfs
def create_wfs(spin, q, k, kpt_c, weight):
wfs = old_ibzwfs._get_wfs(k, spin)
return wfs.morph(
builder.wf_desc,
builder.relpos_ac,
builder.atomdist)
ibzwfs = ibzwfs.create(
ibz=builder.ibz,
ncomponents=old_ibzwfs.ncomponents,
create_wfs_func=create_wfs,
kpt_comm=old_ibzwfs.kpt_comm,
kpt_band_comm=old_ibzwfs.kpt_band_comm,
comm=self.comm)
write_atoms(atoms, builder.initial_magmom_av, builder.grid, log)
ibzwfs.summary(log)
log(density)
log(potential)
log(builder.setups)
log(scf_loop)
log(pot_calc)
return DFT.from_components(
atoms, ibzwfs, density, potential,
builder.setups, scf_loop, pot_calc, log,
params=params, energies=energies,
converge=False)
[docs]
def change_mode(self,
mode: str | dict | Mode,
*,
nbands: int | None = None) -> None:
"""In-place convertion from one mode to another.
**Only LCAO to PW or FD mode implemented!**
"""
from gpaw.dft import Mode
from gpaw.new.lcao.ibzwfs import LCAOIBZWaveFunctions
if not isinstance(self.ibzwfs, LCAOIBZWaveFunctions):
raise ValueError
self.params.mode = Mode.from_param(mode)
builder = self.params.dft_component_builder(
self.atoms, log=self.log, comm=self.comm)
self.scf_loop = builder.create_scf_loop()
self.pot_calc = builder.create_potential_calculator()
if builder.mode == 'pw':
self.density.nct_aX = builder.get_pseudo_core_densities()
self.density.tauct_aX = builder.get_pseudo_core_ked()
self.density = self.density.new(builder.grid,
builder.interpolation_desc)
self.density.normalize(self.pot_calc.charge)
if self.density.nt_sR.xp is np:
self.ibzwfs.kpt_band_comm.broadcast(self.density.nt_sR.data, 0)
self.potential, self.energies, _ = self.pot_calc.calculate(
self.density)
self.ibzwfs = self.ibzwfs.convert_to(
builder.mode,
grid=builder.grid,
pw=builder.wf_desc,
nbands=nbands)
self.results = {}
[docs]
def write_gpw_file(self, filename,
*,
include_wfs: bool = False,
precision: str = 'double',
include_projections: bool = True):
"""Write calculator object to a file.
Parameters
----------
filename:
File to be written.
include_wfs:
Use ``include_wfs=True``
to include wave functions in the file.
precision:
'double' (the default) or 'single'.
include_projections:
Use ``include_projections=False`` to not include
the PAW-projections.
"""
mode = ''
if include_wfs:
mode = 'all'
self.log(f'Writing to {filename} (mode={mode!r})\n')
self.ase_calculator().write(filename,
mode=mode,
precision=precision,
include_projections=include_projections)
self.comm.barrier()
[docs]
@classmethod
def from_gpw_file(cls,
filename: str | Path | IO[str],
log: Logger | str | Path | IO[str] | None = None,
comm=None,
parallel: dict[str, Any] = None,
dtype=None,
force_complex_dtype: bool = False,
object_hooks: dict[str,
Callable[[dict], Any]] | None = None
) -> DFT:
(atoms,
dft,
builder) = read_gpw(filename=filename,
log=log,
comm=comm,
parallel=parallel,
dtype=dtype,
force_complex_dtype=force_complex_dtype,
object_hooks=object_hooks)
return dft
[docs]
def get_state(self):
return DFTState(self.ibzwfs, self.density, self.potential)
@property
def state(self):
warnings.warn('Use of deprecated DFT.state attribute. '
'Use ibzwfs, density and potential attributes instead.')
return self.get_state()
def write_atoms(atoms: Atoms,
magmom_av: Array2D,
grid: UGDesc,
log) -> None:
symbols = atoms.get_chemical_symbols()
rows = []
for a, (x, y, z) in enumerate(atoms.positions):
symbol = symbols[a]
X, Y, Z = magmom_av[a]
rows.append([f'{a}',
symbol,
f'{x:.6f}', f'{y:.6f}', f'{z:.6f}',
f'({X:7.4f}, {Y:7.4f}, {Z:7.4f})'])
log.table(
'\nAtoms',
comment='Å, Bohr magnetons',
header=',symbol,x,y,z,initial magnetic moments'.split(','),
rows=rows)
par = cell_to_cellpar(atoms.cell)
log.table(
'\nUnit cell',
comment='Å',
header='periodic,x,y,z,lengths,angles'.split(','),
rows=[['yes' if p else ' no',
f'{x:.6f}', f'{y:.6f}', f'{z:.6f}',
f'{L:.6f}', f'{A:.6f}']
for p, (x, y, z), L, A in zip(atoms.pbc,
atoms.cell,
par[:3],
par[3:])])
log()
[docs]
class DFTState:
def __init__(self,
ibzwfs: IBZWaveFunctions,
density: Density,
potential: Potential,
energies):
"""State of a Kohn-Sham calculation."""
self.ibzwfs = ibzwfs
self.density = density
self.potential = potential
self.energies = energies