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)[source]#

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)

DFTB+

Contents