DftbPlus

Introduction

DftbPlus is a density-functional based tight-binding code using atom centered orbitals. This interface makes it possible to use DftbPlus as a calculator in ASE. You need to register at DftbPlus site to download the code. Additionally you need Slater-Koster files for the combination of atom types of your system. These can be obtained at dftb.org.

Environment variables

Set environment variables in your configuration file (what is the directory for the Slater-Koster files and what is the name of the executable):

  • 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 ‘geom.out.gen’ contains the input and output geometry and it will be updated during the dftb calculations.

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

All Keywords to the dftb 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 defauls, it is an error if the restart file is missing or broken.
label: str (default ‘dftb’)
Name used for all files. May contain a directory.
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): Gamma-point
  • (n1,n2,n3): Monkhorst-Pack grid
  • (n1,n2,n3,'gamma'): Shifted Monkhorst-Pack grid that includes \(\Gamma\)
  • [(k11,k12,k13),(k21,k22,k23),...]: Explicit list in units of the reciprocal lattice vectors
  • kpts=3.5: \(\vec k\)-point density as in 3.5 \(\vec k\)-points per Å\(^{-1}\).
run_manyDftb_steps: bool (default False)
If True the Dftb calculator is running many steps by its own. If False all the relaxations/ molecular dynamis is done by ASE

Example: Geometry Optimization by ASE

from ase.calculators.dftb import Dftb
from ase.optimize import QuasiNewton
from ase.io import write

from ase.build import molecule
test = molecule('H2O')
test.set_calculator(Dftb(label='h2o',
                         atoms=test,
                         Hamiltonian_MaxAngularMomentum_='',
                         Hamiltonian_MaxAngularMomentum_O='"p"',
                         Hamiltonian_MaxAngularMomentum_H='"s"',
                         ))

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

Example: Geometry Optimization by DFTB

from ase.calculators.dftb import Dftb
from ase.io import write, read

from ase.build import molecule
system = molecule('H2O')
calc = Dftb(label='h2o', atoms=system,
            run_manyDftb_steps=True,
            Driver_='ConjugateGradient',
            Driver_MaxForceComponent='1E-4',
            Driver_MaxSteps=1000,
            Hamiltonian_MaxAngularMomentum_='',
            Hamiltonian_MaxAngularMomentum_O='"p"',
            Hamiltonian_MaxAngularMomentum_H='"s"')
system.set_calculator(calc)
calc.calculate(system)
final = read('geo_end.gen')
write('test.final.xyz', final)

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

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.calculators.dftb import Dftb
from ase.build import molecule
from ase.md.verlet import VelocityVerlet
from ase.md import MDLogger
from ase.units import fs
from ase.io.dftb import read_dftb_velocities, write_dftb_velocities
from ase.io import read, write

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]))
test = o2 + h2_1 + h2_2

# 1fs = 41.3 au
# 1000K = 0.0031668 au
calculator_NVE = Dftb(label='h2o',
                      atoms=test,
                      run_manyDftb_steps=True,
                      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',
    atoms=test,
    run_manyDftb_steps=True,
    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 deg Celcius
    # Driver_Thermostat_Temperature=0.0, # 0 deg Kelvin
    Driver_Thermostat_CouplingStrength=0.01)

write_dftb_velocities(test, 'velocities.txt')
os.system('rm md.log.* md.out* geo_end*xyz')
test.set_calculator(calculator_NVE)
dyn = VelocityVerlet(test, 0.000 * fs)  # fs time step.
dyn.attach(MDLogger(dyn, test, 'md.log.NVE', header=True, stress=False,
                    peratom=False, mode='w'), interval=1)
dyn.run(1)  # run NVE ensemble using DFTB's own driver
test = read('geo_end.gen')
write('test.afterNVE.xyz', test)

read_dftb_velocities(test, filename='geo_end.xyz')
write_dftb_velocities(test, 'velocities.txt')

os.system('mv md.out md.out.NVE')
os.system('mv geo_end.xyz geo_end_NVE.xyz')

test.set_calculator(calculator_NVT)
os.system('rm md.log.NVT')
dyn.attach(MDLogger(dyn, test, 'md.log.NVT', header=True, stress=False,
                    peratom=False, mode='w'), interval=1)
dyn.run(1)  # run NVT ensemble using DFTB's own driver
test = read('geo_end.gen')
read_dftb_velocities(test, filename='geo_end.xyz')

os.system('mv md.out md.out.NVT')
os.system('mv geo_end.xyz geo_end_NVT.xyz')
write_dftb_velocities(test, 'velocities.txt')

# to watch:
#  ase gui geo_end_NVE.xyz geo_end_NVT.xyz