# How to calculate and plot molecular orbitals from a DFT/HF calculation using PySCF❓🤔 | TUTORIAL👨‍💻

In this tutorial, I will show you guys how to calculate and plot molecular orbtials (MO) after density functional theory (DFT) or Hartree-Fock (HF) calculations using the PySCF python package. Once the MO files are generated, we can visualize them using CrysX-3D Viewer or Jmol or Vesta.

Essentially, our end-goal is to be able to plot and visualize something like the following:

## Preliminaries

This tutorial assumes that you already have the PySCF package installed and working.
For details on how to install the PySCF package please refer to my following video:

Furthermore, it is also assumed that you know how to run DFT/HF calculations using PySCF. If not, I recommend you check out this blog post.

## Molecular Orbital plotting

Whenever we run a DFT calculation, we solve Kohn-Sham (KS) equations and we get MO energies and MO coefficients as the result. These MO coefficients are usually stored as a square matrix with each column containing the coefficients of the atomic orbitals for each MO. These can then be used to calculate the values of a particular MO at some points in space. For example, to create the figure shown above (Benzene LUMO), one would calculate the value of the LUMO using the corresponding coefficients on a cubic grid encompassing the Benzene molecule. These values are stored in a file, which can, later on, be read by visualization programs like CrysX-3D Viewer / Jmol / Vesta. Although there exist many file formats for storing such volumetric data, the most popular is the Gaussian CUBE format.

Fortunately, PySCF already comes with the utility to create such cube files.

The following code snippet runs a DFT calculation for a Benzene molecule and then uses the coefficients obtained from that to plot the 1st and 22nd (LUMO) of the Benzene molecule using the `cubegen` function in the `tools` module of `pyscf`. (Note: The 22nd MO corresponds to the LUMO because Benzene (C6H6) has 42 electrons, therefore 21 MOs are occupied, making 21st MO the HOMO and the 22nd MO the LUMO.)

```from pyscf import gto, dft, scf
from pyscf.tools import cubegen

benzene_coordinates = '''
C         -0.65914       -1.21034        3.98683
C          0.73798       -1.21034        4.02059
C         -1.35771       -0.00006        3.96990
C          1.43653       -0.00004        4.03741
C         -0.65915        1.21024        3.98685
C          0.73797        1.21024        4.02061
H         -1.20447       -2.15520        3.97369
H          1.28332       -2.15517        4.03382
H         -2.44839       -0.00006        3.94342
H          2.52722       -0.00004        4.06369
H         -1.20448        2.15509        3.97373
H          1.28330        2.15508        4.03386
'''

mol = gto.Mole()
mol.atom = benzene_coordinates
mol.basis = 'def2-SVP'
mol.build()

mf = dft.RKS(mol)
# mf = scf.RHF(mol) # For Hartree-Fock
mf.kernel()

# 1st MO
cubegen.orbital(mol, 'benzene_mo_1.cub', mf.mo_coeff[:,0])
# 22nd MO (LUMO)
cubegen.orbital(mol, 'benzene_mo_22.cub', mf.mo_coeff[:,21])
```

## Visualization

After running the above code, you will get two files (`benzene_mo_1.cub` and `benzene_mo_22.cub`) in your working directory. We can visualize them by dragging-and-dropping these files over the VESTA or Jmol window or use the File menu.

To make the orbitals look bigger or smaller, you can play with the isosurface values. Isosurfaces are drawn based on the given positive and negative isosurface cutoff values. A smaller absolute cutoff value would result in larger orbitals and vice versa.

You can change the isosurface cutoff value in VESTA as shown below:
Click on ‘Properties’ in the ‘Style’ tab and go to the ‘Isosurfaces’ tab in the new window that opens. Then you can adjust the Isosurface cutoff/level and the color by selecting the isosurface.

While VESTA and Jmol are good for quick visualization and analysis, I would recommend CrysX-3D Viewer for visually appealing and aesthetic visualizations as shown in the starting of this post, especially for manuscripts, posters, presentations, cover pages, etc. Furthermore, CrysX also works on Android smartphones and tablets, so you can download it from Google play here.
Here is how you would do it

## Related

You might also find the following video helpful if you want to plot the molecular orbitals and densities of molecular and periodic systems with TURBOMOLE

Also, if you’re interested in finding out how to visualize cube files with Jmol then check out the following YouTube tutorial

[wpedon id="7041" align="center"]

This site uses Akismet to reduce spam. Learn how your comment data is processed.