How can the units of the "AtomCoordinates" in a Molecule be changed while keeping the resulting MoleculePlot3D consistent?

Question

Mathematica provides the Molecule function that can be visualized with MoleculePlot3D since version 12.0. Before, importing "XYZ" files would produce a Graphics3D object immediately. In this object, the coordinates of the atom were given in picometers, but not as Quantity. Changing the units here necessitated also fiddling with the Spheres and Cylinders used for the plot, which made it rather unfeasible.

Now, a Molecule has "AtomCoordinates" which are given as QuantityArray with the unit "Angstroms". Since ModifyMolecule with "TransformAtomCoordinates" allows applying any function receiving and returning a QuantityArray, I was hopeful that a UnitConvert would work. Sadly, that seems not to be the case. Firstly,

MoleculeModify[
  Molecule["Benzene"], {"TransformAtomCoordinates", 
   UnitConvert[#, "Picometers"] &}]["AtomCoordinates"]
Out[]="QuantityArray[Dimensions:{12,3},Units:Angstroms]

i.e. the units seem unchanged. Plotting the molecule as

Show[MoleculePlot3D[
  MoleculeModify[
   Molecule["Benzene"], {"TransformAtomCoordinates", 
    UnitConvert[#, "Picometers"] &}], PlotTheme -> "BallAndStick"], 
 Axes -> True]

where the correct plot

becomes

reveals that the coordinates are indeed changed in some way, but the theme cannot keep up. I guess the data gets changed, but the unit information itself gets overwritten again after the transformation?

Can anything be done in this case? Different from the old XYZ plot, it seems even harder to access the internals of a Molecule to change how it is plotted, but maybe there is an option that is either undocumented or which I missed.

Answer 1

There are two issues here - first a bug that was introduced in version 13.1 and should be fixed in a future release, and secondly some confusion about the design of Molecule.

By design Molecule will always convert any coordinates input into a QuantityArray with "Angstroms" unit. The following was the behavior from version 12.0 to 13.0:

In[1]:= m1 = 
  MoleculeModify[Molecule["benzene"], "ComputeAtomCoordinates"];
m2 = MoleculeModify[m1, {"TransformAtomCoordinates", UnitConvert[#, "Picometers"] &}];

In[3]:= m1["AtomCoordinates"] == m2["AtomCoordinates"]

Out[3]= True

So the MoleculeModify operation above should be a no-op.

Starting in version 13.1 the input coordinates are no longer converted correctly (only the unit is changed), this is a bug.

The workaround is to make sure that the input AtomCoordinates are always in Angstroms.

If you want to plot the molecule in 3D using picometers I would do something like this

ReplaceAll[
MoleculePlot3D[Molecule["adenosine"], Axes -> True],
 gc_GraphicsComplex :> 
  GeometricTransformation[gc, ScalingTransform[{100, 100, 100}]]
]

Adding an option to MoleculePlot3D to use a different coordinate for the axes is a good suggestion.