Symmetry unique atom coordinates

Jason B saw this before I did, but I think the following reproduces some of the example GAMES input reasonably well.

f[mol_] := Module[{al, out},
  al = AtomList[Molecule[mol], 
    All, {"AtomicNumber", "AtomCoordinates"}];
  out = QuantityMagnitude /@ 
    Flatten /@ al[[First /@ Molecule[mol]["SymmetryEquivalentAtoms"]]];
  out /. z_Integer :> Sequence[ElementData[z, "Abbreviation"], z]
  ]

f["Toluene"]

The general plan for solving this is

1. Generate all symmetry transformations for a given molecule.
2. Apply these transformations to each atom coordinate, giving a list of coordinates for each atom.
3. Group atoms which give the equivalent lists of coordinates are considered equivalent.

Sadly the Wolfram developers don't give the actual transformation functions associated with a given symmetry element via any built-in function. But they do give us enough information in the "SymmetryElements" property to construct these ourselves:

In[26]:= Molecule["methane"]["SymmetryElements"] // pf2

Out[26]= {
    <|
        "Operation" -> "Rotation", "Name" -> Subscript["C", "3"],
        "Degree" -> 3, "UniqueOperationsCount" -> 2, 
  "RotationAxis" -> InfiniteLine[
                {0., 0., 0.},
                {0.9312106494091753, 0.3062515387515941, 0.19762773448891885}
            ]
    |>,
    ........,
    <|
        "Operation" -> "Reflection", "Name" -> "\[Sigma]", 
  "Degree" -> 1, "UniqueOperationsCount" -> 1,
        "SymmetryPlane" -> Hyperplane[
                {-0.6671653488434035, -0.16935533665066543, -0.7254027620919287},
                {0., 0., 0.}
            ]
    |>
 }

By examining the structure of that output, we can write a function to return the transformation from the symmetry element. I like to use KeyValuePattern for easily readable defitions:

symmetryOperation[KeyValuePattern[{"Operation"->"Rotation","Degree"->d_,"RotationAxis"->InfiniteLine[point_,direction_]}]] := RotationTransform[(2 * Pi) / d, direction, point];
symmetryOperation[KeyValuePattern[{"Operation"->"ImproperRotation","Degree"->d_,"RotationAxis"->InfiniteLine[point_,direction_]}]] := ReflectionTransform[direction, point] @* RotationTransform[(2 * Pi) / d, direction, point];
reflectpoint[point_, center_] := point + 2 * (center + -point);
symmetryOperation[KeyValuePattern[{"Operation"->"Inversion","InversionCenter"->Point[center_]}]] := Composition[
    ReflectionTransform[{1, 0, 0}, center], 
    ReflectionTransform[{0, 1, 0}, center], 
    ReflectionTransform[{0, 0, 1}, center]
];
symmetryOperation[KeyValuePattern[{"Operation"->"Reflection","SymmetryPlane"->Hyperplane[normal_,point_]}]] := ReflectionTransform[normal, point]

Now we take write a function to return all symmetry transformations for a molecule, correcting the oversight that Wolfram has made by not including the Identity element:

symmetryTransforms[mol_] := Join[{Identity}, Map[symmetryOperation, mol @ "SymmetryElements"]];

Now wrap it all together with a function to apply each transformation to each atom coordinate, and then gather those that produce the same coordinate list:

symmetryUniqueAtomIndices[mol_, tolerance_:0.1] := Module[
    {
        transforms = symmetryTransforms @ mol,
        points = QuantityMagnitude @ mol @ "AtomCoordinates"
    },
    PrependTo[transforms, Identity];
    GatherBy[Range @ Length @ points,
        Sort[
            DeleteDuplicates[Round[Through[transforms[Part[points, #]]], tolerance]]
        ]&
    ]
]

This uses GatherBy to group equivalent atoms. The important part here is to make a function to canonicalize the transformed coordinates, and for that I'm just rounding the numeric values, deleting duplicates, and then sorting them. There is probably room for improvement in this step.

You can look at the different cyclohexane conformations from this example:

labels = {"planar", "chair", "twist-boat", "boat", "half-boat", "half-chair"};
conformers = AssociationThread[
    labels -> CloudImport[
        CloudObject["https://www.wolframcloud.com/objects/555b1b48-9f89-45ef-a9e2-49c8fe5228b6"],
        "SDF"
    ]
];

Compare the symmetry of the different conformations:

In[10]:= symmetryUniqueAtomIndices /@ conformers

Out[10]= <|"planar" -> {{1, 2, 3, 4, 5, 6}, {7, 8, 9, 10, 11, 12, 13, 
    14, 15, 16, 17, 18}}, 
 "chair" -> {{1, 2, 3, 4, 5, 6}, {7, 9, 12, 13, 16, 18}, {8, 10, 11, 
    14, 15, 17}}, 
 "twist-boat" -> {{1, 4}, {2, 3, 5, 6}, {7, 8, 13, 14}, {9, 12, 15, 
    17}, {10, 11, 16, 18}}, 
 "boat" -> {{1, 4}, {2, 3, 5, 6}, {7, 14}, {8, 13}, {9, 11, 15, 
    18}, {10, 12, 16, 17}}, 
 "half-boat" -> {{1}, {2, 6}, {3, 5}, {4}, {7}, {8}, {9, 18}, {10, 
    17}, {11, 15}, {12, 16}, {13}, {14}}, 
 "half-chair" -> {{1, 4}, {2, 3}, {5, 6}, {7, 13}, {8, 14}, {9, 
    12}, {10, 11}, {15, 17}, {16, 18}}|>

If you want only one atom from each equivalence group, use something like

In[11]:= Map[First] /@ %

Out[11]= <|"planar" -> {1, 7}, "chair" -> {1, 7, 8}, 
 "twist-boat" -> {1, 2, 7, 9, 10}, "boat" -> {1, 2, 7, 8, 9, 10}, 
 "half-boat" -> {1, 2, 3, 4, 7, 8, 9, 10, 11, 12, 13, 14}, 
 "half-chair" -> {1, 2, 5, 7, 8, 9, 10, 15, 16}|>

You can visualize the symmetry groups via something like

MoleculePlot3D[conformers["chair"], 
 symmetryUniqueAtomIndices@conformers["chair"]]

In this image, all atoms of a given color are equivalent under the symmetry operations available. You can see that the hydrogen atoms now fall into two categories, the equatorial (radiating 'out' from the ring) in purple and the axial (with bonds parallel to the main symmetry axis) in blue.