Ctrl+K

Exercises

Exercises#

Dioxetane chemiluminescence#

1,2-Dioxetane is a small molecule that when heated emits a weak light and produces formaldehyde. It is also a model compound for the study of chemi-and bio-luminescence as many such compounds have a similar moiety.

Here we will study the dissociation from dioxetane to 2 molecules of formaldehyde on both ground and excited states at the SA-CASSCF level. 21 snapshots along the reaction path have been computed, so we will only perform a scan along this reaction path. You can download those files here.

import matplotlib.pyplot as plt
import multipsi as mtp
import numpy as np
import py3Dmol
import veloxchem as vlx
Hide code cell output 
* Warning * Environment variable OMP_NUM_THREADS not set.
* Warning * Setting OMP_NUM_THREADS to 6.

# Visualize the reaction path
dioxetane_xyz = ""
for i in range(1, 22):
    f = open("../../data/dioxetane/dioxetane." + str(i) + ".xyz", "r")
    dioxetane_xyz += f.read()

viewer = py3Dmol.view(width=400, height=300)
viewer.addModelsAsFrames(dioxetane_xyz)
viewer.animate({"loop": "forward"})
viewer.setStyle({"stick": {}})
# rotate for a better initial view
viewer.rotate(-90, "y")
viewer.rotate(-90, "x")
viewer.show()

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

As seen above, the reaction is breaking 2 bonds, the O-O and C-C. Knowing this, what is the minimum active space needed to describe the ground state of the reaction?

Answer 4 orbitals, corresponding to the $\sigma$ and $\sigma^*$ of the C-C and O-O bond. Those then become the $\pi$ and $\pi^*$ of the 2 formaldehyde products.

Here we will also compute the first 2 excited states. Can you guess what is the lowest excitation the formaldehyde molecule (for example)? What does that imply about our minimum active space to describe both ground and excited states?

Answer The lowest excitation is likely $n \rightarrow \pi^*$ with $n$ the oxygen lone pair. We thus need to add 2 extra orbitals to describe that excitation, otherwise our active space may not be stable due to competition between strong correlation and excitations.

Try our initial minimal active space and our extended active space on the formaldehyde molecule with 3 states in the state-averaging. Use the orbital viewer before and after the calculation to check the orbitals. Compare the excitation energies.

formaldehyde_str="""
 O     0.00000000     0.00000000     1.08516846
 C     0.00000000     0.00000000    -0.12131123
 H     0.93777476     0.00000000    -0.71598634
 H    -0.93777476     0.00000000    -0.71598634
"""
formaldehyde = vlx.Molecule.read_molecule_string(formaldehyde_str, units='angstrom')
basis = vlx.MolecularBasis.read(formaldehyde,"cc-pvdz")

#Hartree-Fock guess
scf_drv = vlx.ScfRestrictedDriver()
scf_results = scf_drv.compute(formaldehyde,basis)

...

Now that we have settle for our active space, let’s perform the scan over the entire reaction path, still with 3 states.

# 1) Find the orbitals for the product

molecule=vlx.Molecule.read_xyz_file("../../data/dioxetane/dioxetane.21.xyz")
basis = vlx.MolecularBasis.read(molecule,"cc-pvdz")
scf_drv = vlx.ScfRestrictedDriver()
scf_results = scf_drv.compute(molecule,basis)

# Visualize and select the active space
viewer=mtp.OrbitalViewer()
viewer.plot(molecule,basis,scf_drv.mol_orbs)


# Compute a SA-CASSCF with 3 states
...


# 2) Compute the states over the entire reaction path in reverse
Energies=[]
for i in range(21,0,-1):
    molecule=vlx.Molecule.read_xyz_file("../../data/dioxetane/dioxetane."+str(i)+".xyz")

    ...
    
    Energies.append(mcscf_drv.get_energies())

# 4) Verify the final active space (if it has changed, there is probably a warning message in the output of MultiPsi telling you so, you may fix the active space and restart from there)

viewer=mtp.OrbitalViewer()
viewer.plot(molecule,basis,space)


# 5) Plot the results

plt.figure(figsize=(6,4))
plt.title('State energies during dioxetane dissociation')
x = np.array(range(21,0,-1))
Ene = (np.array(Energies)-Energies[-1][0])*627.5
plt.plot(x, Ene[:,0], label='Ground state')
for i in range(2):
    plt.plot(x, Ene[:,i+1], label='State'+str(i+1))
plt.ylim(-80, 60) 
plt.legend()
plt.tight_layout(); plt.show()

From this result, can you guess why this reaction may produce light?

Answer There is a long region where ground and excited state surfaces cross, allowing population of the excited state of formaldehyde which may then emits light.

If you want, you can compare the curve to a Hartree-Fock ground state calculation. The reaction barrier predicted by restricted Hartree-Fock is significantly higher, since the molecule becomes open-shell during a significant portion of the reaction. Also, attempts to compute excitation energies at the TD-HF level have trouble converging.

Contents