Forte Tutorial 5: Visualizing orbitals in Jupyter notebooks

In this tutorial we are going to explore how to generate orbitals and visualize them in Jupyter notebooks using the Python API.

Let's import the psi4 and forte modules, including the forte.utils submodule

Note: in this tutorial we do not call forte.startup()/cleanup() because we are not going to do any computation using forte

We will start by generating SCF orbitals via psi4 using the forte.util

Generating the cube files using the forte.utils.psi4_cubeprop function

Next we use a convenience function that generates cube files. Here is some information on it:

def psi4_cubeprop(wfn, path = '.', orbs = [], nocc = 3, nvir = 3, load = False):
    """
    Run a psi4 cubeprop computation to generate cube files from a given Wavefunction object
    By default this function plots from the HOMO -2 to the LUMO + 2

    Parameters
    ----------
    wfn : psi4Wavefunction
        A psi4 Wavefunction object
    path : str
        The path of the directory that will contain the cube files
    orbs : list
        The list of orbitals to convert to cube files (one based).
    nocc : int
        The number of occupied orbitals
    nvir : int
        The number of virtual orbitals

We want to plot two occupied and two virtual orbitals so we pass the arguments nocc=2 and nvir=2.

After we execute this command, four cube files are generated in the folder ./cubes. We can list them with the os.listdir function:

Reading cube files

To load these cube files we will use the CubeFile class. This class reads a cube file and lets us manipulate them. Let's load the HOMO, which is orbital 15 (the cube files use 1-based indexing)

A CubeFile object cannot be directly visualized. However, we can access all the information about the cube file.

The CubeFile class provides a convenient method (compute_levels) to determine the isocontour levels that contain a given fraction of the electron density. For example, to compute the isovalues that encompass 85% of the density we call compute_levels('mo',0.85), where the argument 'mo' indicates that the cube file stores the values of a molecular orbital (this is necessary, because the levels are computed differently for the density)

This information is also generated and stored by psi4 in the cube file (this is not standard, only psi4 implements this feature)

Rendering cube files

Forte implements a simple 3D rederer based on the pythreejs module. To install it via conda type in the command line: conda install -c conda-forge pythreejs. The renderer is an object of the Pt3JSRendered class (a python class)

Once created, we can add it a CubeFile object and display it. This is a Python widget and it is interactive. Try to interact with it by clicking and dragging the pointer around. You can also zoom the molecule in and out.

The CubeFile class will find the isolevels that correspond to 85% of the density. However, we can select a different value by passing the sumlevel argument.

We can even specify custom levels and colors via the levels and colors parameters

Rendering many cube files

Forte implements a few trick to plot several orbitals all at once. For example, if we pass the argument load=True to the function forte.utils.psi4_cubeprop we can automatically load all the cube file generated

There is also a convenience function (cube_viewer) that allows to plot several cube files al at once. This can be invoked with

This function is customizable, for example we can change the size and select a different font (or even not display the text)

Manipulating cube files

The CubeFile class supports three type of operations:

• scale(double factor): scale all the values on the grid by factor: $$ \phi(\mathbf{r}_i) \leftarrow \mathtt{factor} * \phi(\mathbf{r}_i) $$
• add(CubeFile cube): add to this cube file the grid values stored in cube: $$ \phi(\mathbf{r}_i) \leftarrow \phi(\mathbf{r}_i) + \psi(\mathbf{r}_i) $$
• pointwise_product(CubeFile cube): multiply each point of this cube file with the values stored in cube: $$ \phi(\mathbf{r}_i) \leftarrow \phi(\mathbf{r}_i) * \psi(\mathbf{r}_i) $$

To evaluate the density corresponing to the HOMO we can do just compute the pointwise product with itself

Then we can render it using the compute_levels function to find the appropriate levels. Here we cannot let the renderer determine the levels because it cannot tell that this is a density.

If we want, we can also save a cubefile to disk using the save function