Merge pull request #70 from molssi-seamm/dev

Bugfix and improvements to symmetry handling.

HISTORY.rst

@@ -1,6 +1,12 @@
 =======
 History
 =======
+2023.11.5 -- Bugfix and improved symmetry handling
+    * Fixed bug with symmetry operators containing blanks, e.g. 'x, y, z' rather than
+      'x,y,z'
+    * Added handling of symmetry when get properties of atoms
+    * Added method to lower symmetry to P1/C1
+
 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,

molsystem/atoms.py

@@ -5,6 +5,7 @@
 
 from itertools import zip_longest
 import logging
+import pprint  # noqa: F401
 from typing import Any, Dict, TypeVar
 
 import numpy as np
@@ -772,7 +773,7 @@ def diff(self, other):
 
         return result
 
-    def get_as_dict(self, *args):
+    def get_as_dict(self, *args, asymmetric=False):
         """Return the atom data as a Python dictionary of lists.
 
         Parameters
@@ -785,12 +786,29 @@ def get_as_dict(self, *args):
         dict(str: [])
             A dictionary whose keys are the column names and values as lists
         """
-        rows = self.atoms(*args)
-        columns = [x[0] for x in rows.description]
-        data = {key: [] for key in columns}
-        for row in rows:
-            for key, value in zip(columns, row):
-                data[key].append(value)
+        if asymmetric:
+            rows = self.atoms(*args)
+            columns = [x[0] for x in rows.description]
+            data = {key: [] for key in columns}
+            for row in rows:
+                for key, value in zip(columns, row):
+                    data[key].append(value)
+        else:
+            data = {}
+            properties = [*self.attributes.keys()]
+            for name in properties:
+                if name in ("id", "x", "y", "z", "vx", "vy", "vz", "name"):
+                    continue
+                data[name] = self.get_property(name)
+
+            if "name" in properties:
+                data["name"] = self.get_names()
+
+            xyz = self.get_coordinates()
+            data["x"] = [v[0] for v in xyz]
+            data["y"] = [v[1] for v in xyz]
+            data["z"] = [v[2] for v in xyz]
+            data["id"] = [*range(len(xyz))]
 
         return data
 
@@ -1106,6 +1124,34 @@ def get_n_atoms(self, *args):
         self.cursor.execute(sql, parameters)
         return self.cursor.fetchone()[0]
 
+    def get_property(self, name, asymmetric=False):
+        """Return the property from the atoms, expanded to symmetric atoms.
+
+        Parameters
+        ----------
+        name : str
+            The property (attribute) to return
+
+        asymmetric : bool = False
+            Return just the names for the asymmetric atoms.
+
+        Returns
+        -------
+        [any]
+            The values of the property on the symmetric atoms
+        """
+        if name == "name":
+            return self.get_names(asymmetric=asymmetric)
+
+        data = self.get_column_data(name)
+
+        symmetry = self.configuration.symmetry
+        if asymmetric or symmetry.n_symops == 1:
+            return data
+
+        # Expand to the asymmetric atoms
+        return [data[i] for i in symmetry.atom_to_asymmetric_atom]
+
     def get_velocities(
         self,
         fractionals=True,
@@ -1180,9 +1226,13 @@ def set_coordinates(self, xyz, fractionals=True):
 
         # May need to handle symmetry
         if n_coords != self.n_asymmetric_atoms and n_coords == self.n_atoms:
+            # print("set_coordinates: symmetrizing the coordinates")
+            # pprint.pprint(xyz)
             xyz, error = self.configuration.symmetry.symmetrize_coordinates(
                 xyz, fractionals=fractionals
             )
+            # print("results in")
+            # pprint.pprint(xyz)
 
         x_column = self.get_column("x")
         y_column = self.get_column("y")

molsystem/cif.py

@@ -220,6 +220,7 @@ def to_cif_text(self):
                     lines.append(f"  {names[i]} {names[j]} {r:.4f} {op1} {op2}")
 
         # And that is it!
+        lines.append("")
         return "\n".join(lines)
 
     def from_cif_text(self, text):
@@ -235,7 +236,6 @@ def from_cif_text(self, text):
         None
         """
 
-        print("entering from_cif_text")
         result = ""
         cif = CifFile.ReadCif(io.StringIO(text))
 
@@ -533,18 +533,18 @@ def from_cif_text(self, text):
             if before != after:
                 print("!!!!!!!!!!!!!!!!!!!! SYMOPS CHANGED !!!!!!!!!!!!!!!!!!!!!")
 
-        import pprint
+        # 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("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))
+        # print("\n\nSymmetry info")
+        # pprint.pprint(self.get_symmetry_data(self.symmetry.hall_number))
 
         return result
 
@@ -659,6 +659,7 @@ def to_mmcif_text(self):
                 lines.append(f"MOL1 {nm1} {nm2} {order}")
 
         # And that is it!
+        lines.append("")
         return "\n".join(lines)
 
     def from_mmcif_text(self, text):

molsystem/configuration.py

@@ -679,6 +679,69 @@ def get_symmetry_data(self, hall_number):
         tmp = spglib.get_symmetry_dataset(cell_in, hall_number=hall_number)
         return tmp
 
+    def lower_symmetry(self, other=None):
+        """Lower the symmetry to P1 or C1, optionally from another configuration.
+
+        Parameters
+        ----------
+        other : molsystem._Configuration = None
+            Another configuration to use as the source
+        """
+        if other is None:
+            other = self
+
+        periodicity = other.periodicity
+
+        if periodicity != 0:
+            cell_parameters = other.cell.parameters
+
+        # Get the atom and bond information for low symmetry
+        atom_data = other.atoms.get_as_dict()
+        del atom_data["id"]
+        bond_data = other.bonds.get_as_dict()
+        del bond_data["id"]
+
+        # It is not clear that it makes sense to handle velocities, so drop
+        if "vx" in atom_data:
+            del atom_data["vx"]
+            del atom_data["vy"]
+            del atom_data["vz"]
+
+        self.clear()
+        self.new_atomset()
+        self.new_bondset()
+        if periodicity != 0:
+            self.new_cell()
+            self.cell.parameters = cell_parameters
+        self.new_symmetry()
+        self.periodicity = periodicity
+
+        ids = self.atoms.append(**atom_data)
+
+        iatoms = bond_data["i"]
+        jatoms = bond_data["j"]
+        bond_data["i"] = [ids[i] for i in iatoms]
+        bond_data["j"] = [ids[j] for j in jatoms]
+
+        self.bonds.append(**bond_data)
+
+        # Finally, copy over the charge multiplicity, etc.
+        if self is not other:
+            self.charge = other.charge
+            self.spin_multiplicity = other.spin_multiplicity
+            self.n_active_electrons = other.n_active_electrons
+            self.n_active_orbitals = other.n_active_orbitals
+
+    def new_atomset(self):
+        """Setup a new, empty atomset for this configuration."""
+        self._atomset = self.system_db["atomset"].append(n=1)[0]
+        self.db.execute(
+            "UPDATE configuration SET atomset = ? WHERE id = ?",
+            (self._atomset, self.id),
+        )
+        self.db.commit()
+        return self._atomset
+
     def new_bondset(self):
         """Setup a new, empty bondset for this configuration."""
         self._bondset = self.system_db["bondset"].append(n=1)[0]
@@ -689,6 +752,29 @@ def new_bondset(self):
         self.db.commit()
         return self._bondset
 
+    def new_cell(self):
+        """Setup a new cell for this configuration."""
+        self._cell_id = self.system_db["cell"].append(n=1)[0]
+        sql = "UPDATE configuration SET cell = ? WHERE id = ?"
+        self.db.execute(sql, (self._cell_id, self.id))
+        self.db.commit()
+
+        return self.cell
+
+    def new_symmetry(self):
+        """Create a new symmetry object for this configuration."""
+        self._symmetry_id = None
+
+        table = _Table(self.system_db, "symmetry")
+        self._symmetry_id = table.append()[0]
+        self.cursor.execute(
+            "UPDATE configuration SET symmetry = ? WHERE id = ?",
+            (self._symmetry_id, self.id),
+        )
+        self.db.commit()
+
+        return self.symmetry
+
     def primitive_cell(self, symprec=1.0e-05, spg=True):
         """Find the symmetry of periodic systems and transform to conventional cell.
 

molsystem/symmetry.py

@@ -2,14 +2,14 @@
 
 from fractions import Fraction
 import logging
-import pprint
+import pprint  # noqa: F401
 import re
 
 import numpy as np
 import spglib
 
 logger = logging.getLogger(__name__)
-logger.setLevel("INFO")
+# logger.setLevel("INFO")
 
 
 class _Symmetry(object):
@@ -246,6 +246,7 @@ def symmetry_matrices(self):
                     [0.0, 0.0, 0.0, 1.0],
                 ]
                 for row, op in enumerate(line.split(",")):
+                    op = op.strip()
                     match = regexp.fullmatch(op)
                     if match is None:
                         raise ValueError(
@@ -661,9 +662,9 @@ def _expand(self):
                         f"Mismatch of number of atoms in symmetry: {uvw0.shape[0]} != "
                         f"{n_atoms}"
                     )
-                logger.info("Original coordinates")
-                logger.info(uvw0)
-                logger.info(f"{uvw0.shape=}")
+                logger.debug("Original coordinates")
+                logger.debug(uvw0)
+                logger.debug(f"{uvw0.shape=}")
 
                 logger.debug(f"{n_atoms=}")
                 uvw = np.ndarray((n_atoms, 4))
@@ -689,23 +690,22 @@ def _expand(self):
                 logger.debug(tmp)
 
                 n_sym_atoms = 0
-                print(f"{n_atoms=}")
                 for i in range(n_atoms):
                     values, I1, I2 = np.unique(
                         np.round(tmp[:, :3, i], 4),
                         axis=0,
                         return_index=True,
                         return_inverse=True,
                     )
-                    logger.info(f"{i=}")
-                    pprint.pprint(np.round(tmp[:, :3, i], 4).tolist())
-                    logger.info("values")
-                    pprint.pprint(values.tolist())
+                    logger.debug(f"{i=}")
+                    # pprint.pprint(np.round(tmp[:, :3, i], 4).tolist())
+                    logger.debug("values")
+                    # pprint.pprint(values.tolist())
                     logger.debug(values)
-                    logger.info("I1")
-                    logger.info(I1)
-                    logger.info("I2")
-                    logger.info(I2)
+                    logger.debug("I1")
+                    logger.debug(I1)
+                    logger.debug("I2")
+                    logger.debug(I2)
                     I2 += n_sym_atoms
                     self._atom_generators.append(I1.tolist())
                     self._symop_to_atom.append(I2.tolist())

tests/test_cif.py

@@ -298,7 +298,8 @@ def test_to_mmcif_text(AceticAcid):
 MOL1 C C2 SING
 MOL1 C2 O DOUB
 MOL1 C2 O2 SING
-MOL1 O2 H4 SING"""
+MOL1 O2 H4 SING
+"""
 
     text = AceticAcid.to_mmcif_text()
 
@@ -309,7 +310,7 @@ def test_to_mmcif_text(AceticAcid):
 
 
 def test_to_cif_text(vanadium):
-    """Write a manually created system to mmcif"""
+    """Write a manually created system to cif"""
     correct = """\
 # Generated by MolSSI SEAMM
 data_V
@@ -337,7 +338,8 @@ def test_to_cif_text(vanadium):
  _atom_site_fract_z
  _atom_site_occupancy
   V V1  1  0.000 0.000 0.000  1
-  V V2  1  0.500 0.500 0.500  1"""
+  V V2  1  0.500 0.500 0.500  1
+"""
 
     text = vanadium.to_cif_text()
 
@@ -467,7 +469,8 @@ def test_benzene_to_cif(benzene):
   C1 H1 0.9300 . .
   C2 C3 1.3761 . 5_566
   C2 H2 0.9306 . .
-  C3 H3 0.9296 . ."""
+  C3 H3 0.9296 . .
+"""
 
     # benzene.symmetry.loglevel = "DEBUG"
     # benzene.bonds.loglevel = "DEBUG"