Workflow for linear molecules

molsym/pgdetect/flowchart.py

@@ -19,9 +19,14 @@ def find_point_group(mol):
     paxis = [0,0,0]
     saxis = [0,0,0]
     moit = calcmoit(mol)
-    Ia_mol, Ib_mol, Ic_mol = np.sort(np.linalg.eigh(moit)[0])
+    I_evals, I_evecs = np.linalg.eigh(moit)
+    mask = np.argsort(I_evals)
+    Ia_vec, Ib_vec, Ic_vec = I_evecs[mask]
+    Ia_mol, Ib_mol, Ic_mol = I_evals[mask]
+    #Ia_mol, Ib_mol, Ic_mol = np.sort(np.linalg.eigh(moit)[0])
     # Linear tops
     if np.isclose(Ia_mol, 0.0, atol=mol.tol):
+        paxis = Ia_vec
         if isequivalent(mol, mol.transform(inversion_matrix())):
             pg = "D0h"
         else:

molsym/salcs/projection_op.py

@@ -36,9 +36,9 @@ def eckart_conditions(symtext, translational=True, rotational=True):
         y[3 * i + 1] = smass
         z[3 * i + 2] = smass
         atomx, atomy, atomz = mol.coords[i, 0], mol.coords[i, 1], mol.coords[i, 2]
-        tval0 = atomx * evec[0,0] + atomy * evec[1,0] + atomz * evec[2, 0];
-        tval1 = atomx * evec[0,1] + atomy * evec[1,1] + atomz * evec[2, 1];
-        tval2 = atomx * evec[0,2] + atomy * evec[1,2] + atomz * evec[2, 2];
+        tval0 = atomx * evec[0,0] + atomy * evec[1,0] + atomz * evec[2, 0]
+        tval1 = atomx * evec[0,1] + atomy * evec[1,1] + atomz * evec[2, 1]
+        tval2 = atomx * evec[0,2] + atomy * evec[1,2] + atomz * evec[2, 2]
         rx[3 * i + 0] = (tval1 * evec[0,2] - tval2 * evec[0,1]) * smass
         rx[3 * i + 1] = (tval1 * evec[1,2] - tval2 * evec[1,1]) * smass
         rx[3 * i + 2] = (tval1 * evec[2,2] - tval2 * evec[2,1]) * smass
@@ -77,6 +77,9 @@ def ProjectionOp(symtext, fxn_set, project_Eckart=True):
     numred = len(fxn_set)
     salcs = SALCs(symtext, fxn_set)
     orthogonalize = False
+    #project_Eckart = False
+    if project_Eckart:
+        eckart_cond = eckart_conditions(symtext)
     for ir, irrep in enumerate(symtext.irreps):
         if symtext.pg.is_linear:
             irrmat = None
@@ -94,7 +97,7 @@ def ProjectionOp(symtext, fxn_set, project_Eckart=True):
             # Project out Eckart conditions when constructing SALCs of Cartesian displacements
             if isinstance(fxn_set, CartesianCoordinates) and project_Eckart:
                 orthogonalize = True
-                eckart_cond = eckart_conditions(symtext)
+                #eckart_cond = eckart_conditions(symtext)
                 for i in range(irrep.d):
                     for j in range(irrep.d):
                         if not np.allclose(salc[i,j,:], np.zeros(salc[i,j,:].shape), atol=salcs.tol):

molsym/symtext/general_irrep_mats.py

@@ -57,20 +57,20 @@
     sh = Symel("sigma_h", np.array([0,0,1]), reflection_matrix(z), 0, 0, "sigma_h")
     if pg.is_linear:
         if pg.family == "C":
-            symels = [#Symel("E", None, None, None, None, None), E is included in C 
-                      Symel("C", z, None, None, None, None), 
-                      Symel("sigma_v", None, None, None, None, None)]
+            symels = [Symel("E", None, np.eye(3), None, None, None),
+                      Symel("C", z, np.eye(3), 0, 0, "E"), 
+                      Symel("sigma_v", [1,0,0], reflection_matrix([1,0,0]), 0, 0, "sigma_v")]
             irreps = [Irrep("Sigma^+", None, None, 1), Irrep("Sigma^-", None, None, 1)]
             irreps += [Irrep(greek[i], i, None, 2) for i in range(len(greek))]
             irrep_mats = None
             return symels, irreps, irrep_mats
         elif pg.family == "D":
-            symels = [#Symel("E", None, None, None, None, None), 
-                      Symel("C", z, None, None, None, None), 
-                      Symel("sigma_v", None, None, None, None, None),
-                      #Symel("i", None, None, None, None, None), i is included in S
-                      Symel("S", z, None, None, None, None),
-                      Symel("C_2'", None, None, None, None, None)]
+            symels = [Symel("E", None, np.eye(3), None, None, None), 
+                      Symel("C", z, np.eye(3), None, None, None), 
+                      Symel("sigma_v", [1,0,0], reflection_matrix([1,0,0]), None, None, None),
+                      Symel("i", None, -1*np.eye(3), None, None, None),
+                      Symel("S", z, reflection_matrix(z), None, None, None),
+                      Symel("C_2'", [1,0,0], Cn([1,0,0], 2), None, None, None)]
             irreps = [Irrep("Sigma_g^+", None, None, 1), Irrep("Sigma_g^-", None, None, 1)]
             irreps += [Irrep(greek[i]+"_g", i, None, 2) for i in range(len(greek))]
             irreps += [Irrep("Sigma_u^+", None, None, 1), Irrep("Sigma_u^-", None, None, 1)]

molsym/symtext/symtext.py

@@ -38,7 +38,11 @@
             self.order = len(symels)
         self.irreps = irreps
         self.irrep_mats = irrep_mats
-        self.get_character_table()
+        if self.pg.is_linear:
+            self.get_character_table_linear()
+        else:
+            self.get_character_table()
+
 
     def __len__(self):
         return len(self.symels)
@@ -125,6 +129,18 @@
         for irrep_idx, irrep in enumerate(self.irreps):
             for class_idx, class_name in enumerate(self.classes):
                 self.character_table[irrep_idx,class_idx] = np.trace(self.irrep_mats[irrep.symbol][self.symel_to_class_map.index(class_idx)])
+
+    def get_character_table_linear(self):
+        if self.pg.family == "C":
+            self.classes = ["E", "C", "sigma_v"]
+            self.symel_to_class_map = [0, 1, 2]
+            self.class_orders = [1, 0, 0]
+            self.character_table = None
+        elif self.pg.family == "D":
+            self.classes = ["E", "C", "sigma_v", "i", "S", "C_2'"]
+            self.symel_to_class_map = [0, 1, 2, 3, 4, 5]
+            self.class_orders = [1, 0, 0, 1, 0, 0]
+            self.character_table = None
 
     def direct_product(self, *args):
         """
@@ -158,7 +174,6 @@
         p = np.multiply(dp_vector, self.class_orders)
         return round(p.sum()/(self.order)) > 0
 
-
     def reduction_coefficients(self, rrep_characters):
         """
         Return reduction coefficients for the characters of a representation.

molsym/symtext/symtext_helper.py

@@ -214,7 +214,10 @@
             else:
                 raise Exception(f"Atom {atom} not mapped to another atom under symel i")
         return np.column_stack((amap, ungerade_map))
-    return amap
+    elif pg.family == "C":
+        return np.arange(mol.natoms, dtype=int).reshape((mol.natoms,1))
+    else:
+        raise Exception(f"Linear point group {pg.str} not supported for atom mapping")
 
 def where_you_go(mol, atom, symel):
     """