From f30c5a05ecb966b6e17632298a5e9c98fff96ca2 Mon Sep 17 00:00:00 2001 From: qzhu2017 Date: Wed, 28 Jul 2021 16:35:10 -0700 Subject: [PATCH] add the option to generate xtal by blocks #155 --- pyxtal/__init__.py | 10 +-- pyxtal/block_crystal.py | 151 ++++++++++++++++++++++++++++++++++++ pyxtal/molecular_crystal.py | 6 +- 3 files changed, 159 insertions(+), 8 deletions(-) create mode 100644 pyxtal/block_crystal.py diff --git a/pyxtal/__init__.py b/pyxtal/__init__.py index bd3a1853..e4269e65 100644 --- a/pyxtal/__init__.py +++ b/pyxtal/__init__.py @@ -13,7 +13,7 @@ # PyXtal imports #avoid * from pyxtal.version import __version__ -from pyxtal.molecular_crystal import molecular_crystal +from pyxtal.block_crystal import block_crystal from pyxtal.crystal import random_crystal from pyxtal.symmetry import Group, Wyckoff_position, search_matched_position from pyxtal.operations import apply_ops, SymmOp, get_inverse @@ -235,16 +235,14 @@ def from_random( while True: count += 1 if self.molecular: - if numIons is None: - numIons = [len(Group(group)[0])]*len(species) - - struc = molecular_crystal(dim, + struc = block_crystal(dim, group, species, numIons, factor, thickness = thickness, area = area, + block = block, lattice = lattice, torsions = torsions, sites = sites, @@ -271,7 +269,7 @@ def from_random( break if count >= max_count: - raise RuntimeError("It takes long time to generate the structure, check inputs") + raise RuntimeError("long time to generate structure, check inputs") if quit: self.valid = struc.valid diff --git a/pyxtal/block_crystal.py b/pyxtal/block_crystal.py new file mode 100644 index 00000000..ce3e083f --- /dev/null +++ b/pyxtal/block_crystal.py @@ -0,0 +1,151 @@ +""" +Module for generating crystals based on building blocks +""" + +# Standard Libraries +import numpy as np + +# External Libraries +from pymatgen.core import Molecule + +# PyXtal imports +from pyxtal.molecular_crystal import molecular_crystal as mol_xtal +from pyxtal.molecule import pyxtal_molecule, compare_mol_connectivity, Orientation +from pyxtal.io import search_molecules_in_crystal +from pyxtal.wyckoff_site import mol_site + +def block_crystal(dim, group, molecules, numMols, factor, thickness, + area, block, lattice, torsions, sites, conventional, diag, tm): + + # If block is None, directly generate mol. xtal. + # Otherwise, generate crystal from building block + if block is None: + struc = mol_xtal(dim, + group, + molecules, + numMols, + factor, + thickness = thickness, + area = area, + lattice = lattice, + torsions = torsions, + sites = sites, + conventional = conventional, + diag = diag, + tm = tm) + return struc + + else: + p_mol = pyxtal_molecule(block) + block_mols = search_molecules_in_crystal(p_mol.mol, tol=0.2, once=False) + xtal_mols = [pyxtal_molecule(m, fix=True) for m in molecules] + + orders = [] + for m1 in xtal_mols: + for m2 in block_mols: + if len(m1.mol) == len(m2): + #match, mapping = compare_mol_connectivity(m2, m1.mol) + match, mapping = compare_mol_connectivity(m1.mol, m2) + #print(match, len(m1.mol), len(m2)) + if match: + orders.append([mapping[at] for at in range(len(m2))]) + break + #print(mapping) + #print("smile"); print(m1.mol.to('xyz')) + #print("reference"); print(m2.to('xyz')) + #print(orders[-1]) + #import sys; sys.exit() + + if len(orders) != len(molecules): + raise ValueError("Block is inconsistent with the molecules") + + # rearrange the order of block molecules + numbers = [] + coords = np.zeros([len(p_mol.mol),3]) + count = 0 + for order, m in zip(orders, block_mols): + numbers.extend([m.atomic_numbers[o] for o in order]) + coords[count:count+len(m)] += m.cart_coords[order] + count += len(m) + mol = Molecule(numbers, coords) + #print(mol.to('xyz')); import sys; sys.exit() + + for i in range(10): + struc = mol_xtal(dim, + group, + [mol], + None, + factor, + thickness = thickness, + area = area, + lattice = lattice, + torsions = torsions, + sites = sites, + conventional = conventional, + diag = diag, + tm = tm) + + if struc.valid: + break + + struc.molecules = xtal_mols + struc.numMols = struc.numMols * len(xtal_mols) + + # XYZ from the block_sites + b_site = struc.mol_sites[0] + xyz, _ = b_site._get_coords_and_species(absolute=True, first=True) + + # Create mol sites + ori = Orientation(np.eye(3)) + mol_sites = [] + + count = 0 + for i in range(len(molecules)): + mol = xtal_mols[i] + m_xyz = xyz[count:count+len(mol.mol)] + center = mol.get_center(m_xyz) + mol.reset_positions(m_xyz-center) + #print(mol.mol.to('xyz')) + position = np.dot(center, np.linalg.inv(struc.lattice.matrix)) + position -= np.floor(position) + m_site = mol_site(mol, position, ori, b_site.wp, struc.lattice, struc.diag) + mol_sites.append(m_site) + count += len(mol.mol) + + struc.mol_sites = mol_sites + + return struc + + +if __name__ == "__main__": + from pyxtal import pyxtal + from pyxtal.representation import representation + import pymatgen.analysis.structure_matcher as sm + + smiles = ["C1=C(C=C(C=C1[N+](=O)[O-])[N+](=O)[O-])C(=O)O.smi", + "CC1=CC2=C(C=C1)N3CC4=C(C=CC(=C4)C)N(C2)C3.smi"] + + for i in range(30): + s = pyxtal(molecular=True) + s.from_random(3, 14, smiles, block='xxv') + + rep2 = representation.from_pyxtal(s) + strs = rep2.to_string() + print(i, strs) + + rep3 = representation.from_string(strs, smiles) + s1 = rep3.to_pyxtal() + + pmg1 = s.to_pymatgen() + pmg2 = s1.to_pymatgen() + pmg1.remove_species("H") + pmg2.remove_species("H") + if not sm.StructureMatcher().fit(pmg1, pmg2): + print("Mismatch") + print("S1", s.check_short_distances()) + print("S2", s1.check_short_distances()) + s.to_file('1.cif') + s1.to_file('2.cif') + break + + diff --git a/pyxtal/molecular_crystal.py b/pyxtal/molecular_crystal.py index c3ca1f61..023c3ea6 100644 --- a/pyxtal/molecular_crystal.py +++ b/pyxtal/molecular_crystal.py @@ -63,7 +63,6 @@ def __init__( factor = 1.1, thickness = None, area = None, - block = None, select_high = True, allow_inversion = True, orientations = None, @@ -110,7 +109,10 @@ def __init__( self.symbol = self.group.symbol # Composition - numMols = np.array(numMols) # must convert it to np.array + if numMols is None: + numMols = [len(self.group[0])] * len(molecules) + else: + numMols = np.array(numMols) # must convert it to np.array if not conventional: mul = cellsize(self.group) else: