Practical Examples

This page provides ready-to-use examples with configuration files and commands.

Configuration File Format 

DFTBaby uses INI format for configuration files (dftbaby.cfg):

Sections: [SectionName]
Settings: key = value
Comments: # comment text
Boolean: 0 (False) or 1 (True)

Example 1: Basic Excited State Calculation 

Goal: Calculate the first 10 excited states of a molecule with absorption spectrum.

File Structure 

my_calculation/
├── molecule.xyz
└── dftbaby.cfg

molecule.xyz 

10
Naphthalene
C   0.000000   0.000000   0.000000
C   1.397000   0.000000   0.000000
C   2.093000   1.207000   0.000000
C   1.397000   2.414000   0.000000
C   0.000000   2.414000   0.000000
C  -0.696000   1.207000   0.000000
H  -0.544000  -0.940000   0.000000
H   1.941000  -0.940000   0.000000
H  -0.544000   3.354000   0.000000
H   1.941000   3.354000   0.000000

dftbaby.cfg 

[DFTBaby]
# Ground state settings
charge = 0
multiplicity = 1
scf_conv = 1e-8

# Long-range correction for better CT states
long_range_correction = 1

# Number of excited states
nstates = 10

# Convergence for TD-DFTB
iter_conv = 1e-6

Run Command 

LR_TDDFTB.py molecule.xyz

Expected Output 

Excited State Energies:
S1:  3.245 eV   f=0.0234   (HOMO->LUMO)
S2:  3.876 eV   f=0.1567   (HOMO->LUMO+1)
S3:  4.123 eV   f=0.0012   (HOMO-1->LUMO)
...

Visualization 

Add --graphical=1 for interactive 3D visualization:

LR_TDDFTB.py molecule.xyz --graphical=1

Example 2: Absorption Spectrum Calculation 

Goal: Calculate absorption spectrum with many states and long-range correction.

dftbaby.cfg 

[DFTBaby]
charge = 0
multiplicity = 1
scf_conv = 1e-10

# Essential for accurate spectrum
long_range_correction = 1

# Calculate many states for broad spectrum
nstates = 50

# Tight convergence for accurate energies
iter_conv = 1e-8

# Active space (optional, for large molecules)
# nr_active_occ = 50
# nr_active_virt = 50

Run Command 

LR_TDDFTB.py molecule.xyz > spectrum_output.txt

Post-processing 

Extract data and plot:

grep "eV" spectrum_output.txt > energies.dat

# Plot with Python
python3 << 'EOF'
import matplotlib.pyplot as plt
import numpy as np

# Read energies and oscillator strengths
data = np.loadtxt('energies.dat', usecols=(1,3))
energies = data[:,0]
osc_strengths = data[:,1]

# Create spectrum with Gaussian broadening
x = np.linspace(2.0, 8.0, 1000)
y = np.zeros_like(x)
sigma = 0.2  # eV

for E, f in zip(energies, osc_strengths):
    y += f * np.exp(-(x-E)**2/(2*sigma**2))

plt.plot(x, y)
plt.xlabel('Energy (eV)')
plt.ylabel('Absorption (arb. units)')
plt.title('Absorption Spectrum')
plt.savefig('spectrum.png', dpi=300)
plt.show()
EOF

Example 3: Ground State Geometry Optimization 

Goal: Optimize molecular geometry to the nearest local minimum.

dftbaby.cfg 

[DFTBaby]
charge = 0
scf_conv = 1e-10

[GeometryOptimization]
# Optimize ground state
state = 0

# Coordinate system: 'cartesian' or 'internal'
coord_system = cartesian

# Convergence criteria
grad_tol = 1.0e-5
func_tol = 1.0e-8

# Maximum steps
max_steps = 500

# Optimization method: 'BFGS', 'CG', 'Newton', 'Steepest Descent'
method = BFGS

# Calculate vibrational frequencies after optimization
calc_hessian = 1

Run Command 

GeometryOptimization.py molecule.xyz

Output Files 

optimized.xyz: Optimized geometry
optimization.xyz: Full optimization trajectory
hessian.dat: Hessian matrix (if calc_hessian=1)
vib.molden: Vibrational modes (view with Molden)

Example 4: Excited State Geometry Optimization 

Goal: Find the optimized geometry of the first excited state (S1).

dftbaby.cfg 

[DFTBaby]
charge = 0
scf_conv = 1e-10
long_range_correction = 1
nstates = 5

[GeometryOptimization]
# Optimize S1 (first excited state)
state = 1

coord_system = cartesian
grad_tol = 1.0e-4
func_tol = 1.0e-7
max_steps = 300
method = BFGS

# Don't calculate Hessian for excited states (expensive)
calc_hessian = 0

Run Command 

GeometryOptimization.py molecule.xyz

Use Case: Emission Spectrum 

Optimize ground state (S0):

# Set state=0 in dftbaby.cfg
GeometryOptimization.py molecule.xyz
cp optimized.xyz s0_geometry.xyz

Optimize excited state (S1):

# Set state=1 in dftbaby.cfg
GeometryOptimization.py s0_geometry.xyz
cp optimized.xyz s1_geometry.xyz

Calculate vertical emission:
```
LR_TDDFTB.py s1_geometry.xyz
```
The emission energy is S1→S0 at S1 geometry.

Example 5: Surface Hopping Dynamics 

Goal: Simulate non-adiabatic dynamics with trajectory surface hopping.

File Structure 

dynamics/
├── molecule.xyz
├── dftbaby.cfg
└── dynamics.in

dftbaby.cfg 

[DFTBaby]
charge = 0
scf_conv = 1e-10
long_range_correction = 1

# IMPORTANT: Disable DIIS for dynamics
density_mixer = None

[SurfaceHopping]
# Start in S1 (first excited state)
initial_state = 1

# Calculate 3 excited states
nstates = 3

# Simulation time
nstep = 5000
nuclear_step = 0.1     # 0.1 fs time step

# Run at constant energy (NVE)
dyn_mode = "E"

# Force hop to ground state at conical intersection
switch_to_groundstate = 1

# Decoherence correction (recommended)
decoherence_correction = 1

# Scalar coupling threshold
scalar_coupling_threshold = 0.01

# Output settings
output_step = 10       # Save every 10th step

# Time series (optional)
time_series = ['particle-hole charges current', 'Lambda2 current']

dynamics.in 

Initial conditions with positions (Bohr) and velocities (a.u.).

Format: First the number of atoms and comment line, then atomic positions, followed by velocities in the same order:

 10
 c   0.0000   0.0000   0.0000
 c   2.6400   0.0000   0.0000
 c   3.9500   2.2800   0.0000
 c   2.6400   4.5600   0.0000
 c   0.0000   4.5600   0.0000
 c  -1.3100   2.2800   0.0000
 h  -1.0300  -1.7700   0.0000
 h   3.6700  -1.7700   0.0000
 h  -1.0300   6.3400   0.0000
 h   3.6700   6.3400   0.0000
 0.0001  -0.0002   0.0003
-0.0001   0.0004  -0.0001
 0.0002  -0.0001   0.0002
-0.0003   0.0001  -0.0001
 0.0001   0.0003   0.0002
-0.0002  -0.0002   0.0001
 0.0005   0.0001  -0.0003
-0.0004   0.0002   0.0004
 0.0003  -0.0004   0.0001
-0.0001   0.0003  -0.0002

Run Command 

cd dynamics/
SurfaceHopping.py

Output Files 

dynamics.xyz: Nuclear trajectory
energy_0.dat, energy_1.dat, ...: Energy of each electronic state vs. time
state.dat: Active state vs. time (for hopping analysis)
coeff_0.dat, coeff_1.dat, ...: Quantum populations for each state
particle_hole_charges.xyz: Charge distribution (if requested)

Analysis 

Plot state populations:

python3 << 'EOF'
import matplotlib.pyplot as plt
import numpy as np

# Read active state vs. time
data = np.loadtxt('state.dat')
time_fs = data[:, 0]
state = data[:, 1]

plt.figure(figsize=(10,6))
plt.plot(time_fs, state, linewidth=0.5)
plt.xlabel('Time (fs)')
plt.ylabel('Electronic State')
plt.title('Surface Hopping Trajectory')
plt.yticks([0,1,2,3], ['S0', 'S1', 'S2', 'S3'])
plt.grid(True, alpha=0.3)
plt.savefig('hopping.png', dpi=300)
plt.show()
EOF

Example 6: Ground State Equilibration 

Goal: Equilibrate geometry at finite temperature before excited state dynamics.

dftbaby.cfg 

[DFTBaby]
charge = 0
scf_conv = 1e-8

[SurfaceHopping]
# Stay on ground state
initial_state = 0
nstates = 0           # Don't calculate excited states

# Equilibration time
nstep = 10000
nuclear_step = 0.5    # Larger step for faster equilibration

# Constant temperature (Berendsen thermostat)
dyn_mode = "T"
temp = 300.0
timecoupling = 1.0

output_step = 50

Run Command 

SurfaceHopping.py

Then extract equilibrated geometry:

tail -11 dynamics.xyz | head -10 > equilibrated.xyz

Example 7: QM/MM Calculation for Molecular Crystals 

Goal: Calculate excited states in a periodic molecular crystal using QM/MM.

File Structure 

crystal/
├── crystal.xyz
├── crystal.ff
└── dftbaby.cfg

crystal.ff 

Force field definition with periodic box:

# Lattice vectors (Angstrom)
50.0  0.0  0.0
0.0  50.0  0.0
0.0   0.0  50.0

# Atom types and positions
C  0.0  0.0  0.0
C  1.4  0.0  0.0
H  2.0  0.9  0.0
...

dftbaby.cfg 

[DFTBaby]
charge = 0
scf_conv = 1e-10
long_range_correction = 1

# QM region (Python list syntax)
qmmm_partitioning = "[0,1,2,3,4,5,6,7,8,9]"

# Force field file
periodic_force_field = crystal.ff

# Excited states
nstates = 10

# Active space (important for large QM regions)
nr_active_occ = 50
nr_active_virt = 50

Run Command 

For excited states:

LR_TDDFTB.py crystal.xyz

For geometry optimization:

optimize.py crystal.xyz 0

(The QM/MM setup is read automatically from dftbaby.cfg)

Example 8: Constrained Geometry Optimization 

Goal: Optimize geometry while fixing certain bond lengths or angles.

constraints.txt 

# Fix C-C bond between atoms 0 and 1 to 1.54 Angstrom
bond 0 1 1.54

# Fix angle C-C-C (atoms 0,1,2) to 120 degrees
angle 0 1 2 120.0

# Fix dihedral angle (atoms 0,1,2,3) to 180 degrees
dihedral 0 1 2 3 180.0

dftbaby.cfg 

[DFTBaby]
charge = 0
scf_conv = 1e-10

[GeometryOptimization]
state = 0
coord_system = cartesian
grad_tol = 1.0e-5
max_steps = 500
method = BFGS

Run Command 

GeometryOptimization.py molecule.xyz --constraints=constraints.txt

Example 9: Charged Species Calculation 

Goal: Calculate excited states of a cation or anion.

dftbaby.cfg (Cation)

[DFTBaby]
# Total charge +1
charge = 1

# Doublet (one unpaired electron)
multiplicity = 2

scf_conv = 1e-10
long_range_correction = 1
nstates = 10

Run Command 

LR_TDDFTB.py cation.xyz

dftbaby.cfg (Anion)

[DFTBaby]
# Total charge -1
charge = -1

# Doublet
multiplicity = 2

scf_conv = 1e-10
long_range_correction = 1
nstates = 10

Example 10: High-Accuracy Calculation 

Goal: Maximize accuracy for publication-quality results.

dftbaby.cfg 

[DFTBaby]
charge = 0

# Very tight SCF convergence
scf_conv = 1e-14

# Disable DIIS for better convergence
density_mixer = None

# Long-range correction
long_range_correction = 1

# Dispersion correction (if needed)
dispersion_correction = 1

# Excited states
nstates = 20

# Very tight TD-DFTB convergence
iter_conv = 1e-10

# Increase maximum iterations if needed
max_iter = 200

Run Command 

LR_TDDFTB.py molecule.xyz > high_accuracy.log

Quick Reference: Common Settings 

DFTBaby Section 

Parameter	Default	Description
`charge`	0	Total molecular charge
`multiplicity`	1	Spin multiplicity (2S+1)
`scf_conv`	1e-6	SCF convergence threshold
`long_range_correction`	0	Enable LC-DFTB (0 or 1)
`dispersion_correction`	0	Enable Grimme’s D3 (0 or 1)
`nstates`	5	Number of excited states
`iter_conv`	1e-6	TD-DFTB convergence
`density_mixer`	“DIIS”	Mixing algorithm (None, DIIS, Pulay)

SurfaceHopping Section 

Parameter	Default	Description
`initial_state`	0	Starting electronic state
`nstates`	2	Number of excited states
`nstep`	1000	Number of MD steps
`nuclear_step`	0.1	Time step (fs)
`dyn_mode`	“E”	“E” (NVE), “T” (NVT), “M” (meta)
`temp`	300.0	Temperature (K) for NVT
`decoherence_correction`	0	Enable decoherence (0 or 1)
`output_step`	1	Output frequency

GeometryOptimization Section 

Parameter	Default	Description
`state`	0	Electronic state to optimize
`coord_system`	“cartesian”	“cartesian” or “internal”
`grad_tol`	1e-5	Gradient convergence
`func_tol`	1e-8	Energy convergence
`max_steps`	100000	Maximum optimization steps
`method`	“CG”	BFGS, CG, Newton, Steepest Descent
`calc_hessian`	0	Calculate frequencies (0 or 1)

Tips and Best Practices 

Start Simple
- Begin with default settings
- Add options only when needed
- Test with small molecules first
Convergence Issues
- Increase scf_conv and iter_conv
- Try density_mixer = None
- Check geometry for problems
Long Simulations
- Increase output_step to reduce file size
- Use smaller nstates if high states not needed
- Enable decoherence_correction for surface hopping
Accuracy vs. Speed
- Default settings: good balance
- High accuracy: tighter convergence, more states
- Fast screening: looser convergence, fewer states
Memory Management
- Use nr_active_occ and nr_active_virt for large molecules
- Reduce nstates if memory issues occur
- Consider QM/MM for very large systems

Complete Workflow Example 

Here’s a complete workflow for studying photodynamics:

Step 1: Ground State Optimization 

Create dftbaby.cfg:

[DFTBaby]
charge = 0
long_range_correction = 1

[GeometryOptimization]
state = 0
method = BFGS
calc_hessian = 1

Run:

GeometryOptimization.py molecule.xyz
cp optimized.xyz s0_min.xyz

Step 2: Vertical Excitation 

Update dftbaby.cfg:

[DFTBaby]
charge = 0
long_range_correction = 1
nstates = 20

Run:

LR_TDDFTB.py s0_min.xyz > absorption.log

Step 3: S1 Optimization 

Update dftbaby.cfg:

[GeometryOptimization]
state = 1

Run:

GeometryOptimization.py s0_min.xyz
cp optimized.xyz s1_min.xyz

Step 4: Emission Energy 

Run:

LR_TDDFTB.py s1_min.xyz > emission.log

Step 5: Non-adiabatic Dynamics 

Generate initial conditions:

# Create multiple Wigner samples
python3 << 'EOF'
from DFTB.Dynamics import WignerSampling
ws = WignerSampling('s1_min.xyz')
ws.sample(nsamples=50)
EOF

Update dftbaby.cfg:

[SurfaceHopping]
initial_state = 1
nstates = 3
nstep = 10000
nuclear_step = 0.1
dyn_mode = "E"
decoherence_correction = 1

Run trajectories:

for i in {1..50}; do
    mkdir traj_$i
    cd traj_$i
    cp ../dynamics_$i.in dynamics.in
    cp ../dftbaby.cfg .
    SurfaceHopping.py &
    cd ..
done
wait

Analyze results:

python3 << 'EOF'
import numpy as np
import matplotlib.pyplot as plt
from glob import glob

# Calculate population decay
trajectories = glob('traj_*/state.dat')
# Analysis code here...
EOF

Additional Resources 

Full example files: examples/ directory in repository
Command-line help: program.py --help
Python API: usage_guide.rst
Troubleshooting: usage_guide.rst troubleshooting section