From 473f919e1a2cd83fe0bd642fcd9df150556ef600 Mon Sep 17 00:00:00 2001 From: Qiang Zhu Date: Tue, 5 Dec 2023 20:26:29 -0800 Subject: [PATCH] rewrite the representation class for mol.xtal --- pyxtal/__init__.py | 5 ++-- pyxtal/molecule.py | 8 +++--- pyxtal/representation.py | 11 ++++---- pyxtal/wyckoff_site.py | 58 +++++++++++++++++++++++++--------------- 4 files changed, 50 insertions(+), 32 deletions(-) diff --git a/pyxtal/__init__.py b/pyxtal/__init__.py index 1077acdb..feb5e2e9 100644 --- a/pyxtal/__init__.py +++ b/pyxtal/__init__.py @@ -1313,9 +1313,10 @@ def get_1D_representation(self): Get the 1D representation class for molecular crystals """ if self.molecular: - return representation.from_pyxtal(self) + rep = representation.from_pyxtal(self) else: - return representation_atom.from_pyxtal(self) + rep = representation_atom.from_pyxtal(self) + return rep def transform(self, trans, lattice=None): """ diff --git a/pyxtal/molecule.py b/pyxtal/molecule.py index 3f402035..195be1cf 100644 --- a/pyxtal/molecule.py +++ b/pyxtal/molecule.py @@ -850,7 +850,7 @@ def get_center(self, xyz, geometry=False): def get_principle_axes(self, xyz, rdmt=True): """ - get the principle axis for a rotated xyz + get the principle axis for a rotated xyz, sorted by the moments """ if self.smile is None or len(self.smile)==1 or not rdmt: Inertia = get_inertia_tensor(xyz) @@ -1374,9 +1374,9 @@ def get_orientations_in_wp(self, wp, rtol=1e-2): for i in list_i: orientations_new.append(orientations[i]) - #Check each of the found orientations for consistency with Wyckoff site. - #If consistent, put into an array of valid orientations - #print("======", orientations_new) + # Check each of the found orientations for consistency with Wyckoff site. + # If consistent, put into an array of valid orientations + # print("======", orientations_new) allowed = [] for o in orientations_new: op = o.get_op() diff --git a/pyxtal/representation.py b/pyxtal/representation.py index 3d6b002d..eaff05b5 100644 --- a/pyxtal/representation.py +++ b/pyxtal/representation.py @@ -262,7 +262,8 @@ def from_string(cls, inputs, smiles, composition=None): def to_standard_setting(self): xtal = self.to_pyxtal() - self.x = representation.from_pyxtal(xtal, standard=True).x + rep0 = representation.from_pyxtal(xtal, standard=True) + self.x = rep0.x def to_pyxtal(self, smiles=None, composition=None): """ @@ -328,13 +329,13 @@ def to_pyxtal(self, smiles=None, composition=None): dicts['PBC'] = [1, 1, 1] #dicts['number'] = number dicts['hn'] = struc.group.hall_number - dicts['index'] = 0 + dicts['index'] = v[0] dicts['lattice'] = struc.lattice.matrix dicts['lattice_type'] = ltype - dicts['center'] = v[:3] + dicts['center'] = v[1:4] if smile not in ["Cl-"]: - dicts['orientation'] = np.array(v[3:6]) - dicts['rotor'] = v[6:-1] + dicts['orientation'] = np.array(v[4:7]) + dicts['rotor'] = v[7:-1] dicts['reflect'] = int(v[-1]) site = mol_site.from_1D_dicts(dicts) site.type = i diff --git a/pyxtal/wyckoff_site.py b/pyxtal/wyckoff_site.py index 8561f83f..215836d3 100644 --- a/pyxtal/wyckoff_site.py +++ b/pyxtal/wyckoff_site.py @@ -335,19 +335,35 @@ def __init__(self, mol, position, orientation, wp, lattice=None, stype=0): else: self.lattice = Lattice.from_matrix(lattice) self.PBC = self.wp.PBC - self.mol = mol.mol # A Pymatgen molecule object + #self.mol = mol.mol # A Pymatgen molecule object self.symbols = mol.symbols #[site.specie.value for site in self.mol.sites] - self.numbers = self.mol.atomic_numbers + self.numbers = self.molecule.mol.atomic_numbers self.tols_matrix = mol.tols_matrix self.radius = mol.radius self.type = stype + def update_molecule(self, mol): + self.molecule = mol + self.numbers = mol.mol.atomic_numbers + self.symbols = mol.symbols + self.tols_matrix = mol.tols_matrix + self.radius = mol.radius + + def update_orientation(self, angles): + # QZ: Symmetrize the angle to the compatible orientation first + self.orientation.r = R.from_euler('zxy', angles, degrees=True) + self.orientation.matrix = self.orientation.r.as_matrix() + + def update_lattice(self, lattice): + # QZ: Symmetrize the angle to the compatible orientation first + self.lattice = lattice + def __str__(self): if not hasattr(self.wp, "site_symm"): self.wp.get_site_symmetry() self.angles = self.orientation.r.as_euler('zxy', degrees=True) - formula = self.mol.formula.replace(" ","") + formula = self.molecule.mol.formula.replace(" ","") s = "{:12s} @ [{:7.4f} {:7.4f} {:7.4f}] ".format(formula, *self.position) s += "WP [{:s}] ".format(self.wp.get_label()) s += "Site [{:}]".format(self.wp.site_symm.replace(" ","")) @@ -406,9 +422,9 @@ def encode(self): #if len(xyz)==3: print("encode: \n", self.molecule.mol.cart_coords) rotor = self.molecule.get_torsion_angles(xyz) ori, _, reflect = self.molecule.get_orientation(xyz) - return list(self.position) + list(ori) + rotor + [reflect] + return [self.wp.index] + list(self.position) + list(ori) + rotor + [reflect] else: - return list(self.position) + [0] + return [self.wp.index] + list(self.position) + [0] def to_1D_dicts(self): """ @@ -489,7 +505,7 @@ def _get_coords_and_species(self, absolute=False, PBC=False, first=False, unitce atomic coords: a numpy array of atomic coordinates in the site species: a list of atomic species for the atomic coords """ - coord0 = self.mol.cart_coords.dot(self.orientation.matrix.T) # + coord0 = self.molecule.mol.cart_coords.dot(self.orientation.matrix.T) wp_atomic_sites = [] wp_atomic_coords = None @@ -521,7 +537,6 @@ def _get_coords_and_species(self, absolute=False, PBC=False, first=False, unitce else: wp_atomic_coords = np.append(wp_atomic_coords, tmp, axis=0) wp_atomic_sites.extend(self.symbols) - if first: break @@ -600,7 +615,7 @@ def rotate(self, ax_id=0, ax_vector=None, angle=180): if ax_vector is not None: ax = ax_vector/np.linalg.norm(ax_vector) else: - xyz = self.mol.cart_coords.dot(p.as_matrix().T) + xyz = self.molecule.mol.cart_coords.dot(p.as_matrix().T) ax = self.molecule.get_principle_axes(xyz).T[ax_id] q = R.from_rotvec(ax*rad*angle) @@ -608,10 +623,6 @@ def rotate(self, ax_id=0, ax_vector=None, angle=180): self.orientation.r = o self.orientation.matrix = o.as_matrix() - def update_orientation(self, angles): - self.orientation.r = R.from_euler('zxy', angles, degrees=True) - self.orientation.matrix = self.orientation.r.as_matrix() - #def is_compatible_symmetry(self, tol=0.3): # """ # Check if the molecular symmetry matches the site symmetry @@ -639,7 +650,7 @@ def get_mol_object(self, id=0): Returns: a molecule object """ - coord0 = self.mol.cart_coords.dot(self.orientation.matrix.T) # + coord0 = self.molecule.mol.cart_coords.dot(self.orientation.matrix.T) # # Obtain the center in absolute coords if not hasattr(self.wp, "generators"): self.wp.set_generators() @@ -697,16 +708,16 @@ def update(self, coords, lattice=None, absolute=False, update_mol=True): mol = Molecule(self.symbols, coords-center) #match, _ = compare_mol_connectivity(mol, self.mol, True) - match, _ = compare_mol_connectivity(mol, self.mol) + match, _ = compare_mol_connectivity(mol, self.molecule.mol) if match: #position = np.mean(coords, axis=0).dot(self.lattice.inv_matrix) position = center.dot(self.lattice.inv_matrix) self.position = position - np.floor(position) if update_mol: self.orientation = Orientation(np.eye(3)) - self.mol = mol + self.molecule.mol = mol else: - m1 = pybel.readstring('xyz', self.mol.to('xyz')) + m1 = pybel.readstring('xyz', self.molecule.mol.to('xyz')) m2 = pybel.readstring('xyz', mol.to('xyz')) aligner = openbabel.OBAlign(True, False) aligner.SetRefMol(m1.OBMol) @@ -725,9 +736,9 @@ def update(self, coords, lattice=None, absolute=False, update_mol=True): else: import pickle with open('wrong.pkl', "wb") as f: - pickle.dump([mol, self.mol], f) + pickle.dump([mol, self.molecule.mol], f) mol.to(filename='Wrong.xyz', fmt='xyz') - self.mol.to(filename='Ref.xyz', fmt='xyz') + self.molecule.mol.to(filename='Ref.xyz', fmt='xyz') raise ValueError("molecular connectivity changes! Exit") # todo check if connectivty changed @@ -737,6 +748,10 @@ def _create_matrix(self, center=False, ignore=False): conditions. When multiplied with a set of points, generates additional points in cells adjacent to and diagonal to the original cell + Args: + center: + ignore: + Returns: A numpy array of matrices which can be multiplied by a set of coordinates @@ -784,6 +799,7 @@ def get_distances(self, coord1, coord2, m2=None, center=True, ignore=False): coord2: fractional coordinates of the reference neighbors m2: the length of reference molecule center: whether or not consider the self image for coord2 + ignore: Returns: distance matrix: [m1*m2*pbc, m1, m2] @@ -938,7 +954,7 @@ def get_neighbors_auto(self, factor=1.1, max_d=4.0, ignore_E=True, detail=False, min_ds: list of shortest distances neighs: list of neighboring molecular xyzs """ - mol_center = np.dot(self.position, self.lattice.matrix) + mol_center = np.dot(self.position-np.floor(self.position), self.lattice.matrix) numbers = self.molecule.mol.atomic_numbers coord1, _ = self._get_coords_and_species(first=True, unitcell=True) tm = Tol_matrix(prototype="vdW", factor=factor) @@ -976,7 +992,7 @@ def get_neighbors_auto(self, factor=1.1, max_d=4.0, ignore_E=True, detail=False, n1, n2 = numbers[ids[0][id]], numbers[ids[1][id]] if 1 not in [n1, n2]: pos = coord2[i][ids[1][id]] - mol_center - pairs.append((n2, pos)) + pairs.append((n2, pos))#; print('add self', i, n1, n2, pos, d[i][ids[0][id], ids[1][id]], id, np.linalg.norm(pos)) engs.append(eng[ids[0][id], ids[1][id]]) dists.append(d[i][ids[0][id], ids[1][id]]) else: @@ -1024,7 +1040,7 @@ def get_neighbors_auto(self, factor=1.1, max_d=4.0, ignore_E=True, detail=False, n1, n2 = numbers[ids[0][id]], numbers[ids[1][id]] if 1 not in [n1, n2]: pos = coord2[i][ids[1][id]] - mol_center - pairs.append((n2, pos)) + pairs.append((n2, pos))#; print('add other', i, n1, n2, pos, d[i][ids[0][id], ids[1][id]], np.linalg.norm(pos)) engs.append(eng[ids[0][id], ids[1][id]]) dists.append(d[i][ids[0][id], ids[1][id]])