index | modules | gitlab | page source

# 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.calculators.dftb import Dftb
    from ase.io import write
    from ase.build import molecule
    from ase.optimize import QuasiNewton

    atoms = molecule('H2O')
    calc = Dftb(atoms=atoms,
                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.calculators.dftb import Dftb
    from ase.io import write, read
    from ase.build import molecule

    atoms = molecule('H2O')
    calc = Dftb(atoms=atoms,
                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.io import read, write
    from ase.io.dftb import read_dftb_velocities, write_dftb_velocities
    from ase.calculators.dftb import Dftb
    from ase.build import molecule

    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(atoms=atoms,
                          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(atoms=atoms,
                          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, **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 (<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)