2
$\begingroup$

When plotting a molecule from a xyz data file, Mathematica automatically adds the bonds between the atoms. But in my case it adds too many bonds, not only between nearest neighbours but also between atoms which are further away from each other.

I assume that Mathematica uses a cutoff distance to determine whether to draw a bond or not.

How can I change this cutoff distance?

Here is the content of the .xyz file:

11
MC cluster
O       -0.186388 0.742079 0.766868     
O       0.805206 0.229831 -0.506672     
O       0.737583 0.398046 0.505395  
O       0.620872 -0.517223 0.12745  
O       0.0221594 -0.21901 0.858977     
O       0.128471 -0.483467 -0.732157    
O       -0.796687 -0.122409 -0.555157   
O       -0.0442577 0.206399 -0.0372612  
O       -0.353133 -0.731401 0.0979684   
O       -0.864308 0.0458028 0.456909    
O       -0.0695165 0.451352 -0.98232 

When I increase all coordinates (and thereby the distances) by a factor of 10, Mathematica draws no bonds at all.

$\endgroup$
  • $\begingroup$ Can you give an example *.xyz file where this happens? $\endgroup$ – J. M. will be back soon Mar 6 '16 at 13:29
  • 1
    $\begingroup$ You will find the answer here: mathematica.stackexchange.com/q/24045/12 The options I show there can be directly used in Import, e.g. Import["a.xyz","XYZ", "InferBondsMinDistance" -> 20000] $\endgroup$ – Szabolcs Mar 6 '16 at 15:09
2
$\begingroup$

It does seem possible to change a bond cutoff length using built-in functions, as mentioned by Szabolcs in comments, via something like

Import["ExampleData/caffeine.xyz", "XYZ", 
   "InferBondsMinDistance" -> #] & /@ {10000, 15000, 20000}

Mathematica graphics

But I have no idea what units those are supposed to be. They aren't ångströms (which the "XYZ" files are usually written in), nor are they picometers (which the coordinates are usually converted to before plotting). You can investigate the function Graphics`MoleculePlotDump`InferBonds to learn more, but it doesn't seem all that configurable from the outside.

So why not rebuild the XYZ-plotting function from the ground up to have complete control over the minimum and maximum bond-length cutoffs?

variableBondPlot[xyzString_, maxBondLength_, minBondLength_: 0, 
  atomRadius_: .24, bondRadius_: .08] := 
 Module[{data, colors, coords, spheres, bonds},
  data = ImportString[xyzString, "Table"][[3 ;;]];
  {colors, coords} = {ColorData["Atoms", First@#], Rest@#} & /@ data //
     Transpose;
  spheres = {colors[[#]], Sphere[#, atomRadius]} & /@ 
       Range[Length@data] //
      GatherBy[#, First] & // Flatten // 
    DeleteDuplicates;
  bonds = (UpperTriangularize@DistanceMatrix[Rest /@ data]) // 
     MapIndexed[
       If[minBondLength < # <= maxBondLength, 
         Cylinder[#2]] &, #, {2}] & // Cases[#, _Cylinder, Infinity] &;
  bonds = Module[{i, j, m},
          {i, j} = First@#; m = Length@coords + 1; 
          AppendTo[coords, Mean@coords[[{i, j}]]];
          {{colors[[i]], Cylinder[{i, m}, bondRadius]}, {colors[[j]], 
            Cylinder[{m, j}, bondRadius]}}
          ] & /@ bonds // Flatten[#, 1] & //

      GatherBy[#, First] & // Flatten // DeleteDuplicates;
  Graphics3D[{Specularity[GrayLevel[1], 100], EdgeForm[None], 
    AbsoluteThickness[3], GraphicsComplex[coords, {spheres, bonds}]}, 
   Boxed -> False, Lighting -> "Neutral"]
  ]

The arguments should be self-explanatory. Here it is applied to the OP's data:

variableBondPlot["11
 MC cluster
 O       -0.186388 0.742079 0.766868     
 O       0.805206 0.229831 -0.506672     
 O       0.737583 0.398046 0.505395  
 O       0.620872 -0.517223 0.12745  
 O       0.0221594 -0.21901 0.858977     
 O       0.128471 -0.483467 -0.732157    
 O       -0.796687 -0.122409 -0.555157   
 O       -0.0442577 0.206399 -0.0372612  
 O       -0.353133 -0.731401 0.0979684   
 O       -0.864308 0.0458028 0.456909    
 O       -0.0695165 0.451352 -0.98232",
 .98]

Mathematica graphics

And you can get creative with the bond cutoffs as well:

Manipulate[
 variableBondPlot[Import["ExampleData/caffeine.xyz", "Text"], maxbond,
   minbond],
 {{maxbond, 1.5}, 1.0, 8, .01},
 {{minbond, 0}, 0, maxbond, .01}]

Mathematica graphics

$\endgroup$
2
$\begingroup$

According to Wikipedia and other sources, a bond is usually drawn between two atoms if their distance is less than the sum of their covalent radii. This is what Mathematica does, too, except that it won't draw bonds for molecules that are too close together.

The internally stored radii used for conversion, Graphics`MoleculePlotDump`elementCovalentRadii, do not agree with the radii returned by ElementData[]. I'm not a chemist, so I'm not sure what's up there. The entry for oxygen, elementCovalentRadii[[8]], agrees with Google, namely 73 pm.

So the problem must be in the units of the coordinates of the XYZ file. They should be in Angströms in the XYZ file; Mathematica converts them to picometers for the "VertexCoordinates" of the graphics. They seem only slightly too close. Scaling by 1.06 removes some of the bonds:

rules = ImportString[
   "11
   MC cluster
   O       -0.186388 0.742079 0.766868     
   O       0.805206 0.229831 -0.506672     
   O       0.737583 0.398046 0.505395  
   O       0.620872 -0.517223 0.12745  
   O       0.0221594 -0.21901 0.858977     
   O       0.128471 -0.483467 -0.732157    
   O       -0.796687 -0.122409 -0.555157   
   O       -0.0442577 0.206399 -0.0372612  
   O       -0.353133 -0.731401 0.0979684   
   O       -0.864308 0.0458028 0.456909    
   O       -0.0695165 0.451352 -0.98232 ",
   {"XYZ", "Rules"}]
(*  {"Graphics3D" -> ..., "VertexCoordinates" -> {...}, "VertexTypes" -> {...}}  *)

Export["/tmp/foo.xyz",
  MapAt[1.06 # &, Rest@rules, {1, 2}],  (* scale coords at position {1, 2} of Rest@rules *)
  "Rules"];
Import[%]

Mathematica graphics

At around a scaling factor of 1.65, more bonds begin to disappear. At a factor of 1.77, all bonds are gone. I'll leave the chemical problem of what is the most appropriate solution to the reader.

$\endgroup$
  • 2
    $\begingroup$ You posted an answer here some time ago: mathematica.stackexchange.com/q/24045/12 (I always forget about my own answers) $\endgroup$ – Szabolcs Mar 6 '16 at 15:10
  • $\begingroup$ @Szabolcs Thanks, I totally forgot. I hardly ever use these functions. My fading understanding of first-year chemistry makes me wonder if there's something wrong with the OP's example data, though. Or just what a bond means if the atoms are so close each could be bonded to the others. Perhaps some really interesting physics going on that I don't get at all. $\endgroup$ – Michael E2 Mar 6 '16 at 15:23
  • $\begingroup$ The example data is an atomic cluster. Those are not really oxygen atoms, but simply arbitrary rare gase atoms, which interact via a lennard-jones potential. (I had to choose an atom type because otherwise Mathematica wouldn't plot the data.) When the atoms are cooled down, they form cluster structures. $\endgroup$ – Blub Bla Mar 6 '16 at 15:43
  • $\begingroup$ @BlubBla I wondered if it were something like that, but I know only a little of such things. So do you have particular atoms you wish to join? Look at the answers to the question Szabolcs' linked. They might help. $\endgroup$ – Michael E2 Mar 6 '16 at 15:48
  • 2
    $\begingroup$ @BlubBla Sounds like you're better off importing the coordinates only, and doing the visualization yourself. $\endgroup$ – Szabolcs Mar 6 '16 at 16:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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