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
References
- Linear combination of atomic orbitals – Wikipedia
- Molecular orbital – Wikipedia
- http://paulbourke.net/dataformats/cube/
- https://h5cube-spec.readthedocs.io/en/latest/cubeformat.html
- https://gaussian.com/cubegen/
- https://docs.chemaxon.com/display/docs/gaussian-cube-format.md
- orbitals (uni-muenchen.de)
- Tutorial: Display of Orbitals and Molecular Surfaces (marshall.edu)
- Visualization and Analysis of Quantum Chemical and Molecular Dynamics Data with VMD (ruhr-uni-bochum.de)
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
I’m a physicist specializing in computational material science with a PhD in Physics from Friedrich-Schiller University Jena, Germany. 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.