Quick and dirty electrostatic potential maps

In a previous post I showed how to compute an electrostatic potential map superimposed on the 0.002 isodensity surface of a molecule based on data computed using quantum chemical methods such as RHF/6-31G(d).

Previous posts (such as this one) have demonstrated that the van der Waals surface is a reasonable substitute for the 0.002 isodensity surface. Similarly, the electrostatic potential due to atomic charges can be a reasonable substitute for the electrostatic potential due to the electronic density.

The screencast above shows how to make such electrostatic potential maps using Avogadro and Jmol.

Avogadro uses an empirical method to determine the atomic charges (an integral part of the MMFF force field). It is possible to change the surface using the so-called "Iso Value" but I could not find any documentation on how that actually works. An Iso Value of 0 seems to correspond to the van der Waals sphere surface. It is currently not possible to alter the color range as far as I can see.

Jmol is not able to determine the charges, but the information can be transferred from Avogadro by saving a mol2 file. The electrostatic potential option in the Jmol menu corresponds to the following set of commands (as far as I can determine):

isosurface solvent color range -0.05 0.05 map mep
color isosurface translucent 0.5

and this set of commands can thus be used to control the color range and the nature of the surface. The default surface is a solvent accessible surface, which is slightly larger than the van der Waals surface.

A static and interactive version of the Avogadro and Jmol electrostatic potential map, respectively, can be found here.

Electrostatic potential map made with Avogadro
Click on the picture for an interactive version made with Jmol

Related blogpost: The polarity and solvation option in MolCalc

Steric strain vs electrostatic attraction

Figure 3.5. (a) RHF/6-31G(d) 0.002 au isodensity surface with superimposed electrostatic potential for (a) cis-HO(H)C=C(H)OH and (b) cis-CH3(H)C=C(H)CH3 and. In both cases, the maximum potential value is 0.05 au. (click on the picture for a bigger version).
From Molecular Modeling Basics CRC Press, May 2010.

This figure shows how the electrostatic potential superimposed on the 0.002 au isodensity surface can be used to rationalize why cis-HO(H)C=C(H)OH is more stable than trans-HO(H)C=C(H)OH, while the opposite is true for CH3(H)C=C(H)CH3.

The figure clearly shows the difference in polarity between hydroxyl and methyl groups. For the OH substitutent case the positive (O)H atom is close to the negative O(H) atom in the cis isomer. The CH3 group is non-polar and larger and the overlapping density indicates steric strain.

The figure is made with MacMolPlt as described in a previous post. Below is an interactive version made with Jmol (as described in a previous post). You'll notice that the color scheme is somewhat different, because Jmol maps the electrostatic potential value to a color in a different way than MacMolPlt. However, the general conclusion drawn from both programs is clearly the same.

Figure 3.5. (a) RHF/6-31G(d) 0.002 au isodensity surface with superimposed electrostatic potential for (a) cis-HO(H)C=C(H)OH and (b) cis-CH3(H)C=C(H)CH3 and.
Click on the picture for an interactive version

Electrostatic potential maps reloaded

Here is an interactive version of a figure I described in a previous post. Click on the picture to load it. Remember it's Jmol so you can rotate it and zoom as you like (Mac users: this works best with Safari).

Figure 3.4. The RHF/6-31G(d) electrostatic potential of water.
Click on the picture for an interactive version

The Jmol animation loads a file (fig3-2-2.xyz), which you can access by right-clicking on the Jmol animation once you have loaded it (as described in a previous post). This file also contains all the necessary commands.

Figure 3.4.c is the most common depiction of electrostatic potential maps and the Jmol general syntax for the command for this is

isosurface 0.002 "density.cube.gz" color range -0.05 0.05 "potential.cube.gz"

The screencast below shows how I made the cube files that contain the electron density and electrostatic potential information using MacMolPlt. The first part is identical to the screencast in a previous post on obtaining the electron density cube file.

The program I used to convert the MacMolPlt file to a cube was written by Jonathan Gutow and can be download here.

I start by loading the RHF/6-31G(d) optimized geometry I computed for a previous post, and reorienting. The orientation makes it easier to define the plotting plane for the contour plots.

Wishlist

Christmas is only 363 days away. It's never to early to start working on a new wish list.

[source]

[source]

[source]

[source]

Electrostatic potential maps

Figure 3.4. (a) Contour plot of the electrostatic potential of H₂O. The maximum and minimum contour value is 0.5 Hartrees/electronic charge (au). (b) The corresponding 0.05 au isopotential surface. (c) The electrostatic potential displayed on the 0.002 au isodensity surface of water. The maximum (darkest blue) value corresponds to 0.05 au.
From Molecular Modeling Basics CRC Press, May 2010.

Here is a screencast of how I made the figure:

The files I used were created in a previous post. The various values specified above are determined using trial and error, i.e. I kept fiddling with it until it looked good to me. When using plots like this it is important to specify these values, because using different values can lead to very different looking plots.

Also, different programs use different color scales and color intensity to indicate positive and negative charge, so it is rarely possible to directly compare electrostatic potential maps from different programs.

See this post about using screen capture to get a file with the graphic. This was done for each plot, and the combined in a word processor.

The interactive version of this figure is the subject of a future post.

Electron density doesn't always tell the whole story

Figure 3.3. (a) RHF/6-31G(d) 0.002 au isodensity surface and (b) van der Waals surface of cis-HO(H)C=C(H)OH (click on the picture for a bigger version).
From Molecular Modeling Basics CRC Press, May 2010.

Looking at the electron density or the van der Waals surfaces one would expect steric strain in cis-HO(H)C=C(H)OH. (See this post on how to make such pictures and this post on how to make an interactive version such as the one shown below.)

However, contrary to CH₃(H)C=C(H)CH₃, the cis isomer is more stable than the trans isomer. This is of course due to the electrostatic attraction of the O and H atom in the hydrogen bond. Here the electron density tells only part of the story because this attraction involves the nuclei as well. Here an electrostatic potential map is more useful. How to make such plots will be the subject of future posts.

Figure 3.3. The RHF/6-31G(d) electron density of cis-HO(H)C=C(H)OH.
Click on the picture for an interactive version

Common error messages in GAMESS: charge and multiplicity

  *** CHECK YOUR INPUT CHARGE AND MULTIPLICITY ***
THERE ARE    33 ELECTRONS, WITH CHARGE ICHARG=  0
BUT YOU SELECTED MULTIPLICITY MULT=  1

This error message was produced by the following input file. Can you see what's wrong?

 $contrl runtyp=optimize $end
 $basis gbasis=pm3 $end
 $data
Title
C1
N     7.0    -0.39094     1.95659     0.14008
H     1.0     0.38874     1.60529    -0.40413
H     1.0    -0.08386     2.76975     0.70945
H     1.0    -0.72485     1.22934     0.80007
H     1.0    -1.14035     2.30329    -0.48754
O     8.0    -0.64579     0.16732     2.03360
H     1.0    -0.26212    -0.73042     2.10569
H     1.0    -1.00756     0.26750     2.93979
O     8.0    -1.80535     3.31298    -1.59619
H     1.0    -1.39440     3.81065    -2.33214
H     1.0    -2.74148     3.57968    -1.71559
O     8.0     0.26578     4.05264     1.54485
H     1.0     1.03270     4.27032     2.11226
H     1.0    -0.26135     4.87344     1.64760
 $end

As GAMESS tells you, it is not possible to have singlet state (i.e. a multiplicity of 1) if you have an odd number of electrons. You need to tell GAMESS whether you want to change the multiplicity (for example, $contrl mult=2 scftyp=rohf $end) or add ($contrl icharg=-1 $end) or remove ($contrl icharg=1 $end).

When looking at the molecule it is pretty clear that you want to remove an electron, since ammonium has a positive charge. So the first line in the input should be changed to

 $contrl runtyp=optimize icharg=1 $end

Click on the picture for an interactive version

Notice that you don't have to specify that the positive charge is on ammonium. The quantum mechanics will take care of that.

Electron density and steric strain

Figure 3.2. (a) RHF/6-31G(d) 0.002 au isodensity surface and (b) van der Waals surface of cis-CH₃(H)C=C(H)CH₃ (click on the picture for a bigger version).
From Molecular Modeling Basics CRC Press, May 2010.

This figure shows how the electron density and van der Waals surfaces can be used to visualize steric strain. See this post on how to make such pictures and this post on how to make an interactive version such as the one shown below.

Figure 3.2. The RHF/6-31G(d) electron density of cis-CH₃(H)C=C(H)CH₃.
Click on the picture for an interactive version

Common error messages in GAMESS: error reading IDUM

 **** ERROR READING VARIABLE IDUM     CHECK COLUMN  1
H     1.0    -0.08386     2.76975     0.70945                                
....V....1....V....2....V....3....V....4....V....5....V....6....V....7....V....8

This error message was produced by the following input file. Can you see what's wrong?

 $contrl runtyp=optimize $end
$basis gbasis=pm3 $end
 $data
Title
C1
N     7.0    -0.39094     1.95659     0.14008
H     1.0     0.38874     1.60529    -0.40413
H     1.0    -0.08386     2.76975     0.70945
H     1.0    -0.72485     1.22934     0.80007
H     1.0    -1.14035     2.30329    -0.48754
O     8.0    -0.64579     0.16732     2.03360
H     1.0    -0.26212    -0.73042     2.10569
H     1.0    -1.00756     0.26750     2.93979
O     8.0    -1.80535     3.31298    -1.59619
H     1.0    -1.39440     3.81065    -2.33214
H     1.0    -2.74148     3.57968    -1.71559
O     8.0     0.26578     4.05264     1.54485
H     1.0     1.03270     4.27032     2.11226
H     1.0    -0.26135     4.87344     1.64760
 $end

The answer is the missing space in front of the $basis, and the solution is to add a space in front of it, just like for $contrl and $basis.

Why does GAMESS complain about a line in $data when the problem is with $basis? The reason is that it is possible to specify a basis set for each atom in $data, so GAMESS tries to read this information if no $basis group is found.

The error message would clearly benefit from the line: "Did you mean to specify a $basis group?"

The future of scientific publishing is here today

I recently came across this paper on embedding interactive 3D graphics in pdf files. Unfortunately, I couldn't animate the 3D figures, and I wrote the author (Vlad Vasilyev) who kindly wrote back that:

(1) the publisher (Springer) is a bit behind the times so the interactive figures are found in the supplemental material, but

(2) the ACS is more on the ball, as can be seen in this paper (this is the one I use in the screencast), and

(3) he has made a bunch of tutorials on how to create such diagrams on this page.

This page also contains some pdf files with 3D graphics you can play around with if you don't have access to the papers. Be sure to use Adobe Acrobat Reader 9.

Note that it is also possible to include animations such as vibrational motion (NB: 16 MB file)

Very cool! I haven't tried this myself, but you can be sure I will. Unfortunately, it seems that you need a piece of commercial software (Adobe Acrobat 9 Pro Extended) to do this.

On a related note, check you the interactive figure of interactive figures created by fellow blogger Henry Rzepa (please give the page some time to load).