How to calculate and plot molecular orbitals from a DFT/HF calculation using PySCF❓🤔

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:

Benzene Lowest Unoccupied Molecular Orbital Visualization

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.

Visualizing MO using VESTA by dragging-and-dropping

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.

Changing the isosurface cutoff/level value in VESTA

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

References

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

Get Article as PDF

Manas Sharma

Ph.D. researcher at Friedrich-Schiller University Jena, Germany. I’m a physicist specializing in computational material science. I write efficient codes for simulating light-matter interactions at atomic scales. I like to develop Physics, DFT, and Machine Learning related apps and software from time to time. Can code in most of the popular languages. I like to share my knowledge in Physics and applications using this Blog and a YouTube channel.

manas.bragitoff.com/