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.

Since I am visualizing my work using Mathematica, I import some 3D models of molecules to my notebook. So, the question is, is there any way to modify somehow that imported models inside Mathematica to, for example, label bonds' length? Currently I use .mol2 format, but maybe some others formats provide that ability? Or can Mathematica process that models to Graphics3D?

Thanks.

share|improve this question
add comment

2 Answers 2

up vote 4 down vote accepted

First we import the model:

mol = First@Import["ExampleData/abacavir.mol2"]

It will look like this:

mol1

After that, we can confirm that Mathematica has converted the model into a Graphics3D object by typing mol//InputForm. This answers one of your questions.

As an example of how to edit this model in Mathematica, I will create the object you mentioned, one where the bond lengths are shown.

First we need to know where to put the labels so in order to do that we calculate the midpoints of each bond. To do this we will use information provided by the MOL2 file, EdgeRules and VertexCoordinates. The first tells us which atoms link to which and the latter tells us the position of all the atoms. The positions are given in picometer, and the relative distances between atoms will be the bond length in that unit.

edgeRules = 
  First@Import["ExampleData/abacavir.mol2", "EdgeRules"] /. 
   Rule -> List;
vertexCoordinates = 
  MapIndexed[First@#2 -> # &, 
   First@Import["ExampleData/abacavir.mol2", "VertexCoordinates"]];
paired = edgeRules /. vertexCoordinates;

The above sets up a list paired in which each element is a list of two coordinates, these coordinates correspond to atoms that are bound to each other. This visualization confirms that we can now find the midpoints:

midpoints = #[[2]] - (#[[2]] - #[[1]])/2 & /@ paired;
Graphics3D[{Join[{Green, Sphere[midpoints, 20]}, List @@ mol]}]

mol2

The final step is to write the bond lengths on the corresponding midpoints, which is easy thanks to the paired list.

lengths = Norm[#[[2]] - #[[1]]] & /@ paired;
labels = Text[#[[2]], #[[1]], Background -> Black] & /@ 
   Transpose@{midpoints, lengths};
Graphics3D[Join[{Green, labels}, List @@ mol]]

enter image description here

share|improve this answer
add comment

Anon does a fine job in answering the question of how to annotate a molecular model in Mathematica but fails to address the 800 pound gorilla in the room, which is that Mathematica doesn't import multiple bonds. Taking a look at the wiki page for Abacavir, we see that there are a number of double bonds that don't show up when using Import.

abacavir

Fortunately, the information is there, but we essentially have to build our molecule ourselves.

First@Import["ExampleData/abacavir.mol2", "EdgeTypes"]
(* {"Aromatic", "Aromatic", "Aromatic", "Aromatic", "Aromatic", "Single", ... *)

I have been working on a molecular viewer (it's still in progress) but it has some of the features that the OP is looking for. I apologize in advance for long code (and I'll be happy to remove this post if it is inappropriate to provide this much detail). First, I import the molecule information:

{atomtypes, atomcoords, bonds, bondtypes} = First@Import["ExampleData/abacavir.mol2",
    #] & /@ {"VertexTypes", "VertexCoordinates", "EdgeRules", "EdgeTypes"};
bondtypes[[{2,3,4,5,21,22}]]="Single";

The 2nd command is to address the fact that I haven't yet coded aromatic bonds (they by default show up as double bonds instead of a bond/dashed bond).

The code to build a molecule, including options for highlighting and labeling:

Options[makemolecule] = {hilight -> {}, labels -> {}, 
   labeloffset -> {1, 1, 1}};
makemolecule::badcoord = "2D coordinates suspected, use \
PadRight[#,3]&/@coords"; 
makemolecule[atoms_, coords_, bonds_, bondtypes_, OptionsPattern[]] :=
  Module[{colors, sizes, hb, b1, b2},
  colors = {
    "H" -> RGBColor[0.65, 0.7, 0.7],
    "C" -> RGBColor[0.4, 0.4, 0.4],
    "N" -> RGBColor[0.29, 0.44, 0.89],
    "O" -> RGBColor[0.8, 0.2, 0.2],
    "F" -> RGBColor[0.58, 0.86, 0.41]
    };
  sizes = {
    "C" -> 34, "H" -> 24, "O" -> 31, "N" -> 32, "F" -> 30};
  hb = Mean@coords[[#]] & /@ ({#[[1]], #[[2]]} & /@ bonds);
  b1 = coords[[bonds[[All, 1]]]];
  b2 = coords[[bonds[[All, 2]]]];
  If[(Dimensions@coords)[[2]] != 3,
   Message[makemolecule::badcoord]
   ];
  GraphicsGroup[{EdgeForm[None],
    (* Make the atoms *)
    GraphicsComplex[coords, {
      MapThread[{#1 /. colors, Sphere[#2, #1 /. sizes]} &, {atoms, 
        Range[Length[atoms]]}]}],
    (* Make the bonds *)
    GraphicsComplex[Join[b1, hb], {
      Riffle[atoms[[bonds[[All, 1]]]] /. colors,
       Table[GeometricTransformation[
           Cylinder[{#, Length@b1 + #}, 10],
           TranslationTransform[ k Normalize[hb[[#]]\[Cross]b1[[#]]]]
           ], {k, 
           bondtypes[[#]] /. {"Single" -> {0}, "Double" -> {-15, 15}, 
             "Aromatic" -> {-15, 15}, "Triple" -> {-15, 0, 15}}}] & /@
         Range[Length@b1]
       ]
      }],
    GraphicsComplex[Join[b2, hb], {
      Riffle[atoms[[bonds[[All, 2]]]] /. colors,
       Table[GeometricTransformation[
           Cylinder[{#, Length@b2 + #}, 10],
           TranslationTransform[k Normalize[hb[[#]]\[Cross]b2[[#]]]]
           ], {k, 
           bondtypes[[#]] /. {"Single" -> {0}, "Double" -> {-15, 15}, 
             "Aromatic" -> {-15, 15}, "Triple" -> {-15, 0, 15}}}] & /@ 
        Range[Length@b2]
       ]
      }],

    (* Hilight atoms *)
    If[Length@OptionValue[hilight] > 0,
     GraphicsComplex[
      coords, {Opacity[0.6], Green, 
       Sphere[#, 35] & /@ OptionValue[hilight]}]],
    (* Add labels *)
    If[Length@OptionValue[labels] == Length@atoms,
     GraphicsComplex[
      coords, {Black, 
       MapThread[
        Text[#1, OptionValue[labeloffset] + #2, 
          BaseStyle -> {FontSize -> 12, FontFamily -> "Helvetica", 
            Orange}] &, {OptionValue[labels], coords}]}]]

    }
   ]];

(Yes, only 5 elements in my periodic table at the moment...) I can then view a molecule using:

molview[molecule_] := Graphics3D[{
   Specularity[GrayLevel[1], 100], molecule
   },
  Boxed -> False,
  Lighting -> "Neutral",
  SphericalRegion -> True, Background -> Black]

Putting all that together, I get:

molview@makemolecule[atomtypes, atomcoords, bonds, bondtypes]

abacavir with dbl bonds

And the makemolecule function can accept some rudimentary highlighting and labeling options

molview@makemolecule[atomtypes, atomcoords, bonds, bondtypes, 
    hilight -> {1, 3, 5, 7,     9}, labels -> atomtypes, labeloffset -> {50, 50, 0}]

enter image description here

Now I have to take a closer look at anon's answer to see if there's anything I can incorporate into my function. Hope this is useful.

share|improve this answer
add comment

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.