From f550e152b66cb6c06a888216d98bf6157faabefb Mon Sep 17 00:00:00 2001 From: Justin Bich Date: Fri, 24 Jun 2022 17:22:01 +0200 Subject: [PATCH] Implement a calculator based on dftbplus_ptools --- setup.cfg | 1 + src/dftbplus/calculator.py | 276 +++++++++++++ src/dftbplus/dftbplus.py | 629 ------------------------------ src/dftbplus/io.py | 130 ++++++ src/dftbplus/legacy_calculator.py | 435 +++++++++++++++++++++ 5 files changed, 842 insertions(+), 629 deletions(-) create mode 100644 src/dftbplus/calculator.py delete mode 100644 src/dftbplus/dftbplus.py create mode 100644 src/dftbplus/io.py create mode 100644 src/dftbplus/legacy_calculator.py diff --git a/setup.cfg b/setup.cfg index 97351a0..6ccff6c 100644 --- a/setup.cfg +++ b/setup.cfg @@ -19,6 +19,7 @@ install_requires = numpy ase hsd + dftbplus-ptools tests_require = pytest diff --git a/src/dftbplus/calculator.py b/src/dftbplus/calculator.py new file mode 100644 index 0000000..01f2cd4 --- /dev/null +++ b/src/dftbplus/calculator.py @@ -0,0 +1,276 @@ +#------------------------------------------------------------------------------# +# DFTB+: general package for performing fast atomistic simulations # +# Copyright (C) 2006 - 2022 DFTB+ developers group # +# # +# See the LICENSE file for terms of usage and distribution. # +#------------------------------------------------------------------------------# + +"""This module defines an ASE interface to DFTB+, based on FileIO.""" + + +import os +import copy +import warnings + +import dftbplus_ptools.resultstag +from dftbplus_ptools.hsdinput import Hsdinput +from dftbplus_ase.io import atom_to_geometry +from ase.io import write +from ase.units import Hartree, Bohr +from ase.calculators.calculator import FileIOCalculator + + +class DftbPlus(FileIOCalculator): + """ASE interface to DFTB+. + + A DFTB+ calculator with ase-FileIOCalculator nomenclature. + """ + + implemented_properties = ('energy', 'forces', 'charges', 'stress', + 'occupations', 'fermi_levels', 'fermi_level', + 'eigenvalues', 'dipole') + + + def __init__(self, inp, binary='dftb+', + restart=None, label='dftbplus', atoms=None, **kwargs): + """Initializes a ctypes DFTB+ calculator object. + + Args: + inp (dict): dictionary containing (additional) + initialization options for DFTB+ + binary (str): location of dftb+ binary + restart (str): Prefix for restart file. May contain a directory. + label (str): Name used for all files. + """ + + self._hsdinput = Hsdinput(dictionary=copy.deepcopy(inp)) + self._hamiltonian = self._hsdinput.get_hamiltonian() + self._do_forces = False + self._atoms = None + self._atoms_input = None + + FileIOCalculator.__init__(self, restart=restart, label=label, + atoms=atoms, command=binary, **kwargs) + + # Determine number of spin channels + try: + spinpolkeys = list(self._hsdinput['hamiltonian']\ + [self._hamiltonian]['spinpolarisation'].keys()) + self.nspin = 2 if 'colinear' in spinpolkeys else 1 + except KeyError: + self.nspin = 1 + + + def _write_dftb_in(self, filename): + """Dumps Python input dictionary to HSD format and writes input file. + If necessary, geometry-relevant entries in the input are updated. + If not all max. angular momenta have been specified by user, an + attempt is made to extract remaining ones out of Slater-Koster + files. + + Args: + filename (str): name of input file to be written + + Returns: + (file): containing dictionary in HSD format + """ + geometry = atom_to_geometry(self._atoms_input) + + specieslist = [] + for species in geometry.specienames: + specieslist.append(species.lower()) + + speciesdiff = set(specieslist) - \ + set(self._hsdinput['hamiltonian'][self._hamiltonian] + ['maxangularmomentum'].keys()) + maxangs = self._hsdinput['hamiltonian'][self._hamiltonian]\ + ['maxangularmomentum'] + + path = os.path.join(self.directory, filename) + self._hsdinput.set_filename(path) + self._hsdinput.set_geometry(geometry) + + if speciesdiff and self._hamiltonian == "dftb": + self._hsdinput.set_maxang(maxangs=maxangs, try_reading=speciesdiff) + elif speciesdiff and self._hamiltonian != "dftb": + raise ValueError(f"max. angular momanta for '{speciesdiff}' " + + "missing!") + else: + self._hsdinput.set_maxang(maxangs=maxangs) + + self._hsdinput.write_resultstag() + warnings.warn("Change in input: Because the calculator requires " + + "'WriteResultsTag' to be set to 'Yes', it was set to " + + "'Yes' regardless of its initial value.") + + self._hsdinput.write_hsd() + + + def set(self, **kwargs): + """sets and changes parameters + + Returns: + changed_parameters (dict): containing changed parameters + """ + changed_parameters = FileIOCalculator.set(self, **kwargs) + if changed_parameters: + self.reset() + return changed_parameters + + + def check_state(self, atoms, tol=1e-15): + """Check for any system changes since last calculation. + + Args: + atoms (atoms object): atom to be compared to + tol (float): tolerance for comparison + + Returns: + system_changes (list): containing changed properties + """ + system_changes = FileIOCalculator.check_state(self, atoms, tol=1e-15) + # Ignore unit cell for molecules: + if not atoms.pbc.any() and 'cell' in system_changes: + system_changes.remove('cell') + return system_changes + + + def write_input(self, atoms, properties=None, system_changes=None): + """Writes input files + + Args: + atoms (atoms object): atom to read from + properties (list): list of what needs to be calculated + system_changes (list): List of what has changed since last + calculation + """ + if properties is not None: + if 'forces' in properties or 'stress' in properties: + self._do_forces = True + self._hsdinput.calc_forces() + if 'charges' in properties: + self._hsdinput.calc_charges() + + FileIOCalculator.write_input(self, atoms, properties, system_changes) + + # self._atoms is none until results are read out, + # then it is set to the ones at writing input + self._atoms_input = atoms + self._atoms = None + + self._write_dftb_in('dftb_in.hsd') + + write(os.path.join(self.directory, 'geo_end.gen'), atoms) + + + def read_results(self): + """All results are read from results.tag file. + """ + + self._atoms = self._atoms_input + + path = os.path.join(self.directory, "results.tag") + resultstag = dftbplus_ptools.resultstag.Output(path) + + if resultstag.get_gross_atomic_charges() is not None: + self.results['charges'] = resultstag.get_gross_atomic_charges() + self.results['energy'] = resultstag.get_total_energy() * Hartree + + if resultstag.get_fermi_level() is not None: + self.results['fermi_levels'] = resultstag.get_fermi_level() \ + * Hartree + self.results['fermi_level'] = max(resultstag.get_fermi_level() + * Hartree) + + if resultstag.get_eigenvalues() is not None: + self.results['eigenvalues'] = resultstag.get_eigenvalues() \ + * Hartree + + if self._do_forces: + self.results['forces'] = resultstag.get_forces() * Hartree / Bohr + + if resultstag.get_stress() is not None: + self.results['stress'] = - resultstag.get_stress() * Hartree \ + / Bohr**3 + + if resultstag.get_filling() is not None: + self.results['occupations'] = resultstag.get_filling() + + if resultstag.get_dipole_moments() is not None: + self.results['dipole'] = resultstag.get_dipole_moments() * Bohr + os.remove(os.path.join(self.directory, 'results.tag')) + + + def get_charges(self, atoms): + """Get the calculated charges. This is inhereted to atoms object. (In + get_charges of atoms object: "return self._calc.get_charges(self)") + + Args: + atoms (atoms object): atom to read charges from + """ + + return self.get_property('charges', allow_calculation=True).copy() + + + def get_number_of_spins(self): + """Auxiliary function to extract result: number of spins. + + Returns: + self.nspin (int): number of spins + """ + + return self.nspin + + + def get_eigenvalues(self, kpt=0, spin=0): + """Auxiliary function to extract result: eigenvalues. + + Args: + kpt (int): kpt where eigenvalues are desired + spin (int): spin where eigenvalues are desired + + Returns: + (array): desired eigenvalues + """ + eigenval = self.get_property('eigenvalues', allow_calculation=False) + return eigenval[spin][kpt].copy() + + + def get_fermi_levels(self): + """Auxiliary function to extract result: Fermi level(s). + + Returns: + (array): containing fermi levels + """ + + return self.get_property('fermi_levels', + allow_calculation=False).copy() + + + def get_fermi_level(self): + """Auxiliary function to extract result: result: max. Fermi level. + + Returns: + (np.float64): max. Fermi level + """ + + return self.get_property('fermi_level', allow_calculation=False) + + + def get_occupation_numbers(self, kpt=0, spin=0): + """Auxiliary function to extract result: occupations. + + Args: + kpt (int): kpt where eigenvalues are desired + spin (int): spin where eigenvalues are desired + + Returns: + occs (array): desired occupations + """ + + occs = self.get_property('occupations', allow_calculation=False).copy() + return occs[kpt, spin] + + def get_dipole_moment(self, atoms=None): + return self.get_property('dipole', atoms, + allow_calculation=False).copy() diff --git a/src/dftbplus/dftbplus.py b/src/dftbplus/dftbplus.py deleted file mode 100644 index 03b7cab..0000000 --- a/src/dftbplus/dftbplus.py +++ /dev/null @@ -1,629 +0,0 @@ -#------------------------------------------------------------------------------# -# dftbplus-ase: Interfacing DFTB+ with the Atomic Simulation Environment # -# Copyright (C) 2006 - 2022 DFTB+ developers group # -# # -# See the LICENSE file for terms of usage and distribution. # -#------------------------------------------------------------------------------# - - -'''This module defines an ASE interface to DFTB+, based on FileIO.''' - - -import os -import copy -import numpy as np - -import hsd -from ase.io import write -from ase.units import Hartree, Bohr -from ase.calculators.calculator import FileIOCalculator - - -def read_input(filename): - '''Deserializes HSD input into nested Python dictionaries. - - Args: - filename (str): name of HSD input to read - - Returns: - dftbinp (dict): DFTB+ input as nested dictionaries - ''' - - dictbuilder = hsd.HsdDictBuilder() - parser = hsd.HsdParser(eventhandler=dictbuilder) - with open(os.path.join(os.getcwd(), filename), 'r') as fobj: - parser.feed(fobj) - dftbinp = dictbuilder.hsddict - - return dftbinp - - -def get_dftb_parameters(skdir='./', driver=None, drivermaxforce=None, - drivermaxsteps=None, scc=False, scctol=None, - maxang=None, kpts=None): - '''Generates a suitable Python dictionary for a - calculation with the hamiltonian set to DFTB. - - Args: - skdir (str): path to Slater-Koster files - driver (str): DFTB+ driver for geometry optimization - drivermaxforce (float): max. force component as convergence - criterion of geometry optimization - drivermaxsteps (int): max. number of geometry steps - scc (bool): True for self-consistent calculations - scctol (float): convergence criterion of SCC cycles - maxang (dict): max. angular momentum of atom types - kpts (list/tuple): K-points for Brillouin zone sampling - (periodic structures) - - Returns: - dftbinp (dict): DFTB+ input as nested dictionaries - ''' - - dftbinp = {} - - if driver is not None: - dftbinp['Driver'] = {} - dftbinp['Driver'][str(driver)] = {} - if drivermaxforce is not None: - dftbinp['Driver'][str(driver)]['MaxForceComponent'] = drivermaxforce - if drivermaxsteps is not None: - dftbinp['Driver'][str(driver)]['MaxSteps'] = drivermaxsteps - - dftbinp['Hamiltonian'] = {} - dftbinp['Hamiltonian']['DFTB'] = {} - dftbinp['Hamiltonian']['DFTB']['Scc'] = 'Yes' if scc else 'No' - - if scctol is not None: - dftbinp['Hamiltonian']['DFTB']['SccTolerance'] = scctol - - if maxang is not None: - dftbinp['Hamiltonian']['DFTB']['MaxAngularMomentum'] = maxang - else: - dftbinp['Hamiltonian']['DFTB']['MaxAngularMomentum'] = {} - - if not skdir.endswith('/'): - skdir += '/' - - skfiledict = { - 'Prefix': skdir, - 'Separator': '"-"', - 'Suffix': '".skf"' - } - - dftbinp['Hamiltonian']['DFTB']['SlaterKosterFiles'] = { - 'Type2FileNames': skfiledict - } - - # Handle different K-point formats - # (note: the ability to handle bandpaths has not yet been implemented) - if kpts is not None: - if np.array(kpts).ndim == 1: - # Case: K-points as (gamma-centered) Monkhorst-Pack grid - mp_mesh = kpts - offsets = [0.] * 3 - props = [elem for elem in mp_mesh if isinstance(elem, str)] - props = [elem.lower() for elem in props] - tgamma = 'gamma' in props and len(props) == 1 - if tgamma: - mp_mesh = mp_mesh[:-1] - eps = 1e-10 - for ii in range(3): - offsets[ii] *= mp_mesh[ii] - assert abs(offsets[ii]) < eps \ - or abs(offsets[ii] - 0.5) < eps - if mp_mesh[ii] % 2 == 0: - offsets[ii] += 0.5 - elif not tgamma and len(mp_mesh) != 3: - raise ValueError('Illegal K-Points definition: ' + str(kpts)) - kpts_mp = np.vstack((np.eye(3) * mp_mesh, offsets)) - dftbinp['Hamiltonian']['DFTB']['KPointsAndWeights'] = { - 'SupercellFolding': kpts_mp - } - elif np.array(kpts).ndim == 2: - # Case: single K-points explicitly as (N x 3) or (N x 4) matrix - kpts_coord = np.array(kpts) - if np.shape(kpts_coord)[1] == 4: - kptsweights = kpts_coord - kpts_coord = kpts_coord[:, :-1] - elif np.shape(kpts_coord)[1] == 3: - kptsweights = np.hstack([kpts_coord, [[1.0],] * \ - np.shape(kpts_coord)[0]]) - else: - raise ValueError('Illegal K-Points definition: ' + str(kpts)) - dftbinp['Hamiltonian']['DFTB']['KPointsAndWeights'] = kptsweights - else: - raise ValueError('Illegal K-Points definition: ' + str(kpts)) - - dftbinp['Options'] = { - 'WriteResultsTag': 'Yes' - } - dftbinp['ParserOptions'] = { - 'IgnoreUnprocessedNodes': 'Yes' - } - - dftbinp['Analysis'] = {} - dftbinp['Analysis']['CalculateForces'] = 'No' - - return dftbinp - - -class DftbPlus(FileIOCalculator): - '''ASE interface to DFTB+. - - A DFTB+ calculator with ase-FileIOCalculator nomenclature. - ''' - - implemented_properties = ('energy', 'forces', 'charges', 'stress') - - - def __init__(self, binary='dftb+', parameters=None, - restart=None, ignore_bad_restart_file=False, - label='dftbplus', **kwargs): - '''Initializes a ctypes DFTB+ calculator object. - - Args: - parameters (dict): dictionary containing (additional) - initialization options for DFTB+ - ''' - - self.input = copy.deepcopy(parameters) - - self._do_forces = False - - self._atoms = None - self._pcpot = None - self._lines = None - self._symids = None - self._atoms_input = None - self._mmpositions = None - - self._skdir = self.input['Hamiltonian']['DFTB'] \ - ['SlaterKosterFiles']['Type2FileNames']['Prefix'] - - FileIOCalculator.__init__(self, restart, ignore_bad_restart_file, - label, self._atoms, binary, **kwargs) - - # Determine number of spin channels - try: - spinpolkeys = list(self.input['Hamiltonian']\ - ['DFTB']['SpinPolarisation'].keys()) - self.nspin = 2 if 'colinear' in spinpolkeys else 1 - except KeyError: - self.nspin = 1 - - - def _write_dftb_in(self, filename): - '''Dumps Python input dictionary to HSD format and writes input file. - If necessary, geometry-relevant entries in the input are updated. - If not all max. angular momenta have been specified by user, an - attempt is made to extract remaining ones out of Slater-Koster files. - - Args: - filename (str): name of input file to be written - ''' - - # Setting up the desired, geometry-relevant entries - specieslist, self._symids = self._chemsym2indices() - - geodict = {} - geodict['TypeNames'] = specieslist - typesandcoords = np.empty([len(self._symids), 4], dtype=object) - typesandcoords[:, 0] = self._symids - typesandcoords[:, 1:] = np.array(self._atoms_input.get_positions(), - dtype=float) - geodict['TypesAndCoordinates'] = typesandcoords - geodict['TypesAndCoordinates.attribute'] = 'Angstrom' - - periodic = any(self._atoms_input.get_pbc()) - - if periodic: - geodict['Periodic'] = 'Yes' - geodict['LatticeVectors'] = self._atoms_input.get_cell()[:] - geodict['LatticeVectors.attribute'] = 'Angstrom' - else: - geodict['Periodic'] = 'No' - - self.input['Geometry'] = geodict - - # Setting up max. angular momenta if not specified by user - # extract maxang's out of Slater-Koster files - speciesdiff = set(specieslist) - \ - set(self.input['Hamiltonian']['DFTB']['MaxAngularMomentum'].keys()) - maxangs = self.input['Hamiltonian']['DFTB']['MaxAngularMomentum'] - for species in speciesdiff: - path = os.path.join(self._skdir, '{0}-{0}.skf'.format(species)) - maxang = read_max_angular_momentum(path) - if maxang is None: - msg = 'Error: Could not read max. angular momentum from ' + \ - 'Slater-Koster file "{0}-{0}.skf".'.format(species) + \ - ' Please specify manually!' - raise TypeError(msg) - if maxang not in range(0, 4): - msg = 'The obtained max. angular momentum from ' + \ - 'Slater-Koster file "{0}-{0}.skf"'.format(species) + \ - ' is out of range. Please check!' - raise ValueError(msg) - maxang = '"{}"'.format('spdf'[maxang]) - maxangs[species] = maxang - - with open(filename, 'w') as outfile: - outfile.write(hsd.dumps(self.input)) - - - def _chemsym2indices(self): - '''Assigns indices to the chemical symbols of the Atoms object. - - Returns: - specieslist (list): occurring atom types in Atoms object - symids (list): indices assigned to atom types in specieslist - ''' - - chemsyms = self._atoms_input.get_chemical_symbols() - - symdict = {} - for sym in chemsyms: - if not sym in symdict: - symdict[sym] = len(symdict) + 1 - - specieslist = list(['null'] * len(symdict.keys())) - for sym in symdict: - specieslist[symdict[sym] - 1] = sym - - symids = [] - for sym in chemsyms: - symids.append(symdict[sym]) - - return specieslist, symids - - - def set(self, **kwargs): - changed_parameters = FileIOCalculator.set(self, **kwargs) - if changed_parameters: - self.reset() - return changed_parameters - - - def check_state(self, atoms): - system_changes = FileIOCalculator.check_state(self, atoms) - # Ignore unit cell for molecules: - if not atoms.pbc.any() and 'cell' in system_changes: - system_changes.remove('cell') - if self._pcpot and self._pcpot.mmpositions is not None: - system_changes.append('positions') - return system_changes - - - def write_input(self, atoms, properties=None, system_changes=None): - if properties is not None: - if 'forces' in properties or 'stress' in properties: - self._do_forces = True - self.input['Analysis']['CalculateForces'] = 'Yes' - FileIOCalculator.write_input( - self, atoms, properties, system_changes) - - # self._atoms is none until results are read out, - # then it is set to the ones at writing input - self._atoms_input = atoms - self._atoms = None - - self._write_dftb_in(os.path.join(self.directory, 'dftb_in.hsd')) - - write(os.path.join(self.directory, 'geo_end.gen'), atoms) - - if self._pcpot: - self._pcpot.write_mmcharges('dftb_external_charges.dat') - - - def read_results(self): - ''' All results are read from results.tag file. - It will be destroyed after it is read to avoid - reading it once again after some runtime error. - ''' - - myfile = open(os.path.join(self.directory, 'results.tag'), 'r') - self._lines = myfile.readlines() - myfile.close() - - self._atoms = self._atoms_input - charges, energy = self.read_charges_and_energy() - if charges is not None: - self.results['charges'] = charges - self.results['energy'] = energy - if self._do_forces: - forces = self.read_forces() - self.results['forces'] = forces - self._mmpositions = None - - # stress stuff begins - sstring = 'stress' - have_stress = False - stress = list() - for iline, line in enumerate(self._lines): - if sstring in line: - have_stress = True - start = iline + 1 - end = start + 3 - for i in range(start, end): - cell = [float(x) for x in self._lines[i].split()] - stress.append(cell) - if have_stress: - stress = -np.array(stress) * Hartree / Bohr**3 - self.results['stress'] = stress.flat[[0, 4, 8, 5, 2, 1]] - # stress stuff ends - - # eigenvalues and fermi levels - fermi_levels = self.read_fermi_levels() - if fermi_levels is not None: - self.results['fermi_levels'] = fermi_levels - - eigenvalues = self.read_eigenvalues() - if eigenvalues is not None: - self.results['eigenvalues'] = eigenvalues - - # calculation was carried out with atoms written in write_input - os.remove(os.path.join(self.directory, 'results.tag')) - - - def read_forces(self): - '''Read forces from DFTB+ output file (results.tag).''' - - # Force line indexes - for iline, line in enumerate(self._lines): - fstring = 'forces ' - if line.find(fstring) >= 0: - index_force_begin = iline + 1 - line1 = line.replace(':', ',') - index_force_end = iline + 1 + \ - int(line1.split(',')[-1]) - break - - gradients = [] - for j in range(index_force_begin, index_force_end): - word = self._lines[j].split() - gradients.append([float(word[k]) for k in range(0, 3)]) - - return np.array(gradients) * Hartree / Bohr - - - def read_charges_and_energy(self): - '''Get partial charges on atoms - in case we cannot find charges they are set to None - ''' - - infile = open(os.path.join(self.directory, 'detailed.out'), 'r') - lines = infile.readlines() - infile.close() - - for line in lines: - if line.strip().startswith('Total energy:'): - energy = float(line.split()[2]) * Hartree - break - - qm_charges = [] - for nn, line in enumerate(lines): - if ('Atom' and 'Charge' in line): - chargestart = nn + 1 - break - else: - # print('Warning: did not find DFTB-charges') - # print('This is ok if flag SCC=No') - return None, energy - - lines1 = lines[chargestart:(chargestart + len(self._atoms))] - for line in lines1: - qm_charges.append(float(line.split()[-1])) - - return np.array(qm_charges), energy - - - def get_charges(self, atoms): - '''Get the calculated charges. - This is inhereted to atoms object. - ''' - - if 'charges' in self.results: - return self.results['charges'] - - return None - - - def read_eigenvalues(self): - '''Read Eigenvalues from DFTB+ output file (results.tag). - Unfortunately, the order seems to be scrambled. - ''' - - # Eigenvalue line indexes - index_eig_begin = None - for iline, line in enumerate(self._lines): - fstring = 'eigenvalues ' - if line.find(fstring) >= 0: - index_eig_begin = iline + 1 - line1 = line.replace(':', ',') - ncol, nband, nkpt, nspin = map(int, line1.split(',')[-4:]) - break - else: - return None - - # Take into account that the last row may lack - # columns if nkpt * nspin * nband % ncol != 0 - nrow = int(np.ceil(nkpt * nspin * nband * 1. / ncol)) - index_eig_end = index_eig_begin + nrow - ncol_last = len(self._lines[index_eig_end - 1].split()) - self._lines[index_eig_end - 1] += ' 0.0 ' * (ncol - ncol_last) - - eig = np.loadtxt(self._lines[index_eig_begin:index_eig_end]).flatten() - eig *= Hartree - nn = nkpt * nband - eigenvalues = [eig[i * nn:(i + 1) * nn].reshape((nkpt, nband)) - for i in range(nspin)] - - return eigenvalues - - - def read_fermi_levels(self): - '''Read Fermi level(s) from DFTB+ output file (results.tag).''' - - # Fermi level line indexes - for iline, line in enumerate(self._lines): - fstring = 'fermi_level ' - if line.find(fstring) >= 0: - index_fermi = iline + 1 - break - else: - return None - - fermi_levels = [] - words = self._lines[index_fermi].split() - assert len(words) == 2 - - for word in words: - ee = float(word) - if abs(ee) > 1e-08: - # Without spin polarization, one of the Fermi - # levels is equal to 0.000000000000000E+000 - fermi_levels.append(ee) - - return np.array(fermi_levels) * Hartree - - - def get_number_of_spins(self): - '''Auxiliary function to extract result: number of spins.''' - - return self.nspin - - - # !!!BROKEN!!! - def get_ibz_k_points(self): - '''Auxiliary function to extract: K-points.''' - - return self.kpts_coord.copy() - - - # !!!BROKEN!!! - def get_eigenvalues(self, kpt=0, spin=0): - '''Auxiliary function to extract result: eigenvalues.''' - - return self.results['eigenvalues'][spin][kpt].copy() - - - # !!!BROKEN!!! - def get_fermi_levels(self): - '''Auxiliary function to extract result: Fermi level(s).''' - - return self.results['fermi_levels'].copy() - - - # !!!BROKEN!!! - def get_fermi_level(self): - '''Auxiliary function to extract result: result: max. Fermi level.''' - - return max(self.get_fermi_levels()) - - -def read_max_angular_momentum(path): - '''Read maximum angular momentum from .skf file. - See dftb.org/fileadmin/DFTB/public/misc/slakoformat.pdf - for a detailed description of the Slater-Koster file format. - - Args: - path (str): path to Slater-Koster file - - Returns: - maxang (int): max. angular momentum (None if extraction fails) - ''' - - with open(path, 'r') as skf: - line = skf.readline() - if line[0] == '@': - # Skip additional line for extended format - line = skf.readline() - - # Replace any commas that may appear - # (inconsistency in .skf files) - line = line.replace(',', ' ').split() - - if len(line) == 3: - # max. angular momentum specified: - # extraction possible - maxang = int(line[2]) - 1 - else: - # max. angular momentum not specified - # or wrong format: extraction not possible - maxang = None - - return maxang - - -# def embed(self, mmcharges=None, directory='./'): -# '''Embed atoms in point-charges (mmcharges).''' -# -# self._pcpot = PointChargePotential(mmcharges, self.directory) -# return self._pcpot - - -#class PointChargePotential: -# def __init__(self, mmcharges, directory='./'): -# '''Point-charge potential for DFTB+.''' -# -# self.mmcharges = mmcharges -# self.directory = directory -# self._mmpositions = None -# self.mmforces = None -# -# -# def set_positions(self, mmpositions): -# self._mmpositions = mmpositions -# -# -# def set_charges(self, mmcharges): -# self.mmcharges = mmcharges -# -# -# def write_mmcharges(self, filename='dftb_external_charges.dat'): -# '''mok all -# Write external charges as monopoles for DFTB+. -# ''' -# -# if self.mmcharges is None: -# print('DFTB: Warning: not writing external charges ') -# return -# charge_file = open(os.path.join(self.directory, filename), 'w') -# for [pos, charge] in zip(self._mmpositions, self.mmcharges): -# [x, y, z] = pos -# charge_file.write('%12.6f %12.6f %12.6f %12.6f \n' -# % (x, y, z, charge)) -# charge_file.close() -# -# -# def get_forces(self, calc, get_forces=True): -# '''Returns forces on point charges if the flag get_forces=True.''' -# -# if get_forces: -# return self.read_forces_on_pointcharges() -# else: -# return np.zeros_like(self._mmpositions) -# -# -# def read_forces_on_pointcharges(self): -# '''Read forces from DFTB+ output file (detailed.out).''' -# -# from ase.units import Hartree, Bohr -# infile = open(os.path.join(self.directory, 'detailed.out'), 'r') -# lines = infile.readlines() -# infile.close() -# -# external_forces = [] -# for n, line in enumerate(lines): -# if 'Forces on external charges' in line: -# chargestart = n + 1 -# break -# else: -# raise RuntimeError( -# 'Problem in reading forces on MM external-charges') -# lines1 = lines[chargestart:(chargestart + len(self.mmcharges))] -# for line in lines1: -# external_forces.append( -# [float(i) for i in line.split()]) -# return np.array(external_forces) * Hartree / Bohr - diff --git a/src/dftbplus/io.py b/src/dftbplus/io.py new file mode 100644 index 0000000..c774bdf --- /dev/null +++ b/src/dftbplus/io.py @@ -0,0 +1,130 @@ +#------------------------------------------------------------------------------# +# DFTB+: general package for performing fast atomistic simulations # +# Copyright (C) 2006 - 2022 DFTB+ developers group # +# # +# See the LICENSE file for terms of usage and distribution. # +#------------------------------------------------------------------------------# + +import numpy as np +from ase.atoms import Atoms +from dftbplus_ptools.xyz import Xyz +from dftbplus_ptools.gen import Gen +from dftbplus_ptools.geometry import Geometry +from dftbplus_ptools.hsdinput import Hsdinput + + +def read_dftb(filename): + """Method to read coordinates from the Geometry section + of a DFTB+ input file (typically called "dftb_in.hsd"). + + As described in the DFTB+ manual, this section can be + in a number of different formats. This reader supports + the GEN format, the so-called "explicit" format and the XYZ format. + + The "explicit" format is unique to DFTB+ input files. + The GEN format can also be used in a stand-alone fashion, + as coordinate files with a `.gen` extension. Reading and + writing such files is implemented in `ase.io.gen`. + """ + + try: + xyz = Xyz.fromhsd(filename) + geo = xyz.geometry + except KeyError: + try: + gen = Gen.fromhsd(filename) + geo = gen.geometry + except KeyError: + raise NotImplementedError(f"The Geometry scetion in '{filename}'" + + " is not supported!") + + return geometry_to_atom(geo) + + +def atom_to_geometry(atom): + """converts atom object to dftbplus_ptools.geometry.Geometry object + + Args: + atom (ase atoms object): atom to be converted + + Returns: + (dftbplus_ptools.geometry.Geometry object): converted geometry + """ + chemsyms = atom.get_chemical_symbols() + + symdict = {} + for sym in chemsyms: + if sym not in symdict: + symdict[sym] = len(symdict) + + specienames = list(['null'] * len(symdict.keys())) + for sym, num in symdict.items(): + specienames[num] = sym + + indexes = [] + for sym in chemsyms: + indexes.append(symdict[sym]) + coords = atom.get_positions() + + if any(atom.get_pbc()) and all(atom.get_pbc()): + latvecs = atom.get_cell().array + elif not any(atom.get_pbc()) and not all(atom.get_pbc()): + latvecs = None + else: + raise NotImplementedError("The converted Atoms object contains " + + "1D or 2D PBC's. The geometry object " + + "only supports 3D PBC's!") + + return Geometry(specienames, indexes, coords, latvecs=latvecs) + + +def geometry_to_atom(geometry): + """converts dftbplus_ptools.geometry.Geometry object to atom object + + Args: + geometry (dftbplus_ptools.geometry.Geometry object): geometry to be + converted + + Returns: + (ase atoms object): atom with geometry + """ + specienames = geometry.specienames + indexes = geometry.indexes + atoms_pos = geometry.coords + my_pbc = geometry.periodic + + atom_symbols = [] + for index in indexes: + atom_symbols.append(specienames[index]) + + if my_pbc: + mycell = geometry.latvecs + elif not my_pbc: + mycell = np.zeros((3, 3)) + + return Atoms(positions=atoms_pos, symbols=atom_symbols, cell=mycell, + pbc=my_pbc) + + +def atom_to_hsd(atom, filename="dftb_in.hsd", directory=".", dictionary=None, + get_class=False): + """help function for reading geometry from atom and adding to hsd + + Args: + atom (ase atoms object): atom with geometry + directory (str): directory of file + filename (str): name of file + dictionary (dict): option to change existing dict + get_class (bool): True if dftbplus_ptools Hsdinput class is returned, + False if dict is returned + + Returns: + (dftbplus_ptools Hsdinput class/dict): Hsdinput class/dict with + geometry added + """ + geo = atom_to_geometry(atom) + hsd = Hsdinput(filename, directory, dictionary) + hsd.set_geometry(geo) + if get_class: + return hsd + return hsd.get_hsd() diff --git a/src/dftbplus/legacy_calculator.py b/src/dftbplus/legacy_calculator.py new file mode 100644 index 0000000..4020c2c --- /dev/null +++ b/src/dftbplus/legacy_calculator.py @@ -0,0 +1,435 @@ +#------------------------------------------------------------------------------# +# DFTB+: general package for performing fast atomistic simulations # +# Copyright (C) 2006 - 2022 DFTB+ developers group # +# # +# See the LICENSE file for terms of usage and distribution. # +#------------------------------------------------------------------------------# + +""" This module defines a FileIOCalculator for DFTB+ + +http://www.dftbplus.org/ +http://www.dftb.org/ + +Initial development: markus.kaukonen@iki.fi +""" + +import os +from io import StringIO +import copy +import numpy as np +from ase.calculators.calculator import (FileIOCalculator, kpts2ndarray, + kpts2sizeandoffsets) +import hsd +import dftbplus_ptools.resultstag +from dftbplus_ptools.hsdinput import Hsdinput +from dftbplus_ase.calculator import DftbPlus +from ase.units import Hartree, Bohr + + +class Dftb(DftbPlus): + """legacy calculator based on new dftbplus calculator""" + + discard_results_on_any_change = True + + def __init__(self, restart=None, + ignore_bad_restart_file=FileIOCalculator._deprecated, + label='dftb', atoms=None, kpts=None, + slako_dir=None, + **kwargs): + """ + All keywords for the dftb_in.hsd input file (see the DFTB+ manual) + can be set by ASE. Consider the following input file block: + + >>> Hamiltonian = DFTB { + >>> SCC = Yes + >>> SCCTolerance = 1e-8 + >>> MaxAngularMomentum = { + >>> H = s + >>> O = p + >>> } + >>> } + + This can be generated by the DFTB+ calculator by using the + following settings: + + >>> calc = Dftb(Hamiltonian_='DFTB', # line is included by default + >>> Hamiltonian_SCC='Yes', + >>> Hamiltonian_SCCTolerance=1e-8, + >>> Hamiltonian_MaxAngularMomentum_='', + >>> Hamiltonian_MaxAngularMomentum_H='s', + >>> Hamiltonian_MaxAngularMomentum_O='p') + + In addition to keywords specific to DFTB+, also the following keywords + arguments can be used: + + restart: str + Prefix for restart file. May contain a directory. + Default is None: don't restart. + ignore_bad_restart_file: bool + Ignore broken or missing restart file. By default, it is an + error if the restart file is missing or broken. + label: str (default 'dftb') + Prefix used for the main output file (