Here is my python script
coding: utf-8
“”"
generate_atomic.py
Génère un fichier LAMMPS au format ‘atom_style atomic’ (positions + types + masses).
Dépendances: ase, numpy
Usage: python generate_atomic.py
Place dans le même dossier:
- epon862.sdf
- detda.sdf
“”"
import random
import numpy as np
from ase.build import nanotube
from ase import Atoms
from ase.io import read, write
import sys
---------------- Parameters (edit if needed) ----------------
SEED = 42
n, m = 8, 8
target_length = 51.83
unit_length = 2.49
NUM_EPON = 192
NUM_DETDA = 96
BOX_X = 50.0
BOX_Y = 50.0
BOX_Z = 53.0
MIN_DIST = 1.4 # distance minimale entre atomes (Å)
MAX_ATTEMPTS = 2000
OUTFILE = “resine_cnt.atomic.data”
ATOM_STYLE = “atomic” # atomic → pas de charge / pas de molecule id
---------------- Init ----------------
random.seed(SEED)
np.random.seed(SEED)
Construire CNT
N = round(target_length / unit_length)
cnt = nanotube(n, m, length=N)
print(f"[CNT] atoms: {len(cnt)} | cell z: {cnt.cell[2,2]:.4f} Å")
Lire molécules une seule fois
try:
epon_ref = read(“epon862.sdf”)
detda_ref = read(“detda.sdf”)
except Exception as e:
print(“Erreur: impossible de lire epon862.sdf / detda.sdf”, e)
sys.exit(1)
def copy_and_shift_to_origin(mol):
mol2 = mol.copy()
mol2.translate(-mol2.get_center_of_mass())
return mol2
def random_translation():
return np.array([random.uniform(0, BOX_X),
random.uniform(0, BOX_Y),
random.uniform(0, BOX_Z)])
def too_close(existing_pos, cand_pos, min_dist):
if existing_pos.shape[0] == 0:
return False
d = existing_pos[:, None, :] - cand_pos[None, :, :]
d2 = np.sum(dd, axis=2)
return np.any(d2 < (min_distmin_dist))
packing simple
resin = Atoms()
positions_existing = np.empty((0, 3))
EPON
print(“[Placement] EPON”)
for i in range(NUM_EPON):
mol = copy_and_shift_to_origin(epon_ref)
placed = False
attempts = 0
while not placed and attempts < MAX_ATTEMPTS:
attempts += 1
t = random_translation()
mol_t = mol.copy()
mol_t.translate(t)
pos = mol_t.get_positions()
if not too_close(positions_existing, pos, MIN_DIST):
resin += mol_t
positions_existing = np.vstack((positions_existing, pos))
placed = True
if not placed:
print(f"Warning: EPON #{i+1} not placed after {MAX_ATTEMPTS} attempts. Adding at center (may overlap).")
mol.translate([BOX_X/2, BOX_Y/2, BOX_Z/2])
resin += mol
positions_existing = np.vstack((positions_existing, mol.get_positions()))
DETDA
print(“[Placement] DETDA”)
for i in range(NUM_DETDA):
mol = copy_and_shift_to_origin(detda_ref)
placed = False
attempts = 0
while not placed and attempts < MAX_ATTEMPTS:
attempts += 1
t = random_translation()
mol_t = mol.copy()
mol_t.translate(t)
pos = mol_t.get_positions()
if not too_close(positions_existing, pos, MIN_DIST):
resin += mol_t
positions_existing = np.vstack((positions_existing, pos))
placed = True
if not placed:
print(f"Warning: DETDA #{i+1} not placed after {MAX_ATTEMPTS} attempts. Adding at center (may overlap).")
mol.translate([BOX_X/2, BOX_Y/2, BOX_Z/2])
resin += mol
positions_existing = np.vstack((positions_existing, mol.get_positions()))
print(f"[Resin] atoms: {len(resin)}")
Center CNT and combine
cnt_width = cnt.cell[0,0] if cnt.cell is not None else max(cnt.get_positions()[:,0]) - min(cnt.get_positions()[:,0])
cnt_length = cnt.cell[2,2] if cnt.cell is not None else max(cnt.get_positions()[:,2]) - min(cnt.get_positions()[:,2])
cnt.translate([(BOX_X - cnt_width)/2.0, BOX_Y/2.0, (BOX_Z - cnt_length)/2.0])
system = resin + cnt
system.set_cell([BOX_X, BOX_Y, BOX_Z])
system.set_pbc([True, True, True])
system.center()
print(f"[System] total atoms: {len(system)}")
To write lammps-data atomic we ask ASE to write atom_style atomic
try:
write(OUTFILE, system, format=“lammps-data”, atom_style=ATOM_STYLE)
print(f"Written {OUTFILE} (atom_style={ATOM_STYLE})")
except Exception as e:
print(“Erreur écriture LAMMPS data:”, e)
sys.exit(1)
print(“Done.”)
and here is my input file
in.resine_reax.corrected2.lmp
units real
atom_style atomic
boundary p p p
read_data resine_cnt.atomic.data
Optional: masses ignored by ReaxFF
mass 1 12.011
mass 2 1.008
mass 3 15.999
mass 4 14.007
Soft relaxation
pair_style soft 6.0
pair_coeff * * 1.0
timestep 0.1
fix soft_nve all nve
run 5000
unfix soft_nve
ReaxFF
pair_style reaxff NULL safezone 4.0 mincap 200
pair_coeff * * ffield.reax C H O N
neighbor 2.0 bin
neigh_modify every 10 delay 0 check yes
fix qeq_min all qeq/reax 10 0.0 10.0 1.0e-6 reaxff
reset_timestep 0
min_style sd
minimize 1e-2 1e-4 1000 10000
min_style cg
minimize 1e-4 1e-6 1000 10000
unfix qeq_min
fix qeq all qeq/reax 1 0.0 20.0 1e-6 reaxff
write_data relaxed_after_minim.data
Equilibration NVT
timestep 0.25
velocity all create 298.0 12345 mom yes dist gaussian
fix temp_ber all temp/berendsen 298 298 100
run 80000
unfix temp_ber
Equilibration NPT
fix npt_eq all npt temp 298 298 100 iso 1 1 1000 tchain 1 pchain 1
run 80000
unfix npt_eq
Production
fix npt_prod all npt temp 298 298 100 iso 1 1 1000 tchain 1 pchain 1
run 2000000
write_data final_relaxed.data