7. Computing and visualization of Electrostatic Potential Map over molecular isodensity surface

In the previous post i explained how to compute and visualize an EPM over a Van der Waals surface using nwchem and jmol. However, actually it was found more realistic, theoretically and experimentally, to define the surface of a molecule as an isodensity surface, that is, a surface of equal electronic density, derived from the molecular orbitals influence.

A typically useful isodensity surface for estimating molecular shape and surface, is the one defined by a cutoff of 0.001 electron/Bohr³, that is, the one that encloses the space where the electron density is above such threshold. Such a cutoff encloses more than the 95% of the molecular electronic density.

Taking the case of l-proline again, the nwchem file for generating an electron density grid from the vector files generated in the previous steps would be something like:

dplot
  limitxyz
    -5.0 5.0 100
    -5.0 5.0 100
    -5.0 5.0 100
  gaussian
  output l-proline.den.cube
end
task dplot 

The limitxyz block here is used to define the limits and number of points for each axis. In order to estimate good limits for the grid, you can open the model in jmol and type on the console

jmol$ show boundbox

and substract two unities to the resuling first corner coordinates, and add two unities to the resulting second corner coordinate, and round it.

Once got the density grid output, file, in jmol console you can do:

jmol$ isosurface cutoff 0.001 "l-proline.den.cube" map "l-proline.elp.cube" jmol$ write isosurface "l-proline.den.jvxl"

Compare the result (first figure) with the one obtained previously with VDW surface (second figure):

6. Simple computing and visualization of molecular electron potential map

The most used kind of electrostatic potential map (EPM) is a mapping of the electrostatic potential generated by the charge distribution of the molecule, over an electron isodensity surface, but there is a simple and reasonable approximation to EPM using as mapping isosurface the Van der Waals surface of the atoms, and jmol is able to generate such isosurface.

The second requirement for building an EPM is to provide a volumetric grid of electrostatic potential (ESP) around the molecule. jmol does not generate it, but nwchem does it (see section 4). Using the generated data in section 4, we modify our input file l-proline.nw:

property
  esp
  grid
end
task scf property

save it and run nwchem:

linux$ nwchem l-proline

After process finished, some new files have been generated. One of them, a file named l-proline.elp.cube, contains the ESP volume grid and the optimized molecular geometry in cube format. Worth to mention that, although nwchem was always able to generate ESP grid data, dumping it in cube format is only available since version 6.2, so be sure this, or a recent one, is the version installed in your system¹

So next step is to generate the mol file from the optimized geometry contained in the cube file, and load with jmol:

linux$ babel l-proline.elp.cube l-proline.mol
linux$ jmol l-proline.mol

The final step is to generate the molecular surface coloured with the ESP projected on it. So, from jmol console:

jmol$ isosurface ID myepmsurface molecular map "l-proline.elp.cube"

The molecular keyword is just a minor fix that jmol makes over the pure Van der Waals (vdw) surface, in order to slightly fill the spaces between the Van der Waals spheres, thus generating a surface more similar to the one defined in terms of electronic isodensity.

Because cube files are very large even for small molecules and contains much more data than needed in order to display an isosurface, they are not the best suited for loading into a jmol applet from the web. Instead, you can dump the generated isosurfaces into a file, in jvxl file format (a jmol specific file for storing jmol generated data in more convenient way):

jmol$ write isosurface "l-proline.jvxl"

and load them from the applet, using its name as reference. Below we can see the result in a jmol applet initialized with the script

load http://dl.dropbox.com/u/xxxxxxxx/jmol/l-proline2.mol; isosurface color translucent 0.3 "http://dl.dropbox.com/u/xxxxxxxx/jmol/l-proline.jvxl"; spin on; 

The load command is not needed if you only want to see the isosurface. But you may want to make it slightly translucent in order to view the molecule inside it, as actually this script does:

Notes

[1] Check this document if you need to compile nwchem by yourself.