2
$\begingroup$

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:

enter image description here

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]

enter image description here

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.

$\endgroup$
2

1 Answer 1

2
$\begingroup$

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]

enter image description here

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:

enter image description here

$\endgroup$
3
  • $\begingroup$ The second part is a little hard to follow $\endgroup$ Apr 24, 2023 at 8:22
  • $\begingroup$ @我心永恒 Notice a easier way to rearrange the data with Transpose has been mentioned in this answer under your previous question. If you still feel it difficult to understand, feel free to ask in the comment. (But please be specific when asking. ) $\endgroup$
    – xzczd
    Apr 24, 2023 at 8:31
  • $\begingroup$ Let me delve into it first $\endgroup$ Apr 24, 2023 at 8:35

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.