# Rovibronic states#

The sum-over-states expressions for response functions are expressed in terms of the exact rovibronic eigenstates of the molecular Hamiltonian in the absence of external fields. These states are rarely available in practical calculations, where instead electronic states are employed. However, for a diatomic system, we can determine the rovibrational wave functions and the associated transition moments and energies.

Adopting the adiabatic Born–Oppenheimer approximation, the total molecular rovibronic wave function takes the form of a product

$\Psi_{K, k}(\mathbf{r}, \mathbf{R}) = \Phi_{K;k}(\mathbf{R}) \Psi_K^\mathrm{e}(\mathbf{r}; \mathbf{R})$

where $$K$$ and $$k$$ are indices identifying the electronic and rovibrational levels, respectively. For a diatomic system, the rovibrational wave function is in turn separated into its vibrational, $$\psi$$, and rotational, $$\psi_\mathrm{r}$$, components

$\Phi_{K;k}(\mathbf{R}) = \psi_k(R) \psi_{\mathrm{r}; k}(\theta,\phi)$

where $$k$$ is understood as a compound index uniquely identifying a given rotational level of a given vibrational state.

## Rotational Schrödinger equation#

The Schrödinger equation for diatomic molecule with fixed bond length equals

$- \frac{\hbar^2}{2\mu R_0^2} \left[ \frac{1}{\sin\theta}\frac{\partial}{\partial \theta} \left(\sin\theta \frac{\partial}{\partial \theta} \right) + \frac{1}{\sin^2\theta} \frac{\partial^2}{\partial \phi^2} \right] \psi_\mathrm{r}(\theta,\phi) = E \, \psi_\mathrm{r}(\theta,\phi)$

or, equivalently,

$\frac{\hat{L}^2\psi_\mathrm{r}(\theta,\phi)}{2\mu R_0^2} = E \, \psi_\mathrm{r}(\theta,\phi)$

where $$\hat{L}^2$$ is the square of the total angular moment operator. The angles $$\theta$$ and $$\phi$$ describe the orientation of the diatomic molecule with respect to a laboratory-fixed coordinate system and $$\mu$$ is the reduced mass defined as

$\mu=\dfrac{M_A M_B}{M_A + M_B}$

The eigenfunctions to the total angular momentum operator are known to the the sperical harmonics, $$Y_J^M(\theta, \phi)$$, and the associated rotational energies are given by

$E_J = \frac{\hbar^2 J(J+1)}{2\mu R_0^2} = B J(J+1)$

introducing the rotational constant $$B$$. We note that rotational levels are $$(2J+1)$$-fold degenerate.

Transitions between rotational levels are classified according to

1. Q-branch: when $$\Delta J = 0$$

2. R-branch: when $$\Delta J = +1$$

3. P-branch: when $$\Delta J = -1$$

## Numerov’s method and the vibrational Schrödinger equation#

Numerov’s method can be employed to solve 1D differential equations using finite differences on a grid of points. It applies to a differential equation on the form

$\frac{d^2 y}{dx^2} = -g(x)y(x) + s(x)$

This method will be used to solve the vibrational part of the Schrödinger equation that for the lowest rotational level ($$J=0$$) takes the form

$- \dfrac{\hbar^2}{2\mu} \dfrac{\partial^2\psi(R)}{\partial R^2} +V(R) \psi(R) = E \psi(R)$

For notational convenience, we introduce $$x = (R - R_0)$$ to arrive at

$\psi^{\prime \prime} (x) = - k^{2} (x) \psi (x)$

where

$k^{2} (x) = \dfrac {2\mu} {\hbar^{2}} [ E - V (x) ]$

The values of the wave function in three equidistant points separated by $$h$$ are related as

$\psi(x_{i+1}) = \dfrac{ \psi(x_{i}) \left( 2 + \dfrac {10 h^{2}} {12} k^{2}(x_{i}) \right) - \psi(x_{i-1}) \left( 1 - \dfrac {h^{2}} {12} k^{2}(x_{i-1}) \right) }{ 1 - \dfrac {h^{2}} {12} k^{2}(x_{i+1}) }$

from which solutions are found by the algorithm

1. Guess the energy $$E_\mathrm{guess}$$

2. For $$E_\mathrm{guess}$$, determine the classically allowed interval [$$x_a,x_b$$] inside $$V(x)$$

3. Pick two points $$x_\mathrm{min}$$ and $$x_\mathrm{max}$$ far into the classically forbidden region and impose the wave function to be zero at these points

4. Determine iteratively the wave functions inside the interval $$[x_\mathrm{min}, x_\mathrm{max}]$$:

1. If $$\psi(x_a)<\epsilon$$ and $$\psi(x_b)<\epsilon$$, accept the wave function

2. If not, renew the energy guess

## Numerov’s method in VeloxChem#

import matplotlib.pyplot as plt
import numpy as np
import veloxchem as vlx


We will use HCl as an example.

molecule_string = """
Cl   0.0     0.0   0.0
H    1.274   0.0   0.0
"""

molecule = vlx.Molecule.read_str(molecule_string, units="angstrom")

basis_set_label = "6-31+G"
basis = vlx.MolecularBasis.read(molecule, basis_set_label)

scf_settings = {"conv_thresh": 1.0e-6}
method_settings = {"xcfun": "b3lyp", "grid_level": 4}


The following settings are available for the Numerov solver in VeloxChem.

num_drv = vlx.NumerovDriver()
num_drv.print_keywords()

                     ================================================================================
@numerov
--------------------------------------------------------------------------------
pec_displacements   sequence       PEC scanning range [bohr]
pec_potential       string         potential type for fitting
reduced_mass        float          reduced mass of the molecule
el_transition       boolean        include an electronic transition
final_el_state      integer        final state of the electronic transition
exc_conv_thresh     float          excited state calculation threshold
n_vib_states        integer        number of vibrational states to be resolved
temp                float          temperature in Kelvin
n_rot_states        integer        number of rotational states to be resolved
conv_thresh         float          convergence threshold for vibronic energies
max_iter            integer        max. iteration number per vibronic state
n_margin            float          negative margin of the grid [bohr]
p_margin            float          positive margin of the grid [bohr]
steps_per_au        integer        number of grid points per Bohr radius
================================================================================


We will determine the potential on an interval of

numerov_settings = {
"n_vib_states": 4,
"pec_displacements": "-0.3 - 0.9 (0.2)",
"el_transition": "no",
"pec_potential": "morse",
"n_rot_states": 10,
}
num_drv.update_settings(numerov_settings, scf_settings, method_settings)
num_results = num_drv.compute(molecule, basis)


Potential Energy Curve Driver Setup
=====================================

Number of Geometries          : 7
Wave Function Model           : Spin-Unrestricted Kohn-Sham
SCF Convergece Threshold      : 1.0e-06
DFT Functional                : B3LYP
Molecular Grid Level          : 4


               *** Geometry 1/7:  Initial state converged

               *** Geometry 2/7:  Initial state converged

               *** Geometry 3/7:  Initial state converged

               *** Geometry 4/7:  Initial state converged

               *** Geometry 5/7:  Initial state converged

               *** Geometry 6/7:  Initial state converged

               *** Geometry 7/7:  Initial state converged


Potential Energy Curve Information
----------------------------------
rel. GS energy
Bond distance:     1.1152 angs      0.85345 eV
Bond distance:     1.2211 angs      0.16198 eV
Bond distance:     1.3269 angs      0.00000 eV
Bond distance:     1.4328 angs      0.13267 eV
Bond distance:     1.5386 angs      0.42875 eV
Bond distance:     1.6444 angs      0.80427 eV
Bond distance:     1.7503 angs      1.21606 eV



Numerov Driver Setup
======================

Number of Vibrational States : 4
Energy Convergence Threshold : 1.0e-12
Max. Number of Iterations    : 1000
Grid Points per Bohr Radius  : 500
Margin right of scanned PEC  : 1.50
Margin left of scanned PEC   : 0.50


               *** Vibronic state 0 converged in 120 iterations

               *** Vibronic state 1 converged in 168 iterations

               *** Vibronic state 2 converged in 155 iterations

               *** Vibronic state 3 converged in 174 iterations


*** All vibronic states converged for the initial electronic state


               Rovibronic IR Spectrum
----------------------
Excitation energy    Oscillator strength
P Branch
2613.48 cm-1            3.94260e-09
2593.84 cm-1            5.43630e-09
2574.19 cm-1            5.72716e-09
2554.55 cm-1            5.03997e-09
2534.91 cm-1            3.83494e-09
2515.27 cm-1            2.56641e-09
2495.63 cm-1            1.52520e-09
2475.99 cm-1            8.09790e-10
2456.35 cm-1            3.85657e-10
Q Branch
2633.12 cm-1            5.85360e-08
R Branch
2652.76 cm-1            1.37689e-09
2672.40 cm-1            3.75712e-09
2692.04 cm-1            5.18055e-09
2711.68 cm-1            5.45773e-09
2731.33 cm-1            4.80287e-09
2750.97 cm-1            3.65453e-09
2770.61 cm-1            2.44567e-09
2790.25 cm-1            1.45345e-09
2809.89 cm-1            7.71694e-10
2829.53 cm-1            3.67514e-10



The results are stored in a dictionary.

print(num_results.keys())

dict_keys(['grid', 'potential', 'psi', 'vib_levels', 'excitation_energies', 'oscillator_strengths'])

grid = num_results["grid"]
potential = num_results["potential"]

vib_psi2 = num_results["psi"] ** 2
vib_energies = num_results["vib_levels"]

exc = num_results["excitation_energies"]
osc = num_results["oscillator_strengths"]


The vibrational levels can be plotted against the potential energy curve.

fig, ax1 = plt.subplots()

for energy, psi2 in zip(vib_energies, vib_psi2):
ax1.plot(grid, energy + psi2, "k")

ax = plt.subplot()
ax.plot(grid, potential)
ax.set_ylim(-0.001, 0.05)
ax.set_xlabel("Displacements [Å] respect to the input structure")
ax.set_ylabel("Energy [a.u.]")

plt.grid(True)

plt.show() ## Rovibrational spectra#

Once we have the vibrational wave functions we can calculate the dipole transition moments

$\mu^{fi} = \langle \psi_{f} | \mu(R) | \psi_{i} \rangle$

where $$\mu(R)$$ is the electronic dipole moment at internuclear separation distance $$R$$.

Only the first transition dipole moment (corresponding to the Q-branch) is relevant and the rest are negligible. We will obtain the fine structure by deriving the rotational excited states for the R and P branches.

Applying the rotational selection rules derived before for the rigid rotor, we get the P-branch ($$\Delta J = -1$$)

$E = \hbar \omega_{\text{stretch}} + \left( \dfrac{h^2}{8\pi^2I} \right)((J - 1)J - J(J + 1)) = \hbar \omega_{\text{stretch}} -2\left( \dfrac{h^2}{8\pi^2 I} \right)J$

and the R-branch ($$\Delta J = +1$$)

$E = \hbar \omega_{\text{stretch}} + \left( \dfrac{h^2}{8\pi^2I} \right)((J+1)(J+2) - J(J + 1)) = \hbar \omega_{\text{stretch}} + 2\left( \dfrac{h^2}{8\pi^2I} \right)(J + 1)$

where

$I = \frac{m_a m_b}{m_a + m_b} R^2$

VeloxChem takes care of the spectra and the branches storing the data on the keys: excitation_energies and oscillator_strengths.

As well, computes the dipole moment for every point of the potential energy curve

ax1 = plt.bar(exc["P"], osc["P"], width=5.0, label="P-Branch")
ax2 = plt.bar(exc["R"], osc["R"], width=5.0, label="R-Branch")
plt.xlabel("Excitation energies in cm-1")
plt.ylabel("Oscillator strength")
plt.legend()

plt.show() ## Vibronic transitions#

In this case, we want to consider the transition dipole moments between rovibrational states in different electronic states. The expression for the transition dipole integral will in this case involve the ground, $$|0\rangle$$, and final, $$|F\rangle$$, electronic states.

Expanding the electronic transition moment in a Taylor series about the equilibrium distance of the ground electronic state, we get

$\mu^{F0}(R) = \mu_{F0}(R_0) + \left. \dfrac{\partial \mu^{F0}}{\partial R} \right|_{R_{0}} (R - R_{0}) + \cdots .$

Truncating this expansion to first and second order are known as the Frank–Condon and Herzberg–Teller approximations, respectively. We will here adopt the former.

num_drv = vlx.NumerovDriver()

numerov_settings = {
"n_vib_states": 3,
"pec_displacements": "-0.3 - 0.9 (0.2)",
"el_transition": "yes",
"final_el_state": 2,
"pec_potential": "morse",
}
num_drv.update_settings(numerov_settings, scf_settings, method_settings)
num_results = num_drv.compute(molecule, basis)


Potential Energy Curve Driver Setup
=====================================

Number of Geometries          : 7
Wave Function Model           : Spin-Unrestricted Kohn-Sham
SCF Convergece Threshold      : 1.0e-06
DFT Functional                : B3LYP
Molecular Grid Level          : 4
Targeted Exicted State        : 2
Excited State Threshold       : 1.0e-04


               *** Geometry 1/7:  Initial state converged,  Final state converged

               *** Geometry 2/7:  Initial state converged,  Final state converged

               *** Geometry 3/7:  Initial state converged,  Final state converged

               *** Geometry 4/7:  Initial state converged,  Final state converged

               *** Geometry 5/7:  Initial state converged,  Final state converged

               *** Geometry 6/7:  Initial state converged,  Final state converged

               *** Geometry 7/7:  Initial state converged,  Final state converged


Potential Energy Curve Information
----------------------------------
rel. GS energy  rel. ES energy
Bond distance:     1.1152 angs      0.85345 eV     11.00173 eV
Bond distance:     1.2211 angs      0.16198 eV     10.18858 eV
Bond distance:     1.3269 angs      0.00000 eV      9.77834 eV
Bond distance:     1.4328 angs      0.13267 eV      9.68542 eV
Bond distance:     1.5386 angs      0.42875 eV      9.83600 eV
Bond distance:     1.6444 angs      0.80427 eV     10.11118 eV
Bond distance:     1.7503 angs      1.21606 eV     10.28647 eV



Numerov Driver Setup
======================

Number of Vibrational States : 3
Energy Convergence Threshold : 1.0e-12
Max. Number of Iterations    : 1000
Grid Points per Bohr Radius  : 500
Margin right of scanned PEC  : 1.50
Margin left of scanned PEC   : 0.50


               *** Vibronic state 0 converged in 115 iterations

               *** Vibronic state 1 converged in 168 iterations

               *** Vibronic state 2 converged in 163 iterations


*** All vibronic states converged for the initial electronic state


               *** Vibronic state 0 converged in 99 iterations

               *** Vibronic state 1 converged in 137 iterations

               *** Vibronic state 2 converged in 138 iterations


*** All vibronic states converged for the final electronic state


               UV/vis Absorption/Emission Spectrum
-----------------------------------
Excitation energy    Oscillator strength
Absorption
78438.20 cm-1            2.62965e-01
80567.17 cm-1            5.19146e-02
82599.54 cm-1            3.99141e-03
Emission
78438.20 cm-1            2.62965e-01
75805.08 cm-1            1.01191e-01
73278.86 cm-1            6.49875e-03


print(num_results.keys())

dict_keys(['grid', 'i_potential', 'i_psi', 'i_vib_levels', 'f_potential', 'f_psi', 'f_vib_levels', 'excitation_energies', 'oscillator_strengths'])

grid = num_results["grid"]

dE_elec = np.min(num_drv.pec_energies["f"]) - np.min(num_drv.pec_energies["i"])

gs_potential = num_results["i_potential"]
gs_vib_psi2 = num_results["i_psi"] ** 2
gs_vib_energies = num_results["i_vib_levels"]

es_potential = num_results["f_potential"]
es_vib_psi2 = num_results["f_psi"] ** 2
es_vib_energies = num_results["f_vib_levels"]

exc = num_results["excitation_energies"]
osc = num_results["oscillator_strengths"]

displacements = num_drv.pec_displacements
gs_energies = num_drv.pec_energies["i"]
es_energies = num_drv.pec_energies["f"]

delta_E = np.min(es_energies) - np.min(gs_energies)


The vibrational levels can be then plotted against the potential energy curves:

fig, ax1 = plt.subplots(figsize=(6, 10))

for energy, psi2 in zip(gs_vib_energies, gs_vib_psi2):
ax1.plot(grid, energy + psi2, "k")

for energy, psi2 in zip(es_vib_energies, es_vib_psi2):
ax1.plot(grid, delta_E + energy + psi2, "k")

ax = plt.subplot()
ax.plot(grid, gs_potential)
ax.plot(grid, es_potential + delta_E)
ax.set_xlim(-0.8, 2.0)
ax.set_ylim(-0.01, 0.45)
ax.grid(True)
ax.set_xlabel("Displacements [Å] w.r.t. input structure")
ax.set_ylabel("Energy [a.u.]")

plt.show() Now we can just plot the spectrum by representing the results from the driver:

plt.bar(exc["emission"], osc["emission"], width=500.0, label="emission")
plt.bar(exc["absorption"], osc["absorption"], width=500.0, label="absorption")

plt.legend()
plt.xlabel("Excitation energies in cm-1")
plt.ylabel("Oscillator strength")

plt.show() ## User defined potential#

As an alternative to the automatized determination of potential energy curves and transition moments, these may be instead be determined separately, perhaps with a highly accurate post-Hartree–Fock approach. The method read_pec_data is then used.

As a simple illustration, we will use two harmonic potentials and transition moment that is independent of nuclear displacements.

num_drv = vlx.NumerovDriver()

x = np.arange(-0.7, 2.01, 0.1)

delta_E = 0.2

bond_lengths = x + 2.0
gs_energies = 0.75 * x**2
es_energies = 0.5 * (x - 0.35) ** 2 + delta_E

tmom = [0, 0, 1.0]
trans_moments = [tmom for displacement in x]

num_drv.read_pec_data(bond_lengths, trans_moments, gs_energies, es_energies)

num_drv.set_reduced_mass(1.0)

numerov_settings = {
"n_vib_states": 10,
"el_transition": "yes",
"final_el_state": 1,
"pec_potential": "harmonic",
}
num_drv.update_settings(numerov_settings, scf_settings, method_settings)

num_results = num_drv.compute()


Numerov Driver Setup
======================

Number of Vibrational States : 10
Energy Convergence Threshold : 1.0e-12
Max. Number of Iterations    : 1000
Grid Points per Bohr Radius  : 500
Margin right of scanned PEC  : 1.50
Margin left of scanned PEC   : 0.50


               *** Vibronic state 0 converged in 120 iterations

               *** Vibronic state 1 converged in 194 iterations

               *** Vibronic state 2 converged in 193 iterations

               *** Vibronic state 3 converged in 190 iterations

               *** Vibronic state 4 converged in 191 iterations

               *** Vibronic state 5 converged in 187 iterations

               *** Vibronic state 6 converged in 188 iterations

               *** Vibronic state 7 converged in 194 iterations

               *** Vibronic state 8 converged in 195 iterations

               *** Vibronic state 9 converged in 200 iterations


*** All vibronic states converged for the initial electronic state


               *** Vibronic state 0 converged in 103 iterations

               *** Vibronic state 1 converged in 172 iterations

               *** Vibronic state 2 converged in 167 iterations

               *** Vibronic state 3 converged in 170 iterations

               *** Vibronic state 4 converged in 170 iterations

               *** Vibronic state 5 converged in 167 iterations

               *** Vibronic state 6 converged in 172 iterations

               *** Vibronic state 7 converged in 173 iterations

               *** Vibronic state 8 converged in 171 iterations

               *** Vibronic state 9 converged in 167 iterations


*** All vibronic states converged for the final electronic state


               UV/vis Absorption/Emission Spectrum
-----------------------------------
Excitation energy    Oscillator strength
Absorption
43314.88 cm-1            5.67445e-04
46035.10 cm-1            3.61293e-03
48755.33 cm-1            1.10780e-02
51475.56 cm-1            2.17810e-02
54195.79 cm-1            3.08405e-02
56916.02 cm-1            3.34778e-02
59636.25 cm-1            2.89534e-02
62356.48 cm-1            2.04635e-02
65076.70 cm-1            1.20261e-02
67796.93 cm-1            5.94669e-03
Emission
43314.88 cm-1            5.67445e-04
39983.29 cm-1            2.56212e-03
36651.70 cm-1            5.98373e-03
33320.12 cm-1            9.60179e-03
29988.53 cm-1            1.18635e-02
26656.95 cm-1            1.19878e-02
23325.36 cm-1            1.02682e-02
19993.77 cm-1            7.61992e-03
16662.19 cm-1            4.95672e-03
13330.60 cm-1            2.83167e-03



Retrieving results from the calculation.

grid = num_results["grid"]

gs_potential = num_results["i_potential"]
gs_psi2 = num_results["i_psi"] ** 2
gs_vib_energies = num_results["i_vib_levels"]

es_potential = num_results["f_potential"]
es_psi2 = num_results["f_psi"] ** 2
es_vib_energies = num_results["f_vib_levels"]

exc = num_results["excitation_energies"]
osc = num_results["oscillator_strengths"]


Plotting potentials and vibrational wave functions.

fig, ax1 = plt.subplots(figsize=(6, 10))

for energy, psi2 in zip(gs_vib_energies, gs_psi2):
ax1.plot(grid, energy + psi2, "k")

for energy, psi2 in zip(es_vib_energies, es_psi2):
ax1.plot(grid, delta_E + energy + psi2, "k")

ax = plt.subplot()

ax.plot(grid, gs_potential)
ax.plot(grid, es_potential + delta_E)

ax.set_xlim(-1.5, 2.2)
ax.set_ylim(-0.01, 0.35)

ax.grid(True)

plt.show() Plotting the vibronic absorption and emission spectra.

plt.bar(exc["absorption"], osc["absorption"], width=1000.0, label="Absorption")
plt.bar(exc["emission"], osc["emission"], width=1000.0, label="Emission")

plt.legend()
plt.xlabel(r"Excitation energies [cm$^{-1}$]")
plt.ylabel("Oscillator strength")

plt.show() 