add the option to generate xtal by blocks #155

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

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
+
+

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: