Obtain Total Energy - Distance Curves for different orbitals in Hydrogen H2 for use in General Chemistry Teaching

Dear PSI4 community,

I am not an expert in quantum chemical calculations. My objective is to adapt certain textbook plots to use them in teaching “General Chemistry”.
In particular, I would like to use PSI4 to calculate the total energy-distance curves for the hydrogen molecule H2. I succeed in obtaining the curve for the ground state (see the attached code which has been put together with the help of the available PSI4 tutorials). I can also directly obtain the energies of the molecular orbitals (not shown in the code), but what is needed is the total energy at a certain distance under the assumption that the electrons occupy higher orbitals.

Any suggestions on how I could obtain these values are greatly appreciated.

Thank you very much.

ps. The plot I would like to adapt is taken from Atkins / Friedmann (Molecular Quantum Mechanics).

import numpy as np
import pandas as pd

memory 600 mb

molecule h {

set reference uhf
h_energy, h_wfn = energy("ccsd(t)/cc-pVDZ", base_type="cp", return_wfn=True)

molecule h2 {
	H 1 R

Rvals = np.linspace(0.15,3.5,50)
set reference rhf

rows = []
table = []

for R in Rvals:
    h2.R = R
    h2_energy, h2wfn = energy("ccsd(t)/cc-pVDZ", base_type="cp", return_wfn = True)
    bonding = psi_hartree2kJmol * (h2_energy  - 2*h_energy)

df = pd.DataFrame(table, columns = ["bonding"], index = rows)

First, you can simplify ccsd(t)/cc-pVDZ to ccsd/cc-pVDZ and remove the base_type keyword. base_type is only useful for computations with multiple non-interacting fragments, and the (t) is only meaningful for systems with more than two electrons.

Technical note: wavefunctions are physical. We make the Born-Oppenheimer approximation to separate out the electronic and nuclear wavefunctions, and then make the Hartree-Fock approximation on the electronic wavefunction to get orbitals. Orbitals are (usually) an excellent approximation to the ground state wavefunction, but they tend to become poorer approximations as the energy levels increase. For this reason, “under the assumption that the electrons occupy higher orbitals” is less accurate than “for excited electronic states.”

Now for the how-to. Because H2/cc-pVDZ is a small system, we’ll get our roots via Full Configuration Interaction.

molecule {
symmetry c1
H 1 1.0 

set num_roots 11


This is the core of “get lots of roots, ignoring symmetry”. Your graph calls for 10 (remembering that there are two solutions for every pi state), so I grabbed 11. To get more than that, you need a larger basis set. It isn’t clear to me why Psi does that.