Skip to main content
`Ctrl`+`K`

# DFTB+

## Introduction

DFTB+ is a density-functional based tight-binding code using atom-centered orbitals. This interface makes it possible to use DFTB+ as a calculator in ASE. You need Slater-Koster files for the combination of atom types of your system. These can be obtained at dftb.org.

## Environment variables

The default command that ASE will use to start DFTB+ is `dftb+ > PREFIX.out`. You can change this command by setting the `ASE_DFTB_COMMAND` environment variable, e.g.:

    $ export ASE_DFTB_COMMAND="/path/to/dftb+ > PREFIX.out"

For compatibility, also the old `DFTB_COMMAND` variable can be used, and the resulting command will be `$DFTB_COMMAND > PREFIX.out`. Before each DFTB+ calculation, also make sure that the `DFTB_PREFIX` variable points to the directory where the Slater-Koster files are kept, e.g.:

    $ export DFTB_PREFIX=/path/to/mio-0-1/

## Parameters

As a FileIOCalculator, `ase.calculators.dftb.Dftb` writes input files, runs DFTB+, and extracts the required information from the resulting output. The input files are `dftb_in.hsd` (the calculation settings), `geo_end.gen` (the initial geometry) and `dftb_external_charges.dat` (the external point charges in case of electrostatic QM/MM embedding).

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',  # this 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)

## Examples

Below are three examples of how to use the DFTB+ calculator. The required SKF files (i.e. `H-H.skf`, `H-O.skf`, `O-H.skf` and `O-O.skf`) can for example be chosen from the mio-1-1 parameter set.

### Geometry optimization by ASE

    from ase.build import molecule
    from ase.calculators.dftb import Dftb
    from ase.io import write
    from ase.optimize import QuasiNewton

    atoms = molecule('H2O')
    calc = Dftb(
        label='h2o',
        Hamiltonian_MaxAngularMomentum_='',
        Hamiltonian_MaxAngularMomentum_O='p',
        Hamiltonian_MaxAngularMomentum_H='s',
    )
    atoms.calc = calc

    dyn = QuasiNewton(atoms, trajectory='test.traj')
    dyn.run(fmax=0.01)
    write('final.xyz', atoms)

### Geometry optimization by DFTB+

    from ase.build import molecule
    from ase.calculators.dftb import Dftb
    from ase.io import read, write

    atoms = molecule('H2O')
    calc = Dftb(
        label='h2o',
        Driver_='ConjugateGradient',
        Driver_MaxForceComponent=1e-4,
        Driver_MaxSteps=1000,
        Hamiltonian_MaxAngularMomentum_='',
        Hamiltonian_MaxAngularMomentum_O='p',
        Hamiltonian_MaxAngularMomentum_H='s',
    )

    atoms.calc = calc
    calc.calculate(atoms)

    # The 'geo_end.gen' file written by the ASE calculator
    # (containing the initial geometry), has been overwritten
    # by DFTB+ and now contains the final, optimized geometry.
    final = read('geo_end.gen')
    write('final.xyz', final)

### NVE-MD followed by NVT-MD (both by DFTB+)

Note that this example is unphysical for at least two reasons:

> - no spin polarization for the oxygen molecule

> - the Berendsen coupling is too strong (0.01 here should be 0.0001)

    # fun collision of:  2 H2 + O2 -> 2 H2O
    import os

    from ase.build import molecule
    from ase.calculators.dftb import Dftb
    from ase.io import read, write
    from ase.io.dftb import read_dftb_velocities, write_dftb_velocities

    o2 = molecule('O2')
    h2_1 = molecule('H2')
    h2_2 = molecule('H2')
    o2.translate([0, 0.01, 0])
    h2_1.translate([0, 0, 3])
    h2_1.euler_rotate(center='COP', theta=90)
    h2_2.translate([0, 0, -3])
    h2_2.euler_rotate(center='COP', theta=90)
    o2.set_velocities(([0, 0, 0], [0, 0, 0]))
    h2_1.set_velocities(([0, 0, -3.00], [0, 0, -3.000]))
    h2_2.set_velocities(([0, 0, 3.000], [0, 0, 3.000]))
    atoms = o2 + h2_1 + h2_2

    # 1fs = 41.3 au
    # 1000K = 0.0031668 au
    calculator_NVE = Dftb(
        label='h2o',
        Hamiltonian_MaxAngularMomentum_='',
        Hamiltonian_MaxAngularMomentum_O='p',
        Hamiltonian_MaxAngularMomentum_H='s',
        Driver_='VelocityVerlet',
        Driver_MDRestartFrequency=10,
        Driver_Velocities_='',
        Driver_Velocities_empty='<<+ "velocities.txt"',
        Driver_Steps=1000,
        Driver_KeepStationary='Yes',
        Driver_TimeStep=4.13,
        Driver_Thermostat_='None',
        Driver_Thermostat_empty='',
    )

    # 1fs = 41.3 au
    # 1000K = 0.0031668 au
    calculator_NVT = Dftb(
        label='h2o',
        Hamiltonian_MaxAngularMomentum_='',
        Hamiltonian_MaxAngularMomentum_O='p',
        Hamiltonian_MaxAngularMomentum_H='s',
        Driver_='VelocityVerlet',
        Driver_MDRestartFrequency=5,
        Driver_Velocities_='',
        Driver_Velocities_empty='<<+ "velocities.txt"',
        Driver_Steps=500,
        Driver_KeepStationary='Yes',
        Driver_TimeStep=8.26,
        Driver_Thermostat_='Berendsen',
        Driver_Thermostat_Temperature=0.00339845142,  # 800 degC
        Driver_Thermostat_CouplingStrength=0.01,
    )

    write_dftb_velocities(atoms, 'velocities.txt')

    atoms.calc = calculator_NVE
    atoms.get_potential_energy()  # run NVE ensemble using DFTB+'s own driver
    atoms = read('geo_end.gen')
    write('after_NVE.xyz', atoms)

    read_dftb_velocities(atoms, filename='geo_end.xyz')
    write_dftb_velocities(atoms, 'velocities.txt')
    os.system('mv geo_end.xyz geo_end_NVE.xyz')

    atoms.calc = calculator_NVT
    atoms.get_potential_energy()  # run NVT ensemble using DFTB+'s own driver
    atoms = read('geo_end.gen')
    write('after_NVT.xyz', atoms)

    read_dftb_velocities(atoms, filename='geo_end.xyz')
    write_dftb_velocities(atoms, 'velocities.txt')
    os.system('mv geo_end.xyz geo_end_NVT.xyz')

    # to watch:
    #  ase gui geo_end_NVE.xyz geo_end_NVT.xyz

## DFTB+ calculator class

class ase.calculators.dftb.Dftb(restart=None, ignore_bad_restart_file=<object object>, label='dftb', atoms=None, kpts=None, slako_dir=None, command=None, profile=None, with_tags=False, **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)

  with_tags: bool (default False)

    Whether to use atom tags to distinguish between different chemical species. When enabled, atoms with different tags are treated as distinct species, even if they have the same chemical symbol. For example, hydrogen atoms in a nanoparticle (tag 1) and hydrogen atoms in a solvent (tag 2) would be represented as “H1” and “H2”, requiring separate Slater-Koster files: “H1-H1.skf”, “H1-H2.skf”, and “H2-H2.skf”.

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

  pcpot: PointCharge object

    An external point charge potential (for QM/MM calculations)