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