add search cp-plane function #196

doc/conf.py

@@ -39,6 +39,7 @@
 # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
 # ones.
 extensions = [
+    "sphinx_rtd_theme",
     "sphinx.ext.autodoc",
     "sphinx.ext.todo",
     "sphinx.ext.coverage",

pyxtal/__init__.py

@@ -1140,15 +1140,23 @@ def get_num_torsions(self):
             N_torsions += len(s.molecule.torsionlist)
         return N_torsions
 
-    def to_ase(self, resort=True):
+    def to_ase(self, resort=True, center_only=False):
         """
         Export to ase Atoms object.
         """
         if self.valid:
             if self.dim > 0:
                 lattice = self.lattice.copy()
                 if self.molecular:
-                    coords, species = self._get_coords_and_species(True)
+                    if center_only:
+                        coords, species = [], []
+                        for site in self.mol_sites:
+                            _coords = site.wp.apply_ops(site.position)
+                            _coords -= np.floor(_coords)
+                            coords.extend(_coords.dot(self.lattice.matrix))
+                            species.extend([site.type+1] * site.wp.multiplicity)
+                    else:
+                        coords, species = self._get_coords_and_species(True)
                     latt, coords = lattice.add_vacuum(coords, frac=False, PBC=self.PBC)
                     atoms = Atoms(species, positions=coords, cell=latt, pbc=self.PBC)
                 else:

pyxtal/plane.py

@@ -0,0 +1,202 @@
+"""
+crystal plane class
+"""
+from pyxtal import pyxtal
+import numpy as np
+from math import gcd
+import itertools
+
+def has_reduction(hkl):
+    h, k, l = hkl
+    gcf = gcd(h, gcd(k, l))
+    if gcf > 1:
+        # like [2, 2, 0]
+        return True
+    elif hkl[np.nonzero(np.array(hkl))[0][0]] < 0:
+        # like [-2, 0, 0]
+        return True
+    return False
+
+def reduced_hkl(hkl):
+    h, k, l = hkl
+    gcf = gcd(h, gcd(k, l))
+    if gcf > 1:
+        return [int(h/gcf), int(k/gcf), int(l/gcf)], gcf
+    else:
+        return [h, k, l], 1
+
+def get_structure_factor(atoms, hkl):
+    return structure_factor(atoms.get_scaled_positions(), hkl)
+
+def structure_factor(pos, hkl, total=True):
+    coords = np.dot(pos, hkl)
+    F = np.exp(-2*np.pi*(1j)*coords)
+    if total:
+        return F.sum()
+    else:
+        return F
+
+def get_dspacing(inv_matrix, hkl):
+    return 1/np.linalg.norm(inv_matrix.dot(np.array(hkl)))
+
+
+
+class planes():
+    """
+    This is a database class to process crystal data
+
+    Args:
+        db_name: *.db format from ase database
+    """
+
+    def __init__(self, extent=6, d_min=1.0):
+        self.d_min = d_min
+        self.extent = extent
+        self.set_planes()
+        #self.set_xtal()
+
+    def set_xtal(self, xtal):
+        self.xtal = xtal
+        self.atoms = xtal.to_ase(center_only=True)
+        self.cell_reciprocal = xtal.lattice.inv_matrix
+
+    def set_planes(self):
+        planes = list(itertools.product(range(-self.extent, self.extent+1), repeat=3))
+        planes = [hkl for hkl in planes if hkl != (0, 0, 0) and not has_reduction(hkl)]
+        self.planes = planes
+
+    def search_close_packing_planes(self, N_max=10):
+        """
+        Search for the close-packed molecular plane for a given crystal
+        """
+        cp_planes = []
+        for hkl in self.planes:
+            dspacing = get_dspacing(self.cell_reciprocal, hkl)
+            for n in range(1, N_max):
+                if dspacing/n > self.d_min:
+                    hkl1 = n*np.array(hkl)
+                    F = get_structure_factor(self.atoms, hkl1)
+                    if np.abs(F) >= len(self.atoms) * 0.5: 
+                        # Scan the plane with high density
+                        plane = self.get_separation(hkl1)
+                        if plane is not None:
+                            cp_planes.append(plane)
+                        break
+                else:
+                    # early termination for very short ranged hkl
+                    break
+        if len(cp_planes) > 0:
+            cp_planes = sorted(cp_planes, key=lambda x: -x[-1][0])
+        return cp_planes
+
+    def get_separation(self, hkl, tol=-0.1, slab_factor=1.5):
+        """
+        Compute the separation for the given hkl plane
+    
+        Args:
+            - hkl: three indices
+            - slab factor
+        """
+        hkl_reduced, hkl_factor = reduced_hkl(hkl)
+        d_spacing = get_dspacing(self.cell_reciprocal, hkl_reduced)
+
+        # Go over each molecule to compute the molecular slab
+        slabs = []
+        for mol_site in self.xtal.mol_sites:
+            N_atoms = len(mol_site.numbers)
+            center_frac = mol_site.position # fractional
+            xyz, species = mol_site._get_coords_and_species(unitcell=True)
+            for id_mol, op in enumerate(mol_site.wp.ops):
+                start, end = id_mol*N_atoms, (id_mol+1)*N_atoms
+                coords = xyz[start:end, :] # frac
+                center_frac = op.operate(mol_site.position)
+                center_frac -= np.floor(center_frac)
+                coords -= center_frac # place the center to (0, 0, 0)
+                center_hkl = np.dot(center_frac, hkl_reduced)
+                center_hkl -= np.floor(center_hkl)
+                coords_hkl = np.dot(coords, hkl_reduced)
+
+                # Terminate only if the molecular slab width is s
+                width = coords_hkl.max() - coords_hkl.min()
+                if width > slab_factor/hkl_factor:
+                    return None
+                lower, upper = center_hkl+coords_hkl.min(), center_hkl+coords_hkl.max()
+                slabs.append([center_hkl, lower, upper])
+
+        groups = self.group_slabs(slabs, 0.5/hkl_factor)
+        separations = self.find_unique_separations(groups, d_spacing)
+        return (hkl, d_spacing, separations)
+
+    def group_slabs(self, slabs, tol):
+        groups = [] 
+        for slab in slabs:
+            new = True
+            center, lower, upper = slab
+            for group in groups:
+                if group[1] <= center <= group[2]:
+                    new = False
+                else:
+                    dist = center-group[0]
+                    shift = np.round(dist)
+                    if abs(dist-shift) < tol:
+                        new = False
+                        lower -= shift
+                        upper -= shift
+
+                if not new:
+                    if lower < group[1]:
+                        group[1] = lower
+                    if upper > group[2]:
+                        group[2] = upper
+                    center = (group[1] + group[2])/2
+                    break
+            if new:
+                groups.append([center, lower, upper])
+
+        return sorted(groups)
+
+    def find_unique_separations(self, groups, d):
+        groups.append(groups[0])
+        groups[-1] = [g+1 for g in groups[-1]]
+        separations = []
+        for i in range(len(groups)-1):
+            separations.append(d*(groups[i+1][1] - groups[i][2])) 
+        separations = np.unique(-1*np.array(separations).round(decimals=3))
+        return -np.sort(separations)
+
+    def gather(self, planes):
+        for _plane in planes:
+            (hkl, _, separations) = _plane
+            for separation in separations:
+                p = plane(hkl, self.cell_reciprocal, separation)
+                print(p)
+
+class plane():
+    """
+    This simplest possible plane object 
+    """
+
+    def __init__(self, hkl, cell_reciprocal, separation=None):
+        self.indices = hkl
+        self.cell_reciprocal = cell_reciprocal
+        self.simple_indices, self.factor = reduced_hkl(hkl)
+        self.d_spacing = get_dspacing(self.cell_reciprocal, self.indices)
+        self.separation = separation
+
+    def __str__(self):
+        s = "({:2d} {:2d} {:2d})".format(*self.indices)
+        s += " Spacing: {:6.3f}".format(self.d_spacing)
+        s += " Separation: {:6.3f}".format(self.separation)
+
+        return s
+
+if __name__ == "__main__":
+    from pyxtal.db import database
+
+    db = database('pyxtal/database/mech.db') 
+    for code in db.codes[:2]:
+        print('\n', code)
+        p = planes()
+        p.set_xtal(db.get_pyxtal(code))
+        cp_planes = p.search_close_packing_planes()
+        p.gather(cp_planes)

pyxtal/symmetry.py

@@ -1416,7 +1416,7 @@ def from_symops_wo_group(ops):
         _, spg_num = get_symmetry_from_ops(ops)
         wp = Wyckoff_position.from_group_and_index(spg_num, 0)
         if isinstance(ops[0], str):
-            ops = [SymmOp.from_xyz_string(op) for op in ops]
+            ops = [SymmOp.from_xyz_str(op) for op in ops]
         wp.ops = ops
         match_spg, match_hm = wp.update()
         #print("match_spg", match_spg, "match_hall", match_hm)
@@ -1437,7 +1437,7 @@ def from_symops(ops, G):
 
         """
         if isinstance(ops[0], str):
-            ops = [SymmOp.from_xyz_string(op) for op in ops]
+            ops = [SymmOp.from_xyz_str(op) for op in ops]
 
         for wp in G:
             if wp.has_equivalent_ops(ops):
@@ -1699,12 +1699,12 @@ def process_ops(self):
         opss = [self.ops]
         if self.number in [5, 12] and self.index > 0:
             # replace y with y+1/2
-            op2 = SymmOp.from_xyz_string('x, y+1/2, z')
+            op2 = SymmOp.from_xyz_str('x, y+1/2, z')
             ops = [op2*op for op in self.ops]
             opss.append(ops)
 
         if self.number in [13] and self.index > 0:
-            op2 = SymmOp.from_xyz_string('x, -y, z')
+            op2 = SymmOp.from_xyz_str('x, -y, z')
             ops = [op2*op for op in self.ops]
             opss.append(ops)
 
@@ -2559,7 +2559,7 @@ def letter_from_index(index, group, dim=3):
     elif dim == 0:
         checko = True
     if checko is True:
-        if len(group[-1]) == 1 and group[-1][0] == SymmOp.from_xyz_string("0,0,0"):
+        if len(group[-1]) == 1 and group[-1][0] == SymmOp.from_xyz_str("0,0,0"):
             # o comes before a
             letters1 = "o" + letters
     length = len(group)
@@ -2586,7 +2586,7 @@ def index_from_letter(letter, group, dim=3):
     elif dim == 0:
         checko = True
     if checko is True:
-        if len(group[-1]) == 1 and group[-1][0] == SymmOp.from_xyz_string("0,0,0"):
+        if len(group[-1]) == 1 and group[-1][0] == SymmOp.from_xyz_str("0,0,0"):
             # o comes before a
             letters1 = "o" + letters
     length = len(group)
@@ -3066,7 +3066,7 @@ def get_wyckoffs(num, organized=False, dim=3):
             if dim == 0:
                 wyckoffs[-1].append(SymmOp(y))
             else:
-                wyckoffs[-1].append(SymmOp.from_xyz_string(y))
+                wyckoffs[-1].append(SymmOp.from_xyz_str(y))
     if organized:
         wyckoffs_organized = [[]]  # 2D Array of WP's organized by multiplicity
         old = len(wyckoffs[0])
@@ -3117,7 +3117,7 @@ def get_wyckoff_symmetry(num, dim=3):
                 if dim == 0:
                     op = SymmOp(z)
                 else:
-                    op = SymmOp.from_xyz_string(z)
+                    op = SymmOp.from_xyz_str(z)
                 symmetry[-1][-1].append(op)
     return symmetry
 
@@ -3156,7 +3156,7 @@ def get_generators(num, dim=3):
         # Loop over ops
         for y in x:
             if dim > 0:
-                op = SymmOp.from_xyz_string(y)
+                op = SymmOp.from_xyz_str(y)
             else:
                 op = SymmOp(y)
             generators[-1].append(op)
@@ -3738,7 +3738,7 @@ def get_symmetry_from_ops(ops, tol=1e-5):
     from spglib import get_hall_number_from_symmetry
 
     if isinstance(ops[0], str):
-        ops = [SymmOp.from_xyz_string(op) for op in ops]
+        ops = [SymmOp.from_xyz_str(op) for op in ops]
     rot = [op.rotation_matrix for op in ops]
     tran = [op.translation_vector for op in ops]
     hall_number = get_hall_number_from_symmetry(rot, tran, tol)