DFTBaby Wiki
Brought to you by:
quack-quack
Menu
▾
▴
Home
Tight-binding DFT is a semiempirical electron structure method based on DFT.
Features
- Electronic Prametrization
- Pseudo-orbitals can be generated from atomic DFT calculations
using non-hybrid functionals from the libxc-library - Slater-Koster files for the hamiltonian, the overlap, and dipoles
can be generated from the pseudo-orbitals
- Pseudo-orbitals can be generated from atomic DFT calculations
- Repulsive Potentials
- scaled repulsive potentials of Hotbit
- DFTB
- long-range correction
- Grimme's dispersion correction (experimental)
- TD-DFTB
- excited states with and without long-range correction
- analytical excited state gradients
- QM/MM
- QM part treated with TD-DFTB, MM part with UFF of Gaussian 09
- QM/MM calculations of molecular crystals (with periodic boundary conditions) using
rudimentary implementation of the DREIDING force field
- Analysis
- cube files for molecular orbitals and transition and difference densities of excited states
Installation
You can find the latest version of DFTB on its Sourceforge page
After downloading the release archive the DFTB module can be installed locally by typing
$ tar -xvf DFTBaby-#.#.#.tar.gz
$ cd DFTBaby-#.#.#/
$ python setup.py install --user
You will find more information in the README.txt located in the top folder.
If the installation succeeds you should be able to execute
$ DFTB2.py --help
All programs in the DFTBaby package come with a --help option that shows all the available command line options
with a short explanation. The options can be passed on the command line or can be placed in a configuration file
called dftbaby.cfg in the current directory.
Alternatively, DFTBaby can be installed by adding the folder DFTBaby-#### to the PYTHONPATH
environment variable. The subfolders DFTB/ , DFTB/Analyse , DFTB/Modeling and DFTB/Formats
should be added to the PATH environment variable. This setup has the advantage that you
can modify the code and use it without reinstalling it repeatedly.
The script 'dftbaby_env.sh' contains the necessary lines of code to set those environment
variables and can be sourced from the .bashrc profile.
If you are using environment modules you can use the file 'dftbaby.modulefile' as a template.
In some cases you will have to compile the extensions for the architecture of your machine
with f2py:
cd DFTB/extensions/
make
Programs
DFTBaby provides the module DFTB that can be imported from python. In addition
the following programs (among other scripts) can be called directly from the command line:
- DFTB2.py - only ground state calculation, MOs
- LR_TDDFTB.py - ground and excited state calculations, gradients
- optimize.py - optimizes geometries on ground and excited states
All programs have a help function that shows all available command line options, e.g.:
$ LR_TDDFTB.py --help
Graphical Analysis Tool
The results of a TD-DFTB calculation can be visualized graphically. To start the graphical user interface
at the end of a calculation it is enough to add the option --graphical=1 when calling LR_TDDFTB.py.
The graphical interface requires 'Mayavi' which can be obtained from http://docs.enthought.com/mayavi/mayavi/
or by installing Enthought's Canopy.
Limitations and Bugs
The following known limitations should be kept in mind when using (lc)-TD-DFTB:
-
DFTB uses a minimal basis set of valence orbitals. Being numerically exact
atomic orbitals obtained from a DFT calculations, the pseudo-orbitals are much better
than for instance, STO-3G basis functions. However, Rydberg states, which are frequent
in organic conjugated molecules, cannot be described in this way. -
The Coulomb interaction between density fluctuations around the neutral reference density
is calculated using the monompole approximation. In reality, the charge fluctuations around an atom
are not spherically symmetric. In my experience the monopole approximation leads to occupied Kohn-Sham orbitals that are too high in energy. Excitations from these orbitals to unoccupied ones produce spurious low-lying
states, that are dark.
Adding the interaction between Mulliken monopoles and dipoles should in principle fix this problem, however
the Mulliken dipoles are too large and lead to SCF convergence problems.
Not all features work as expected and probably there are still some bugs. If you have any doubt, I am happy to try answer your questions. Please write to alexander.humeniuk@gmail.com.
Usage
Below the usage of the program is illustrated by a few simple calculations.
Ground state calculations
The program DFTB2.py constructs the semiempirical DFTB2 hamiltonian from Slater-Koster files
and equilibrates the partial charges self-consistently. The only command line argument expected
is the name of a file with the molecular geometry in the xyz-format.
Create the file 'benzene.xyz' with the following geometry:
12
benzene
H 0.000000 2.484212 0.000000
H 0.000000 -2.484212 0.000000
H 2.151390 1.242106 0.000000
H -2.151390 -1.242106 0.000000
H -2.151390 1.242106 0.000000
H 2.151390 -1.242106 0.000000
C 0.000000 1.396792 0.000000
C 0.000000 -1.396792 0.000000
C 1.209657 0.698396 0.000000
C -1.209657 -0.698396 0.000000
C -1.209657 0.698396 0.000000
C 1.209657 -0.698396 0.000000
Now execute the DFTB2.py program
$ DFTB2.py benzene.xyz | tee /tmp/dftb.out
For each SCF iteration the program will print the Mulliken charges and energies of the molecular orbitals
...
*******************
* Iteration: 37 *
*******************
Convergence
===========
relative change: 8.417431e-11
Electronic Energies
===================
band structure energy E_bs: -13.0260674 hartree -354.4574798 eV
Coulomb energy E_coulomb : 0.0008605 hartree 0.0234163 eV
long range HF-x E_x : -1.0064646 hartree -27.3873081 eV
total electronic energy : -14.0316715 hartree -381.8213716 eV
repulsion energy : 0.8358768 hartree 22.7453747 eV
total energy : -13.1957947 hartree -359.0759969 eV
HOMO-LUMO gap : 0.3818691 hartree 10.3911927 eV
Orbital Energies
================
0: -0.5341218 hartree -14.5342008 eV (2.000 e)
1: -0.5328525 hartree -14.4996595 eV (2.000 e)
2: -0.5328096 hartree -14.4984937 eV (2.000 e)
3: -0.5263937 hartree -14.3239083 eV (2.000 e)
4: -0.5263200 hartree -14.3219007 eV (2.000 e)
5: -0.5225932 hartree -14.2204908 eV (2.000 e)
6: -0.4904442 hartree -13.3456711 eV (2.000 e)
7: -0.4775354 hartree -12.9944059 eV (2.000 e)
8: -0.4741544 hartree -12.9024028 eV (2.000 e)
9: -0.4588021 hartree -12.4846469 eV (2.000 e)
10: -0.4587841 hartree -12.4841559 eV (2.000 e)
11: -0.4255102 hartree -11.5787278 eV (2.000 e)
12: -0.4255102 hartree -11.5787256 eV (2.000 e)
13: -0.3848252 hartree -10.4716322 eV (2.000 e)
14: -0.3848233 hartree -10.4715806 eV (2.000 e)
15: -0.0029542 hartree -0.0803879 eV (0.000 e)
16: -0.0029522 hartree -0.0803333 eV (0.000 e)
17: 0.1380796 hartree 3.7573385 eV (0.000 e)
18: 0.2853296 hartree 7.7642163 eV (0.000 e)
19: 0.3334102 hartree 9.0725563 eV (0.000 e)
20: 0.3334376 hartree 9.0733023 eV (0.000 e)
21: 0.4609968 hartree 12.5443657 eV (0.000 e)
22: 0.4775578 hartree 12.9950144 eV (0.000 e)
23: 0.4776308 hartree 12.9970018 eV (0.000 e)
24: 0.7182194 hartree 19.5437533 eV (0.000 e)
25: 0.7184974 hartree 19.5513173 eV (0.000 e)
26: 0.7488216 hartree 20.3764803 eV (0.000 e)
27: 0.7492221 hartree 20.3873784 eV (0.000 e)
28: 0.9148878 hartree 24.8953737 eV (0.000 e)
29: 1.0419351 hartree 28.3525094 eV (0.000 e)
Mulliken Charges (Partial Charges)
NOTE: charges are in units of e-!
=========================================================================
h-0 : 0.968 (-0.032)
h-1 : 0.968 (-0.032)
h-2 : 0.968 (-0.032)
h-3 : 0.968 (-0.032)
h-4 : 0.968 (-0.032)
h-5 : 0.968 (-0.032)
c-6 : 4.032 (+0.032)
c-7 : 4.032 (+0.032)
c-8 : 4.032 (+0.032)
c-9 : 4.032 (+0.032)
c-10 : 4.032 (+0.032)
c-11 : 4.032 (+0.032)
=========================================================================
sum of charges : +30.000
sum of partial charges: +0.000
!!! CONVERGED after 38 iterations (relative change = 8.42e-11 < 1.00e-10 = threshold)!!!
TOTAL ENERGY OF STRUCTURE 0 : -13.1957947 hartree -359.0759969 eV
It is important to note that the counting of the MOs starts at 0, therefore the HOMO is the orbital #14 and not #15 as one
would expect for a molecule with 30 valence electrons. Also note that Mulliken charges are given in units of the electronic charge,
so a partial Mulliken charge of -0.032 on the hydrogens means that they are actually slightly positively charged as compared
to the carbon atoms.
If you wish to also see the MO coefficients in each iteration, you should run
$ DFTB2.py benzene.xyz --verbose=2 | tee /tmp/dftb.out
...
| MO10 | MO11 | MO12 | MO13 | MO14 | MO15 | MO16 | MO17 | MO18 | MO19 |
------------------------------------------------------------------------------------------------------------------------------------
en. (a.u.)| -0.4588 | -0.4255 | -0.4255 | -0.3848 | -0.3848 | -0.0030 | -0.0030 | 0.1381 | 0.2853 | 0.3334 |
occ. | 2.0000 | 2.0000 | 2.0000 | 2.0000 | 2.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 |
------------------------------------------------------------------------------------------------------------------------------------
h-0 1s | -0.0611 | 0.2747 | 0.0728 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | -0.4241 | -0.5476 |
h-1 1s | 0.0611 | 0.2747 | 0.0728 | -0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | -0.4241 | 0.5476 |
h-2 1s | 0.1668 | -0.2004 | 0.2014 | 0.0000 | -0.0000 | -0.0000 | 0.0000 | -0.0000 | -0.4241 | -0.3997 |
h-3 1s | -0.1668 | -0.2004 | 0.2014 | 0.0000 | -0.0000 | -0.0000 | 0.0000 | 0.0000 | -0.4241 | 0.3997 |
h-4 1s | -0.2271 | -0.0743 | -0.2742 | -0.0000 | -0.0000 | 0.0000 | -0.0000 | 0.0000 | -0.4242 | -0.1481 |
h-5 1s | 0.2271 | -0.0743 | -0.2742 | -0.0000 | 0.0000 | 0.0000 | -0.0000 | -0.0000 | -0.4242 | 0.1481 |
c-6 2s | 0.0357 | -0.0023 | -0.0006 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.2200 | 0.1739 |
c-6 2py | 0.0937 | -0.2550 | -0.0675 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | -0.2732 | -0.4701 |
c-6 2pz | -0.0000 | 0.0000 | -0.0000 | -0.1448 | -0.5192 | 0.1560 | 0.6130 | 0.4889 | -0.0000 | -0.0000 |
c-6 2px | 0.1267 | -0.0807 | 0.3050 | 0.0000 | -0.0000 | 0.0000 | -0.0000 | 0.0000 | -0.0000 | -0.0432 |
c-7 2s | -0.0357 | -0.0023 | -0.0006 | -0.0000 | -0.0000 | -0.0000 | 0.0000 | -0.0000 | 0.2200 | -0.1739 |
c-7 2py | 0.0937 | 0.2550 | 0.0675 | -0.0000 | 0.0000 | 0.0000 | -0.0000 | 0.0000 | 0.2732 | -0.4701 |
c-7 2pz | -0.0000 | -0.0000 | 0.0000 | 0.1448 | 0.5192 | 0.1560 | 0.6130 | -0.4889 | -0.0000 | 0.0000 |
c-7 2px | 0.1267 | 0.0807 | -0.3050 | -0.0000 | 0.0000 | -0.0000 | -0.0000 | -0.0000 | 0.0000 | -0.0432 |
c-8 2s | -0.0978 | 0.0017 | -0.0017 | 0.0000 | -0.0000 | -0.0000 | 0.0000 | -0.0000 | 0.2202 | 0.1270 |
c-8 2py | -0.2082 | 0.2867 | 0.0991 | -0.0000 | 0.0000 | 0.0000 | -0.0000 | -0.0000 | -0.1367 | -0.2752 |
c-8 2pz | -0.0000 | 0.0000 | -0.0000 | 0.3772 | -0.3851 | 0.4528 | -0.4416 | -0.4889 | 0.0000 | 0.0000 |
c-8 2px | -0.1754 | 0.0493 | -0.2732 | -0.0000 | 0.0000 | 0.0000 | -0.0000 | -0.0000 | -0.2365 | -0.2372 |
c-9 2s | 0.0978 | 0.0017 | -0.0017 | 0.0000 | -0.0000 | -0.0000 | -0.0000 | 0.0000 | 0.2202 | -0.1270 |
c-9 2py | -0.2082 | -0.2867 | -0.0991 | -0.0000 | 0.0000 | -0.0000 | 0.0000 | -0.0000 | 0.1367 | -0.2752 |
c-9 2pz | -0.0000 | -0.0000 | 0.0000 | -0.3772 | 0.3851 | 0.4528 | -0.4416 | 0.4889 | 0.0000 | -0.0000 |
c-9 2px | -0.1754 | -0.0493 | 0.2732 | 0.0000 | -0.0000 | -0.0000 | -0.0000 | -0.0000 | 0.2365 | -0.2372 |
c-10 2s | 0.1339 | 0.0006 | 0.0023 | -0.0000 | -0.0000 | 0.0000 | -0.0000 | 0.0000 | 0.2204 | 0.0471 |
c-10 2py | 0.2039 | 0.2981 | 0.0557 | 0.0000 | 0.0000 | 0.0000 | -0.0000 | 0.0000 | -0.1365 | -0.2047 |
c-10 2pz | -0.0000 | -0.0000 | 0.0000 | -0.5221 | -0.1342 | -0.6088 | -0.1714 | -0.4889 | 0.0000 | 0.0000 |
c-10 2px | -0.2855 | 0.0924 | -0.2617 | -0.0000 | -0.0000 | 0.0000 | -0.0000 | -0.0000 | 0.2365 | 0.0286 |
c-11 2s | -0.1339 | 0.0006 | 0.0023 | -0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.2204 | -0.0471 |
c-11 2py | 0.2039 | -0.2981 | -0.0558 | 0.0000 | -0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.1365 | -0.2047 |
c-11 2pz | -0.0000 | 0.0000 | -0.0000 | 0.5221 | 0.1342 | -0.6088 | -0.1714 | 0.4889 | 0.0000 | -0.0000 |
c-11 2px | -0.2855 | -0.0924 | 0.2617 | 0.0000 | -0.0000 | -0.0000 | 0.0000 | -0.0000 | -0.2365 | 0.0286 |
...
Since for larger molecules the output can be huge it should be redirected to a file in the scratch or tmp directory.
To check the progress in the SCF-convergence type
$ grep "relative change" /tmp/dftb.out
Visualizing Orbitals
The converged Kohn-Sham orbitals can be written to a file that can be visualized using Molden.
$ DFTB2.py benzene.xyz --molden_file=benzene.molden --verbose=0
$ molden benzene.molden
Since DFTB uses numerical atomic orbitals, while Molden requires Gaussian-type orbitals, the orbitals are converted into the STO-6G basis. Also note that molden starts counting orbitals from 1!
To visualize the exact orbitals, you should write the orbitals of interest to cube-files:
The command
$ DFTB2.py benzene.xyz --cubedir=. --points_per_bohr=5.0 --verbose=0
will create the cube files 'benzene_HOMO.cube' and 'benzene_LUMO.cube'. For an explanation on how to write out orbitals other than the frontier orbitals, add the --help option:
$ DFTB2.py benzene.xyz --cubedir=. --points_per_bohr=5.0 --verbose=0 --help
and search for the relevant section
Cube-Files:
--cubedir=CUBEDIR directory where cube files are stored [default: none]
--cube_orbitals=CUBE_ORBITALS
list indeces of molecular orbitals which should be
saved to a cube file, e.g. '[0,1,2]' (if None, the
HOMO-1, HOMO, LUMO and LUMO+1 are saved) [default: []]
--points_per_bohr=POINTS_PER_BOHR
resolution of grid for cube files [default: 2.0]
--cube_threshold=CUBE_THRESHOLD
atomic orbitals with coefficients below
threshold*(largest coefficient) are neglected when
evaluating MO amplitudes on a cube grid., [default:
0.001]
To avoid writing long sequences of command line arguments, you can put all options into a file 'dftbaby.cfg':
[DFTBaby]
cubedir=.
points_per_bohr=5.0
verbose=0
Now executing
$ DFTB2.py benzene.xyz
in the same folder where 'dftbaby.cfg' is located will also create the desired cube files.
SCF-Convergence Issues
The most common problem encountered is that the SCF calculation does not converge. To solve this issue you should
-
check that the input geometry is reasonable. A missing hydrogen or a fragmented molecule can be the culprits. Also if you happen to hit a conical intersection between the ground and an excited state, the SCF calculation can hop back and forth between two solutions.
-
inspect the HOMO-LUMO gap which is printed at each iteration. If the HOMO-LUMO gap drops below a threshold, the virtual orbitals are shifted up which usually resolves the issue automatically.
- Keep in mind, that convergence issues point to a serious problem with using the DFTB method for the particular molecule. Even is you can enforce convergence by using some combination of level-shifting, DIIS and linear-mixing of the densities at the early stages of the SCF-cycle (see the section 'SCF-Convergence' in the help menue), this solution is probably not reasonable.
Closed Shell Ions
To compute the electronic structure of a closed-shell ion (cation or anion) the charge has to be set inside the xyz-file on the comment line.
For the Ag5+ cation the xyz-file could be:
5
charge=+1
Ag -5.106663 2.600666 -0.719950
Ag -3.294714 3.534874 0.875566
Ag -3.035389 1.106012 -0.096583
Ag -2.540480 -1.332515 0.751398
Ag -1.196625 -0.380320 -1.246348
Excited states
Excited state calculations are performed with the program 'LR_TDDFTB.py', which first performs an SCF-calculation and then solves the linear response equations of TD-DFTB. By default, the full TD-DFTB matrix is constructed and diagonalized. While this is feasible for small molecules, for larger ones you should use the Davidson algorithm which iteratively determines the lowest excited states only. The iterative diagonalization is requested by adding the option --nstates with the number of states.
To compute the lowest 6 states of benzene type:
$ LR_TDDFTB.py benzene.xyz --nstates=6 --diag_conv=1.0e-8
At the end of the output you should find a table with a summary of the properties of the excited states:
Excited States
==============
Spn exc. en. /hartree exc. en. / eV exc. en.y / nm osc. strength dominant contrib. en_KS / eV transition dipole [TDx TDy TDz] /\ diagn.
1 S : 0.1696006 hartree 4.6150703 eV 268.6508244 nm 0.0000000 14->15(+7.069e-01),13->16(-7.069e-01) 10.3912 [-0.0001003 -0.0000575 -0.0000000] 0.8111
2 S : 0.1713666 hartree 4.6631237 eV 265.8823805 nm 0.0000000 13->15(+7.022e-01),14->16(+7.022e-01) 10.3912 [+0.0000068 -0.0000062 -0.0000000] 0.9997
3 S : 0.2140494 hartree 5.8245820 eV 212.8637613 nm 0.0013724 12->15(-7.039e-01),11->16(-6.998e-01) 11.4983 [+0.0000000 +0.0000000 +0.0980699] 0.9160
4 S : 0.2152690 hartree 5.8577711 eV 211.6577128 nm 0.0000000 12->16(+5.548e-01),11->15(+5.502e-01) 11.4984 [+0.0000000 +0.0000000 -0.0000478] 0.8752
5 S : 0.2152692 hartree 5.8577753 eV 211.6575605 nm 0.0000000 11->16(-5.540e-01),12->15(+5.509e-01) 11.4984 [-0.0000000 -0.0000000 -0.0000580] 0.8920
6 S : 0.2165301 hartree 5.8920866 eV 210.4250174 nm 0.0000000 11->15(+7.006e-01),12->16(-6.946e-01) 11.4983 [-0.0000000 -0.0000000 -0.0001795] 0.8508
-
The first columns indicate the spin (S = Singlet, T = Triplet) and symmetry (if available)
-
The next columns give the excitation energy in different units (Hartree, eV and nm) followed by the oscillator strength.
-
The column 'dominant contrib.' lists the dominant single excitations present in an excited state with their phases. Note that orbitals are counted from 0, again. So 14->15 is a HOMO->LUMO excitation.
-
The column 'en_KS' lists the energy difference between the Kohn-Sham orbital energies of the dominant excitation. If this value is equal to the excitation energy, the coupling matrix in TD-DFTB did not shift the energy of this excitation.
-
The column 'transition dipole' shows the transtion dipole vector.
-
The last column contains a metric for detecting charge transfer states, a value close to 0 (< 0.3) indicates a charge transfer state, while a value close to 1.0 is indicative of a valence excitation.
You can save the excitation energies and oscillator strengths as a table to a file with the option --spectrum_file
$ LR_TDDFTB.py benzene.xyz --nstates=6 --spectrum_file=benzene.spec
Symmetry
The symmetry of excited states can be determined by looking at the transformation properties of the transition densities and comparing the sign changes with the character tables. This only works for non-degenerate representations. The following symmetry groups can be detected: C1, Cs, C2v, C2h, D2h, D3h, D4h, D6h, D2d.
Using symmetry does not speed up the calculation, quite the opposite, also the molecule will be reoriented.
$ LR_TDDFTB.py benzene.xyz --nstates=6 --diag_conv=1.0e-8 --use_symmetry=1
...
Excited States
==============
Spn D6h exc. en. /hartree exc. en. / eV exc. en.y / nm osc. strength dominant contrib. en_KS / eV transition dipole [TDx TDy TDz] /\ diagn.
1 S B1u: 0.1696006 hartree 4.6150703 eV 268.6508244 nm 0.0000000 14->15(-7.069e-01),13->16(-7.069e-01) 10.3912 [+0.0000000 -0.0001003 +0.0000575] 0.8111
2 S B2u: 0.1713666 hartree 4.6631237 eV 265.8823805 nm 0.0000000 13->15(+7.022e-01),14->16(-7.022e-01) 10.3912 [-0.0000000 +0.0000068 +0.0000097] 0.9997
3 S A2u: 0.2140494 hartree 5.8245820 eV 212.8637613 nm 0.0013725 12->15(-7.039e-01),11->16(-6.998e-01) 11.4983 [+0.0980732 -0.0000000 +0.0000000] 0.9160
4 S ?: 0.2152690 hartree 5.8577711 eV 211.6577128 nm 0.0000000 12->16(-5.548e-01),11->15(-5.502e-01) 11.4984 [+0.0000478 -0.0000000 +0.0000000] 0.8752
5 S ?: 0.2152692 hartree 5.8577753 eV 211.6575605 nm 0.0000000 11->16(-5.540e-01),12->15(+5.509e-01) 11.4984 [-0.0000580 +0.0000000 -0.0000000] 0.8920
6 S A1u: 0.2165301 hartree 5.8920866 eV 210.4250174 nm 0.0000000 11->15(-7.006e-01),12->16(+6.946e-01) 11.4983 [+0.0001795 -0.0000000 +0.0000000] 0.8508
Note how the symmetry of the degenerate 4th and 5th states could not be determined.
Active Space
For large molecules or if computations for similar geometries have to be performed many times as in an MD simulation,
a significant speed-up can be achieved by restricting the excitations to a space of active orbitals.
This does not change results, if the lowest excited states are contained
in a subspace of excitations around the HOMO-LUMO gap.
As an extreme example consider S1 state of benzene.
$ LR_TDDFTB.py benzene.xyz --nstates=1
The section "Active Space" in the output informs you that you could have limited the active
space to only 4 frontier orbitals:
Active Space
============
State 1: contains excitations from window occ 13 - virt 16
The suggested active space for describing the lowest 1 states consists of
the highest 2 occupied and the lowest 2 virtual orbitals.
This active space accounts for 0.9999900 of the total probability for each state.
You should check that the active space contains all neighbouring (almost) degenerate orbitals.
Dimension of active space: 4
Dimension of full space: 225
In subsequent calculations you can then restrict the active space to excitations from the highest 2 occupied to the lowest 2 unoccupied orbitals:
$ LR_TDDFTB.py benzene.xyz --nstates=1 --nr_active_occ=2 --nr_active_virt=2
It is important to check that the active space is large enough, otherwise the results are completely arbitrary.
Visualizing Transition densities
Transition densities can be saved to cube files.
The following command saves the transition densities for the lowest 3 states as cube files in the current directory:
$ LR_TDDFTB.py benzene.xyz --nstates=6 --cubedir=. --cube_states="[1,2,3]"
Analytic gradients
Analytic gradients can be calculated using the auxilary functional method developed by Furche.
The following command computes the gradient on the S1 state:
LR_TDDFTB.py benzene.xyz --nstates=6 --gradient_state=1 --gradient_file=grad.xyz
The xyz-file 'grad.xyz' contains the forces acting on each atom in atomic units:
12
Gradient of total energy of state 1 units=au
h 0.000003714972337 0.011424857772237 0.000000000000000
h -0.000003714977431 -0.011424864937982 -0.000000000000000
h 0.009909015691574 0.005697788230926 -0.000000000000000
h -0.009909015900698 -0.005697788853081 0.000000000000000
h -0.009901920585340 0.005721201254311 -0.000000000000000
h 0.009901920805706 -0.005721201871464 0.000000000000000
c -0.000013020265626 -0.023010988607407 -0.000000000000003
c 0.000013020297149 0.023011018910976 0.000000000000003
c -0.019966502699402 -0.011502269134064 0.000000000000000
c 0.019966505175863 0.011502258167657 -0.000000000000000
c 0.019976498416076 -0.011548286500111 -0.000000000000002
c -0.019976500930213 0.011548275491433 0.000000000000002
Using DFTB with Python
Being written in python, you can integrate the DFTB module in other python programs.
# import module for excited state calculations
from DFTB.LR_TDDFTB import LR_TDDFTB
# geometry of a N2 molecule
atomlist = [(7,(0.0, 0.0, 0.0)), # N atom at position 0,0,0 bohr
(7,(0.0, 0.0, 1.0))] # N atom at position 0,0,1 bohr
tddftb = LR_TDDFTB(atomlist)
tddftb.setGeometry(atomlist, charge=0)
options={"nstates": 2} # 2 excited states
tddftb.getEnergies(**options) # perform calculation
print "Excited state energies / Hartree"
print tddftb.getExcEnergies()
For a more complicated example, have a look at the script 'DFTB/optimize.py'
Parallel Runs
Parts of the DFTB program are parallelized using OpenMP and Intel's Math Kernel Library (MKL).
On a shared memory machine, you can set the number of processors via the environment variables OMP_NUM_THREADS and MKL_NUM_THREADS,
which should equal the number of processors requested in the PBS submit script:
# request 4 processors on a single node
#PBS -l nodes=1:ppn=4
...
export OMP_NUM_THREADS=4
export MKL_NUM_THREADS=4
Requesting more that 4 processors does not speed up the calculations considerably.
In particular when running MD simulations on the cluster, the number of processors
should be specified explicitly using these two environment variables,
since by default all available processors would be used.