Implement permutation handling for geometric symmetry at API level

cosymlib/__init__.py

@@ -68,6 +68,47 @@ def _get_table_format(labels, molecules_names, data):
     return txt
 
 
+def _get_table_format_chir(labels, molecules_names, data, axis):
+
+    txt = 'Structure'
+    max_len_name = 12
+    names = []
+    for name in molecules_names:
+        if len(name) > max_len_name:
+            names.append(name[:(max_len_name - 1)])
+        elif len(name) > 0:
+            names.append(name)
+
+    n = max_len_name - 3
+    for label in labels:
+        n += len(label)
+        txt += '{}'.format(label.rjust(n))
+        n = 11 - len(label)
+    txt += '\n\n'
+
+    for idx, name in enumerate(names):
+        txt += '{:{width}}'.format(name + ',', width=max_len_name)
+        n = 11
+        for idn, label in enumerate(labels):
+            txt += '{:>{width}.{prec}f},'.format(data[idn][idx], width=n, prec=3)
+            n = 10
+        txt += '\n'
+
+        txt += '\n'
+
+        for idn, c_lab in enumerate(['x', 'y', 'z']):
+            txt += ' ' * (max_len_name-1) + c_lab + ' '
+
+            for c in np.array(axis).T[idn][0]:
+                txt += '{:>{width}.{prec}f},'.format(c, width=n, prec=3)
+            n = 10
+            txt += '\n'
+
+        txt += '\n'
+    txt += '\n'
+    return txt
+
+
 def _get_axis_info(molecule, group, axis, axis2, center):
 
     txt = ''
@@ -319,7 +360,7 @@ def print_geometric_measure_info(self, label, multi=1, central_atom=0, center=No
 
         output.write(txt)
 
-    def print_geometric_symmetry_measure(self, label, central_atom=0, center=None, output=sys.stdout):
+    def print_geometric_symmetry_measure(self, label, central_atom=0, center=None, permutation=None, output=sys.stdout):
         """
         Prints geometric symmetry measure in format
 
@@ -328,26 +369,32 @@ def print_geometric_symmetry_measure(self, label, central_atom=0, center=None, o
         :param central_atom: Position of the central atom
         :type central_atom: int
         :param center: Center of symmetry in Cartesian coordinates. If None center is optimized
-        :type center: list
+        :type center: list, tuple
+        :param permutation: Define permutation
+        :type permutation: list, tuple
         :param output: Display hook
         :type output: hook
         """
         kwargs = _get_symmetry_arguments(locals())
 
         txt = 'Evaluating symmetry operation : {}\n \n'.format(label)
+
+        txt += '              CSM      permutation\n'
         for idx, molecule in enumerate(self._molecules):
             csm = molecule.geometry.get_symmetry_measure(**kwargs)
+            permutation = molecule.geometry.get_symmetry_permutation(**kwargs)
+
             max_name = len(max(molecule.name, key=len))
             txt += '{} '.format(molecule.name)
             if max_name < 9:
                 n = 18 - len(molecule.name)
             else:
                 n = 9 + max_name - len(molecule.name)
-            txt += '{:{width}.{prec}f}\n'.format(csm, width=n, prec=3)
+            txt += '{:{width}.{prec}f}'.format(csm, width=n, prec=3) + '   {}\n'.format(permutation)
 
         output.write(txt)
 
-    def print_symmetry_nearest_structure(self, label, central_atom=0, center=None, output=sys.stdout):
+    def print_symmetry_nearest_structure(self, label, central_atom=0, center=None, permutation=None, output=sys.stdout):
         """
         Prints the nearest structure to ideal symmetric structure
 
@@ -357,6 +404,8 @@ def print_symmetry_nearest_structure(self, label, central_atom=0, center=None, o
         :type central_atom: int
         :param center: Center of symmetry in Cartesian coordinates. If None center is optimized
         :type center: int
+        :param permutation: Define permutation
+        :type permutation: list, tuple
         :param output: Display hook
         :type output: hook
         """
@@ -366,15 +415,28 @@ def print_symmetry_nearest_structure(self, label, central_atom=0, center=None, o
             geometry = molecule.geometry.get_symmetry_nearest_structure(**kwargs)
             output.write(file_io.get_file_xyz_txt(geometry))
 
-    def print_chirality_measure(self, order=1, central_atom=0, center=None, output=sys.stdout):
+    def print_chirality_measure(self, order=1, central_atom=0, center=None, permutation=None, output=sys.stdout):
+        """
+        Prints the chirality measure
 
+        :param order: Order of the chirality measure (1: Cs, 2:Ci, n:S_n)
+        :type order: int
+        :param central_atom: Position of the central atom
+        :type central_atom: int
+        :param center: Center of symmetry in Cartesian coordinates. If None center is optimized
+        :type center: int
+        :param permutation: Define permutation
+        :type permutation: list, tuple
+        :param output: Display hook
+        :type output: hook
+        """
         if order < 1:
             raise Exception('Chirality order should be greater than 0')
 
         if order == 1:
             reference = ['Cs']
-        elif order == 2:
-            reference = ['Ci']
+        # elif order == 2:
+        #     reference = ['Ci']
         elif int(order) > 8:
             raise ValueError('Maximum value available for Sn chirality measure is 8')
         else:
@@ -386,12 +448,18 @@ def print_chirality_measure(self, order=1, central_atom=0, center=None, output=s
         reference = [r.lower() for r in reference]
 
         csm_list = []
+        axis_list = []
         for label in reference:
-            csm_list.append([geometry.get_symmetry_measure(label, central_atom=central_atom, center=center)
+            csm_list.append([geometry.get_symmetry_measure(label, central_atom=central_atom,
+                                                           center=center, permutation=permutation)
+                             for geometry in self.get_geometries()])
+            axis_list.append([geometry.get_symmetry_optimum_axis(label, central_atom=central_atom,
+                                                                 center=center, permutation=permutation)
                              for geometry in self.get_geometries()])
 
         molecules_names = [molecule.name for molecule in self._molecules]
-        txt = _get_table_format(reference, molecules_names, csm_list)
+        txt = 'Chirality measure and symmetry axis\n\n'
+        txt += _get_table_format_chir(reference, molecules_names, csm_list, axis_list)
 
         output.write(txt)
 

cosymlib/file_io/__init__.py

@@ -555,7 +555,10 @@ def get_connectivity_from_file(file_name):
 
     connect = []
     for line in file_txt.split('\n'):
-            connect.append([int(num) for num in line.split()])
+            try:
+                connect.append([int(num) for num in line.split()])
+            except ValueError:
+                raise Exception('Error reading connectivity file')
 
     connectivity = []
     for atom in connect:
@@ -565,6 +568,18 @@ def get_connectivity_from_file(file_name):
     return connectivity
 
 
+def get_permutation_from_file(file_name):
+    permutation = []
+    with open(file_name, 'r') as f:
+        for line in f:
+            try:
+                permutation += np.array(line.strip().split(), dtype=int).tolist()
+            except ValueError:
+                raise Exception('Error reading permutation file')
+
+    return permutation
+
+
 # Get OUPUT files
 def get_file_xyz_txt(structure):
 

cosymlib/molecule/geometry/__init__.py

@@ -91,6 +91,12 @@ def __init__(self,
     def __len__(self):
         return len(self._positions)
 
+    def __str__(self):
+        txt = '{}\n\n'.format(len(self._positions))
+        for s, p in zip(self._symbols, self._positions):
+            txt += '{:3}'.format(s) + '{:10.5f} {:10.5f} {:10.5f}\n'.format(*p)
+        return txt
+
     @property
     def name(self):
         return self._name
@@ -194,7 +200,7 @@ def get_generalized_coordinate(self, shape_label1, shape_label2, central_atom=0)
 
     # Symmetry methods
     @set_parameters
-    def get_symmetry_measure(self, label, central_atom=0, center=None, multi=1):
+    def get_symmetry_measure(self, label, central_atom=0, center=None, multi=1, permutation=None):
         """
         Get the symmetry measure
 
@@ -204,14 +210,35 @@ def get_symmetry_measure(self, label, central_atom=0, center=None, multi=1):
         :type central_atom: int
         :param center: center of the measure in Cartesian coordinates
         :type center: list
+        :param permutation: define permutation
+        :type permutation: list, tuple
         :return: The symmetry measure
         :rtype: float
         """
 
         return self._symmetry.measure(label)
 
     @set_parameters
-    def get_symmetry_nearest_structure(self, label, central_atom=0, center=None,  multi=1):
+    def get_symmetry_permutation(self, label, central_atom=0, center=None, multi=1, permutation=None):
+        """
+        Get the optimum atoms permutation for geometrical symmetry measures
+
+        :param label: Symmetry point group
+        :type label: str
+        :param central_atom: central atom position (0 if no central atom)
+        :type central_atom: int
+        :param center: center of the measure in Cartesian coordinates
+        :type center: list
+        :param permutation: define permutation
+        :type permutation: list, tuple
+        :return: The symmetry measure
+        :rtype: float
+        """
+
+        return self._symmetry.optimum_permutation(label)
+
+    @set_parameters
+    def get_symmetry_nearest_structure(self, label, central_atom=0, center=None, multi=1, permutation=None):
         """
         Returns the nearest ideal structure
 
@@ -221,6 +248,8 @@ def get_symmetry_nearest_structure(self, label, central_atom=0, center=None,  mu
         :type central_atom: int
         :param center: center of the measure in Cartesian coordinates
         :type center: list
+        :param permutation: Define permutation
+        :type permutation: list, tuple
         :return: The structure
         :rtype: Geometry
         """
@@ -230,7 +259,7 @@ def get_symmetry_nearest_structure(self, label, central_atom=0, center=None,  mu
         # return self._symmetry.nearest_structure(label)
 
     @set_parameters
-    def get_symmetry_optimum_axis(self, label, central_atom=0, center=None):
+    def get_symmetry_optimum_axis(self, label, central_atom=0, permutation=None, center=None):
         return self._symmetry.optimum_axis(label)
 
     def get_pointgroup(self, tol=0.01):

cosymlib/symmetry/__init__.py

@@ -6,14 +6,15 @@
 import numpy as np
 
 
-def _get_key_symgroup(label, center, central_atom, connectivity, multi, connect_thresh):
+def _get_key_symgroup(label, center, central_atom, connectivity, multi, connect_thresh, permutation):
     group_key = label.lower()
     center_key = ' '.join(['{:10.8f}'.format(n) for n in center]) if center is not None else None
     connectivity_key = np.array2string(np.array(connectivity), precision=10) if connectivity is not None else None
     multi_key = int(multi)
     central_atom_key = int(central_atom)
     connect_thresh_key = '{:10.8f}'.format(connect_thresh)
-    return group_key, center_key, central_atom_key, connectivity_key, multi_key, connect_thresh_key
+    permutation_key = tuple(permutation) if permutation is not None else None
+    return group_key, center_key, central_atom_key, connectivity_key, multi_key, connect_thresh_key, permutation_key
 
 
 def _get_key_wfnsym(group, vector_axis1, vector_axis2, center, alpha_occupancy, beta_occupancy):
@@ -54,6 +55,7 @@ def __init__(self,
                  multi=1,
                  axis=None,
                  axis2=None,
+                 permutation=None,
                  ):
 
         try:
@@ -79,6 +81,7 @@ def __init__(self,
         self._multi = multi
         self._axis = axis
         self._axis2 = axis2
+        self._permutation = permutation
         self._results = {}
 
     # Modifier methods
@@ -109,7 +112,7 @@ def _get_symgroup_results(self, group):
 
         # Crude calculation call methods
         key = _get_key_symgroup(group, self._center, self._central_atom, self._connectivity, self._multi,
-                                self._connect_thresh)
+                                self._connect_thresh, self._permutation)
         if key not in self._results:
             self._results[key] = Symgroupy(self._coordinates,
                                            group=group,
@@ -119,8 +122,13 @@ def _get_symgroup_results(self, group):
                                            center=self._center,
                                            connectivity=self._connectivity,
                                            connect_thresh=self._connect_thresh,
-                                           permutation=None)
+                                           permutation=self._permutation)
 
+        permu = self._results[key].optimum_permutation
+        key_2 = _get_key_symgroup(group, self._center, self._central_atom, self._connectivity, self._multi,
+                                self._connect_thresh, permu)
+
+        self._results[key_2] = self._results[key]
         return self._results[key]
 
     def _get_wfnsym_results(self, group):
@@ -157,7 +165,7 @@ def _get_wfnsym_results(self, group):
         elif isinstance(self._electronic_structure, ProtoElectronicStructure):
 
             key = _get_key_symgroup(group, self._center, self._central_atom, self._connectivity, self._multi,
-                                    self._connect_thresh)
+                                    self._connect_thresh, self._permutation)
             if self._center is None:
                 coord_a = np.array(self._coordinates)
                 sym_m = []

requirements.txt

@@ -5,6 +5,6 @@
 numpy>=1.16.6
 matplotlib
 PyYAML>=5.1
-symgroupy
-wfnsympy
+symgroupy>=0.5.9
+wfnsympy>=0.3.5
 huckelpy>=0.2