Academic Articles & Supplements

Plotting Complex Orbitals from Cube File Data

Rulin Feng, Xiaojuan Yu, Jochen Autschbach

Author

Rulin Feng, Xiaojuan Yu, Jochen Autschbach

Title

Plotting Complex Orbitals from Cube File Data

Description

Supplemental notebook to "Spin–Orbit Natural Transition Orbitals and Spin-Forbidden Transitions"

Plotting Complex Orbitals from Cube File Data

Given two cube file, one with the real part, one with the imaginary part, of a complex molecular orbital, this notebook will plot an isosurface of the absolute value of the orbital and color the surface with the complex phase.

If you use this notebook, or code based on it, for your own research, please cite the reference below and provide the URL from where you obtained this notebook.

In[]:=

SetDirectory[NotebookDirectory[]];

Some default settings. Change as needed.

In[]:=

opa=0.7;(*defaultopacity*)spec=50;(*defaultspecularity*)(*pluss=ColorData["RedBlueTones",0.15];minuss=ColorData["RedBlueTones",0.85];*)(*libertyBlue=RGBColor[0,51/255,153/255];myOrange=RGBColor[1.0,0.588,0.](*Sameas#FF9600usedinadfGUI*);*)myGreen=ColorData["HTML"]["Green"];(*myYellow=ColorData["HTML"]["Yellow"];*)vp=ViewPoint/.Options[Graphics3D];vv=ViewVertical/.Options[Graphics3D];

Define the cube file on the next line, without the file name extension
.cube. You can also define the file name extension, to accommodate .cub vs .cube

In[]:=

cubefileRe="2022-02-471gnt6_Supplements/orbital-re";cubefileIm="2022-02-471gnt6_Supplements/orbital-im";ext="cube";

The next set of commands import the cube file data and visualize the molecule that’s defined in the cube file. Available PlotThemes are
HeavyAtom
BallAndStick
SpaceFilling
Tubes
WireFrame

In[]:=

cube=cubefileRe<>"."<>ext;tmp=Import[cube,"Cube"];Export["tmp.xyz",tmp];mol=Import["tmp.xyz","InferBondTypes"False];volRe=Import[cube,"VolumetricData"];cube=cubefileIm<>"."<>ext;volIm=Import[cube,"VolumetricData"];molplot=MoleculePlot3D[mol,PlotTheme"Tubes"]

Out[]=

Check that data ranges in the Re and Im cube file are the same. Otherwise, we definitely have a problem with incompatible grids for the two cube files.

In[]:=

If[volRe["DataRange"]≠volIm["DataRange"],Print[Style["WARNING! Incompatible Data ranges. DO NOT CONTINUE unless you know exactly what you're doing",Large,Bold,Red]],Print[Style["Data ranges are compatible. Continuing.",Bold,myGreen]]]

Data ranges are compatible. Continuing.

Now we create the data for the absolute value and phase of the orbital, and grid specs

In[]:=

volComplex=Normal[volRe[["Data",1]]]+INormal[volIm[["Data",1]]];volRange=volRe["DataRange"];volAbs=Abs[volComplex];volArg=Arg[volComplex]/(2Pi);n=Table[Length[volAbs[[i]]],{i,1,3}];min=Table[volRange[[i,1]],{i,1,3}];max=Table[volRange[[i,2]],{i,1,3}];step=Table[(max[[i]]-min[[i]])/(n[[i]]-1),{i,1,3}];

next: should be the step sizes shown in the cube file, in case you want to cross check that

In[]:=

step/0.529177

Out[]=

{0.16,0.16,0.16}

create interpolation function for the array with the phases (this step can be a bit slow)
The idea comes from
https://mathematica.stackexchange.com/questions/105623/color-a-single-contour-from-a-listcontourplot3d-using-values-from-a-list
(the second option from that posting uses ListSliceDensityPlot3D, but there seems to be no easy option to set the opacity of the surface, and the plot takes longer to make)
Note that Mathematica 12.3 gives a warning(?) with that example

Function

:Too many parameters in {x,y,z} to be filled from Function[{x,y,z},Hue[func[x,y,z]]][0.3].

when you execute the example from that web page, but the result looks exactly like the one on StackExchange. We get the same warning below, but the plots definitely look OK.

In[]:=

cfunc=Interpolation[Flatten[Table[{min[[1]]+(x-1)step[[1]],min[[2]]+(y-1)step[[2]],min[[3]]+(z-1)step[[3]],volArg[[z,y,x]]},{x,1,n[[1]]},{y,1,n[[2]]},{z,1,n[[3]]}],2]]

Out[]=

InterpolatingFunction

Domain: {{-4.23,4.15},{-4.23,4.15},{-4.23,4.15}}

Output: scalar

Data not saved. Save now



In[]:=

cfunc[1,2,3]

Out[]=

0.250242

Define the desired contour (iso-surface) values on the next line, then
create the isosurfaces and combine with the molecule plot.

In[]:=

contours=0.03;combo=Show[ListContourPlot3D[volAbs,DataRange->volRange,BoxRatiosAutomatic,Contours{contours},ContourStyle{Directive[Opacity[opa],Specularity[White,spec]]},ColorFunction->Function[{x,y,z},Hue[cfunc[x,y,z]]],ColorFunctionScaling->False],molplot,Lighting"Neutral",BoxedFalse,BoxStyleDirective[Thickness[0.01],Gray],AxesFalse,MeshTrue,AxesLabelNone,TicksNone,PlotRangeAll,ViewPointDynamic[vp],ViewVerticalDynamic[vv]]

Function

:Too many parameters in {x,y,z} to be filled from Function[{x,y,z},Hue[cfunc[x,y,z]]][0.03].

Out[]=

Rotate unless you are happy with the display, then export. After the Export, white space should be removed from the plot, e.g. with an ImageMagick tool:
mogrify -trim <filename>

In[]:=

Export["orbitalplot.png",combo,"PNG",ImageResolution300];

Cite this as: Rulin Feng, Xiaojuan Yu, Jochen Autschbach, "Plotting Complex Orbitals from Cube File Data" from the Notebook Archive (2022), https://notebookarchive.org/2022-02-471gnt6