DFTB+

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

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

pcpot: PointCharge object

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

Contents