Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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:

enter image description here

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

enter image description here

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.

share|improve this question

2 Answers 2

up vote 8 down vote accepted

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.

share|improve this answer
    
This doesn't seem to be working with Version 10. Any ideas? –  paw Sep 12 at 7:47
    
@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. –  Szabolcs Sep 12 at 11:44
    
Ok, I'll try my best and add an answer if I find a solution. –  paw Sep 12 at 11:49

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

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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