Symmetry unique atom coordinates

Question

One thing I love about Mathematica is how easily I can go from the name of a molecule to estimated coordinates of its atoms, with a command like

AtomList[Molecule[Entity["Chemical", "Toluene"]], All, {"AtomicNumber", "AtomCoordinates"}]

(although, oddly enough, "AtomCoordinates" does not appear in the "AtomList" documentation)

I can also easily get the point group:

Molecule[Entity["Chemical", "Toluene"]]["PointGroup"]

This is exciting because this is exactly the input I need to run GAMESS and do quantum chemistry calculations (starting with a geometry optimization, of course, since JM has informed me that these coordinates are heuristic guesses).

But, really, this is not exactly the input that I need: what I really need are coordinates of only the symmetry-unique atoms.

I don't suppose there's a way to get coordinates of symmetry-unique atoms, which I can use for GAMESS input? I know there's some functions related to point group symmetry, but I haven't thought of how to do it.

To clarify, I'm looking for an answer that matches the point group given by the PointGroup property. For example, the code above will give a point group of D3d for cyclohexane, which corresponds to the chair conformation. So there should be three equivalence classes: carbons, equatorial hydrogens, and axial hydrogens, since axial and equatorial can't be transformed into each other by the symmetry elements in D3d.

Answer 1

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.

Answer 2

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