Plotting the molecular isosurface electrostatic potential map

Question

After performing Gaussian calculation, cube files can be obtained. The distributions are the potential cube and esp cube, both of which contain volume data. What I want to do now is to project the potential cube onto the esp cube file, which is the "molecular isosurface electrostatic potential map".

"molecular isosurface electrostatic potential map" like this:

This image was created using the VESTA crystallography software. Although it is quite good, it cannot meet some of my specific requirements.

Below is a part of the code I am working on, with comments.

(* === 1. Define file path === *)
path1 = "http://pic.wxyh.top/C3H4O2-esp.cube";
path2 = "http://pic.wxyh.top/C3H4O2-Potential.cube";

(* === 2. Import volumetric data === *)
espCube = Import[path1, "VolumetricData"]; (* Import esp.cube file *)
potentialCube = Import[path2, "VolumetricData"]; (* Import Potential.cube file *)

(* === 3. Get data information === *)
espData = espCube["Data"]; (* Extract data from esp.cube file *)
potentialData = potentialCube["Data"]; (* Extract data from Potential.cube file *)
dataRange = espCube["DataRange"]; (* Get data range *)

(* === 4. Process data === *)
espDataNormal = Normal[espData]; (* Normalize espData *)
potentialDataNormal = Normal[potentialData]; (* Normalize potentialData *)

(* One of the processing steps for isosurface electrostatic potential maps that calculates the electrostatic interaction energy for each ion *)
projectedData = MapThread[Times, {espDataNormal, potentialDataNormal}, 3]

projectedDataNormal = Normal[projectedData]

The data obtained in the end is of this type. I know that ListContourPlot3D can be used to plot the graph later, but I don't know how to process this data. I'm not sure how to convert it into a three-dimensional coordinate format. The data is too long and I'm feeling a bit overwhelmed.

Answer 1

As mentioned in the comment above, the problem is strongly related to

How can I convert a 3D contour plot into a region for use as a custom slice in ListSliceDensityPlot3D?

And you don't need any advanced technique to obtain the desired plot:

mole = Import[path1, "Molecule"];

moleplot = MoleculePlot3D@mole
(* Definition of espData, etc. are the same as yours. *)
surface = ListContourPlot3D[espData[[1]], Contours -> 0.03 {-1, 1}, 
   DataRange -> dataRange] // DiscretizeGraphics

potentialPlot = 
 ListSliceDensityPlot3D[potentialData[[1]], surface, DataRange -> dataRange, 
  ColorFunction -> ({Opacity[0.8], Blend["M10DefaultDensityGradient", #1]} &)]

Show[moleplot, potentialPlot]

Oh, perhaps the M10DefaultDensityGradient should be regarded as advanced, which is discussed in e.g.

What's the default colormap (or color scheme) used in mathematica?

But it's not necessary.

---

It's possible to achieve the goal with ListContourPlot3D only, but this solution is a bit advanced:

potefunc = 
  ListInterpolation[Normal@Transpose[potentialData[[1]], {3, 2, 1}], dataRange];

help[opt___] := 
 ListContourPlot3D[espData[[1]], Contours -> 0.03 {-1, 1}, DataRange -> dataRange, 
  opt]

range = 
 potefunc @@@ 
   First@Cases[help[], GraphicsComplex[p : {{_, _, _} ..}, __] :> p, Infinity] // 
  MinMax
(* {-0.0615143, 0.366115} *)

moleplot~Show~
 help[ColorFunctionScaling -> False, 
  ColorFunction -> 
   Function[{x, y, z, f}, {Opacity[0.8], 
     Blend["M10DefaultDensityGradient", Rescale[potefunc[x, y, z], range]]}], 
  Mesh -> None, PlotRange -> All]

The advantage of this solution is, the obtained surface is smoother: