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=False, label='dftb', atoms=None, kpts=None, slako_dir=None, **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:

>>> 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)