Skip to main content
`Ctrl`+`K`

#

# Source code for ase.calculators.dftb

    # fmt: off

    """ This module defines a FileIOCalculator for DFTB+

    http://www.dftbplus.org/
    http://www.dftb.org/

    Initial development: markus.kaukonen@iki.fi
    """

    import os

    import numpy as np

    from ase.calculators.calculator import (
        BadConfiguration,
        FileIOCalculator,
        kpts2ndarray,
        kpts2sizeandoffsets,
    )
    from ase.units import Bohr, Hartree

    class Dftb(FileIOCalculator):
        implemented_properties = ['energy', 'forces', 'charges',
                                  'stress', 'dipole']
        discard_results_on_any_change = True

        fileio_rules = FileIOCalculator.ruleset(
            configspec=dict(skt_path=None),
            stdout_name='{prefix}.out')

        def __init__(self, restart=None,
                     ignore_bad_restart_file=FileIOCalculator._deprecated,
                     label='dftb', atoms=None, kpts=None,
                     slako_dir=None,
                     command=None,
                     profile=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:

            >>> from ase.calculators.dftb import Dftb
            >>>
            >>> 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 (<label>.out).
            atoms: Atoms object (default None)
                Optional Atoms object to which the calculator will be
                attached. When restarting, atoms will get its positions and
                unit-cell updated from file.
            kpts: (default None)
                Brillouin zone sampling:

                * ``(1,1,1)`` or ``None``: Gamma-point only
                * ``(n1,n2,n3)``: Monkhorst-Pack grid
                * ``dict``: Interpreted as a path in the Brillouin zone if
                  it contains the 'path_' keyword. Otherwise it is converted
                  into a Monkhorst-Pack grid using
                  ``ase.calculators.calculator.kpts2sizeandoffsets``
                * ``[(k11,k12,k13),(k21,k22,k23),...]``: Explicit (Nkpts x 3)
                  array of k-points in units of the reciprocal lattice vectors
                  (each with equal weight)

            Additional attribute to be set by the embed() method:

            pcpot: PointCharge object
                An external point charge potential (for QM/MM calculations)
            """

            if command is None:
                if 'DFTB_COMMAND' in self.cfg:
                    command = self.cfg['DFTB_COMMAND'] + ' > PREFIX.out'

            if command is None and profile is None:
                try:
                    profile = self.load_argv_profile(self.cfg, 'dftb')
                except BadConfiguration:
                    pass

            if command is None:
                command = 'dftb+ > PREFIX.out'

            if slako_dir is None:
                if profile is not None:
                    slako_dir = profile.configvars.get('skt_path')

                if slako_dir is None:
                    slako_dir = self.cfg.get('DFTB_PREFIX', './')
                if not slako_dir.endswith('/'):
                    slako_dir += '/'

            self.slako_dir = slako_dir
            if kwargs.get('Hamiltonian_', 'DFTB') == 'DFTB':
                self.default_parameters = dict(
                    Hamiltonian_='DFTB',
                    Hamiltonian_SlaterKosterFiles_='Type2FileNames',
                    Hamiltonian_SlaterKosterFiles_Prefix=self.slako_dir,
                    Hamiltonian_SlaterKosterFiles_Separator='"-"',
                    Hamiltonian_SlaterKosterFiles_Suffix='".skf"',
                    Hamiltonian_MaxAngularMomentum_='',
                    Options_='',
                    Options_WriteResultsTag='Yes',
                    ParserOptions_='',
                    ParserOptions_ParserVersion=1,
                    ParserOptions_IgnoreUnprocessedNodes='Yes')
            else:
                self.default_parameters = dict(
                    Options_='',
                    Options_WriteResultsTag='Yes',
                    ParserOptions_='',
                    ParserOptions_ParserVersion=1,
                    ParserOptions_IgnoreUnprocessedNodes='Yes')

            self.pcpot = None
            self.lines = None
            self.atoms = None
            self.atoms_input = None
            self.do_forces = False

            super().__init__(restart, ignore_bad_restart_file,
                             label, atoms, command=command,
                             profile=profile, **kwargs)

            # Determine number of spin channels
            try:
                entry = kwargs['Hamiltonian_SpinPolarisation']
                spinpol = 'colinear' in entry.lower()
            except KeyError:
                spinpol = False
            self.nspin = 2 if spinpol else 1

            # kpoint stuff by ase
            self.kpts = kpts
            self.kpts_coord = None

            if self.kpts is not None:
                initkey = 'Hamiltonian_KPointsAndWeights'
                mp_mesh = None
                offsets = None

                if isinstance(self.kpts, dict):
                    if 'path' in self.kpts:
                        # kpts is path in Brillouin zone
                        self.parameters[initkey + '_'] = 'Klines '
                        self.kpts_coord = kpts2ndarray(self.kpts, atoms=atoms)
                    else:
                        # kpts is (implicit) definition of
                        # Monkhorst-Pack grid
                        self.parameters[initkey + '_'] = 'SupercellFolding '
                        mp_mesh, offsets = kpts2sizeandoffsets(atoms=atoms,
                                                               **self.kpts)
                elif np.array(self.kpts).ndim == 1:
                    # kpts is Monkhorst-Pack grid
                    self.parameters[initkey + '_'] = 'SupercellFolding '
                    mp_mesh = self.kpts
                    offsets = [0.] * 3
                elif np.array(self.kpts).ndim == 2:
                    # kpts is (N x 3) list/array of k-point coordinates
                    # each will be given equal weight
                    self.parameters[initkey + '_'] = ''
                    self.kpts_coord = np.array(self.kpts)
                else:
                    raise ValueError('Illegal kpts definition:' + str(self.kpts))

                if mp_mesh is not None:
                    eps = 1e-10
                    for i in range(3):
                        key = initkey + '_empty%03d' % i
                        val = [mp_mesh[i] if j == i else 0 for j in range(3)]
                        self.parameters[key] = ' '.join(map(str, val))
                        offsets[i] *= mp_mesh[i]
                        assert abs(offsets[i]) < eps or abs(offsets[i] - 0.5) < eps
                        # DFTB+ uses a different offset convention, where
                        # the k-point mesh is already Gamma-centered prior
                        # to the addition of any offsets
                        if mp_mesh[i] % 2 == 0:
                            offsets[i] += 0.5
                    key = initkey + '_empty%03d' % 3
                    self.parameters[key] = ' '.join(map(str, offsets))

                elif self.kpts_coord is not None:
                    for i, c in enumerate(self.kpts_coord):
                        key = initkey + '_empty%09d' % i
                        c_str = ' '.join(map(str, c))
                        if 'Klines' in self.parameters[initkey + '_']:
                            c_str = '1 ' + c_str
                        else:
                            c_str += ' 1.0'
                        self.parameters[key] = c_str

        def write_dftb_in(self, outfile):
            """ Write the innput file for the dftb+ calculation.
                Geometry is taken always from the file 'geo_end.gen'.
            """

            outfile.write('Geometry = GenFormat { \n')
            outfile.write('    <<< "geo_end.gen" \n')
            outfile.write('} \n')
            outfile.write(' \n')

            params = self.parameters.copy()

            s = 'Hamiltonian_MaxAngularMomentum_'
            for key in params:
                if key.startswith(s) and len(key) > len(s):
                    break
            else:
                if params.get('Hamiltonian_', 'DFTB') == 'DFTB':
                    # User didn't specify max angular mometa.  Get them from
                    # the .skf files:
                    symbols = set(self.atoms.get_chemical_symbols())
                    for symbol in symbols:
                        path = os.path.join(self.slako_dir,
                                            '{0}-{0}.skf'.format(symbol))
                        l = read_max_angular_momentum(path)
                        params[s + symbol] = '"{}"'.format('spdf'[l])

            if self.do_forces:
                params['Analysis_'] = ''
                params['Analysis_CalculateForces'] = 'Yes'

            # --------MAIN KEYWORDS-------
            previous_key = 'dummy_'
            myspace = ' '
            for key, value in sorted(params.items()):
                current_depth = key.rstrip('_').count('_')
                previous_depth = previous_key.rstrip('_').count('_')
                for my_backslash in reversed(
                        range(previous_depth - current_depth)):
                    outfile.write(3 * (1 + my_backslash) * myspace + '} \n')
                outfile.write(3 * current_depth * myspace)
                if key.endswith('_') and len(value) > 0:
                    outfile.write(key.rstrip('_').rsplit('_')[-1] +
                                  ' = ' + str(value) + '{ \n')
                elif (key.endswith('_') and (len(value) == 0)
                      and current_depth == 0):  # E.g. 'Options {'
                    outfile.write(key.rstrip('_').rsplit('_')[-1] +
                                  ' ' + str(value) + '{ \n')
                elif (key.endswith('_') and (len(value) == 0)
                      and current_depth > 0):  # E.g. 'Hamiltonian_Max... = {'
                    outfile.write(key.rstrip('_').rsplit('_')[-1] +
                                  ' = ' + str(value) + '{ \n')
                elif key.count('_empty') == 1:
                    outfile.write(str(value) + ' \n')
                elif ((key == 'Hamiltonian_ReadInitialCharges') and
                      (str(value).upper() == 'YES')):
                    f1 = os.path.isfile(self.directory + os.sep + 'charges.dat')
                    f2 = os.path.isfile(self.directory + os.sep + 'charges.bin')
                    if not (f1 or f2):
                        print('charges.dat or .bin not found, switching off guess')
                        value = 'No'
                    outfile.write(key.rsplit('_')[-1] + ' = ' + str(value) + ' \n')
                else:
                    outfile.write(key.rsplit('_')[-1] + ' = ' + str(value) + ' \n')
                if self.pcpot is not None and ('DFTB' in str(value)):
                    outfile.write('   ElectricField = { \n')
                    outfile.write('      PointCharges = { \n')
                    outfile.write(
                        '         CoordsAndCharges [Angstrom] = DirectRead { \n')
                    outfile.write('            Records = ' +
                                  str(len(self.pcpot.mmcharges)) + ' \n')
                    outfile.write(
                        '            File = "dftb_external_charges.dat" \n')
                    outfile.write('         } \n')
                    outfile.write('      } \n')
                    outfile.write('   } \n')
                previous_key = key
            current_depth = key.rstrip('_').count('_')
            for my_backslash in reversed(range(current_depth)):
                outfile.write(3 * my_backslash * myspace + '} \n')

        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):
            from ase.io import write
            if properties is not None:
                if 'forces' in properties or 'stress' in properties:
                    self.do_forces = True
            FileIOCalculator.write_input(
                self, atoms, properties, system_changes)
            with open(os.path.join(self.directory, 'dftb_in.hsd'), 'w') as fd:
                self.write_dftb_in(fd)
            write(os.path.join(self.directory, 'geo_end.gen'), atoms,
                  parallel=False)
            # 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
            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 """

            with open(os.path.join(self.directory, 'results.tag')) as fd:
                self.lines = fd.readlines()

            self.atoms = self.atoms_input
            charges, energy, dipole = self.read_charges_energy_dipole()
            if charges is not None:
                self.results['charges'] = charges
            self.results['energy'] = energy
            if dipole is not None:
                self.results['dipole'] = dipole
            if self.do_forces:
                forces = self.read_forces()
                self.results['forces'] = forces
            self.mmpositions = None

            # stress stuff begins
            sstring = 'stress'
            have_stress = False
            stress = []
            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)."""
            from ase.units import Bohr, Hartree

            # Initialise the indices so their scope
            # reaches outside of the for loop
            index_force_begin = -1
            index_force_end = -1

            # 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(3)])

            return np.array(gradients) * Hartree / Bohr

        def read_charges_energy_dipole(self):
            """Get partial charges on atoms
                in case we cannot find charges they are set to None
            """
            with open(os.path.join(self.directory, 'detailed.out')) as fd:
                lines = fd.readlines()

            for line in lines:
                if line.strip().startswith('Total energy:'):
                    energy = float(line.split()[2]) * Hartree
                    break

            qm_charges = []
            for n, line in enumerate(lines):
                if ('Atom' and 'Charge' in line):
                    chargestart = n + 1
                    break
            else:
                # print('Warning: did not find DFTB-charges')
                # print('This is ok if flag SCC=No')
                return None, energy, None

            lines1 = lines[chargestart:(chargestart + len(self.atoms))]
            for line in lines1:
                qm_charges.append(float(line.split()[-1]))

            dipole = None
            for line in lines:
                if 'Dipole moment:' in line and 'au' in line:
                    line = line.replace("Dipole moment:", "").replace("au", "")
                    dipole = np.array(line.split(), dtype=float) * Bohr

            return np.array(qm_charges), energy, dipole

        def get_charges(self, atoms):
            """ Get the calculated charges
            this is inhereted to atoms object """
            if 'charges' in self.results:
                return self.results['charges']
            else:
                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())
            # XXX dirty fix
            self.lines[index_eig_end - 1] = (
                self.lines[index_eig_end - 1].strip()
                + ' 0.0 ' * (ncol - ncol_last))

            eig = np.loadtxt(self.lines[index_eig_begin:index_eig_end]).flatten()
            eig *= Hartree
            N = nkpt * nband
            eigenvalues = [eig[i * N:(i + 1) * N].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) in [1, 2], 'Expected either 1 or 2 Fermi levels'

            for word in words:
                e = float(word)
                # In non-spin-polarized calculations with DFTB+ v17.1,
                # two Fermi levels are given, with the second one being 0,
                # but we don't want to add that one to the list
                if abs(e) > 1e-8:
                    fermi_levels.append(e)

            return np.array(fermi_levels) * Hartree

        def get_ibz_k_points(self):
            return self.kpts_coord.copy()

        def get_number_of_spins(self):
            return self.nspin

        def get_eigenvalues(self, kpt=0, spin=0):
            return self.results['eigenvalues'][spin][kpt].copy()

        def get_fermi_levels(self):
            return self.results['fermi_levels'].copy()

        def get_fermi_level(self):
            return max(self.get_fermi_levels())

        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):
            """ mok all
            write external charges as monopoles for dftb+.

            """
            if self.mmcharges is None:
                print("DFTB: Warning: not writing exernal charges ")
                return
            with open(os.path.join(self.directory, filename), 'w') as charge_file:
                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))

        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 (results.tag)."""
            from ase.units import Bohr, Hartree
            with open(os.path.join(self.directory, 'detailed.out')) as fd:
                lines = fd.readlines()

            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

    def read_max_angular_momentum(path):
        """Read maximum angular momentum from .skf file.

        See dftb.org for A detailed description of the Slater-Koster file format.
        """
        with open(path) as fd:
            line = fd.readline()
            if line[0] == '@':
                # Extended format
                fd.readline()
                l = 3
                pos = 9
            else:
                # Simple format:
                l = 2
                pos = 7

            # Sometimes there ar commas, sometimes not:
            line = fd.readline().replace(',', ' ')

            occs = [float(f) for f in line.split()[pos:pos + l + 1]]
            for f in occs:
                if f > 0.0:
                    return l
                l -= 1