3
$\begingroup$

The following is a description of a typical cube file

enter image description here

Based on the above instructions, I wrote the following mathematica code:

(* === 1. File Name === *)

(* The name and location of the file to be read is specified here. *)

fileName = "https://dl-web.dropbox.com/s/aikd8iaxxvntifv/C3H4O2-Potential.cube";


(* === 2. Read All Lines === *)

(* Read all lines from the specified file. *)

lines = ReadList[fileName, String];



(* === 3. Comment Line === *)

(* Determine the line number for the comment line. *)

commentLine = 2;


(* === 4. Atom Count === *)

(* Determine the number of atoms and the line numbers where the atom coordinates are found. *)

atomNum = (lines[[3]] // StringSplit // ToExpression)[[1]];



(* === 5. Grid Size === *)

(* Determine the size of the grid and the line numbers where the grid dimensions are found. *)

gridLines = lines[[4 ;; 6]];
grid = ToExpression /@ StringSplit /@ gridLines;
gridSize = grid[[All, 1]];



(* === 6. Data Lines === *)

(* Determine the line numbers where the potential energy data starts and ends. *)

dataStart = 2 + 1 + 3 + atomNum + 1;
dataEnd = Length[lines];
dataLines = lines[[dataStart ;; dataEnd]];
data = Flatten[ToExpression /@ StringSplit /@ dataLines];


(* === 7. Potential Energy Data === *)

(* Reshape the potential energy data into a 3D array. *)

potentialData = ArrayReshape[data, gridSize];



(* === 8. Contour Plot === *)

(* Create a 3D contour plot of the potential energy data. *)

ListContourPlot3D@potentialData

enter image description here

The result was terrible, and then I used a different code:

minValue = Min[potentialData];
maxValue = Max[potentialData];

(*生成等值面级别列表*)
contourLevels = Range[minValue, maxValue, (maxValue - minValue)/20];

(*使用等值面级别列表绘制等值面图*)
ListContourPlot3D[potentialData, Contours -> contourLevels, 
 PlotRange -> All]

enter image description here

It was still very bad, and then I changed my strategy

ListDensityPlot3D[potentialData, ColorFunction -> "TemperatureMap", 
 PlotRange -> All]

enter image description here

Looks much better, and much worse than gaussview produced

enter image description here

What can I do to match his performance, or even surpass his?

$\endgroup$
8
  • 3
    $\begingroup$ In the last figure you clearly see 2 contours. Why don't you do the same. You should also plot 2 contours not a huge number of them. $\endgroup$
    – yarchik
    Apr 20, 2023 at 8:21
  • $\begingroup$ I don't know which contour line to use at the moment@yarchik $\endgroup$ Apr 20, 2023 at 8:32
  • 1
    $\begingroup$ Just experiment. I am sure that gaussview does not produce the above figure with default settings. And, by the way, this is not a Mathematica question. $\endgroup$
    – yarchik
    Apr 20, 2023 at 8:36
  • $\begingroup$ But I'm using Gaussian by default@yarchik $\endgroup$ Apr 20, 2023 at 8:38
  • $\begingroup$ You use Gaussian for performing the calculations and gaussview for plotting. In gaussview one can change the contour values. And in MA you should do the same, just use Contours->{-0.1, 0.1} in ListContourPlot3D. $\endgroup$
    – yarchik
    Apr 20, 2023 at 8:46

1 Answer 1

7
$\begingroup$

The first issue you are having is this:

In[262]:= ToExpression["1.06897E-03"]

Out[262]= -0.0942383

It's unfortunate, but that string is interpreted as a number times the constant E minus another number. I made the following adjustment to your code:

data = Flatten[ImportString[StringRiffle[dataLines, "\n"], "Table"]];
potentialData = Transpose[ArrayReshape[data, gridSize], 1 <-> 3]

and it plots nicely. I remembered the Transpose trick from the last time I wrote code like this.

In recent Mathematica versions you can get the plot directly through Import:

Import["/Users/jasonb/Downloads/C3H4O2-Potential.cube", "Graphics3D"]

enter image description here

$\endgroup$
8
  • $\begingroup$ Using the default does work, but it's different from GaussView itself $\endgroup$ Apr 20, 2023 at 13:57
  • $\begingroup$ My manual purpose was to do the following molecular isosurface electrostatic potential diagram, mathematica does not have this function by default $\endgroup$ Apr 20, 2023 at 13:59
  • $\begingroup$ You can also use Import to get the underlying data and plot it using any custom styling - see this example. Choosing the right values for the isosurfaces can be a trial-and-error process $\endgroup$
    – Jason B.
    Apr 20, 2023 at 14:04
  • $\begingroup$ I've seen this, and it doesn't seem to have what I need $\endgroup$ Apr 20, 2023 at 14:05
  • $\begingroup$ you can look this picture sobereva.com/images/443/bar.jpg $\endgroup$ Apr 20, 2023 at 14:07

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