APS 2026 Short Course Tutorial 3: Excited-state geometry optimization using time-dependent density functional theory (TDDFT)

This tutorial outlines how to perform excited-state geometry optimization using the excitation energy and analytical nuclear forces obtained from time-dependent density functional theory (TDDFT). Using the formaldehyde molecule as an example, we will relax the atomic geometry of the \(^1A''\) excited state by combining DFT calculations using the pw.x code in Quantum ESPRESSO, and TDDFT calculations using the wbse.x code in WEST. More details about the implementation of TDDFT in WEST, for spin-conserving and spin-flip excitations with analytical forces and hybrid functionals, may be found in Jin et al., J. Chem. Theory Comput. 19, 8689 (2023).

A TDDFT excited-state geometry optimization is performed in four steps:

  1. Perform a ground-state DFT calculation.

  2. Compute spin-conserving or spin-flip excitation energy and excited-state forces for a target excited state using TDDFT.

  3. Update the atomic positions using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm.

  4. Check for convergence. If not converged, return to step 1 with the updated geometry.

These steps are managed by the WESTpy Python package.

1. Prepare input files of starting geometry

Download the following files to your working directory:

[1]:

%%bash
wget -N -q http://www.quantum-simulation.org/potentials/sg15_oncv/upf/C_ONCV_PBE-1.2.upf
wget -N -q http://www.quantum-simulation.org/potentials/sg15_oncv/upf/H_ONCV_PBE-1.2.upf
wget -N -q http://www.quantum-simulation.org/potentials/sg15_oncv/upf/O_ONCV_PBE-1.2.upf
wget -N -q https://west-code.org/doc/training/formaldehyde/es/pw.in -O pw.in
wget -N -q https://west-code.org/doc/training/formaldehyde/es/wbse.in -O wbse.in

The pw.in file contains the initial atomic coordinates for the excited state geometry optimization:

[2]:

%%bash
cat pw.in

&CONTROL
calculation = 'scf'
pseudo_dir = './'
tprnfor = .true.
/
&SYSTEM
ibrav = 0
nat = 4
ntyp = 3
ecutwfc = 50
nbnd = 10
/
&ELECTRONS
conv_thr = 1.d-8
diago_full_acc = .true.
startingpot = 'file'
startingwfc = 'file'
/
ATOMIC_SPECIES
C  12.0107  C_ONCV_PBE-1.2.upf
H   1.0079  H_ONCV_PBE-1.2.upf
O  16.00    O_ONCV_PBE-1.2.upf
K_POINTS gamma
CELL_PARAMETERS angstrom
15 0 0
0 15 0
0 0 15
ATOMIC_POSITIONS crystal
C  0.46600000  0.500000000  0.500000000
H  0.42652967  0.436863407  0.500000000
H  0.42652967  0.563136597  0.500000000
O  0.54644441  0.500000000  0.520000000

The z coordinate of the oxygen atom is displaced from its ground-state equilibrium (0.5 \(\rightarrow\) 0.52) to establish an initial out-of-plane angle. The self-consistent field convergence threshold conv_thr is tightened to obtain more accurate forces.

The wbse.in file is the input for the wbse.x code that computes excited-state energy and forces using TDDFT:

[3]:

%%bash
cat wbse.in

input_west:
  outdir: ./

wbse_init_control:
  solver: TDDFT

wbse_control:
  wbse_calculation: D
  n_liouville_eigen: 1
  trev_liouville: 0.00000001
  trev_liouville_rel: 0.000001
  spin_excitation: S
  l_forces: True
  forces_state: 1

The keyword l_forces: True turns on the calculation of analytical forces, for the first excited state (\(^1A''\)) as specified by forces_state: 1.

2. Perform excited-state geometry optimization

The workflow of carrying out an excited-state geometry optimization is managed by the WESTpy Python package, which automatically runs the pw.x code to compute the ground-state energy and forces, then runs the wbse.x code to compute the excited-state energy and forces and updates the excited-state geometry.

We initialize a bfgs_iter object by specifying the full commands to run pw.x and wbse.x, and the list of pseudopotentials. Then we call the solve method to perform the BFGS geometry optimization of the excited state.

[ ]:

from westpy import *

run_pw = "mpirun -n 2 pw.x"
run_wbse = "mpirun -n 2 wbse.x"
pp = [f"{element}_ONCV_PBE-1.2.upf" for element in ["C", "H", "O"]]

bfgs = bfgs_iter(run_pw=run_pw, run_wbse=run_wbse, pp=pp)
bfgs.solve()

The output would look like the following:

BFGS Geometry Relaxation

Running pw.x ...
Running wbse.x ...

number of scf cycles    =   1
number of bfgs steps    =   0
energy new              =     -45.2645563135 Ry
total force             =       0.1895441688 Ry/Bohr
new trust radius        =       0.1434817128 Bohr

Running pw.x ...
Running wbse.x ...

number of scf cycles    =   2
number of bfgs steps    =   1
energy old              =     -45.2645563135 Ry
energy new              =     -45.2661820794 Ry
total force             =       0.1428684377 Ry/Bohr
CASE: energy_new < energy_old
new trust radius        =       0.0600584257 Bohr

...

Running pw.x ...
Running wbse.x ...

number of scf cycles    =  15
number of bfgs steps    =  14
energy old              =     -45.2741518651 Ry
energy new              =     -45.2741969921 Ry
total force             =       0.0010256371 Ry/Bohr

Final geometry saved to final_geo

bfgs converged in  15 scf cycles and  14 bfgs steps
(criteria: energy < 1.00e-04 Ry, force < 1.00e-03 Ry/Bohr)

End of BFGS Geometry Optimization

Final energy =     -45.2741969921 Ry

If the reader does NOT have the computational resources to run the calculation, the files needed for the next step can be downloaded as:

[4]:

%%bash
mkdir -p final_geo
mkdir -p final_geo/pwscf.save
mkdir -p final_geo/west.wbse.save
wget -N -q https://west-code.org/doc/training/formaldehyde/es/pwscf.xml -O final_geo/pwscf.save/data-file-schema.xml
wget -N -q https://west-code.org/doc/training/formaldehyde/es/wbse.json -O final_geo/west.wbse.save/wbse.json

The total energy of the optimized excited state can be decomposed into two components. The first component is the ground-state total energy at this optimized excited-state geometry, which can be extracted from the final_geo/pwscf.save/data-file-schema.xml file:

[5]:

from xml.etree import ElementTree as ET

Hartree2eV = 27.2114

fname = "final_geo/pwscf.save/data-file-schema.xml"

xmlData = ET.parse(fname)
root = xmlData.getroot()
es_etot1 = float(root.find("output").find("total_energy").find("etot").text) * Hartree2eV

The second component is the excitation energy at the optimized excited-state geometry, which can be retrieved from the final_geo/west.wbse.save/wbse.json file:

[6]:

import json

fname = "final_geo/west.wbse.save/wbse.json"

with open(fname, "r") as f:
    j = json.load(f)

which_es = int(j["input"]["wbse_control"]["forces_state"])
es_etot2 = float(j["exec"]["davitr"][-1]["ev"][which_es - 1]) * Hartree2eV / 2.0

The total energy of the excited state is therefore:

[7]:

es_etot = es_etot1 + es_etot2
print(f"Excited state total energy: {es_etot:.6f} eV")

Excited state total energy: -615.987142 eV

The adiabatic excitation energy is the difference between the total energies of the excited and ground states at their respective optimized geometries:

[8]:

# Ground state total energy obtained from a separate DFT geometry relaxation
gs_etot = -619.532480

# Adiabatic excitation energy
aee = es_etot - gs_etot
print(f"Adiabatic excitation energy: {aee:.6f} eV")

Adiabatic excitation energy: 3.545338 eV

This value is consistent with the reference value of 3.530 eV.