Charge transfer#
Let us study the ethylene dimer as a simple and concrete example of the issue with charge-transfer (CT) excitations in TDDFT. We first create an xyz-file for an H-aggregate stacked along the z-axis with an easy way to alter the intermolecular separation distance.
Show code cell source
dimer_xyz = """12
C 0.67759997 0.00000000 -zcoord
C -0.67759997 0.00000000 -zcoord
H 1.21655197 0.92414474 -zcoord
H 1.21655197 -0.92414474 -zcoord
H -1.21655197 -0.92414474 -zcoord
H -1.21655197 0.92414474 -zcoord
C 0.67759997 0.00000000 zcoord
C -0.67759997 0.00000000 zcoord
H 1.21655197 0.92414474 zcoord
H 1.21655197 -0.92414474 zcoord
H -1.21655197 -0.92414474 zcoord
H -1.21655197 0.92414474 zcoord"""
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
This system is also found in the section on selection rules, where the Frenkel excitonic states, \(| \sigma \rangle\) and \(| \gamma \rangle\), and the charge-transfer states, \(| \delta \rangle\) and \(| \rho \rangle\), were introduced. They are illustrated in the figure below together with the canonical and delocalized molecular orbitals (MOs) that are involved in the associated electronic excitation processes.
The nature of the excited states is revealed by plotting the excitation energies as a function of dimer separation distance, see figure below. The excitation energies for the excitonic states are only weakly dependent on the distance whereas those of the CT states show a strong dependence.
Moreover, it is seen that the excitation energies of the CT (but not excitonic) states strongly depend on the amount of exact Hartree–Fock (HF) exchange in the exchange–correlation (XC) functional. We are here adopting the B3LYP and BHANDHLYP functionals with some 20% and 50% of HF exchange, respectively.
The excitation energies of the CT states are severely underestimated and unrealistic at the B3LYP level of theory.
Show code cell source
import matplotlib.pyplot as plt
import numpy as np
d = 2 * np.arange(2.0, 11.0, 2.0)
plt.plot(
d,
[8.7460, 9.3348, 9.6384, 9.7892, 9.8796],
"o-",
alpha=0.5,
color="green",
label=r"BHANDHLYP, $\rho$",
)
plt.plot(
d,
[8.4311, 9.3348, 9.6384, 9.7892, 9.8796],
"o-",
alpha=0.5,
color="blue",
label=r"BHANDHLYP, $\delta$",
)
plt.plot(
d,
[7.5915, 8.1239, 8.1364, 8.1393, 8.1403],
"o-",
alpha=0.5,
color="black",
label=r"BHANDHLYP, $\sigma$",
)
plt.plot(
d,
[8.2471, 8.1582, 8.1462, 8.1434, 8.1424],
"o-",
alpha=0.5,
color="red",
label=r"BHANDHLYP, $\gamma$",
)
plt.plot(
d,
[8.1910, 8.1955, 8.2069, 8.2095, 8.2104],
"^-",
alpha=0.5,
color="black",
label=r"B3LYP, $\sigma$",
)
plt.plot(
d,
[8.3491, 8.2266, 8.2157, 8.2132, 8.2123],
"^-",
alpha=0.5,
color="red",
label=r"B3LYP, $\gamma$",
)
plt.plot(
d,
[6.8032, 7.1899, 7.3113, 7.3716, 7.4077],
"^-",
alpha=0.5,
color="green",
label=r"B3LYP, $\rho$",
)
plt.plot(
d,
[6.6507, 7.1899, 7.3113, 7.3716, 7.4077],
"^-",
color="blue",
alpha=0.5,
label=r"B3LYP, $\delta$",
)
plt.legend(fontsize=6)
plt.xlabel("Intermolecular separation (Å)")
plt.ylabel("Excitation energy (eV)")
plt.setp(plt.gca(), ylim=[6.5, 11])
plt.show()
The underlying data for the figure above has been calculated with use of the VeloxChem program and the code cells found below. You can create all data by activating the iterations over distances and functionals in the code cell.
Here we will run the calculation for a single dimer separation of 20 Å and use the well-behaved BHANDHLYP functional.
import veloxchem as vlx
scf_drv = vlx.ScfRestrictedDriver()
scf_drv.ostream.mute()
lreig_drv = vlx.LinearResponseEigenSolver()
lreig_drv.ostream.mute()
lreig_drv.nstates = 6
# for xcfun in ["b3lyp", "bhandhlyp"]:
for xcfun in ["bhandhlyp"]:
print(f"XC functional: {xcfun:16s}")
scf_drv.xcfun = xcfun
# for zcoord in np.arange(2.0, 11.0, 2.0):
for zcoord in np.arange(10.0, 11.0, 2.0):
print(f"Z coordinate: {zcoord:8.2f}")
molecule = vlx.Molecule.read_xyz_string(
dimer_xyz.replace("zcoord", str(zcoord))
)
basis = vlx.MolecularBasis.read(molecule, "6-31g", ostream=None)
scf_results = scf_drv.compute(molecule, basis)
# lreig_drv.detach_attach = True
# lreig_drv.filename = "dimer_20A"
lreig_results = lreig_drv.compute(molecule, basis, scf_results)
print(
f" {'E (eV)':8s}{'mu_x':8s}{'mu_y':8s}{'mu_z':8s}{'m_x':8s}{'m_y':8s}{'m_z':8s}{'f':8s}{'R':8s}"
)
print(72 * "-")
for E, e, m, f, R in zip(
lreig_results["eigenvalues"],
lreig_results["electric_transition_dipoles"],
lreig_results["magnetic_transition_dipoles"],
lreig_results["oscillator_strengths"],
lreig_results["rotatory_strengths"],
):
print(
f"{E*27.2114:8.4f}{e[0]:8.4f}{e[1]:8.4f}{e[2]:8.4f}{m[0]:8.4f}{m[1]:8.4f}{m[2]:8.4f}{f:8.4f}{R:8.4f}"
)
XC functional: bhandhlyp
Z coordinate: 10.00
E (eV) mu_x mu_y mu_z m_x m_y m_z f R
------------------------------------------------------------------------
8.1403 0.0000 0.0000 -0.0000 0.0000 -4.9598 0.0000 0.0000 -0.0000
8.1424 -2.0065 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.8031 0.0000
9.5360 -0.0000 -0.0000 0.0214 0.0000 0.0000 -0.0000 0.0001 -0.0000
9.5360 0.0000 -0.0000 -0.0011 0.0000 0.0000 -0.0000 0.0000 0.0000
9.8796 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000
9.8796 -0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000
With data restricted to a single molecular configuration and functional, it can be less than obvious to determine the character of electronic transitions.
It is a common practice to turn to the dominant elements in the excitation vectors.
for state, exc_str in enumerate(lreig_results['excitation_details']):
print("State:", state + 1)
for orb_str in exc_str:
print(orb_str)
State: 1
HOMO -> LUMO 0.6330
HOMO-1 -> LUMO+1 0.6330
HOMO -> LUMO+1 0.3065
HOMO-1 -> LUMO -0.3065
State: 2
HOMO-1 -> LUMO+1 0.7032
HOMO -> LUMO -0.7032
State: 3
HOMO -> LUMO+2 -0.7265
HOMO-1 -> LUMO+3 0.6540
State: 4
HOMO-1 -> LUMO+3 -0.7208
HOMO -> LUMO+2 -0.6476
State: 5
HOMO-1 -> LUMO 0.9495
HOMO-1 -> LUMO+1 0.2275
HOMO -> LUMO 0.2121
State: 6
HOMO -> LUMO+1 0.9504
HOMO -> LUMO -0.2238
HOMO-1 -> LUMO+1 -0.2082
At this separation distance, excited states number 1, 2, 5, and 6 are all formed from \(\pi\pi^*\)-excitations, and from these coefficients alone, the character of the electronic transitions is not clear.
The next level of analysis is reached by turning to the natural transition orbitals (NTOs) or detachment/attachment densities. Below we will use the Valet module to perform a subgroup division of the dimer system into monomers and give a visual presentation of the electronic transitions based on the detachment/attachment densities.
States 5 and 6 amply reveal their CT character in the figure illustrations.
from valet import transition_analysis_utils as tau
subgroup_name = ["Ethylene 1", "Ethylene 2"]
atom_subgroup_map = [0] * 6 + [1] * 6
print(atom_subgroup_map)
[0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1]
def plot_transition_diagram(state_index, filename, subgroup_name, atom_subgroup_map):
folder = "../../data/tddft/" # Path to where the cube files were saved
density_type_plus = "detach" # or hole
density_type_minus = "attach" # or electron/particle
hole_cube = "%s_S%d_%s.cube" % (folder + filename, state_index, density_type_plus)
particle_cube = "%s_S%d_%s.cube" % (
folder + filename,
state_index,
density_type_minus,
)
# We have to load the detachment and attachment densities and compute the atomic charges
transition = tau.load_transition(hole_cube, particle_cube)
segment_array = tau.compute_atomic_charges([transition])
subgroup_info = tau.SubgroupInfo()
# Determine subgroup charges
subgroup_info.set_subgroups(subgroup_name, atom_subgroup_map)
tau.compute_subgroup_charges(transition, subgroup_info)
diagram_title = "Ethylene dimer: State %d" % (state_index)
tau.create_diagram(subgroup_info, title=diagram_title, save_plot=False)
plot_transition_diagram(1, "dimer_20A", subgroup_name, atom_subgroup_map)
plot_transition_diagram(2, "dimer_20A", subgroup_name, atom_subgroup_map)
plot_transition_diagram(5, "dimer_20A", subgroup_name, atom_subgroup_map)
plot_transition_diagram(6, "dimer_20A", subgroup_name, atom_subgroup_map)
