this question originates from an attempt to visualize molecular electrostatic potentials, but hopefully should be of general interest.

Suppose that you have two sets of data, S1 and S2. Each set consists of a certain number of 3D points (x,y and z coordinates) each one characterised by a certain (scalar) value. The 3D points are the same for the two datasets - the only difference is the scalar values associated with said points.

What I would like to do is to build a 3D isosurface from one of these two datasets (S1, by selecting a specific isovalue) and then map the second dataset (S2) onto this particular isosurface. In other words, to "color" the isosurface according to the S2.

Within a more "chemical" point of view, I would like to map the electrostatic potential (S2) onto a specific isosurface of the electronic density (S1). This can be done by a number of visualization tools (here you'll find more or less what should be obtained in this particular case), but I would really like to take advantage of Mathematica to play around with 3D countours and the like.

In principle, there exist a Wolfram Demonstration Project about this. However, there is very little Mathematica used in that example, all the tricky bits are taken care of by external software - including the mapping of the electrostatic potential onto the isosurface of the electronic density.

In practice, the starting point consists into two .cube files, which can be read into Mathematica and plot via a CubePlot function (stolen from here) which in turn takes advantage of ListContourPlot3D.

I attach the two relevant .cube files (S1 and S2) and a (rather embarassing...) notebook which goes as far as reading S1 and S2 from the two .cube files and plotting the S1 isosurface of interest. In principle it should be possible (he said confidently...) to select a specific isovalue of S1 via a specific value of Contours (this I have done) while coloring the surface via the S2 data... but that's where I'm presently stuck.

Any help would be much appreciated indeed - happy to discuss any of this in greater detail!


I am afraid this is not the end of the story.

Playing around a bit I realised that the distribution of the electrostatic potential (the color map) for the molecule in the example above was not symmetric as it should have been. Easy to miss that for a very small molecule, but for larger ones the issue becomes immediately clear.

This is the correct result (obtained by using VMD) for a slightly larger molecule:

enter image description here

While this is what I get using Mathematica (here you can find the notebook and the relavant files):

enter image description here

Obviously there is something very wrong going on in here. I have played around with a number of things. I believe the problem lies somewhere in the ColorFunction bit:

ColorFunction -> Function[{x, y, z}, ColorData["ThermometerColors", Rescale[esp[x, y, z], espRange]]]

which in turn relies on the interpolation of the electrostatic potential:

dataRangeESP = Map[MinMax, {espxg, espyg, espzg}]; esp = ListInterpolation[Transpose[espdata, {3, 2, 1}], dataRangeESP];

The isosurfaces for the electron density (and for the electrostatic potential itself, I have tried...) are represented correctly - it is "just" the coloring of said isosurfaces with the electrostatic potential that for some reason gives the wrong result.

Any help would be very much appreciated indeed!

  • $\begingroup$ I think what you need to do is leverage the ColorFunction option of ListContourPlot3D. You can make an interpolating function from the second cube file, and base the color function on that. $\endgroup$
    – Jason B.
    Jan 4, 2018 at 2:04
  • $\begingroup$ Hello - and many thanks for this. Your approach is basically what I have tried (with very little success) in the attached notebook: CubePlot[{cubedata, xg, yg, zg, xyz}, Contours -> {0.04}, ColorFunction -> (Function[{espxg, espyg, espzg, espdata}, Hue[espdata]][Rescale[#, {-0.3, 0.01}]] &), ColorFunctionScaling -> False], where cubedata and espdata refer to S1 and S2 respectively. No luck up to now. Could yuo perhaps elaborate on how to make an interpolating function from the second cube file, and base the color function on that ? Many thanks $\endgroup$
    – gcs
    Jan 4, 2018 at 9:14
  • $\begingroup$ @Jason B. In here an updated version of the notebook, where I managed to color a specific isosurface of S1 using a silly interpolation of a function of the x,y, and z of S2. Your comment and this post helped a lot. However, I still cannot map the values of S2 onto the S1 isosurface... any hint? Thanks! $\endgroup$
    – gcs
    Jan 4, 2018 at 17:01
  • $\begingroup$ give this a try, don't have time now for a proper answer, pastebin.com/raw/zWdbDTFn $\endgroup$
    – Jason B.
    Jan 4, 2018 at 18:39
  • $\begingroup$ @Jason B. Terrific - thanks! I am going to post an answer in a minute... $\endgroup$
    – gcs
    Jan 6, 2018 at 9:29

1 Answer 1


The solution takes advantage of Jason B.'s suggestions - hats off to him!

The full notebook is available here.

In summary, the two .cube files (containing the electron density and the electrostatic potential - on the same 3D grid) are read via a ReadCube function (stolen from here).

A particular isosurface of the electron density is then chosen via the Contours option of ListContourPlot3D. The electrostatic potential is mapped onto said surface via the ColorFunction directive of ListContourPlot3D.

Here follows the relevant snippet:

dataRange = Map[MinMax, {xg, yg, zg}]; (* the 3D range *)
espRange = {-0.03, 0.01}; (* the range of values for the electrostatic potential*)
esp = ListInterpolation[Transpose[espdata, {3, 2, 1}], dataRange]; (* Interpolate the electrostatic potential so that we can use it to color the electron density isosurface via the ColorFunction directive *)
vp = ViewPoint -> {1.3, -2.4, 2.0};
labels = Directive[FontSize -> 24, Black, FontFamily -> "Helvetica"];
densityPlot = 
Legended[ListContourPlot3D[Transpose[cubedata, {3, 2, 1}], 
   DataRange -> dataRange, Contours -> {0.001}, PlotRange -> Automatic,
   BoxRatios -> Automatic, Axes -> False, vp, Boxed -> False, 
   MeshStyle -> Directive[Thickness[0.0075], Blue]
   , ColorFunction -> 
    Function[{x, y, z}, 
     ColorData["NeonColors", Rescale[esp[x, y, z], espRange]]], 
   ColorFunctionScaling -> False, MeshFunctions -> {esp[#, #2, #3] &},
Mesh -> {{-0.03, -0.02, -0.01, 0.00, 0.01}}], 
BarLegend[{"NeonColors", espRange}, LabelStyle -> labels]]

And here is the final result:

enter image description here

To my (admittedly quite limited...) knowledge, this is an incredibly quick way to generate this sort of electrostatic potential maps - which look so much better if compared to the vast majority of what you get out of dedicated visualization software.

  • $\begingroup$ This is great! Thanks for sharing your result $\endgroup$
    – Jason B.
    Jan 6, 2018 at 16:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.