import os, numpy as np, time, shutil, subprocess, math, glob

CONFIG = { "LX": 60.0, "LY": 60.0, "LZ": 80.0, "PORE_RADIUS": 8.0, "TEMP": 300.0, "WALL_EPSILON": 0.07, "WALL_SIGMA": 3.55, "STEPS_EQ": 50000, "STEPS_PROD": 5000000, "RESTART_FREQ": 500000, "THERMO_FREQ": 5000, "DUMP_FREQ": 100000, "BASE_SEED": 12345, "OMP_NUM_THREADS": 4, "MEMBRANE_GAP": 10.0, "MIN_CC": 1.42, "MIN_OO": 2.80, "MIN_OC": 2.80, "MIN_HX": 1.50, "CELL_SIZE": 3.0, "MAX_FAIL": 800, "TARGET_WATER": 1200, "REPLICA": 1 }

WORK_DIR = "/kaggle/working/Nanopore_5pores" LOCAL_DIR = "/kaggle/working/sim" os.makedirs(WORK_DIR, exist_ok=True) os.makedirs(LOCAL_DIR, exist_ok=True) os.environ['OMP_NUM_THREADS'] = str(CONFIG['OMP_NUM_THREADS'])

def setup_lammps(): lmp_path = "/kaggle/working/lammps-static/bin/lmp" if not os.path.exists(lmp_path): print("📥 تحميل LAMMPS Static Binary...") subprocess.run([ "wget", "-q", "-O", "/kaggle/working/lammps.tar.gz", "https://download.lammps.org/static/lammps-linux-x86_64-30Mar2026.tar.gz" ], check=True) subprocess.run([ "tar", "-xzf", "/kaggle/working/lammps.tar.gz", "-C", "/kaggle/working/" ], check=True) os.chmod(lmp_path, 0o755) result = subprocess.run([lmp_path, "-h"], capture_output=True, text=True) if "LAMMPS" in result.stdout: print("✅ LAMMPS يعمل – يمكنك إيقاف نوت الـ CPU الآن") else: print("❌ فشل LAMMPS") print(result.stderr[:300]) raise RuntimeError("LAMMPS لم يعمل") return lmp_path

class NeighborGrid2D: def init(self, Lx, Ly, cs): self.Lx, self.Ly, self.cs = Lx, Ly, cs self.nx = max(1, int(Lx/cs)) self.ny = max(1, int(Ly/cs)) self.cells = {} def _idx(self, x, y): return (int((x % self.Lx) / self.Lx * self.nx) % self.nx, int((y % self.Ly) / self.Ly * self.ny) % self.ny) def insert(self, x, y): self.cells.setdefault(self._idx(x, y), []).append((x, y)) def any_within(self, x, y, cutoff): c2 = cutoff * cutoff ix0, iy0 = self._idx(x, y) for di in (-1, 0, 1): for dj in (-1, 0, 1): for (ax, ay) in self.cells.get(((ix0+di) % self.nx, (iy0+dj) % self.ny), []): dx = x - ax; dy = y - ay dx -= self.Lx * round(dx / self.Lx) dy -= self.Ly * round(dy / self.Ly) if dxdx + dydy < c2: return True return False

class NeighborGrid3D: def init(self, Lx, Ly, Lz, cs): self.Lx, self.Ly, self.Lz, self.cs = Lx, Ly, Lz, cs self.nx = max(1, int(Lx/cs)) self.ny = max(1, int(Ly/cs)) self.nz = max(1, int(Lz/cs)) self.cells = {} def _idx(self, x, y, z): return (int((x % self.Lx) / self.Lx * self.nx) % self.nx, int((y % self.Ly) / self.Ly * self.ny) % self.ny, int(((z + self.Lz0.5) / self.Lz) * self.nz) % self.nz) def insert(self, x, y, z): self.cells.setdefault(self._idx(x, y, z), []).append((x, y, z)) def any_within(self, x, y, z, cutoff): c2 = cutoff * cutoff ix0, iy0, iz0 = self._idx(x, y, z) for di in (-1, 0, 1): for dj in (-1, 0, 1): for dk in (-1, 0, 1): for (ax, ay, az) in self.cells.get(((ix0+di) % self.nx, (iy0+dj) % self.ny, (iz0+dk) % self.nz), []): dx = x-ax; dy = y-ay; dz = z-az dx -= self.Lxround(dx/self.Lx) dy -= self.Lyround(dy/self.Ly) dz -= self.Lzround(dz/self.Lz) if dxdx+dydy+dz*dz < c2: return True return False

class Builder: def init(self, cfg): self.cfg = cfg def graphene_pore(self): a, r, hg = 2.46, self.cfg['PORE_RADIUS'], self.cfg['MEMBRANE_GAP']/2.0 centers = [np.array([15,15]), np.array([15,45]), np.array([45,15]), np.array([45,45]), np.array([30,30])] atoms, grid = [], NeighborGrid2D(self.cfg['LX'], self.cfg['LY'], self.cfg['CELL_SIZE']) nx = int(self.cfg['LX']/a)+2 ny = int(self.cfg['LY']/(anp.sqrt(3)/2))+2 for i in range(-1, nx): for j in range(-1, ny): xA = ia+(j%2)a/2 yA = ja*np.sqrt(3)/2 for (x, y) in [(xA, yA), (xA, yA+a/2)]: x %= self.cfg['LX']; y %= self.cfg['LY'] if any(np.linalg.norm([x,y]-c) < r for c in centers): continue if grid.any_within(x, y, self.cfg['MIN_CC']): continue grid.insert(x, y) atoms.append((x, y, hg, 3)) atoms.append((x, y, -hg, 3)) print(f" ✅ كربون: {len(atoms)} ذرة") return atoms

def pack_water(self, carbon, seed):
    rng = np.random.RandomState(seed)
    Lx, Ly, Lz = self.cfg['LX'], self.cfg['LY'], self.cfg['LZ']
    hg = self.cfg['MEMBRANE_GAP']/2.0
    grid = NeighborGrid3D(Lx, Ly, Lz, self.cfg['CELL_SIZE'])
    for (cx, cy, cz, _) in carbon:
        grid.insert(cx, cy, cz)
    O, H1, H2 = [], [], []
    added, fail = 0, 0
    target = self.cfg['TARGET_WATER']
    while fail < self.cfg['MAX_FAIL'] and added < target:
        zlo, zhi = (hg+2.5, Lz/2-1.5) if rng.randint(2)==0 else (-Lz/2+1.5, -hg-2.5)
        ox = rng.rand()*Lx; oy = rng.rand()*Ly; oz = rng.uniform(zlo, zhi)
        if grid.any_within(ox, oy, oz, self.cfg['MIN_OC']): fail+=1; continue
        d = rng.randn(3); d /= max(np.linalg.norm(d), 1e-10)
        perp = np.cross(d, [0,0,1])
        if np.linalg.norm(perp) < 1e-6: perp = np.cross(d, [1,0,0])
        perp /= np.linalg.norm(perp)
        R_OH, th = 0.9572, 109.47*math.pi/180.0
        h1 = (ox+R_OH*d[0], oy+R_OH*d[1], oz+R_OH*d[2])
        h2d = math.cos(th)*d+math.sin(th)*perp
        h2 = (ox+R_OH*h2d[0], oy+R_OH*h2d[1], oz+R_OH*h2d[2])
        if (grid.any_within(h1[0],h1[1],h1[2],self.cfg['MIN_HX']) or
            grid.any_within(h2[0],h2[1],h2[2],self.cfg['MIN_HX'])): fail+=1; continue
        grid.insert(ox,oy,oz); grid.insert(*h1); grid.insert(*h2)
        O.append((ox,oy,oz)); H1.append(h1); H2.append(h2)
        added+=1; fail=0
        if added%200==0: print(f"   💧 {added}/{target}")
    print(f"   🎉 {len(O)} جزيء ماء")
    return O, H1, H2

class Guardian: def init(self, cfg): self.cfg = cfg def validate(self, carbon, O, H1, H2): Lx, Ly, Lz = self.cfg['LX'], self.cfg['LY'], self.cfg['LZ'] cs, hg = self.cfg['CELL_SIZE'], self.cfg['MEMBRANE_GAP']/2.0 for z in (hg, -hg): g2 = NeighborGrid2D(Lx, Ly, cs) for (cx, cy, cz, _) in carbon: if abs(cz-z) > 0.1: continue if g2.any_within(cx, cy, self.cfg['MIN_CC']): print("❌ تداخل C-C"); return False g2.insert(cx, cy) gc3 = NeighborGrid3D(Lx, Ly, Lz, cs) for (cx, cy, cz, _) in carbon: gc3.insert(cx, cy, cz) for (ox, oy, oz) in O: if gc3.any_within(ox, oy, oz, self.cfg['MIN_OC']): print("❌ تداخل O-C"); return False go = NeighborGrid3D(Lx, Ly, Lz, cs) for (ox, oy, oz) in O: if go.any_within(ox, oy, oz, self.cfg['MIN_OO']): print("❌ تداخل O-O"); return False go.insert(ox, oy, oz) print(" ✅ النظام آمن"); return True

def write_data(self, carbon, O, H1, H2, fname):
    nw, nc = len(O), len(carbon)
    Lx, Ly, Lz = self.cfg['LX'], self.cfg['LY'], self.cfg['LZ']
    with open(fname, 'w') as f:
        f.write(f"Nanopore 5pores\n{nw*3+nc} atoms\n{nw*2} bonds\n{nw} angles\n")
        f.write("3 atom types\n1 bond types\n1 angle types\n\n")
        f.write(f"0.0 {Lx:.6f} xlo xhi\n0.0 {Ly:.6f} ylo yhi\n")
        f.write(f"{-Lz/2:.6f} {Lz/2:.6f} zlo zhi\n\nMasses\n\n1 15.9994\n2 1.00794\n3 12.01\n\nAtoms\n\n")
        aid = 1
        for i in range(nw):
            o, h1, h2 = O[i], H1[i], H2[i]
            f.write(f"{aid} {i+1} 1 -0.8476 {o[0]:.4f} {o[1]:.4f} {o[2]:.4f}\n"); aid+=1
            f.write(f"{aid} {i+1} 2  0.4238 {h1[0]:.4f} {h1[1]:.4f} {h1[2]:.4f}\n"); aid+=1
            f.write(f"{aid} {i+1} 2  0.4238 {h2[0]:.4f} {h2[1]:.4f} {h2[2]:.4f}\n"); aid+=1
        for (cx, cy, cz, _) in carbon:
            f.write(f"{aid} {nw+1} 3 0.0 {cx:.4f} {cy:.4f} {cz:.4f}\n"); aid+=1
        f.write("\nBonds\n\n")
        for i in range(nw):
            b = i*3+1
            f.write(f"{i*2+1} 1 {b} {b+1}\n{i*2+2} 1 {b} {b+2}\n")
        f.write("\nAngles\n\n")
        for i in range(nw):
            b = i*3+1
            f.write(f"{i+1} 1 {b+1} {b} {b+2}\n")

def find_latest_restart(case_dir): files = glob.glob(os.path.join(case_dir, "restart.A")) +
glob.glob(os.path.join(case_dir, "restart.B")) return max(files, key=os.path.getmtime) if files else None

def write_lammps_input(data_path, restart_path, cfg, case_dir): in_file = os.path.join(LOCAL_DIR, "run.in") use_restart = restart_path and os.path.exists(restart_path) with open(in_file, 'w') as f: if use_restart: f.write(f"read_restart {restart_path}\n") f.write("pair_style lj/cut/coul/long 10.0\n") f.write("kspace_style pppm 1.0e-4\n") f.write("pair_modify mix arithmetic\n") f.write("group water type 1 2\n") f.write("group carbon type 3\n") f.write("fix freeze carbon setforce 0.0 0.0 0.0\n") f.write("fix shake_pr water shake 0.0001 20 0 b 1 a 1\n") f.write(f"fix nvt_pr water nvt temp {cfg['TEMP']} {cfg['TEMP']} 500.0\n") f.write("timestep 2.0\n") else: f.write(f"""units real atom_style full boundary p p p read_data {data_path}

pair_style lj/cut/coul/long 10.0 pair_coeff 1 1 0.1553 3.166 pair_coeff 2 2 0.0 0.0 pair_coeff 3 3 {cfg['WALL_EPSILON']} {cfg['WALL_SIGMA']} pair_coeff 1 3 0.104 3.36 pair_modify mix arithmetic

bond_style harmonic bond_coeff 1 1000.0 1.0 angle_style harmonic angle_coeff 1 100.0 109.47 kspace_style pppm 1.0e-4

group water type 1 2 group carbon type 3 fix freeze carbon setforce 0.0 0.0 0.0

fix shake_eq water shake 0.0001 20 0 b 1 a 1 velocity water create {cfg['TEMP']} {cfg['BASE_SEED']} dist gaussian fix nvt_eq water nvt temp {cfg['TEMP']} {cfg['TEMP']} 500.0 timestep 1.0 thermo {cfg['THERMO_FREQ']} run {cfg['STEPS_EQ']} unfix nvt_eq unfix shake_eq

fix shake_pr water shake 0.0001 20 0 b 1 a 1 fix nvt_pr water nvt temp {cfg['TEMP']} {cfg['TEMP']} 500.0 timestep 2.0 """) f.write(f"thermo {cfg['THERMO_FREQ']}\n") f.write("thermo_style custom step temp press density etotal\n") f.write(f"dump mydump water custom {cfg['DUMP_FREQ']} {os.path.join(LOCAL_DIR,'traj.lammpstrj')} id type x y z\n") restart_A = os.path.join(case_dir, "restart.A") restart_B = os.path.join(case_dir, "restart.B") f.write(f"restart {cfg['RESTART_FREQ']} {restart_A} {restart_B}\n") f.write(f"run {cfg['STEPS_PROD']} upto\n") f.write(f"write_restart {os.path.join(case_dir,'final.restart')}\n") return in_file

cfg = CONFIG case_dir = os.path.join(WORK_DIR, f"pore5_r{cfg['REPLICA']}") os.makedirs(case_dir, exist_ok=True)

lmp = setup_lammps()

restart_to_use = find_latest_restart(case_dir) if restart_to_use: print(f"♻️ استئناف من: {restart_to_use}") inp_file = write_lammps_input(None, restart_to_use, cfg, case_dir) else: print("🆕 بناء نظام جديد...") builder = Builder(cfg) carbon = builder.graphene_pore() O, H1, H2 = builder.pack_water(carbon, cfg['BASE_SEED']) guard = Guardian(cfg) if not guard.validate(carbon, O, H1, H2): raise RuntimeError("❌ فشل فحص النظام") data_file = os.path.join(LOCAL_DIR, "data.5pores") guard.write_data(carbon, O, H1, H2, data_file) inp_file = write_lammps_input(data_file, None, cfg, case_dir)

log_path = os.path.join(LOCAL_DIR, "log.txt") print("🚀 تشغيل LAMMPS...") t0 = time.time()

with open(log_path, 'w') as log_f: proc = subprocess.Popen( [lmp, "-in", inp_file], stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True ) for line in proc.stdout: print(line, end='', flush=True) log_f.write(line) proc.wait() ret = proc.returncode

dt = time.time() - t0 if ret != 0: print(f"❌ فشل LAMMPS بعد {dt:.0f} ثانية") else: print(f"✅ اكتملت في {dt/60:.1f} دقيقة") for fn in ["traj.lammpstrj", "log.txt"]: src = os.path.join(LOCAL_DIR, fn) dst = os.path.join(case_dir, fn) if os.path.exists(src): if os.path.exists(dst): os.remove(dst) shutil.move(src, dst) print(f"📁 النتائج: {case_dir}")# ============================================================

import os, numpy as np, time, shutil, subprocess, math, glob

CONFIG = { "LX": 60.0, "LY": 60.0, "LZ": 80.0, "PORE_RADIUS": 8.0, "TEMP": 300.0, "WALL_EPSILON": 0.07, "WALL_SIGMA": 3.55, "STEPS_EQ": 50000, "STEPS_PROD": 5000000, "RESTART_FREQ": 500000, "THERMO_FREQ": 5000, "DUMP_FREQ": 100000, "BASE_SEED": 12345, "OMP_NUM_THREADS": 4, "MEMBRANE_GAP": 10.0, "MIN_CC": 1.42, "MIN_OO": 2.80, "MIN_OC": 2.80, "MIN_HX": 1.50, "CELL_SIZE": 3.0, "MAX_FAIL": 800, "TARGET_WATER": 1200, "REPLICA": 1 }

WORK_DIR = "/kaggle/working/Nanopore_5pores" LOCAL_DIR = "/kaggle/working/sim" os.makedirs(WORK_DIR, exist_ok=True) os.makedirs(LOCAL_DIR, exist_ok=True) os.environ['OMP_NUM_THREADS'] = str(CONFIG['OMP_NUM_THREADS'])

def setup_lammps(): lmp_path = "/kaggle/working/lammps-static/bin/lmp" if not os.path.exists(lmp_path): print("📥 تحميل LAMMPS Static Binary...") subprocess.run([ "wget", "-q", "-O", "/kaggle/working/lammps.tar.gz", "https://download.lammps.org/static/lammps-linux-x86_64-30Mar2026.tar.gz" ], check=True) subprocess.run([ "tar", "-xzf", "/kaggle/working/lammps.tar.gz", "-C", "/kaggle/working/" ], check=True) os.chmod(lmp_path, 0o755) result = subprocess.run([lmp_path, "-h"], capture_output=True, text=True) if "LAMMPS" in result.stdout: print("✅ LAMMPS يعمل – يمكنك إيقاف نوت الـ CPU الآن") else: print("❌ فشل LAMMPS") print(result.stderr[:300]) raise RuntimeError("LAMMPS لم يعمل") return lmp_path

class NeighborGrid2D: def init(self, Lx, Ly, cs): self.Lx, self.Ly, self.cs = Lx, Ly, cs self.nx = max(1, int(Lx/cs)) self.ny = max(1, int(Ly/cs)) self.cells = {} def _idx(self, x, y): return (int((x % self.Lx) / self.Lx * self.nx) % self.nx, int((y % self.Ly) / self.Ly * self.ny) % self.ny) def insert(self, x, y): self.cells.setdefault(self._idx(x, y), []).append((x, y)) def any_within(self, x, y, cutoff): c2 = cutoff * cutoff ix0, iy0 = self._idx(x, y) for di in (-1, 0, 1): for dj in (-1, 0, 1): for (ax, ay) in self.cells.get(((ix0+di) % self.nx, (iy0+dj) % self.ny), []): dx = x - ax; dy = y - ay dx -= self.Lx * round(dx / self.Lx) dy -= self.Ly * round(dy / self.Ly) if dxdx + dydy < c2: return True return False

class NeighborGrid3D: def init(self, Lx, Ly, Lz, cs): self.Lx, self.Ly, self.Lz, self.cs = Lx, Ly, Lz, cs self.nx = max(1, int(Lx/cs)) self.ny = max(1, int(Ly/cs)) self.nz = max(1, int(Lz/cs)) self.cells = {} def _idx(self, x, y, z): return (int((x % self.Lx) / self.Lx * self.nx) % self.nx, int((y % self.Ly) / self.Ly * self.ny) % self.ny, int(((z + self.Lz0.5) / self.Lz) * self.nz) % self.nz) def insert(self, x, y, z): self.cells.setdefault(self._idx(x, y, z), []).append((x, y, z)) def any_within(self, x, y, z, cutoff): c2 = cutoff * cutoff ix0, iy0, iz0 = self._idx(x, y, z) for di in (-1, 0, 1): for dj in (-1, 0, 1): for dk in (-1, 0, 1): for (ax, ay, az) in self.cells.get(((ix0+di) % self.nx, (iy0+dj) % self.ny, (iz0+dk) % self.nz), []): dx = x-ax; dy = y-ay; dz = z-az dx -= self.Lxround(dx/self.Lx) dy -= self.Lyround(dy/self.Ly) dz -= self.Lzround(dz/self.Lz) if dxdx+dydy+dz*dz < c2: return True return False

class Builder: def init(self, cfg): self.cfg = cfg def graphene_pore(self): a, r, hg = 2.46, self.cfg['PORE_RADIUS'], self.cfg['MEMBRANE_GAP']/2.0 centers = [np.array([15,15]), np.array([15,45]), np.array([45,15]), np.array([45,45]), np.array([30,30])] atoms, grid = [], NeighborGrid2D(self.cfg['LX'], self.cfg['LY'], self.cfg['CELL_SIZE']) nx = int(self.cfg['LX']/a)+2 ny = int(self.cfg['LY']/(anp.sqrt(3)/2))+2 for i in range(-1, nx): for j in range(-1, ny): xA = ia+(j%2)a/2 yA = ja*np.sqrt(3)/2 for (x, y) in [(xA, yA), (xA, yA+a/2)]: x %= self.cfg['LX']; y %= self.cfg['LY'] if any(np.linalg.norm([x,y]-c) < r for c in centers): continue if grid.any_within(x, y, self.cfg['MIN_CC']): continue grid.insert(x, y) atoms.append((x, y, hg, 3)) atoms.append((x, y, -hg, 3)) print(f" ✅ كربون: {len(atoms)} ذرة") return atoms

def pack_water(self, carbon, seed):
    rng = np.random.RandomState(seed)
    Lx, Ly, Lz = self.cfg['LX'], self.cfg['LY'], self.cfg['LZ']
    hg = self.cfg['MEMBRANE_GAP']/2.0
    grid = NeighborGrid3D(Lx, Ly, Lz, self.cfg['CELL_SIZE'])
    for (cx, cy, cz, _) in carbon:
        grid.insert(cx, cy, cz)
    O, H1, H2 = [], [], []
    added, fail = 0, 0
    target = self.cfg['TARGET_WATER']
    while fail < self.cfg['MAX_FAIL'] and added < target:
        zlo, zhi = (hg+2.5, Lz/2-1.5) if rng.randint(2)==0 else (-Lz/2+1.5, -hg-2.5)
        ox = rng.rand()*Lx; oy = rng.rand()*Ly; oz = rng.uniform(zlo, zhi)
        if grid.any_within(ox, oy, oz, self.cfg['MIN_OC']): fail+=1; continue
        d = rng.randn(3); d /= max(np.linalg.norm(d), 1e-10)
        perp = np.cross(d, [0,0,1])
        if np.linalg.norm(perp) < 1e-6: perp = np.cross(d, [1,0,0])
        perp /= np.linalg.norm(perp)
        R_OH, th = 0.9572, 109.47*math.pi/180.0
        h1 = (ox+R_OH*d[0], oy+R_OH*d[1], oz+R_OH*d[2])
        h2d = math.cos(th)*d+math.sin(th)*perp
        h2 = (ox+R_OH*h2d[0], oy+R_OH*h2d[1], oz+R_OH*h2d[2])
        if (grid.any_within(h1[0],h1[1],h1[2],self.cfg['MIN_HX']) or
            grid.any_within(h2[0],h2[1],h2[2],self.cfg['MIN_HX'])): fail+=1; continue
        grid.insert(ox,oy,oz); grid.insert(*h1); grid.insert(*h2)
        O.append((ox,oy,oz)); H1.append(h1); H2.append(h2)
        added+=1; fail=0
        if added%200==0: print(f"   💧 {added}/{target}")
    print(f"   🎉 {len(O)} جزيء ماء")
    return O, H1, H2

class Guardian: def init(self, cfg): self.cfg = cfg def validate(self, carbon, O, H1, H2): Lx, Ly, Lz = self.cfg['LX'], self.cfg['LY'], self.cfg['LZ'] cs, hg = self.cfg['CELL_SIZE'], self.cfg['MEMBRANE_GAP']/2.0 for z in (hg, -hg): g2 = NeighborGrid2D(Lx, Ly, cs) for (cx, cy, cz, _) in carbon: if abs(cz-z) > 0.1: continue if g2.any_within(cx, cy, self.cfg['MIN_CC']): print("❌ تداخل C-C"); return False g2.insert(cx, cy) gc3 = NeighborGrid3D(Lx, Ly, Lz, cs) for (cx, cy, cz, _) in carbon: gc3.insert(cx, cy, cz) for (ox, oy, oz) in O: if gc3.any_within(ox, oy, oz, self.cfg['MIN_OC']): print("❌ تداخل O-C"); return False go = NeighborGrid3D(Lx, Ly, Lz, cs) for (ox, oy, oz) in O: if go.any_within(ox, oy, oz, self.cfg['MIN_OO']): print("❌ تداخل O-O"); return False go.insert(ox, oy, oz) print(" ✅ النظام آمن"); return True

def write_data(self, carbon, O, H1, H2, fname):
    nw, nc = len(O), len(carbon)
    Lx, Ly, Lz = self.cfg['LX'], self.cfg['LY'], self.cfg['LZ']
    with open(fname, 'w') as f:
        f.write(f"Nanopore 5pores\n{nw*3+nc} atoms\n{nw*2} bonds\n{nw} angles\n")
        f.write("3 atom types\n1 bond types\n1 angle types\n\n")
        f.write(f"0.0 {Lx:.6f} xlo xhi\n0.0 {Ly:.6f} ylo yhi\n")
        f.write(f"{-Lz/2:.6f} {Lz/2:.6f} zlo zhi\n\nMasses\n\n1 15.9994\n2 1.00794\n3 12.01\n\nAtoms\n\n")
        aid = 1
        for i in range(nw):
            o, h1, h2 = O[i], H1[i], H2[i]
            f.write(f"{aid} {i+1} 1 -0.8476 {o[0]:.4f} {o[1]:.4f} {o[2]:.4f}\n"); aid+=1
            f.write(f"{aid} {i+1} 2  0.4238 {h1[0]:.4f} {h1[1]:.4f} {h1[2]:.4f}\n"); aid+=1
            f.write(f"{aid} {i+1} 2  0.4238 {h2[0]:.4f} {h2[1]:.4f} {h2[2]:.4f}\n"); aid+=1
        for (cx, cy, cz, _) in carbon:
            f.write(f"{aid} {nw+1} 3 0.0 {cx:.4f} {cy:.4f} {cz:.4f}\n"); aid+=1
        f.write("\nBonds\n\n")
        for i in range(nw):
            b = i*3+1
            f.write(f"{i*2+1} 1 {b} {b+1}\n{i*2+2} 1 {b} {b+2}\n")
        f.write("\nAngles\n\n")
        for i in range(nw):
            b = i*3+1
            f.write(f"{i+1} 1 {b+1} {b} {b+2}\n")

def find_latest_restart(case_dir): files = glob.glob(os.path.join(case_dir, "restart.A")) +
glob.glob(os.path.join(case_dir, "restart.B")) return max(files, key=os.path.getmtime) if files else None

def write_lammps_input(data_path, restart_path, cfg, case_dir): in_file = os.path.join(LOCAL_DIR, "run.in") use_restart = restart_path and os.path.exists(restart_path) with open(in_file, 'w') as f: if use_restart: f.write(f"read_restart {restart_path}\n") f.write("pair_style lj/cut/coul/long 10.0\n") f.write("kspace_style pppm 1.0e-4\n") f.write("pair_modify mix arithmetic\n") f.write("group water type 1 2\n") f.write("group carbon type 3\n") f.write("fix freeze carbon setforce 0.0 0.0 0.0\n") f.write("fix shake_pr water shake 0.0001 20 0 b 1 a 1\n") f.write(f"fix nvt_pr water nvt temp {cfg['TEMP']} {cfg['TEMP']} 500.0\n") f.write("timestep 2.0\n") else: f.write(f"""units real atom_style full boundary p p p read_data {data_path}

pair_style lj/cut/coul/long 10.0 pair_coeff 1 1 0.1553 3.166 pair_coeff 2 2 0.0 0.0 pair_coeff 3 3 {cfg['WALL_EPSILON']} {cfg['WALL_SIGMA']} pair_coeff 1 3 0.104 3.36 pair_modify mix arithmetic

bond_style harmonic bond_coeff 1 1000.0 1.0 angle_style harmonic angle_coeff 1 100.0 109.47 kspace_style pppm 1.0e-4

group water type 1 2 group carbon type 3 fix freeze carbon setforce 0.0 0.0 0.0

fix shake_eq water shake 0.0001 20 0 b 1 a 1 velocity water create {cfg['TEMP']} {cfg['BASE_SEED']} dist gaussian fix nvt_eq water nvt temp {cfg['TEMP']} {cfg['TEMP']} 500.0 timestep 1.0 thermo {cfg['THERMO_FREQ']} run {cfg['STEPS_EQ']} unfix nvt_eq unfix shake_eq

fix shake_pr water shake 0.0001 20 0 b 1 a 1 fix nvt_pr water nvt temp {cfg['TEMP']} {cfg['TEMP']} 500.0 timestep 2.0 """) f.write(f"thermo {cfg['THERMO_FREQ']}\n") f.write("thermo_style custom step temp press density etotal\n") f.write(f"dump mydump water custom {cfg['DUMP_FREQ']} {os.path.join(LOCAL_DIR,'traj.lammpstrj')} id type x y z\n") restart_A = os.path.join(case_dir, "restart.A") restart_B = os.path.join(case_dir, "restart.B") f.write(f"restart {cfg['RESTART_FREQ']} {restart_A} {restart_B}\n") f.write(f"run {cfg['STEPS_PROD']} upto\n") f.write(f"write_restart {os.path.join(case_dir,'final.restart')}\n") return in_file

cfg = CONFIG case_dir = os.path.join(WORK_DIR, f"pore5_r{cfg['REPLICA']}") os.makedirs(case_dir, exist_ok=True)

lmp = setup_lammps()

restart_to_use = find_latest_restart(case_dir) if restart_to_use: print(f"♻️ استئناف من: {restart_to_use}") inp_file = write_lammps_input(None, restart_to_use, cfg, case_dir) else: print("🆕 بناء نظام جديد...") builder = Builder(cfg) carbon = builder.graphene_pore() O, H1, H2 = builder.pack_water(carbon, cfg['BASE_SEED']) guard = Guardian(cfg) if not guard.validate(carbon, O, H1, H2): raise RuntimeError("❌ فشل فحص النظام") data_file = os.path.join(LOCAL_DIR, "data.5pores") guard.write_data(carbon, O, H1, H2, data_file) inp_file = write_lammps_input(data_file, None, cfg, case_dir)

log_path = os.path.join(LOCAL_DIR, "log.txt") print("🚀 تشغيل LAMMPS...") t0 = time.time()

with open(log_path, 'w') as log_f: proc = subprocess.Popen( [lmp, "-in", inp_file], stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True ) for line in proc.stdout: print(line, end='', flush=True) log_f.write(line) proc.wait() ret = proc.returncode

dt = time.time() - t0 if ret != 0: print(f"❌ فشل LAMMPS بعد {dt:.0f} ثانية") else: print(f"✅ اكتملت في {dt/60:.1f} دقيقة") for fn in ["traj.lammpstrj", "log.txt"]: src = os.path.join(LOCAL_DIR, fn) dst = os.path.join(case_dir, fn) if os.path.exists(src): if os.path.exists(dst): os.remove(dst) shutil.move(src, dst) print(f"📁 النتائج: {case_dir}")