Minor bug fixes, file cleanup

molsym/molecule.py

@@ -2,6 +2,7 @@
 import qcelemental as qcel
 from dataclasses import dataclass
 from copy import deepcopy
+import sys
 global_tol = 1e-6
 
 @dataclass
@@ -69,6 +70,22 @@ def from_schema(cls, schema):
         for (idx, symb) in enumerate(atoms):
             masses[idx] = qcel.periodictable.to_mass(symb)
         return cls(atoms, coords, masses)
+
+    @classmethod
+    def from_psi4_molecule(cls, mol):
+        """
+        Class method for constructing a Molecule from a QCSchema object
+
+        :param schema: Schema dictionary to be converted to Molecule
+        :type schema: dict
+        :rtype: molsym.Molecule
+        """
+        if "psi4" not in sys.modules:
+            raise ImportError("Psi4 is required to use this function")
+        atoms = [mol.symbol(i) for i in range(mol.natom())]
+        coords = mol.geometry().to_array()
+        masses = [mol.mass(i) for i in range(mol.natom())]
+        return cls(atoms, coords, masses)
 
     @classmethod
     def from_file(cls, fn):

molsym/salcs/function_set.py

@@ -9,7 +9,6 @@ def __init__(self, symtext, fxn_list) -> None:
         # fxn_list: List of coordinate objects. Concrete classes will deal with operations on the coordinates
         self.fxns = fxn_list
         self.symtext = symtext
-        #self.fxn_map, self.phase_map = self.get_fxn_map()
         self.fxn_map = self.get_fxn_map()
         self.SE_fxns = self.get_symmetry_equiv_functions()
 

molsym/salcs/internal_coordinates.py

@@ -20,22 +20,12 @@ def operate_on_ic(self, ic_idx, symop):
         :return: New internal coordinate index and phase
         :rtype: int, float
         """
-        symbol = self.symtext.symels[symop].symbol
-        ic_type = self.ic_types[ic_idx][0]
-        if symbol[0] == "i" or symbol[0] == "S" or symbol[0] == "s": # s is for sigma
-            # Are we doing this twice??? TODO here
-            if ic_type in ["D", "O", "L"]: 
-                self.phase_map[ic_idx, symop] = -1.0
-            if ic_type == "D" or ic_type == "O": 
-                self.phase_map[ic_idx, symop] = -1.0
-            elif ic_type == "L":
-                self.phase_map[ic_idx, symop] = -1.0
         mapped_ic = []
         for atom in self.ic_list[ic_idx]:
             atom2 = int(self.symtext.atom_map[atom, symop])
             mapped_ic.append(atom2)
-        index, phase = self.ic_index(mapped_ic) 
-        return index, phase 
+        index, phase = self.ic_index(mapped_ic, symop, ic_idx)
+        return index, phase
 
     def get_fxn_map(self):
         """
@@ -71,35 +61,55 @@ def get_symmetry_equiv_functions(self):
                 reduced_seics = list(set(seics))
                 done += reduced_seics
                 SEICs.append(reduced_seics)
+
         return SEICs
 
-    def ic_index(self, ic):
+    def ic_index(self, ic, symop, ic_idx):
         """
         Permutes internal coordinate indices and tracks phase to avoid redundantly defined coordinates.
         Example, the angle between atoms 1,2,3 is the same as 3,2,1.
 
         :type ic: List[int]
         :rtype: (int, int)
         """
+        symbol = self.symtext.symels[symop].symbol
+        ic_type = self.ic_types[ic_idx][0]
+        phase_from_op = 1.0
+        if symbol[0] == "i" or symbol[0] == "S" or symbol[0] == "s": # s is for sigma
+            if ic_type in ["D", "O", "L"]:
+                phase_from_op = -1.0
         if len(ic) > 3:
             ic2 = deepcopy(ic)
             ic2[2], ic2[3] = ic[3], ic[2]
-            for f, funky in enumerate(self.ic_list):
+            for c, coord in enumerate(self.ic_list):
                 # TODO here
-                if ic == funky:
-                    return f, 1
-                elif list(reversed(ic)) == funky:
-                    return f, 1
-                elif ic2 == funky:
-                    return f, -1
-                elif list(reversed(ic2)) == funky:
-                    return f, -1
+                if ic == coord:
+                    return c, 1 * phase_from_op
+                elif list(reversed(ic)) == coord:
+                    return c, 1 * phase_from_op
+                elif ic2 == coord:
+                    return c, -1 * phase_from_op
+                elif list(reversed(ic2)) == coord:
+                    return c, -1 * phase_from_op
+
+        elif len(ic) == 3:
+            #first, loop through IC list, only after the list is exhausted, loop through reverse indices
+            for c, coord in enumerate(self.ic_list):
+                if ic == coord:
+                    return c, 1
+                elif list(reversed(ic)) == coord:
+                    return c, 1
+            ic2 = deepcopy(ic)
+            #for c, coord in enumerate(self.ic_list):
+            #    if len(coord) == 3:
+            #        if ic[0] == coord[0] and ic[1] == coord[1]:
+            #            return c, -1
         else:
-            for f, funky in enumerate(self.ic_list):
-                if ic == funky:
-                    return f, 1
-                elif list(reversed(ic)) == funky:
-                    return f, 1
+            for c, coord in enumerate(self.ic_list):
+                if ic == coord:
+                    return c, 1
+                elif list(reversed(ic)) == coord:
+                    return c, 1
 
     def span(self, SE_fxn):
         chorker_loaf = np.zeros(self.fxn_map.shape)
@@ -130,4 +140,4 @@ def special_function(self, salc, coord, sidx, irrmat):
         ic2 = self.fxn_map[coord, sidx]
         p = self.phase_map[coord, sidx]
         salc[:,:,ic2] += (irrmat[sidx, :, :]) * p
-        return salc
+        return salc

molsym/salcs/salc.py

@@ -205,10 +205,11 @@ def finish_building(self, orthogonalize=False, remove_complexity=False):
                         B[:,col] /= np.linalg.norm(B[:,col])
                     for idx, salc in enumerate(self.salcs_by_irrep[irrep_idx]):
                         self.salcs[salc].coeffs = B[:,idx]
+                    #raise Exception("BEANS")
                 else:
                     n_pf_sets = round(len(self.salcs_by_irrep[irrep_idx]) / irrep.d)
-                    #B1 = self.basis_transformation_matrix[:,:irrep.d]
-                    B1 = self.basis_transformation_matrix[:,self.salcs_by_irrep[irrep_idx]]
+                    salcs_in_this_irrep = self.salcs_by_irrep[irrep_idx]
+                    B1 = self.basis_transformation_matrix[:,salcs_in_this_irrep]
                     # Gram-Schmidt orthogonalize columns of B1
                     trans_mat = np.eye(n_pf_sets)
                     for col_idx in range(1, n_pf_sets):
@@ -218,13 +219,13 @@ def finish_building(self, orthogonalize=False, remove_complexity=False):
                             trans_mat[:,col_idx] -= proj * trans_mat[:,gs_idx]
                         nrm = np.linalg.norm(B1[:,col_idx])
                         B1[:,col_idx] /= nrm
+                        self.salcs[salcs_in_this_irrep[col_idx]].coeffs = B1[:,col_idx]
                         trans_mat[:,col_idx] /= nrm
+                    B1 = self.basis_transformation_matrix[:,salcs_in_this_irrep]
                     # Transform other partner function sets according to the Gram-Schmidt orthogonalization of B1
-                    for pf_idx in range(irrep.d):
-                        pfxn_set = [pf_idx*irrep.d + i for i in range(n_pf_sets)]
-                        Bi = self.basis_transformation_matrix[:,pfxn_set]
+                    for pf_idx in range(1,irrep.d):
+                        pfxn_set = [pf_idx*n_pf_sets + i for i in range(n_pf_sets)]
+                        Bi = self.basis_transformation_matrix[:,[salcs_in_this_irrep[idx] for idx in pfxn_set]]
                         Bi_trans = Bi @ trans_mat
-                        for Bidx,salc_idx in enumerate(pfxn_set):
-                            self.salcs[salc_idx].coeffs = Bi_trans[:,Bidx] / np.linalg.norm(Bi_trans[:,Bidx])
-
-
+                        for Bidx, salc_idx in enumerate(pfxn_set):
+                            self.salcs[salcs_in_this_irrep[salc_idx]].coeffs = Bi_trans[:,Bidx] / np.linalg.norm(Bi_trans[:,Bidx])

test/test_projection_op.py

@@ -43,6 +43,7 @@
                          ["A_1", "A_1", "E", "E", "E", "E"],
                          ["A'","A'","A'","A'","A'","A'","A'","A'","A''","A''","A''","A''"],
                          ["A_1", "A_1", "E", "E", "T_2", "T_2", "T_2", "T_2", "T_2", "T_2"]]
+
 @pytest.mark.parametrize("i", [i for i in range(len(fns))])
 def test_internal_coordinate_SALCs(i):
     mol = molsym.Molecule.from_file("test/xyz/"+fns[i]+".xyz")