diff --git a/molsym/pgdetect/flowchart.py b/molsym/pgdetect/flowchart.py index 9d6ce67..62c6bf7 100755 --- a/molsym/pgdetect/flowchart.py +++ b/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: diff --git a/molsym/salcs/projection_op.py b/molsym/salcs/projection_op.py index baaeda5..334ad11 100755 --- a/molsym/salcs/projection_op.py +++ b/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): diff --git a/molsym/symtext/general_irrep_mats.py b/molsym/symtext/general_irrep_mats.py index ae0b63d..3408d78 100755 --- a/molsym/symtext/general_irrep_mats.py +++ b/molsym/symtext/general_irrep_mats.py @@ -57,20 +57,20 @@ def pg_to_symels(PG): 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)] diff --git a/molsym/symtext/symtext.py b/molsym/symtext/symtext.py index 26345b0..77da564 100755 --- a/molsym/symtext/symtext.py +++ b/molsym/symtext/symtext.py @@ -38,7 +38,11 @@ def __init__(self, mol, rotate_to_std, reverse_rotate, pg, symels, atom_map, mul 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 @@ def get_character_table(self): 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 @@ def contains_symmetric_irrep(self, dp_vector): 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. diff --git a/molsym/symtext/symtext_helper.py b/molsym/symtext/symtext_helper.py index 5da19a6..0098d18 100755 --- a/molsym/symtext/symtext_helper.py +++ b/molsym/symtext/symtext_helper.py @@ -214,7 +214,10 @@ def get_linear_atom_mapping(mol, pg): 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): """