Added support for customizing bulk regions and improved documentation

Author: jpqiu

docs/structures/README.md

@@ -30,4 +30,52 @@ stc.find_surf_idx(element='Ti', tolerance=1, dsur='dw', check_cross_boundary=Tru
 [4, 9, 34, 39, 64, 69, 94, 99, 124, 129, 154, 159, 184, 189, 214, 219]
 # bottom surface
 [0, 5, 30, 35, 60, 65, 90, 95, 120, 125, 150, 155, 180, 185, 210, 215]
-```
+```
+
+
+## Generate an Interface Containing Ions or Organic Molecules in Distinct Regions
+
+In many interfacial simulations, it is necessary to include ions or organic molecules—especially those that exhibit chemical adsorption. The generate_components_symmetry_box() method allows precise control over the spatial distribution of different species in the system.
+
+For example, the image below shows an interface model constructed with ions and organics distributed in adsorption and bulk regions:
+![components](./figures/slab-components.png)
+
+
+```python
+#
+from ase.build import fcc111
+from ase.io import write
+from ectoolkits.structures.slab import Slab
+
+# Build a 6×6×6 Cu(111) slab with orthogonal unit cell and 1.5 Å vacuum
+slab = fcc111('Cu', size=(6, 6, 6), a=3.57, vacuum=1.5, orthogonal=True)
+slab = Slab(slab)
+
+# Define water box height
+water_box_len = 30
+
+# Generate symmetric water box with:
+# - [Cl-] in the adsorption region (1 molecule)
+# - (CCOCCO), [K+], and [Cl-] in the bulk region (1, 3, and 3 molecules respectively)
+# Adsorption zone is 0.5, bulk region spans from 6 Å (bottom) to 10 Å (top)
+slab.generate_components_symmetry_box(
+    water_box_len,
+    ads_components_smiles=['[Cl-]'], ads_components_counts=[1],
+    bulk_components_smiles=['CCOCCO', '[K+]', '[Cl-]'], bulk_components_counts=[1, 3, 3],
+    ads_limit=0.5, bulk_dw=6, bulk_up=10
+)
+
+# Identify top and bottom surface atom indices
+top_surface_idx = slab.find_surf_idx(element='Cu', tolerance=0.1, dsur='up')
+bottom_surface_idx = slab.find_surf_idx(element='Cu', tolerance=0.1, dsur='dw')
+
+# Generate the final interface structure
+interface = slab.generate_interface(
+    water_box_len=water_box_len,
+    top_surface_idx=top_surface_idx,
+    bottom_surface_idx=bottom_surface_idx
+)
+
+# Write to XYZ file
+write('coord.xyz', interface)
+```

docs/structures/figures/slab-components.png

@@ -351,29 +351,40 @@ def generate_components_symmetry_box(self,
                                             ads_components_counts: list,
                                             bulk_components_smiles: list,
                                             bulk_components_counts: list,
-                                            ads_limit = 0.5
+                                            ads_limit = 0.5,
+                                            bulk_dw= None,
+                                            bulk_up= None
                                             ):
         """
-        Generate water + ion/organic components box in a slab-aligned box
+		Generate a slab-aligned box with symmetric distribution of water, adsorbed molecules, and bulk components.
+
         Args:
             water_box_len (float): Height of the water box along z-direction.
             ads_components_smiles (list): SMILES list of adsorbed components.
             ads_components_counts (list): Corresponding counts of adsorbed components.
             bulk_components_smiles (list): SMILES list of bulk components.
             bulk_components_counts (list): Corresponding counts of bulk components.
-            ads_limit (float): Maximum z-direction distance for adsorbates from the surface.  
-                For large molecules, this value may need to be increased. (default = 0.5 Å)
+            ads_limit (float): Vertical range (in Å) above and below the surface for placing adsorbed species. Default = 0.5 Å.
+            bulk_dw (float, optional): Lower boundary (Å) of the symmetric bulk region (below box center).
+            bulk_up (float, optional): Upper boundary (Å) of the symmetric bulk region (below box center).
+                Must satisfy: 0 < bulk_dw < bulk_up < water_box_len / 2.
+                If set, creates symmetric regions for bulk placement.
+
 
         Returns:
             None. Generates a Packmol input file and runs Packmol to generate `watbox.xyz`.
             
         Example:
             >>> from ectoolkits.structures.tools import generate_components_symmetry_box
             >>> generate_components_symmetry_box(
+            ...     water_box_len = 20
             ...     ads_components_smiles=["[Cl-]"],
             ...     ads_components_counts=[2],
             ...     bulk_components_smiles=["[K+]", "[Cl-]"],
-            ...     bulk_components_counts=[3, 3]  
+            ...     bulk_components_counts=[3, 3],
+            ...     ads_limit=0.5,
+            ...     bulk_dw = 3,
+            ...     bulk_up = 6
             )
 
         
@@ -383,31 +394,40 @@ def generate_components_symmetry_box(self,
 
         assert len(ads_components_smiles) == len(ads_components_counts), "Length mismatch in adsorbed components"
         assert len(bulk_components_smiles) == len(bulk_components_counts), "Length mismatch in bulk components" 
+        custom_bulk = False
+        if bulk_dw is not None and bulk_up is not None:
+            if bulk_dw > 0 and bulk_up < water_box_len / 2:
+                custom_bulk = True
+            else:
+                raise ValueError(f"Invalid custom bulk region: bulk_dw={bulk_dw}, bulk_up={bulk_up}. "
+                 f"Ensure 0 < bulk_dw < bulk_up < water_box_len/2.")
+
         cell = self.get_cell()
         cell_a = cell[0]
         cell_b = cell[1]
+
+		# === Water molecule estimation ===
         header = "-"
         logger.info(header * 50)
         logger.info("Now Generate Water Box")
         space_per_water = 9.86 ** 3 / 32
-        init_wat_num = (np.linalg.norm(np.cross(cell_a, cell_b))
-            * water_box_len) / space_per_water
-        init_wat_num = int(init_wat_num)
+        init_wat_num = int((np.linalg.norm(np.cross(cell_a, cell_b))
+            * water_box_len) / space_per_water)
 
         exclude_water = get_exclude_water_count(
             ads_components_smiles, ads_components_counts,
             bulk_components_smiles, bulk_components_counts
         )
 
         wat_num = max(0, init_wat_num - 2 * exclude_water)
-        wat_num = int(wat_num)
         logger.info(f"final number: {wat_num}")
         logger.info("Read Water Box Length: {0:03f} A".format(water_box_len))
         logger.info("Predict Water Number: {0}".format(wat_num))
 
         n_vec_a, d1_a, d2_a, n_vec_b, d1_b, d2_b = get_plane_eq(cell_a, cell_b)
         logger.info("Calculate Plane Equation")
 
+        # === Create directory and input ===
         if os.path.exists('gen_water'):
             logger.info("found gen_water direcotry, now remove it")
             shutil.rmtree('gen_water')
@@ -422,24 +442,46 @@ def generate_components_symmetry_box(self,
                 "filetype xyz",
                 "output watbox.xyz\n"
             ]
+            # === Water ===
             water_xyz = os.path.join("gen_water", "water.xyz")
             smiles2xyz('O',water_xyz)
             lines.append(get_packmol_inp("water.xyz", wat_num, n_vec_a, n_vec_b, d1_a, d2_a, d1_b, d2_b, 0, water_box_len, radius=1))
 
+            # === Adsorbates ===
             for i, (smiles, count) in enumerate(zip(ads_components_smiles, ads_components_counts)):
                 xyz_file = f"ads_component_{i}.xyz"
                 xyz_path = os.path.join("gen_water", xyz_file)
                 smiles2xyz(smiles, xyz_path)
+                # bottom surface
                 lines.append(get_packmol_inp(xyz_file, count, n_vec_a, n_vec_b, d1_a, d2_a, d1_b, d2_b, 0, ads_limit, radius=2))
+                # top surface
                 lines.append(get_packmol_inp(xyz_file, count, n_vec_a, n_vec_b, d1_a, d2_a, d1_b, d2_b, water_box_len - ads_limit, water_box_len, radius=2))
 
+            # === Bulk Components ===
             for i, (smiles, count) in enumerate(zip(bulk_components_smiles, bulk_components_counts)):
                 xyz_file = f"bulk_component_{i}.xyz"
                 xyz_path = os.path.join("gen_water", xyz_file)
                 smiles2xyz(smiles, xyz_path)
-                lines.append(get_packmol_inp(xyz_file, count, n_vec_a, n_vec_b, d1_a, d2_a, d1_b, d2_b, water_box_len / 4, water_box_len / 2, radius=2))
-                lines.append(get_packmol_inp(xyz_file, count, n_vec_a, n_vec_b, d1_a, d2_a, d1_b, d2_b, water_box_len / 2, 3 * water_box_len / 4, radius=2))
+                if custom_bulk:
+                    # lower layer
+                    lines.append(get_packmol_inp(xyz_file, count, n_vec_a, n_vec_b,
+                                                d1_a, d2_a, d1_b, d2_b,
+                                                bulk_dw, bulk_up, radius=2))
+                    # symmetric upper layer
+                    lines.append(get_packmol_inp(xyz_file, count, n_vec_a, n_vec_b,
+                                                d1_a, d2_a, d1_b, d2_b,
+                                                water_box_len - bulk_up, water_box_len - bulk_dw, radius=2))
+                else:
+                    # default symmetric placement in central box
+                    lines.append(get_packmol_inp(xyz_file, count, n_vec_a, n_vec_b,
+                                                d1_a, d2_a, d1_b, d2_b,
+                                                water_box_len / 4, water_box_len / 2, radius=2))
+                    lines.append(get_packmol_inp(xyz_file, count, n_vec_a, n_vec_b,
+                                                d1_a, d2_a, d1_b, d2_b,
+                                                water_box_len / 2, 3 * water_box_len / 4, radius=2))
+        
             f.write("\n".join(lines))
+        # === Run Packmol ===
         os.chdir("./gen_water")
         os.system("packmol < gen_wat_box.inp")
         os.chdir("../")