10
$\begingroup$

Mathematica's ChemicalData has a lot of useful information, but many compounds only have flat 2D structure diagrams rather than the cool 3D "MoleculePlot".

I thought I could fake it using GraphPlot3D

GraphPlot3D[ChemicalData["VanadiumIVOxide", "EdgeRules"], 
 EdgeRenderingFunction -> (Cylinder[#1, .05] &), 
 VertexRenderingFunction -> ({MapIndexed[
      First@#2 -> ColorData["Atoms"][#1] &, 
      ChemicalData["VanadiumIVOxide", "VertexTypes"]], 
     Sphere[#1, .15]} &), Boxed -> False]

This is close, but no cigar:

enter image description here

The colouring of the nodes doesn't carry across from the VertexRenderingFunction option, which is probably a lighting issue.

Does anyone have any suggestions for replicating the "MoleculePlot" look more closely?

Improve this question
$\endgroup$
10
  • $\begingroup$ I get Missing["NotAvailable"] from ChemicalData["VanadiumIVOxide", "EdgeRules"] on v10.4.1 $\endgroup$
    – Kuba ♦
    May 6, 2016 at 7:05
  • $\begingroup$ "VanadiumVOxide" works though. $\endgroup$
    – Kuba ♦
    May 6, 2016 at 7:06
  • $\begingroup$ @Kuba - just use a compound that has these available. Odd that I have it (on v9) and you don't. I'll check v10.4 when I get home. $\endgroup$
    – Verbeia
    May 6, 2016 at 7:07
  • 1
    $\begingroup$ These answers are all so awesome that I'm embarrassed that I actually asked this for completely flippant reasons related to a draft TV script that will probably never come to light. I've ended up accepting the one that would look prettiest on TV. $\endgroup$
    – Verbeia
    May 7, 2016 at 4:59
  • 1
    $\begingroup$ I'm glad you clarified the rationale for the question. It needs to be stated, however, that the accepted answer does not produce the correct chemical structure, just in case future visitors happen upon this question for less theatrical purposes. $\endgroup$
    – bobthechemist
    May 7, 2016 at 11:26

6 Answers 6

Reset to default
10
$\begingroup$

There is a problem with syntax in VertexRenderingFunction, can't explain more because I don't know what was the goal there.

With[{
  atoms = ChemicalData["Valeraldehyde", "VertexTypes"] 
  }
  ,
  GraphPlot3D[ ChemicalData["Valeraldehyde", "EdgeRules"],
      EdgeRenderingFunction   -> (
         { Specularity[White, 100], Cylinder[#1, .05] }&
      ),
      VertexRenderingFunction -> Function[{pos, name, nr}, 
         {ColorData["Atoms"][atoms[[name]]], Specularity[White, 100], 
          Sphere[pos, .15]
      }],
      Boxed                   -> False
  ]
 ]

enter image description here

Graph3D may be of more use since it supports GraphLayout:

With[{
  atoms = ChemicalData["Valeraldehyde", "VertexTypes"]}
 ,
 Graph3D[
  ChemicalData["Valeraldehyde", "EdgeRules"],
  EdgeShapeFunction -> ({Specularity[White, 100], Cylinder[#1, .1]} &),
  VertexShapeFunction -> 
   Function[{pos, name, nr}, {ColorData["Atoms"][atoms[[name]]], 
     Specularity[White, 100], Sphere[pos, .25]}]
  ,
  GraphLayout -> "SpringEmbedding",

  Boxed -> False,
  Lighting -> "Neutral"
  ]
 ]

enter image description here

Improve this answer
$\endgroup$
3
  • $\begingroup$ I would say that SpringEmbedding is closer to the real structure than SpringElectricalEmbedding $\endgroup$
    – Lucas
    May 6, 2016 at 8:42
  • $\begingroup$ @Lucas Yep, thanks :) $\endgroup$
    – Kuba ♦
    May 6, 2016 at 9:04
  • 2
    $\begingroup$ Using Method -> {"SpringEmbedding", "InferentialDistance" -> 3.} on GraphPlot3D[] yields a reasonable-looking (in the sense of VSEPR) structure for pentanal. $\endgroup$
    – J. M.'s eventual burnout ♦
    May 7, 2016 at 1:13
16
$\begingroup$

Edit Introduction of the Molecule symbol in v12 has broken the package used in this answer. It can be fixed by converting the StructuredArray of atom coordinates into coordinate lists in picometers. I should fix this in the package; however, at the moment I have changed the code in this answer to reflect this hack.

The short (and probably correct) answer is to download Mercury and get the CIF file from the Crystallograph Open Database. That way you can focus on your science and less on the, well fitting a square peg in a round hole.

If you really want that square peg to fit in that round hole, however.

  • First, recognize that Mathematica does not yet get a passing grade in chemical structure. To do so, it needs to understand that molecules do come in three dimensions and atoms can be multiply bonded to one another. A solution needs to address these two issues.

  • second, you need good data. Since the correct molecular coordinates are rarely provided by ChemicalData, we need another source. SDF, MOL, XYZ (and possibly CIFs, although they are much more complicated) contain the type of information we need.

  • Third, with good data in hand, creating graphics complexes of molecules turns out to be not that difficult, and we can use curated data to help us out with colors and atom sizes.

I've created a proof-of-concept molecule viewer package on github. I won't repost the code here, just some highlights. Please note that the code is sloppy, and was originally written to help me make an NMR interpretation instructional video. It wasn't meant for public consumption but I think it gets the point across, nonetheless.

Grab the package (after looking at it) from within Mathematica and import an example chemical structure.

(* View the file first, so you know what you are installing*)<< \
"https://raw.githubusercontent.com/bobthechemist/molviewer/master/\
Molviewer.m"
(* Load sample file *)
f = "https://raw.githubusercontent.com/bobthechemist/molviewer/master/\
623-70-1.sdf";
Import[f, "SDF"]

enter image description here

Note that the default display of a molecule is to completely ignore multiple bonds, although if 3-dimensional atom coordinates are provided, they are used. Here is the same file rendered using the MolViewer package:

data = (First@Import[f, {"SDF", #}] & /@ {"VertexTypes", "EdgeRules", 
    "EdgeTypes", "VertexCoordinates"});
data = data /. {x_StructuredArray :> 
 QuantityMagnitude@UnitConvert[x, "Picometers"]};
molview@makemolecule[Sequence @@ data]

enter image description here

The package includes a couple nice features to add highlights and labels to atoms, if you so desire.

molview@makemolecule[Sequence @@ data, labels -> Range@18, 
  hilight -> {1, 3, 5, 7, 9}, labeloffset -> {25, -50, 0}]

enter image description here

Turning our attention to the molecule of interest, I am using data from here as reference. It's behind a paywall, but you can get the relevant data from the crytallography open database. Let's assume a little magic happens and we are able to convert the CIF into a XYZ file (and the Mg from the paper mysteriously disappears).

v2o5a = ImportString["15
     9012219
     V       1.716593    1.719889    0.000000
     O       1.686578    0.142009   -0.003217
     O       3.685999    2.216045   -0.011437
     O       0.000000    2.198951   -0.011794
     V      -1.716593    1.719889    0.000000
     O      -1.686578    0.142009   -0.003217
     O      -3.685999    2.216045   -0.011437
     V      -4.055407    2.663111   -1.787000
     V      -4.055407    2.663111    1.787000
     V       4.055407    2.663111   -1.787000
     V       4.055407    2.663111    1.787000
     O      -2.086001    2.166955   -1.798437
     O      -2.086001    2.166955    1.775563
     O       2.086001    2.166955   -1.798437
     O       2.086001    2.166955    1.775563", {"XYZ", #}] & /@ \
{"VertexTypes", "VertexCoordinates"};

We can now view this file, but since the XYZ format does not have connectivity information, we will only print the atoms by passing empty lists to the bond-related arguments.

molview@makemolecule[v2o5a[[1]], {}, {}, 
  QuantityMagnitude@UnitConvert[v2o5a[[2]], "Picometers"]]

enter image description here

As with SDF files, Mathematica has a hard time with XYZ files resulting in some pretty interesting connectivity. It seems to do an OK job with this particular file; however caveat emptor.

enter image description here

The molviewer allows you to make connections, which might be useful (as it is in this case) to show the square pyramidal structure of the unit cell.

molview@makemolecule[
  v2o5a[[1]], {5 -> 6, 5 -> 4, 5 -> 13, 5 -> 12, 5 -> 7, 12 -> 4, 
   4 -> 13, 13 -> 7, 7 -> 12, 12 -> 6, 4 -> 6, 13 -> 6, 7 -> 6}, 
  ConstantArray["Single", 14], 
  QuantityMagnitude@UnitConvert[v2o5a[[2]], "Picometers"]]

enter image description here

I don't know if this makes the case for using Mathematica as a molecule viewer, but you can do cool things with it. For example, I know of no molecular viewing software that allows you to easily create a periodic table based on atomic size:

elements = 
  DeleteCases[Molviewer`Private`sizerules[][[All, 1]], 
   x_ /; ElementData[x, "Block"] == "f"];
positions = {100 ElementData[#, "Group"], 
     100 ElementData[#, "Period"], 0} & /@ elements;
molview@makemolecule[elements, {}, {}, positions]

enter image description here

Improve this answer
$\endgroup$
2
  • $\begingroup$ "atoms can be multiply bonded to one another." - and this is not even considering things like diborane. :) In any case, I'm aware of how complicated the following suggestion would be to implement, but it would be cool to have multiple bonds represented like in this depiction of ethylene. $\endgroup$
    – J. M.'s eventual burnout ♦
    May 6, 2016 at 18:17
  • $\begingroup$ @J.M. That would be interesting indeed. Most "proper" structure visualization software is able to always have the multiple bond "facing" the viewer, regardless of how the molecule is rotated. It is something I would like to incorporate but is beyond my programming skills. $\endgroup$
    – bobthechemist
    May 6, 2016 at 22:02
8
$\begingroup$

The problem with the GraphPlot is that you often get something that doesn't really resemble the molecule, it just has the right connectivity. First let's look at the example that is available in version 10,

ChemicalData["VanadiumVOxide"]

enter image description here

Obviously there is some data about atom positions, even if only in the 2D drawing. For the record, the 3D structure would be available via

ChemicalData["VanadiumVOxide", "AtomPositions"]
(* Missing["NotAvailable"] *)

So we can grab the 2D positions of the atoms in that layout, append a 0 for the z coordinate, and then make a molecule plot using functions we aren't supposed to have access to

ChemicalData["H2O", "MoleculePlot"];
data = ChemicalData["VanadiumVOxide", #] & /@ {"VertexTypes", 
    "EdgeRules", "VertexCoordinates"};
data[[3]] = Append[#, 0] & /@ data[[3]];
Graphics`MoleculePlotDump`iMoleculePlot3D[data, 
 ViewPoint -> {0, 0, 5}]

enter image description here

That first line, ChemicalData["H2O", "MoleculePlot"]; is required to initialize the iMoleculePlot3D function. If anyone can figure out what package to load with Needs I would love to hear it.

This works in version 9 as well,

ChemicalData["H2O", "MoleculePlot"];
data = ChemicalData["VanadiumIVOxide", #] & /@ {"VertexTypes", 
    "EdgeRules", "VertexCoordinates"};
data[[3]] = Append[#, 0] & /@ data[[3]];
Graphics`MoleculePlotDump`iMoleculePlot3D[data, 
 ViewPoint -> {0, 0, 5}]

enter image description here

Improve this answer
$\endgroup$
5
$\begingroup$

If you have OpenBabel installed, you can use its built-in structure optimization methods to generate 3D structures from SMILES structure format provided by Mathematica:

Import["!obabel -:\"" <> ChemicalData["Valeraldehyde", "SMILES"] <> 
  "\" -o xyz --gen3d --conformer --nconf 50 --score energy --weighted", "XYZ"]

optimized structure

Unfortunately it does not know how to optimize the vanadium oxide since it is not a molecule.

Improve this answer
$\endgroup$
3
  • $\begingroup$ This is what I was trying to tell Verbeia in chat; I feel this is probably best done by dedicated software instead of trying to do the entire business in Mathematica. $\endgroup$
    – J. M.'s eventual burnout ♦
    May 6, 2016 at 12:18
  • $\begingroup$ @J.M. True - if you were doing actual science. But this is the point that I confess that it was an idea for a friend's sci-fi TV script, where they needed to show the chemical onscreen and look pensive. :) $\endgroup$
    – Verbeia
    May 7, 2016 at 4:35
  • 1
    $\begingroup$ @Verbeia Please don't encourage "TV-science" :P Unless it's a parody $\endgroup$
    – shrx
    May 7, 2016 at 7:48
3
$\begingroup$

Three years after this question was posted, Wolfram introduced, in MMA 12 (2019), an experimental function (MoleculePlot3D) that widens its molecular drawing capabilities.

When it comes to molecular drawing in MMA, one can divide substances into two categories:

1) Those for which ChemicalData contains plotting information.

This is now the case for valderaldehde:

ChemicalData["Valeraldehyde", "MoleculePlot"]

enter image description here

ChemicalData["Valeraldehyde", "SpaceFillingMoleculePlot"]

enter image description here

This is also the case for the "SDF" test molecule in bobthechemist's answer, which is ethyl trans-2-butenoate. MMA doesn't recognize its IUPAC-like name:

ChemicalData["Ethyltrans2butenoate", "MoleculePlot"]

enter image description here

But does recognize its CAS number:

ChemicalData["CAS623-70-1", "MoleculePlot"]

enter image description here

2) Those for which ChemicalData does not contain plotting information

For these, one can import the SMILES information from elsewhere. Here it plots the structure for maitotoxin $(\text{C}_{164} \text{H}_{256} \text{O}_{68} \text{S}_2 \text{Na}_2)$. It doesn't know where to put the two sodium ions, since that information isn't in the SMILES string, but it does get the relative sizes right (radius sodium ion: radius oxygen atom = 102 pm : 60 pm = 1.7):

MoleculePlot3D[
 ImportString[
   "S(=O)(=O)(O[C@H]1[C@@H](O[C@@H]2C[C@@H]3O[C@@H]4[C@H]([C@H]5O[C@H]\
([C@@H]([C@H]([C@@H]5O[C@H]4C[C@H]3O[C@H]2[C@H]1O)O)O)[C@H]1[C@@H]([C@\
H]([C@H]2[C@@H](O1)C[C@H]([C@@H](O2)[C@@H](C[C@H](C[C@H]2[C@@H]([C@H](\
[C@H]1[C@@H](O2)C[C@H]([C@@H](O1)[C@H]1[C@@H](C[C@@]2(O[C@@H]3C[C@@]4(\
O[C@@H]5C[C@@]6(O[C@@H]7CC[C@H](O[C@H]7C[C@H]6O[C@]5(CC[C@H]4O[C@]3([\
C@H]([C@H]2O1)O)C)C)[C@@]1(O[C@@]2(C[C@@H]3O[C@@H]4C[C@@H]5O[C@@H]6C\\\
C=C/[C@@H]7O[C@@]8(C[C@@]9(O[C@@]%10(CC[C@H](O[C@H]%10C[C@H]9O[C@H]8C[\
C@H]7O[C@H]6C[C@]5(O[C@]4(CC[C@]3(O[C@H]2C[C@H]1O)C)C)C)[C@H]([C@@H](\
C[C@H]([C@H](CC=C)C)C)O)O)C)C)C)C)C)C)C)C)O)O)O)O)O)O)O)O)O)O)C[C@@H](\
[C@H]([C@@H]1[C@@H](C[C@H]2O[C@]3(C[C@H]4O[C@]5(C[C@H]([C@H]6O[C@H]([\
C@@H]([C@@H]([C@@H]6O[C@@H]5C[C@@H]4O[C@@H]3[C@H]([C@@H]2O1)O)O)O)[C@\
H](C)[C@@H]([C@@H](CC[C@@H]([C@@H]([C@@H](C[C@H](C(\\C(=C\\CO)\\C)=C)\
O)C)O)OS(=O)(=O)[O-])C)O)O)C)C)O)O)O)[O-].[Na+].[Na+]", 
   "SMILES"][[1]]]

enter image description here

MoleculePlot3D[
 ImportString[
   "S(=O)(=O)(O[C@H]1[C@@H](O[C@@H]2C[C@@H]3O[C@@H]4[C@H]([C@H]5O[C@H]\
([C@@H]([C@H]([C@@H]5O[C@H]4C[C@H]3O[C@H]2[C@H]1O)O)O)[C@H]1[C@@H]([C@\
H]([C@H]2[C@@H](O1)C[C@H]([C@@H](O2)[C@@H](C[C@H](C[C@H]2[C@@H]([C@H](\
[C@H]1[C@@H](O2)C[C@H]([C@@H](O1)[C@H]1[C@@H](C[C@@]2(O[C@@H]3C[C@@]4(\
O[C@@H]5C[C@@]6(O[C@@H]7CC[C@H](O[C@H]7C[C@H]6O[C@]5(CC[C@H]4O[C@]3([\
C@H]([C@H]2O1)O)C)C)[C@@]1(O[C@@]2(C[C@@H]3O[C@@H]4C[C@@H]5O[C@@H]6C\\\
C=C/[C@@H]7O[C@@]8(C[C@@]9(O[C@@]%10(CC[C@H](O[C@H]%10C[C@H]9O[C@H]8C[\
C@H]7O[C@H]6C[C@]5(O[C@]4(CC[C@]3(O[C@H]2C[C@H]1O)C)C)C)[C@H]([C@@H](\
C[C@H]([C@H](CC=C)C)C)O)O)C)C)C)C)C)C)C)C)O)O)O)O)O)O)O)O)O)O)C[C@@H](\
[C@H]([C@@H]1[C@@H](C[C@H]2O[C@]3(C[C@H]4O[C@]5(C[C@H]([C@H]6O[C@H]([\
C@@H]([C@@H]([C@@H]6O[C@@H]5C[C@@H]4O[C@@H]3[C@H]([C@@H]2O1)O)O)O)[C@\
H](C)[C@@H]([C@@H](CC[C@@H]([C@@H]([C@@H](C[C@H](C(\\C(=C\\CO)\\C)=C)\
O)C)O)OS(=O)(=O)[O-])C)O)O)C)C)O)O)O)[O-].[Na+].[Na+]", 
   "SMILES"][[1]], PlotTheme -> "SpaceFilling"]

enter image description here

Finally, a note on vanadium(IV) oxide and vanadium(V) oxide:

You can get MMA to produce drawings for these from its SMILES ChemicalData info., but the structures it outputs don't make sense, since the SMILES info. MMA has is for them is only that of the individual repeating units. This is not surprising, since SMILES is designed for individual molecules. As a consequence, MMA is treating the substances like molecular covalent compounds, when in fact they are ionic—i.e., they don't exist as individual molecules, but rather as repeating units within extended crystalline structures. Consider, for instance, vanadium(V) oxide:

MoleculePlot3D[Molecule[ChemicalData["VanadiumVOxide", "SMILES"]]]

enter image description here

Note also that the depicted double covalent bonds don't make sense for this compound. It's not surprising MMA does this, since

Instead, a better ball-and-stick representation would look something like this, except extended into three dimensions (this is for a monolayer). Source: https://en.wikipedia.org/wiki/Vanadium(V)_oxide:

enter image description here

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ MoleculePlot3D[Molecule["ethyl (2E)-but-2-enoate"]] should work fine. $\endgroup$
    – J. M.'s eventual burnout ♦
    May 23, 2020 at 17:28
2
$\begingroup$

Another possibility with GraphPlot3D (I also had to choose a different oxide because of missing data):

{edges, elements} = 
  ChemicalData["VanadiumVOxide", #] & /@ {"EdgeRules", "VertexTypes"};
coords = GraphPlot3D[edges][[1, 1, 1]];
colorRules = Thread[coords -> elements];
GraphPlot3D[edges, 
 VertexRenderingFunction -> ({ColorData["Atoms"][#1 /. colorRules], 
     Sphere[#1, .15]} &), 
 EdgeRenderingFunction -> (Cylinder[#1, .05] &), Boxed -> False]

colored atoms

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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