Change bond cutoff distance when plotting .xyz files

Question

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.

Answer 1

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}

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]

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}]

Answer 2

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[%]

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.