mathematica.blogoverflow.comAccess mathematicablogoverflowcom Mathematica - Stack

 current community chat Mathematica Mathematica Meta more communities Explore other Stack Exchange communities on stackexchange.com Stack Exchange This page is an archive of this blog, and is presented for historical purposes only. Plotting electronic orbitals using Mathematica 2013-09-27 by Jason B . 8 comments As a chemist it is often useful to plot electronic orbitals.  These are used to describe the wave function of electrons in atoms or molecules.  Typically, these are output from electronic structure software in the form of a cube file , first developed by Gaussian .  These files contain volumetric data for a given orbital on a three-dimensional grid. There exist many applications to visualize cube files, such as VMD  or GaussView , but I wanted to take advantage of Mathematica ‘s  capability to easily combine graphics, as well as the ability to automate the process in order to efficiently create frames for a movie . First off, we need a function to extract the data from the cube file.  In the process, we will create the text for an XYZ file , a format also developed by Gaussian.  The function OutForm is used here to mimic the printf function found in other programming languages. OutForm[num_?NumericQ, width_Integer, ndig_Integer, OptionsPattern[]] := Module[{mant, exp, val}, {mant, exp} = MantissaExponent[num]; mant = ToString[NumberForm[mant, {ndig, ndig}]]; exp = If[Sign[exp] == -1, "-", "+"] <> IntegerString[exp, 10, 2]; val = mant <> "E" <> exp; StringJoin@PadLeft[Characters[val], width, " "] ]; ReadCube[cubeFileName_?StringQ] := Module[ {moltxt, nAtoms, lowerCorner, nx, ny, nz, xstep, ystep, zstep, atoms, desc1, desc2, xyzText, cubeDat, xgrid, ygrid, zgrid, dummy1, dummy2, atomicNumber, atomx, atomy, atomz, tmpString, headerTxt,bohr2angstrom}, bohr2angstrom = 0.529177249; moltxt = OpenRead[cubeFileName]; desc1 = Read[moltxt, String]; desc2 = Read[moltxt, String]; lowerCorner = {0, 0, 0}; {nAtoms, lowerCorner[[1]], lowerCorner[[2]], lowerCorner[[3]]} = Read[moltxt, String] // ImportString[#, "Table"][[1]] &; xyzText = ToString[nAtoms] <> "\n"; xyzText = xyzText <> desc1 <> desc2 <> "\n"; {nx, xstep, dummy1, dummy2} = Read[moltxt, String] // ImportString[#, "Table"][[1]] &; {ny, dummy1, ystep, dummy2} = Read[moltxt, String] // ImportString[#, "Table"][[1]] &; {nz, dummy1, dummy2, zstep} = Read[moltxt, String] // ImportString[#, "Table"][[1]] &; Do[ {atomicNumber, dummy1, atomx, atomy, atomz} = Read[moltxt, String] // ImportString[#, "Table"][[1]] &; xyzText = If[Sign[lowerCorner[[1]]] == 1, xyzText <> ElementData[atomicNumber, "Abbreviation"] <> OutForm[atomx, 17, 7] <> OutForm[atomy, 17, 7] <> OutForm[atomz, 17, 7] <> "\n", xyzText <> ElementData[atomicNumber, "Abbreviation"] <> OutForm[bohr2angstrom atomx, 17, 7] <> OutForm[bohr2angstrom atomy, 17, 7] <> OutForm[bohr2angstrom atomz, 17, 7] <> "\n"]; , {nAtoms}]; cubeDat = Partition[Partition[ReadList[moltxt, Number, nx ny nz], nz], ny]; Close[moltxt]; moltxt = OpenRead[cubeFileName]; headerTxt = Read[moltxt, Table[String, {2 + 4 + nAtoms}]]; Close[moltxt]; headerTxt = StringJoin@Riffle[headerTxt, "\n"]; xgrid = Range[lowerCorner[[1]], lowerCorner[[1]] + xstep (nx - 1), xstep]; ygrid = Range[lowerCorner[[2]], lowerCorner[[2]] + ystep (ny - 1), ystep]; zgrid = Range[lowerCorner[[3]], lowerCorner[[3]] + zstep (nz - 1), zstep]; {cubeDat, xgrid, ygrid, zgrid, xyzText, headerTxt} ]; If you need to create a cube file, then the following function can be used: WriteCube[cubeFileName_?StringQ, headerTxt_?StringQ, cubeData_] := Module[{stream}, stream = OpenWrite[cubeFileName, FormatType -> FortranForm]; WriteString[stream, headerTxt, "\n"]; Map[WriteString[stream, ##, "\n"] & @@ Riffle[ScientificForm[#, {3, 4}, NumberFormat -> (Row[{#1, "E", If[#3 == "", "+00", #3], "\t"}] &), NumberPadding -> {"", "0"}, NumberSigns -> {"-", " "}] & /@ #, "\n", {7, -1, 7}] &, cubeData, {2}]; Close[stream];] Next we need the function to plot the orbital, CubePlot[{cub_, xg_, yg_, zg_, xyz_}, plotopts : OptionsPattern[]] := Module[{xyzplot, bohr2picometer, datarange3D, pr}, bohr2picometer = 52.9177249; datarange3D = bohr2picometer {{xg[[1]], xg[[-1]]}, {yg[[1]], yg[[-1]]}, {zg[[1]], zg[[-1]]}}; xyzplot = ImportString[xyz, "XYZ"]; Show[xyzplot, ListContourPlot3D[Transpose[cub, {3, 2, 1}], Evaluate[FilterRules[{plotopts}, Options[ListContourPlot3D]]], Contours -> {-.02, .02}, ContourStyle -> {Blue, Red}, DataRange -> datarange3D, MeshStyle -> Gray, Lighting -> {{"Ambient", White}}], Evaluate[ FilterRules[{plotopts}, {ViewPoint, ViewVertical, ImageSize}]]] ]; Let’s look at an example. Download the file cys-MO35.cube  and place it in your home directory (or anywhere in your $Path ). Then, read in the cube file with: {cubedata,xg,yg,zg,xyz,header}= ReadCube["cys-MO35.cube"]; and plot it via CubePlot[{cubedata, xg, yg, zg, xyz}] When I want to create a movie file, I want all the images to have exactly the same ViewAngle , ViewPoint , and ViewCenter .  When you give these options to CubePlot , it feeds them directly to the Show function vp = {ViewCenter -> {0.5, 0.5, 0.5}, ViewPoint -> {1.072, 0.665, -3.13}, ViewVertical -> {0.443, 0.2477, 1.527}}; CubePlot[{cubedata, xg, yg, zg, xyz}, vp] Finally, you can also give any options that normally go to ListContourPlot3D CubePlot[{cubedata, xg, yg, zg, xyz}, vp, ContourStyle -> {Texture[ExampleData[{"ColorTexture", "Vavona"}]], Texture[ExampleData[{"ColorTexture", "Amboyna"}]]}, Contours -> {-.015, .015}] Many thanks to Daniel Healion for the ReadCube and WriteCube functions. Filed under chemistry , graphics Wolfram Technology Conference 2012 2012-10-29 by Yves Klett . 2 comments The Wolfram Technology Conference took place from 2012-10-17 to 2012-10-19 in Champaign, IL. This is a loose collection of whatever interesting/entertaining stuff I came across during the conference. Note that I had to sign a non-disclosure agreement (NDA) to attend sessions on future Mathematica releases and upcoming Wolfram technology products, so there will be no infos on on that (make of that what you will). Some general info for those who are not familiar with the Tech Conference: There are several types of talks, and you can easily cram your schedule bumper to bumper: WRI overview talks (presenting some topic, e.g. “ Mathematica Connectivity” or “Image Processing” in general) WRI in-depth talks (taking on a specific area, e.g. “Manipulate Secrets”) Hands-on workshops (big this year: SystemModeler ) User talks that cover the wide range of Mathematica application in education, science, industry and entertainment. See the 2011 schedule and the 2012 videos to get an impression. The Champaign Hilton Garden Inn is a decent venue, the only thing to be aware of is the sometimes lethal airconditioning, so be sure to bring warm clothes. The focus of most presentations is on Mathematica , but since it is the Wolfram Technology Conference, things like Wolfram Alpha, SystemModeler and other technologies start to feature more and more prominently. There is an exceptional density of Wolfram developers (and thus WRI competence) and you can set up meetings with them and other participants with the online conference system (crowdvine) or just try to taser and drag them into the coffee room downstairs (there is also a merchandise and book store offering conference discounts). The program  is quite packed with three to four tracks and other things like meet-ups, lunch round tables on certain subjects and similar. On all evenings there is some kind social event either on-site or in Champaign at large. The crowd is a very easygoing one with a quite high proportion of regulars. 2012-10-15: Landing at ORD. 2012-10-16: Yawn. The good thing about flying west is that getting up early (really early) is easy. Staying awake will be ...