10
$\begingroup$

Mathematica has the ability to plot a molecule using the data contained in an XYZ file. This is a simple text file, of which this is an example.

The molecule is plotted using the Import command. Entering this:

Import["https://gist.github.com/anonymous/8fde18207e8abc7b770c/raw/8491d99e7fec7cc989fe38548319e09d6910ba02/form_MO_28", "XYZ"] 

gives this as the output:

a molecule

I've been able to put this together with orbital plots, with results like this:

molecule + orbital plot

But I can only show other plots alongside the molecule plot, I can't modify the molecule. For example, I'd like to make the tubes forming the bonds smaller in diameter. As it is, the plots don't show if there is a single, double, or triple bond.

Somewhere in Mathematica, there must be a function which takes as its argument the XYZ text and returns the molecule plot. If I could see that function, I could modify it accordingly.

Any clue on how to access this function? The documentation for this feature is minimal.

$\endgroup$

2 Answers 2

9
$\begingroup$

To answer the question about accessing the function that does the plotting: the hints are here and in the SystemFiles/Formats/XYZ directory.

In[20]:= stream = OpenRead["ExampleData/caffeine.xyz"] Out[20]= InputStream["ExampleData/caffeine.xyz", 194] In[22]:= data = System`Convert`XYZDump`ImportXYZ[stream] Out[22]= {"VertexTypes" -> {"H", "N", "C", "N", "C", "C", "C", "N", "O", "C", "O", "H", "C", "H", "H", "H", "H", "H", "C", "C", "N", "H", "H", "H"}, "VertexCoordinates" -> {{-338.041, -112.724, 57.3304}, {96.683, -107.374, -81.9823}, {5.67293, 85.2719, 39.2316}, {-137.517, -102.122, -5.70552}, {-126.15, 25.9071, 52.3414}, {-30.6834, -168.363, -71.6934}, {113.942, 18.7412, -27.009}, {56.0263, 208.391, 82.5159}, {-49.268, -281.806, -120.947}, {-263.281, -173.04, \ -0.60953}, {-223.013, 79.8862, 108.997}, {254.97, 297.35, 62.2959}, {205.274, -173.609, -149.313}, {-248.077, -272.695, 48.8263}, {-300.89, -190.253, -104.98}, {291.761, -184.815, \ -78.5787}, {237.879, -112.119, -237.437}, {171.899, -274.899, \ -184.392}, {-15.1845, 309.7, 153.483}, {189.341, 211.812, 41.9319}, {228.613, 99.6844, -24.403}, {-16.8703, 404.366, 93.0109}, {35.3532, 329.791, 251.777}, {-120.745, 275.376, 172.03}}} 

This showed the format the next function expects:

In[23]:= ImportExport`MoleculePlot3D[data] Out[23]= --Graphics-- 

To get more hints about the available rendering options, check the XYZ format documentation and take a look at the option definition for Graphics`MoleculePlotDump`iMoleculePlot3D:

Options[iMoleculePlot3D] := Flatten[{"Rendering" -> Automatic, "AtomScaling" -> Automatic, ColorRules -> Automatic, "DrawBonds" -> Automatic, "DrawAtoms" -> Automatic, "InferBonds" -> Automatic, "InferBondsMinDistance" -> 40, "InferBondsTolerance" -> 25, "AtomSizeRules" -> Automatic, "BondSize" -> Automatic, "WireframeBonds" -> Automatic, "ImageSizeScaled" -> False, Boxed -> False, ViewPoint -> Automatic, ColorFunction -> Automatic, Lighting -> "Neutral", Sequence @@ Options["Graphics3D"]}] 

If you want to do more spelunking, take a look at the Spelunking package.

$\endgroup$
3
  • $\begingroup$ This doesn't seem to be working with Version 10. Any ideas? $\endgroup$ Commented Sep 12, 2014 at 7:47
  • $\begingroup$ @paw Sorry, not right away, and I'll be away from this site for a while starting today ... Try using Information or Names to see if those symbols changed context in v10. $\endgroup$ Commented Sep 12, 2014 at 11:44
  • $\begingroup$ Ok, I'll try my best and add an answer if I find a solution. $\endgroup$ Commented Sep 12, 2014 at 11:49
4
$\begingroup$

To add a little to Szabolcs's answer, one can get the bonds with

Graphics`MoleculePlotDump`ElementsToBonds 

Example:

mol = Import["ExampleData/Caffeine.xyz", "Rules"]; Graphics`MoleculePlotDump`ElementsToBonds[mol] (* {1 -> 10, 2 -> 7, 2 -> 6, 2 -> 13, 3 -> 8, 3 -> 7, 3 -> 5, 4 -> 5, 4 -> 6, 4 -> 10, 5 -> 11, 6 -> 9, 7 -> 21, 8 -> 20, 8 -> 19, 10 -> 15, 10 -> 14, 12 -> 20, 13 -> 18, 13 -> 16, 13 -> 17, 19 -> 23, 19 -> 22, 19 -> 24, 20 -> 21} *) 

We can turn it into a Graph:

g = Graph[ Range @ Length["VertexTypes" /. mol], UndirectedEdge @@@ Graphics`MoleculePlotDump`ElementsToBonds[mol] ]; 

We can construct a 3D plot of the molecule with GraphPlot3D:

With[{types = "VertexTypes" /. mol}, GraphPlot3D[ g, VertexCoordinateRules -> MapThread[Rule, {Range@Length@#, #} &["VertexCoordinates" /. mol]], VertexRenderingFunction -> ({Specularity[GrayLevel[1], 100], ColorData["Atoms"]@types[[#2]], Sphere[#1, Scaled[0.025]]} &), EdgeRenderingFunction -> ({Specularity[GrayLevel[1], 100], EdgeForm[], ColorData["Atoms"]@types[[First@#2]], Cylinder[{First[#1], Mean[#1]}, Scaled[0.015]], ColorData["Atoms"]@types[[Last@#2]], Cylinder[{Last[#1], Mean[#1]}, Scaled[0.015]]} &) ] ] 

Mathematica graphics

$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.