Plotting the electrostatic potential with PyMOL

It has always been my aim to automatise the plotting of densities and orbitals. You can see my previous efforts in the context of VMD and Jmol in some previous posts, and you can find the associated scripts in the TheoDORE distribution. Let's turn to PyMOL now. The nice thing about PyMOL is that it can be scripted with python, which means that it is easy to add functionality in an integrated way. I started creating a toolbox for using PyMOL with quantum chemistry programs: qc_pymol, which you can find on github.

Let's say we have a cube file called es_1_diff.cube and we want to draw isosurface at isovalues of -0.01 and 0.01, then using qc_pymol we just type into the pymol console:

show_dens es_1_diff.cube, -0.01 0.01, cyan orange

And we get the following picture of the density. In this case, this is the difference density of an nπ* state with respect to the ground state (cyan is where the density is taken away, orange where it is added).

We can also map the electrostatic potential (ESP) onto the density if we have both as cube files. In this case, the command is

map_esp es_1_dens.cube, es_1_dens_esp.cube

Here, you can see that we have positive (blue) potential at the oxygen where the electron is taken away, and negative (yellow) potential at the other oxygen and at the carbon atom where electron density is added above.

We can do the same thing for the next three states and get the following combinations of difference densities and ESPs. We get the same trends as before depletion (cyan) in the difference density corresponds to positive ESP (blue), and addition (orange) corresponds to negative ESP (orange).