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)
orNone
: 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 usingase.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)
orNone
: Gamma-point only(n1,n2,n3)
: Monkhorst-Pack griddict
: Interpreted as a path in the Brillouin zone if it contains the ‘path’ keyword. Otherwise it is converted into a Monkhorst-Pack grid usingase.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)