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)