Coordenadas de átomos únicos de simetría

Question

Una cosa que me encanta de Mathematica es la facilidad con la que puedo pasar del nombre de una molécula a las coordenadas estimadas de sus átomos, con un comando como

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

(aunque, curiosamente, "AtomCoordinates" no aparece en la documentación de "AtomList")

También puedo obtener fácilmente el grupo de puntos:

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

Esto es emocionante porque esta es exactamente la entrada que necesito para ejecutar GAMESS y hacer cálculos de química cuántica (comenzando con una optimización de geometría, por supuesto, ya que JM me ha informado que estas coordenadas son conjeturas heurísticas ).

Pero, en realidad, esta no es exactamente la entrada que necesito: lo que realmente necesito son coordenadas de solo los átomos de simetría única.

Supongo que no hay una manera de obtener coordenadas de átomos de simetría única, que puedo usar para la entrada de GAMESS. Sé que hay algunas funciones relacionadas con la simetría de grupos de puntos , pero no he pensado en cómo hacerlo.

Para aclarar, estoy buscando una respuesta que coincida con el grupo de puntos dado por la propiedad PointGroup. Por ejemplo, el código anterior dará un grupo de puntos de D3d para ciclohexano, que corresponde a la conformación de la silla. Por lo tanto, debería haber tres clases de equivalencia: carbonos, hidrógenos ecuatoriales e hidrógenos axiales, ya que los elementos de simetría en D3d no pueden transformarse entre sí y axiales.

Answer 1

El plan general para resolver esto es

1. Genere todas las transformaciones de simetría para una molécula determinada.
2. Aplique estas transformaciones a las coordenadas de cada átomo, dando una lista de coordenadas para cada átomo.
3. Los átomos de grupo que dan listas equivalentes de coordenadas se consideran equivalentes.

Lamentablemente, los desarrolladores de Wolfram no brindan las funciones de transformación reales asociadas con un elemento de simetría dado a través de ninguna función incorporada. Pero nos dan suficiente información en la "SymmetryElements"propiedad para construirlos nosotros mismos:

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

Al examinar la estructura de esa salida, podemos escribir una función para devolver la transformación del elemento de simetría. Me gusta usar KeyValuePatternpara definiciones fáciles de leer:

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]

Ahora tomamos escribir una función para devolver todas las transformaciones de simetría para una molécula, corrigiendo el descuido que Wolfram ha hecho al no incluir el elemento Identity:

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

Ahora envuélvalo todo junto con una función para aplicar cada transformación a cada coordenada de átomo y luego reúna las que producen la misma lista de coordenadas:

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

Esto se usa GatherBypara agrupar átomos equivalentes. La parte importante aquí es hacer una función para canonicalizar las coordenadas transformadas, y para eso solo estoy redondeando los valores numéricos, eliminando duplicados y luego ordenándolos. Probablemente haya margen de mejora en este paso.

Puede ver las diferentes conformaciones de ciclohexano de este ejemplo :

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 la simetría de las diferentes conformaciones:

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

Si solo desea un átomo de cada grupo de equivalencia, use algo como

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

Puede visualizar los grupos de simetría a través de algo como

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

En esta imagen, todos los átomos de un color dado son equivalentes bajo las operaciones de simetría disponibles. Puede ver que los átomos de hidrógeno ahora se dividen en dos categorías, el ecuatorial (que irradia 'hacia fuera' del anillo) en púrpura y el axial (con enlaces paralelos al eje de simetría principal) en azul.

Answer 2

Jason B vio esto antes que yo, pero creo que lo siguiente reproduce algunos de los ejemplos de entrada de JUEGOS razonablemente bien.

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