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.

Improve this question
$\endgroup$
2
  • $\begingroup$ Stongly related: mathematica.stackexchange.com/q/105642/1871 $\endgroup$ Commented Apr 24, 2023 at 4:27
  • $\begingroup$ Now the problem is dealing with the data, I don't know how to deal with it, too long data@xzczd $\endgroup$ Commented Apr 24, 2023 at 6:08

1 Answer 1

Reset to default
3
$\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

Improve this answer
$\endgroup$
3
  • $\begingroup$ The second part is a little hard to follow $\endgroup$ Commented 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$ Commented Apr 24, 2023 at 8:31
  • $\begingroup$ Let me delve into it first $\endgroup$ Commented Apr 24, 2023 at 8:35

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.