# Symmetry unique atom coordinates

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.

• By symmetry-unique, would you take 3D geometry into account? Or would you consider all the hydrogen atoms in cyclohexane to be equivalent? The molecule property "SymmetryEquivalentAtoms" is in the latter category - it only considers topology, and will give a nested list of atom indices. – Jason B. Sep 19 '20 at 1:36
• Take geometry into account. First of all because I assume that would be necessary for GAMESS input. And second of all because when I say "symmetry unique" here, I mean equivalent under the point group given by "PointGroup" property. So, it seems that "SymmetryEquivalentAtoms" is not the answer. By the way it seems misleadingly named, since its natural to assume the "Symmetry" referred to is the one given by the "PointGroup" property. The axial and equatorial hydrogens aren't equivalent in the D3d symmetry are they? That's what "PointGroup" gives for cyclohexane. – user54038 Sep 19 '20 at 19:39
• I just edited the question to clarify this. – user54038 Sep 19 '20 at 20:07
• I understand your point about the property being misnamed. I was originally calling it "TopologicallyEquivalentAtom" but we changed it to the current name at some point. I think it should be renamed, and the property you describe here could be called "SymmetryUniqueAtoms". – Jason B. Sep 19 '20 at 22:59
• Yeah, although getting equivalence classes is nice. I'll also note that SymmetryEquivalentAtoms is listed under GeometricProperties, not GraphProperties. I don't meant to complain--I love the new molecule functionality and as you can see I've been using it a lot. – user54038 Sep 20 '20 at 1:43

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:= Molecule["methane"]["SymmetryElements"] // pf2

Out= {
<|
"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"};
labels -> CloudImport[
CloudObject["https://www.wolframcloud.com/objects/555b1b48-9f89-45ef-a9e2-49c8fe5228b6"],
"SDF"
]
];


Compare the symmetry of the different conformations:

In:= symmetryUniqueAtomIndices /@ conformers

Out= <|"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:= Map[First] /@ %

Out= <|"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.

• This seems like it should work if Molecule[x]["SymmetryElements"] is closed under composition. Do you know if it is? I'm worried about a pair of atoms that are equivalent but it takes two symmetry operations to turn one into the other. In that case I don't think symmetryUniqueAtomIndices would work, unless those two symmetry operations were also returned as a single, composed symmetry operation. I haven't been able to test this because I actually haven't been able to think of a molecule where it's relevant. – user54038 Oct 22 '20 at 13:22
• Actually, this answer told me everything I needed to know, which is how to actually use the symmetry elements to check equivalence of atoms. I'm still not 100% sure it's correct in every case, but I can take it from here. Accepting. – user54038 Oct 29 '20 at 2:21

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"] • This works for many of the GAMESS document example molecules so long as they are recognized by Molecule. Molecule["Methylnitrene"], for example, does not work. – bobthechemist Sep 19 '20 at 2:02
• In case you don't get a notification, I started a comment here and it got too long, I put it in chat. – Jason B. Sep 19 '20 at 2:16
• Comments above based on the unedited answer. See this chat for reference. – bobthechemist Sep 19 '20 at 11:18