Merge pull request #68 from molssi-seamm/dev

Support for InChI improved, RMSD and PubChem added...

HISTORY.rst

@@ -1,6 +1,17 @@
 =======
 History
 =======
+2023.10.30 -- Support for InChI improved, RMSD and PubChem added...
+    * Adds support for aligning structures and calculating RMSD
+    * Adds support for working directly with PubChem to get structures, IUPAC names,
+      etc.
+    * Improves support for InChI, working around issues in both OpenBabel and RDKit.
+    * Added substantial new functionality for spacegroups and primitive cell handling,
+      but still not complete.
+
+2023.9.20 -- Better support for primitive cells and spacegroups
+    * Added getting the spacegroup from the symmetry operators
+    * Fixed updating the coordinates from the primitive cell
 
 2023.9.5 -- Support for velocities of atoms.
 

devtools/conda-envs/test_env.yaml

@@ -16,6 +16,7 @@ dependencies:
   - mendeleev
   - openbabel
   - pathvalidate
+  - pubchempy
   - pycifrw
   - rdkit
   - seekpath

molsystem/__init__.py

@@ -9,8 +9,8 @@
 # the molsystem package.
 
 import molsystem.elements  # noqa: F401
-from molsystem.system_db import SystemDB  # noqa: F401
-from molsystem.cell import Cell  # noqa: F401
+from .system_db import SystemDB  # noqa: F401
+from .cell import Cell  # noqa: F401
 from .properties import standard_properties, add_properties_from_file  # noqa: F401
 
 # Handle versioneer

molsystem/align.py

@@ -0,0 +1,137 @@
+# -*- coding: utf-8 -*-
+
+"""Alignment methods for configurations"""
+
+import logging
+
+try:
+    import openbabel  # noqa: F401
+    import openbabel.openbabel as ob
+except ModuleNotFoundError:
+    print(
+        "Please install openbabel using conda:\n"
+        "     conda install -c conda-forge openbabel"
+    )
+    raise
+
+logger = logging.getLogger(__name__)
+
+
+class AlignMixin:
+    """A mixin for handling alignment of configuration."""
+
+    def RMSD(self, other, include_h=True, symmetry=False):
+        """Compute the RMSD between configurations.
+
+        Parameters
+        ----------
+        other : Configuration or iterable of Configurations
+            A single configuration or e.g. list of configurations to match
+
+        include_h : bool = True
+            Whether to include hydrogen atoms in the RMSD
+
+        symmetry : bool = False
+            Whether to detect symmetric flips. Note this requires a lot of memory for
+            larger systems.
+
+        Returns
+        -------
+        float or [float]
+            The RMS between the current configuration and the target(s).
+        """
+        from .configuration import _Configuration
+
+        align = ob.OBAlign(include_h, symmetry)
+        align.SetRefMol(self.to_OBMol())
+
+        if isinstance(other, _Configuration):
+            align.SetTargetMol(other.to_OBMol())
+            align.Align()
+            return align.GetRMSD()
+        elif isinstance(other, ob.OBMol):
+            align.SetTargetMol(other)
+            align.Align()
+            return align.GetRMSD()
+
+        result = []
+        try:
+            for target in other:
+                if isinstance(target, _Configuration):
+                    align.SetTargetMol(target.to_OBMol())
+                    align.Align()
+                    result.append(align.GetRMSD())
+                elif isinstance(target, ob.OBMol):
+                    align.SetTargetMol(target)
+                    align.Align()
+                    result.append(align.GetRMSD())
+                else:
+                    raise TypeError(f"RMSD can't handle {type(target)}.")
+        except TypeError:
+            raise
+        except Exception as e:
+            logger.error(f"Configuration.RMSD {e}")
+            raise
+        return result
+
+    def align(self, other, include_h=True, symmetry=False):
+        """Align configurations.
+
+        Parameters
+        ----------
+        other : Configuration or iterable of Configurations
+            A single configuration or e.g. list of configurations to match
+
+        include_h : bool = True
+            Whether to include hydrogen atoms in the RMSD
+
+        symmetry : bool = False
+            Whether to detect symmetric flips. Note this requires a lot of memory for
+            larger systems.
+
+        Returns
+        -------
+        float or [float]
+            The RMS between the current configuration and the target(s).
+        """
+        from .configuration import _Configuration
+
+        align = ob.OBAlign(include_h, symmetry)
+        align.SetRefMol(self.to_OBMol())
+
+        if isinstance(other, _Configuration):
+            mol = other.to_OBMol()
+            align.SetTargetMol(mol)
+            align.Align()
+            align.UpdateCoords(mol)
+            other.coordinates_from_OBMol(mol)
+            return align.GetRMSD()
+        elif isinstance(other, ob.OBMol):
+            align.SetTargetMol(other)
+            align.Align()
+            align.UpdateCoords(other)
+            return align.GetRMSD()
+
+        result = []
+        try:
+            for target in other:
+                if isinstance(target, _Configuration):
+                    mol = target.to_OBMol()
+                    align.SetTargetMol(mol)
+                    align.Align()
+                    align.UpdateCoords(mol)
+                    target.coordinates_from_OBMol(mol)
+                    result.append(align.GetRMSD())
+                elif isinstance(target, ob.OBMol):
+                    align.SetTargetMol(target)
+                    align.Align()
+                    align.UpdateCoords(target)
+                    result.append(align.GetRMSD())
+                else:
+                    raise TypeError(f"RMSD can't handle {type(target)}.")
+        except TypeError:
+            raise
+        except Exception as e:
+            logger.error(f"Configuration.RMSD {e}")
+            raise
+        return result

molsystem/bonds.py

@@ -33,8 +33,6 @@ def __init__(self, configuration):
 
         self._system = None
 
-        self._bondset = self._configuration.bondset
-
         super().__init__(configuration._system_db, "bond")
 
     def __enter__(self):
@@ -78,7 +76,7 @@ def bondorders(self):
     @property
     def bondset(self):
         """The bondset for these bonds."""
-        return self._bondset
+        return self._configuration.bondset
 
     @property
     def bonds_for_asymmetric_bonds(self):

molsystem/cif.py

@@ -133,7 +133,11 @@ def to_cif_text(self):
                 lines.append("space_group_name_H-M_full   'P 1'")
             else:
                 spgname_system = symmetry.spacegroup_names_to_system[spgname]
-                lines.append(f"_space_group_{spgname_system}   '{spgname}'")
+                lines.append(f"_space_group_{spgname_system}       '{spgname}'")
+                hall = symmetry.hall_symbol
+                lines.append(f"_symmetry_space_group_name_Hall '{hall}'")
+                it_number = symmetry.IT_number
+                lines.append(f"_space_group_IT_number          {it_number}")
             lines.append(f"_cell_length_a   {a}")
             lines.append(f"_cell_length_b   {b}")
             lines.append(f"_cell_length_c   {c}")
@@ -231,6 +235,7 @@ def from_cif_text(self, text):
         None
         """
 
+        print("entering from_cif_text")
         result = ""
         cif = CifFile.ReadCif(io.StringIO(text))
 
@@ -264,12 +269,15 @@ def from_cif_text(self, text):
 
             # Where is the symmetry info?
             spgname = None
+            used_symops = False
             if "_space_group_symop" + dot + "operation_xyz" in data_block:
                 symdata = "_space_group_symop" + dot + "operation_xyz"
                 self.symmetry.symops = data_block[symdata]
+                used_symops = True
             elif "_symmetry_equiv" + dot + "pos_as_xyz" in data_block:
                 symdata = "_symmetry_equiv" + dot + "pos_as_xyz"
                 self.symmetry.symops = data_block[symdata]
+                used_symops = True
             else:
                 for section in (
                     "_space_group" + dot + "name_Hall",
@@ -516,6 +524,28 @@ def from_cif_text(self, text):
                 symop2=symop2s,
             )
 
+        # If used symops, need to find the spacegroup name(s)
+        if used_symops:
+            before = self.symmetry.symops
+            international_name = self.symmetry.find_spacegroup_from_operators()
+            self.symmetry.update_group(international_name)
+            after = self.symmetry.symops
+            if before != after:
+                print("!!!!!!!!!!!!!!!!!!!! SYMOPS CHANGED !!!!!!!!!!!!!!!!!!!!!")
+
+        import pprint
+
+        print("Reading CIF file")
+        print(f"{self.symmetry.hall_number=}")
+        print(f"{self.symmetry.hall_symbol=}")
+        print(f"{self.symmetry.group=}")
+        pprint.pprint(self.atoms.get_coordinates(asymmetric=True))
+        print("all atoms")
+        pprint.pprint(self.atoms.get_coordinates(asymmetric=False))
+
+        print("\n\nSymmetry info")
+        pprint.pprint(self.get_symmetry_data(self.symmetry.hall_number))
+
         return result
 
     def to_mmcif_text(self):