Friday 13 September 2019

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).


Unknown said...

Thank you for such a great tool!

How do you change your esp map surface colors such that they correspond to the updated ramp colors? Your default ramp is [-0.4, 0.4], and I found that I need a smaller range. I found a way to update the ramp (ramp_update e_pot_color, range=[-0.1,0,0.1], but not the surface colors.

Felix said...

Hi, nice to hear that you found this. What I did is shown here

I think the commands you are looking for are
stored.ramp_vals = [-0.35,-0.08,-0.02, 0, 0.02, 0.08, 0.35]
stored.ramp_cols = ["green", "marine", "cyan", "white", "yellow", "tv_red", "magenta"]

I wrote this inside the script. But I guess it would also work entering it directly from the command line in pymol (but then the syntax may be a bit different)


Xinglong said...

Hi Felix,

Nice post! For pymol, are we able to plot the density and esp difference between the ground state and the excited state with dens.cube and esp.cube from each state? In this case, do we need to preprocess them to find the dens_diff.cube and esp_diff.cube that captures the differences of the two states?

In your tutorial, how did you generate es_1_diff.cube?

Thank you very much!

Best wishes,

Felix said...


Nice to hear that you are interested in this!

This specific job was run in Q-Chem following Example 10.8 from here
This produces all the files as they are needed.

Otherwise, yes, the density and ESP are both additive quantities. You can get the difference density by subtracting the cube files of excited state and ground state density. And you can get the ESP of the difference density by subtracting their ESPs.

BTW: These plots come from this paper

Xinglong said...

Thanks Felix for your very quick response! Just wondering if PyMol can allow us to load the individual cube files and then do the subtraction directly from the loaded cube files?


Felix said...

Hi Xinglong,

I have a little python Utility for these purposes, here

Then you'd need a bit of python along the lines of

from theodore import lib_util
dens1 = lib_util.cube_file('dens1.cube')
dens2 = lib_util.cube_file('dens2.cube')
dens1.lin_comb(dens2, 1., -1., 'diff.cube').write()

Otherwise, you'd probably find the same kind of utilities in pymol or VMD. But I couldn't tell you exactly how to do it.


Xinglong said...

Thanks Felix,

I will take a look and try it to see it goes!