Visualizing Frontier Orbitals with Psi4 and Python

Frontier Orbitals means the Highest Occupied Molecular Orbital (HOMO) and the Lowest Occupied Orbital (LUMO). Most Organic Chemists are mostly only interested in these orbitals because these are the point of interactions where an electron with high energy could jump into the next orbital. They are the dominant interactions of molecules that would most likely with each other. An example of a HOMO/LUMO diagram is shown in Figure 1.

Figure 1: Example of a HOMO/LUMO Molecular Orbital

However, we know the electron is a wave and there is an electron cloud that is produced in these orbitals. Visualizing these clouds gives us some idea into how the electron behaves and a memorable shape we could perhaps infer some things. To do this let’s go back to our psi4 methanol z-matrix, and as found in previous blog post. Let’s get right into the code:

import psi4

psi4.set_memory('1000mb')
psi4.core.set_num_threads(1)

zmatrix = '''\
    H11
    O11  H11  0.9600
    C11  O11  1.4000  H11  108.0000
    H12  C11  1.1000  O11  112.0000  H11  -60.0000
    H13  C11  1.1000  O11  112.0000  H11   60.0000
    H14  C11  1.1000  O11  111.8699  H11 -180.0000
'''

universe = psi4.geometry(zmatrix)
universe.update_geometry()
universe.print_out()

energy, wave_function = psi4.optimize(
  'hf/6-31g*',
   return_wfn =True,
  molecule=universe
)

Okay so that setups the molecule with optimized geometry. We would want it to return the wavefunction and we will use that to get our HOMO/LUMO.

So we want two things, the number of alpha electrons and the number of molecular orbitals because we need to determine which orbitals are the HOMO/LUMO pair. We use 2 because we adhere to Pauli Exclusion Principle.

number_of_occupied = 2
number_of_virutal = 2

alpha_electrons = wave_function.nalpha()
molecule = wave_function.nmo()

min_orb = max(1, alpha_electrons + 1 - number_of_occupied)
max_orb = min(molecule, alpha_electrons + number_of_virutal)
orbitals = [k for k in range(min_orb, max_orb + 1)]

Now, we get psi4 to produce those orbital files which comes in the form of a cube we also want the electron density to generate the electron cloud. We pass in our initial configurations and the tasks to be performed by psi4, meaning we want to generate the cloud in density form. We pass in which orbitals we would like to generate and then cubeprop initiates the generation of the files.

  cubeprop_tasks = []
  cubeprop_tasks.append('DENSITY')
  cubeprop_tasks.append('ORBITALS')

  psi4.set_options(
    {
      'scf_type': 'df',
      'g_convergence': 'gau_tight',
      'freeze_core': 'true',
      'CUBEPROP_TASKS': cubeprop_tasks,
      'CUBEPROP_FILEPATH': './',
      'CUBEPROP_ORBITALS': orbitals,
    }
  )

  psi4.cubeprop(wave_function)

So our files might look like this where we have the density and the cube file:

Da.cube
Psi_a_9_2-A.cube

Well it would be nice to visualize what this looks like and stick to python to remain consistency. We will be using these two packages which takes the psi4 points to grid and plots it in plotly.

pip install moly
pip install opt_einsum

Once that is installed, to add the cube file. You can add a lot of others too:

import moly

fig = moly.Figure()
fig.add_cube(
    'Psi_a_9_2-A.cube',
    iso=0.03,
    colorscale="rdbu",
    opacity=0.1
)

fig.show()

And we get something like this: