DftbPlus

Introduction

DftbPlus is a density-functional baed tight-binding code using atom-centered orbitals. This interface makes it possible to use DftbPlus 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

Before each DftbPlus calculation, please ensure that the following environment variables are set. DFTB_PREFIX should point to the directory where the Slater-Koster files are kept. DFTB_COMMAND should contain the path of the DftbPlus executable. This variable need not be set if a dftb+ executable can be found in the PATH.

  • bash:

    $ DFTB_PREFIX=/my_disk/my_name/lib/Dftb+sk/mio-0-1/  (an example)
    $ DFTB_COMMAND=~/bin/DFTB+/dftb+_s081217.i686-linux  (an example)
    
  • csh/tcsh:

    $ setenv DFTB_PREFIX /my_disk/my_name/lib/Dftb+sk/mio-0-1/  (an example)
    $ setenv DFTB_COMMAND ~/bin/DFTB+/dftb+_s081217.i686-linux   (an example)
    

DftbPlus Calculator (a FileIOCalculator)

The file ‘geo_end.gen’ contains the input and output geometry and it will be updated during the DftbPlus calculations.

All keywords to the DftbPlus calculator can be set by ASE.

Parameters

restart: str (default None)

If restart == None it is assumed that a new input file ‘dftb_hsd.in’ will be written by ase using default keywords and the ones given by the user.

If restart != None it is assumed that keywords are in file ‘restart’

ignore_bad_restart_file: bool (default False)

Ignore broken or missing restart file. By default, an error is raised if restart!=None and 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 the restart 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)

Example: 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.set_calculator(calc)

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

Example: Geometry Optimization by DftbPlus

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.set_calculator(calc)

calc.calculate(atoms)
final = read('geo_end.gen')
write('final.xyz', final)

Example: NVE md followed by NVT md (both by DftbPlus)

This is unphysical because of at least two reasons:

  • oxygen does not have spin here

  • 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.set_calculator(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.set_calculator(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