This tutorial can be downloaded link.

Intro Tutorial 12: 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).

We will use the formaldehyde molecule as an example to determine the optimized atomic geometry of the \(^1A''\) excited state using TDDFT with the Perdew–Burke–Ernzerhof (PBE) functional and calculate the adiabatic excitation energy.

In the first part of the tutorial, we review the relaxation of atomic geometry in the electronic ground state. We will conduct DFT calculations using the pw.x executable from the Quantum ESPRESSO package. In the second part, we relax the atomic geometry of the electronic excited state by combining DFT calculations - using the pw.x executable - and TDDFT calculations using the wbse.x executable from the WEST code.

More details about the implementation of TDDFT in WEST, for spin-conserving or spin-flip excitations with analytical forces and hybrid functionals, may be found in Y. Jin, V. Yu, M. Govoni, A. Xu, and G. Galli, J. Chem. Theory Comput. 19, 8689 (2023).

Ground state geometry optimization

We first review how to perform the ground state geometry optimization calculation using DFT. 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/gs/pw.in

We can inspect the pw.in file, the input for the pw.x code:

[2]:

%%bash
cat pw.in

&CONTROL
calculation = 'relax'
pseudo_dir = './'
/
&SYSTEM
ibrav = 0
nat = 4
ntyp = 3
ecutwfc = 50
nbnd = 10
/
&ELECTRONS
diago_full_acc = .true.
/
&IONS
ion_dynamics = 'bfgs'
/
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.500000000

We run pw.x on 2 cores.

[ ]:

%%bash
mpirun -n 2 pw.x -i pw.in > pw.out

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

[3]:

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

We extract the total energy at the optimized geometry from the file pwscf.save/data-file-schema.xml.

[4]:

from xml.etree import ElementTree as ET

Hartree2eV = 27.2114

fname = 'pwscf.save/data-file-schema.xml'

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

print(f'Ground state total energy: {gs_etot:.6f} eV')

Ground state total energy: -619.532480 eV

We can also calculate the bond lengths and angles of the formaldehyde molecule based on the optimized geometry, and compare the results with published values. First, we extract the atomic coordinates.

[5]:

import numpy as np

Bohr2Angstrom = 0.529177

# extract atomic coordinates
atomic_structure = root.find('input/atomic_structure')
alat = float(atomic_structure.attrib['alat'])
atoms = atomic_structure.findall('atomic_positions/atom')
gs_geo = [np.array(atom.text.split(), dtype=np.float64) * Bohr2Angstrom for atom in atoms]

We calculate the bond lengths.

[6]:

gs_bl_c_h1 = np.linalg.norm(gs_geo[1] - gs_geo[0])
gs_bl_c_h2 = np.linalg.norm(gs_geo[2] - gs_geo[0])
gs_bl_c_o = np.linalg.norm(gs_geo[3] - gs_geo[0])

print(f'C-H1 bond length: {gs_bl_c_h1:.4f} Å')
print(f'C-H2 bond length: {gs_bl_c_h2:.4f} Å')
print(f'C-O  bond length: {gs_bl_c_o:.4f} Å')

C-H1 bond length: 1.1169 Å
C-H2 bond length: 1.1169 Å
C-O  bond length: 1.2067 Å

The calculated bond lengths agree well with the reference values of 1.118 Å and 1.211 Å for the C-H and C-O bonds, respectively.

We calculate the angle formed between the two C-H bonds.

[7]:

ang = np.dot(gs_geo[1] - gs_geo[0], gs_geo[2] - gs_geo[0])
ang /= np.linalg.norm(gs_geo[1] - gs_geo[0])
ang /= np.linalg.norm(gs_geo[2] - gs_geo[0])
gs_ang_h1_c_h2 = np.arccos(ang) * 180 / np.pi

print(f'H1-C-H2 bond angle: {gs_ang_h1_c_h2:.4f} degrees')

H1-C-H2 bond angle: 115.9762 degrees

This is consistent with the reference value of 116.1 degrees.

We calculate the out-of-plane angle \(\Phi\).

[8]:

n_hch = np.cross(gs_geo[1] - gs_geo[0], gs_geo[2] - gs_geo[0])
ang = np.dot(n_hch, gs_geo[3] - gs_geo[0])
ang /= np.linalg.norm(n_hch)
ang /= np.linalg.norm(gs_geo[3] - gs_geo[0])
gs_ang_oop = np.arccos(ang) * 180 / np.pi - 90

print(f'Out-of-plane angle: {gs_ang_oop:.4f} degrees')

Out-of-plane angle: 0.0000 degrees

This is also consistent with the reference value of 0.0 degrees.

Excited state geometry optimization

Now we explain how to perform an excited state geometry optimization calculation. Download the following files to your working directory:

[9]:

%%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
wget -N -q https://west-code.org/doc/training/formaldehyde/es/wbse.in

We can inspect the pw.in file, which contains the initial atomic coordinates for the excited state geometry optimization:

[10]:

%%bash
cat pw.in

&CONTROL
calculation = 'relax'
pseudo_dir = './'
/
&SYSTEM
ibrav = 0
nat = 4
ntyp = 3
ecutwfc = 50
nbnd = 10
/
&ELECTRONS
diago_full_acc = .true.
/
&IONS
ion_dynamics = 'bfgs'
/
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.500000000

We made a minor adjustment to the z coordinate of the oxygen atom to establish an initial out-of-plane angle.

We can examine the wbse.in file, the input file for the wbse.x executable that computes the excited state energy and forces through TDDFT:

[11]:

%%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. The keyword forces_state: 1 instructs the code to compute the forces of the first excited state, \(^1A''\).

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

We create a bfgs_iter object by specifying the 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.2645433679 Ry
total force             =       0.1896335500 Ry/Bohr
new trust radius        =       0.1435687300 Bohr

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

number of scf cycles    =   2
number of bfgs steps    =   1
energy old              =     -45.2645433679 Ry
energy new              =     -45.2661898561 Ry
total force             =       0.1429407362 Ry/Bohr
CASE: energy_new < energy_old
new trust radius        =       0.0600494533 Bohr

...

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

number of scf cycles    =  15
number of bfgs steps    =  14
energy old              =     -45.2741434562 Ry
energy new              =     -45.2741629611 Ry
total force             =       0.0009790657 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.2741629611 Ry

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

[12]:

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

[13]:

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:

[14]:

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:

[15]:

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

Excited state total energy: -615.986679 eV

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

[16]:

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

Adiabatic excitation energy: 3.545801 eV

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

We can also compute the bond lengths and angles at the optimized excited state geometry. First, we extract the atomic coordinates.

[17]:

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

xmlData = ET.parse(fname)
root = xmlData.getroot()
atomic_structure = root.find('input/atomic_structure')
alat = float(atomic_structure.attrib['alat'])
atoms = atomic_structure.findall('atomic_positions/atom')
es_geo = [np.array(atom.text.split(), dtype=np.float64) * Bohr2Angstrom for atom in atoms]

We calculate the bond lengths.

[18]:

es_bl_c_h1 = np.linalg.norm(es_geo[1] - es_geo[0])
es_bl_c_h2 = np.linalg.norm(es_geo[2] - es_geo[0])
es_bl_c_o = np.linalg.norm(es_geo[3] - es_geo[0])

print(f'C-H1 bond length: {es_bl_c_h1:.4f} Å')
print(f'C-H2 bond length: {es_bl_c_h2:.4f} Å')
print(f'C-O  bond length: {es_bl_c_o:.4f} Å')

C-H1 bond length: 1.1020 Å
C-H2 bond length: 1.1015 Å
C-O  bond length: 1.3111 Å

The calculated bond lengths agree well with the reference values of 1.103 Å and 1.308 Å for the C-H and C-O bonds, respectively.

We calculate the angle formed between the two C-H bonds.

[19]:

ang = np.dot(es_geo[1] - es_geo[0], es_geo[2] - es_geo[0])
ang /= np.linalg.norm(es_geo[1] - es_geo[0])
ang /= np.linalg.norm(es_geo[2] - es_geo[0])
es_ang_h1_c_h2 = np.arccos(ang) * 180 / np.pi

print(f'H1-C-H2 bond angle: {es_ang_h1_c_h2:.4f} degrees')

H1-C-H2 bond angle: 117.1526 degrees

This is consistent with the reference value of 116.8 degrees.

We calculate the out-of-plane angle \(\Phi\).

[20]:

n_hch = np.cross(es_geo[1] - es_geo[0], es_geo[2] - es_geo[0])
ang = np.dot(n_hch, es_geo[3] - es_geo[0])
ang /= np.linalg.norm(n_hch)
ang /= np.linalg.norm(es_geo[3] - es_geo[0])
es_ang_oop = np.arccos(ang) * 180 / np.pi - 90

print(f'Out-of-plane angle: {es_ang_oop:.4f} degrees')

Out-of-plane angle: 31.1072 degrees

This is also consistent with the reference value of 30.0 degrees.

Finally, we compile the results into a table for a comprehensive view.

[21]:

import pandas as pd

data = {
    'State': ['$^1A_1$, this work', '$^1A_1$, ref.', '$^1A^{\prime\prime}$, this work', '$^1A^{\prime\prime}$, ref.'],
    'C-O bond length (Å)': [gs_bl_c_o, 1.211, es_bl_c_o, 1.308],
    'C-H bond length (Å)': [gs_bl_c_h1, 1.118, es_bl_c_h1, 1.103],
    'H-C-H angle (degree)': [gs_ang_h1_c_h2, 116.1, es_ang_h1_c_h2, 116.8],
    'Out-of-plane angle (degree)': [gs_ang_oop, 0.0, es_ang_oop, 30.0],
    'Adiabatic excitation energy (eV)': ['N/A', 'N/A', aee, 3.53]
}

df = pd.DataFrame(data)
pd.options.display.float_format = '{:.3f}'.format
df

[21]:
State C-O bond length (Å) C-H bond length (Å) H-C-H angle (degree) Out-of-plane angle (degree) Adiabatic excitation energy (eV)
0 $^1A_1$, this work 1.207 1.117 115.976 0.000 N/A
1 $^1A_1$, ref. 1.211 1.118 116.100 0.000 N/A
2 $^1A^{\prime\prime}$, this work 1.311 1.102 117.153 31.107 3.546
3 $^1A^{\prime\prime}$, ref. 1.308 1.103 116.800 30.000 3.530