{ "cells": [ { "cell_type": "markdown", "id": "da77c3f6", "metadata": {}, "source": [ "# Resin Blend from Molecular Dynamics\n", "\n", "This Notebook builds a small polyisoprene (PI) / hydrogenated dicyclopentadiene (H-DCPD) resin blend inspired by Rissanou et al., *Polymers* **18**, 594 (2026).\n", "\n", "The paper used 100 PI chains with 232 H-DCPD resin molecules for the high-resin PI blend. Here we scale the same composition down to 20 PI chains and 46 H-DCPD molecules, giving a tutorial-size system of about 10,000 atoms.\n", "\n", "The goal is not to reproduce the statistical quality of the production paper. Instead, this notebook follows the same workflow as the GUI tutorial: pack the blend, compress the low-density starting structure, anneal the material, run a short production trajectory, and analyze resin-rich structure.\n" ] }, { "cell_type": "markdown", "id": "0cc45866", "metadata": {}, "source": [ "## Imports and user settings\n", "\n", "The production trajectory is intentionally short for a tutorial; increase the step counts and sampling for publication-quality statistics.\n" ] }, { "cell_type": "code", "execution_count": 1, "id": "53dc1b9c", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PLAMS working folder: /Users/nonofrio/Downloads/Bridgestone/polymer_resin_tutorial/final/PolymerResin_workdir\n" ] } ], "source": [ "from __future__ import annotations\n", "\n", "from pathlib import Path\n", "from typing import Iterable\n", "import warnings\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "from scm.base import Units\n", "from scm.input_classes import AMS\n", "from scm.plams import AMSJob, Molecule, Settings, finish, init, packmol, view\n", "\n", "WORKDIR = Path(\"PolymerResin_workdir\")\n", "STRUCTURE_DIR = Path(\".\")\n", "SYSTEM = {\n", " \"polymer\": STRUCTURE_DIR / \"PI_30mer.xyz\",\n", " \"resin\": STRUCTURE_DIR / \"HDCPD_trimer.xyz\",\n", "}\n", "\n", "N_PI_CHAINS = 20\n", "N_RESIN_MOLECULES = 46\n", "PACKING_DENSITY_G_CM3 = 0.05\n", "PACKMOL_SEED = 2026\n", "\n", "UFF_RELAX_MAX_ITERATIONS = 1000\n", "PRESHRINK_STEPS = 20_000\n", "PRESHRINK_TARGET_DENSITY_G_CM3 = 1.0\n", "PRESHRINK_LAMBDA = 1.5\n", "\n", "REAXFF_FORCE_FIELD = \"dispersion/CHONSSi-lg.ff\"\n", "TEMPERATURE_K = 413.0\n", "PRESSURE_PA = 101325.0\n", "REAXFF_TIMESTEP_FS = 0.5\n", "\n", "ANNEALING_STEPS = 25_000\n", "ANNEALING_SAMPLING_FREQ = 500\n", "MAX_ANNEALING_ITERATIONS = 7\n", "DENSITY_CONVERGENCE_PERCENT = 0.5\n", "\n", "PRODUCTION_STEPS = 50_000\n", "PRODUCTION_SAMPLING_FREQ = 100\n", "\n", "RDF_PAIRS = ((\"resin\", \"resin\"), (\"resin\", \"polymer\"), (\"polymer\", \"polymer\"))\n", "RDF_RMAX = 40.0\n", "RDF_DR = 0.1\n", "ANALYSIS_START_FRACTION = 0.5\n", "DENSITY_MAP_BINS = 80\n", "\n", "init(folder=str(WORKDIR))\n" ] }, { "cell_type": "markdown", "id": "6e34ea0a", "metadata": {}, "source": [ "## General helper functions" ] }, { "cell_type": "code", "execution_count": 2, "id": "68e7b653", "metadata": {}, "outputs": [], "source": [ "KG_M3_TO_G_CM3 = 1.0 / 1000.0\n", "AU_DENSITY_TO_G_CM3 = (\n", " Units.conversion_factor(\"au\", \"kg\")\n", " / Units.conversion_factor(\"bohr\", \"m\") ** 3\n", " * KG_M3_TO_G_CM3\n", ")\n", "\n", "\n", "def validate_settings(settings: Settings) -> Settings:\n", " AMS.from_settings(settings)\n", " return settings\n", "\n", "\n", "def molecule_density_g_cm3(molecule: Molecule) -> float:\n", " return molecule.get_density() * KG_M3_TO_G_CM3\n", "\n", "\n", "def density_trace_g_cm3(job: AMSJob) -> np.ndarray:\n", " densities_au = job.results.get_history_property(\"Density\", \"MDHistory\")\n", " return np.asarray(densities_au, dtype=float) * AU_DENSITY_TO_G_CM3\n", "\n", "\n", "def time_trace_ps(job: AMSJob) -> np.ndarray:\n", " times_fs = job.results.get_history_property(\"Time\", \"MDHistory\")\n", " return np.asarray(times_fs, dtype=float) / 1000.0\n", "\n", "\n", "def last_third_average(values: Iterable[float]) -> float:\n", " values = np.asarray(list(values), dtype=float)\n", " if len(values) == 0:\n", " raise ValueError(\"Cannot average an empty trace\")\n", " return float(np.mean(values[-max(1, len(values) // 3):]))\n", "\n", "\n", "def guess_bonds(molecule: Molecule) -> Molecule:\n", " if hasattr(molecule, \"guess_bonds\"):\n", " molecule.guess_bonds()\n", " return molecule\n", "\n", "\n", "def run_or_load(molecule: Molecule, settings: Settings, name: str) -> tuple[AMSJob, Molecule]:\n", " job_path = WORKDIR / name\n", " if job_path.exists():\n", " print(f\"Loading {name}\")\n", " job = AMSJob.load_external(job_path)\n", " else:\n", " print(f\"Running {name}\")\n", " job = AMSJob(molecule=molecule, settings=settings, name=name)\n", " job.run()\n", "\n", " output_molecule = job.results.get_main_molecule()\n", " print(f\"{name}: density {molecule_density_g_cm3(output_molecule):.3f} g/cm^3\")\n", " return job, output_molecule\n" ] }, { "cell_type": "markdown", "id": "91c41187", "metadata": {}, "source": [ "## AMS settings builders\n", "\n", "We will perform UFF geometry optimization, UFF pre-shrink, ReaxFF geometry optimization, ReaxFF NPT annealing, and ReaxFF NVT production.\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "03ca7567", "metadata": {}, "outputs": [], "source": [ "def uff_relax_settings() -> Settings:\n", " settings = Settings()\n", " settings.input.ams.task = \"GeometryOptimization\"\n", " settings.input.ams.GeometryOptimization.MaxIterations = UFF_RELAX_MAX_ITERATIONS\n", " settings.input.ams.GeometryOptimization.PretendConverged = \"Yes\"\n", " settings.input.ForceField.Type = \"UFF\"\n", " return validate_settings(settings)\n", "\n", "\n", "def uff_shrink_settings(target_lattice: tuple[float, float, float]) -> Settings:\n", " a, b, c = target_lattice\n", " settings = Settings()\n", " settings.input.ams.task = \"MolecularDynamics\"\n", " settings.input.ams.MolecularDynamics.NSteps = PRESHRINK_STEPS\n", " settings.input.ams.MolecularDynamics.TimeStep = 0.25\n", " settings.input.ams.MolecularDynamics.Trajectory.SamplingFreq = 1000\n", " settings.input.ams.MolecularDynamics.InitialVelocities.Temperature = 5.0\n", " settings.input.ams.MolecularDynamics.Thermostat.Type = \"Berendsen\"\n", " settings.input.ams.MolecularDynamics.Thermostat.Temperature = 5.0\n", " settings.input.ams.MolecularDynamics.Thermostat.Tau = 10.0\n", " settings.input.ams.MolecularDynamics.Deformation.TargetLattice._1 = f\"{a:.3f} 0 0\"\n", " settings.input.ams.MolecularDynamics.Deformation.TargetLattice._2 = f\"0 {b:.3f} 0\"\n", " settings.input.ams.MolecularDynamics.Deformation.TargetLattice._3 = f\"0 0 {c:.3f}\"\n", " settings.input.ForceField.Type = \"UFF\"\n", " return validate_settings(settings)\n", "\n", "\n", "def reaxff_relax_settings() -> Settings:\n", " settings = Settings()\n", " settings.input.ams.task = \"GeometryOptimization\"\n", " settings.input.ams.GeometryOptimization.MaxIterations = 1000\n", " settings.input.ams.GeometryOptimization.PretendConverged = \"Yes\"\n", " settings.input.ReaxFF.ForceField = REAXFF_FORCE_FIELD\n", " return validate_settings(settings)\n", "\n", "\n", "def reaxff_md_settings(nsteps: int, sampling_freq: int, npt: bool) -> Settings:\n", " settings = Settings()\n", " settings.input.ams.task = \"MolecularDynamics\"\n", " settings.input.ams.MolecularDynamics.NSteps = nsteps\n", " settings.input.ams.MolecularDynamics.TimeStep = REAXFF_TIMESTEP_FS\n", " settings.input.ams.MolecularDynamics.Trajectory.SamplingFreq = sampling_freq\n", " settings.input.ams.MolecularDynamics.InitialVelocities.Temperature = TEMPERATURE_K\n", " settings.input.ams.MolecularDynamics.Thermostat.Type = \"Berendsen\"\n", " settings.input.ams.MolecularDynamics.Thermostat.Temperature = TEMPERATURE_K\n", " settings.input.ams.MolecularDynamics.Thermostat.Tau = 100.0\n", " if npt:\n", " settings.input.ams.MolecularDynamics.Barostat.Type = \"Berendsen\"\n", " settings.input.ams.MolecularDynamics.Barostat.Pressure = PRESSURE_PA\n", " settings.input.ams.MolecularDynamics.Barostat.Tau = 500.0\n", " settings.input.ReaxFF.ForceField = REAXFF_FORCE_FIELD\n", " return validate_settings(settings)\n" ] }, { "cell_type": "markdown", "id": "76a44bfd", "metadata": {}, "source": [ "## Build the packed PI/H-DCPD blend\n", "\n", "Packmol assigns region names to the atoms. We use `polymer` for the PI chains and `resin` for the H-DCPD molecules; those same labels are used later for the RDFs and the density map.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "b4e8af17", "metadata": {}, "outputs": [], "source": [ "pi_chain = guess_bonds(Molecule(SYSTEM[\"polymer\"]))\n", "resin = guess_bonds(Molecule(SYSTEM[\"resin\"]))\n", "\n", "blend, pack_details = packmol(\n", " molecules=[pi_chain, resin],\n", " n_molecules=[N_PI_CHAINS, N_RESIN_MOLECULES],\n", " density=PACKING_DENSITY_G_CM3,\n", " region_names=[\"polymer\", \"resin\"],\n", " tolerance=2.0,\n", " seed=PACKMOL_SEED,\n", " return_details=True,\n", ")\n", "\n", "blend.write(\"polymer_resin_initial.in\")\n", "print(f\"Atoms: {len(blend)}\")\n", "print(f\"Initial density: {molecule_density_g_cm3(blend):.3f} g/cm^3\")\n", "view(blend, direction=\"tilt_z\", guess_bonds=True)\n" ] }, { "cell_type": "markdown", "id": "053b3510", "metadata": {}, "source": [ "## Equilibrate the polymer bulk\n", "\n", "First we relax the low-density packed structure with UFF. Then we pre-shrink it into a pseudo-2D periodic cell. The target volume is computed from the current density and the target density, and the cell dimensions are chosen from `a = b = lambda*c`. Finally, we relax the shrunken structure with ReaxFF before starting molecular dynamics.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "d23c5984", "metadata": {}, "outputs": [], "source": [ "def preshrink_lattice(molecule: Molecule) -> tuple[float, float, float, float, float]:\n", " current_density = molecule_density_g_cm3(molecule)\n", " current_volume = molecule.unit_cell_volume()\n", " target_volume = current_density * current_volume / PRESHRINK_TARGET_DENSITY_G_CM3\n", " new_a = (target_volume * PRESHRINK_LAMBDA) ** (1.0 / 3.0)\n", " new_b = new_a\n", " new_c = new_a / PRESHRINK_LAMBDA\n", " return current_density, target_volume, new_a, new_b, new_c\n", "\n", "\n", "uff_relax_job, relaxed_blend = run_or_load(\n", " blend,\n", " uff_relax_settings(),\n", " \"UFF_InitialRelax\",\n", ")\n", "\n", "current_density, target_volume, new_a, new_b, new_c = preshrink_lattice(relaxed_blend)\n", "print(f\"Current density: {current_density:.3f} g/cm^3\")\n", "print(f\"Target density: {PRESHRINK_TARGET_DENSITY_G_CM3:.3f} g/cm^3\")\n", "print(f\"Target volume: {target_volume:.1f} A^3\")\n", "print(f\"Target lattice: a={new_a:.3f}, b={new_b:.3f}, c={new_c:.3f} A\")\n", "\n", "uff_shrink_job, shrunken_blend = run_or_load(\n", " relaxed_blend,\n", " uff_shrink_settings((new_a, new_b, new_c)),\n", " \"UFF_Shrink\",\n", ")\n", "\n", "reaxff_relax_job, reaxff_relaxed_blend = run_or_load(\n", " shrunken_blend,\n", " reaxff_relax_settings(),\n", " \"ReaxFF_InitialRelax\",\n", ")\n" ] }, { "cell_type": "markdown", "id": "fb667d6d", "metadata": {}, "source": [ "## Anneal at 413 K\n", "\n", "The annealing stage uses ReaxFF NPT molecular dynamics at 413 K. The loop repeats the annealing job from the final structure of the previous run until the last-third average density changes by less than about 0.5 percent, or until the maximum number of tutorial iterations is reached.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "3970b9ce", "metadata": {}, "outputs": [], "source": [ "annealing_settings = reaxff_md_settings(\n", " ANNEALING_STEPS,\n", " ANNEALING_SAMPLING_FREQ,\n", " npt=True,\n", ")\n", "annealing_jobs: list[AMSJob] = []\n", "density_summary: list[dict[str, float | str]] = []\n", "\n", "current_molecule = reaxff_relaxed_blend\n", "previous_last_third_density: float | None = None\n", "\n", "for iteration in range(1, MAX_ANNEALING_ITERATIONS + 1):\n", " job_name = f\"ReaxFF_Anneal_{iteration}\"\n", " initial_density = molecule_density_g_cm3(current_molecule)\n", " job, current_molecule = run_or_load(current_molecule, annealing_settings, job_name)\n", " annealing_jobs.append(job)\n", "\n", " densities = density_trace_g_cm3(job)\n", " last_third_density = last_third_average(densities)\n", " if previous_last_third_density is None:\n", " density_change = np.nan\n", " else:\n", " density_change = abs(last_third_density / previous_last_third_density - 1.0) * 100.0\n", "\n", " density_summary.append(\n", " {\n", " \"Annealing step\": iteration,\n", " \"Initial density (g/cm^3)\": initial_density,\n", " \"Final density (g/cm^3)\": float(densities[-1]),\n", " \"Last-third average density (g/cm^3)\": last_third_density,\n", " \"Change (%)\": density_change,\n", " }\n", " )\n", " print(density_summary[-1])\n", "\n", " if previous_last_third_density is not None and density_change < DENSITY_CONVERGENCE_PERCENT:\n", " break\n", " previous_last_third_density = last_third_density\n", "\n", "equilibrated_blend = current_molecule\n" ] }, { "cell_type": "code", "execution_count": null, "id": "5591bda6", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(6.0, 3.2))\n", "for job in annealing_jobs:\n", " ax.plot(time_trace_ps(job), density_trace_g_cm3(job), label=job.name)\n", "ax.set_xlabel(\"Time (ps)\")\n", "ax.set_ylabel(\"Density (g/cm^3)\")\n", "ax.legend(frameon=False)\n", "fig.tight_layout()\n", "fig.savefig(\"density_per_annealing_step.png\")\n" ] }, { "cell_type": "markdown", "id": "5c1ee4bf", "metadata": {}, "source": [ "## Run the production trajectory\n", "\n", "The production run uses ReaxFF NVT molecular dynamics at 413 K. This short trajectory is enough for the tutorial analysis, but longer production runs and independent replicas are recommended for smoother RDFs and better statistics.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "4e27ffc6", "metadata": {}, "outputs": [], "source": [ "production_job, production_blend = run_or_load(\n", " equilibrated_blend,\n", " reaxff_md_settings(PRODUCTION_STEPS, PRODUCTION_SAMPLING_FREQ, npt=False),\n", " \"ReaxFF_Production\",\n", ")\n", "\n", "print(f\"Final production density: {molecule_density_g_cm3(production_blend):.3f} g/cm^3\")\n" ] }, { "cell_type": "markdown", "id": "b6835b99", "metadata": {}, "source": [ "## Center-of-mass RDFs\n", "\n", "To analyze resin miscibility in the polymer matrix, we compute center-of-mass RDFs for resin-resin, resin-polymer, and polymer-polymer pairs. Molecules are identified from the `polymer` and `resin` region labels, using the atom counts of the original input molecules to split each region into molecules.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "02b26d20", "metadata": {}, "outputs": [], "source": [ "def atom_region(atom) -> str | None:\n", " properties = atom.properties\n", " return properties.get(\"region\") if hasattr(properties, \"get\") else getattr(properties, \"region\", None)\n", "\n", "\n", "def region_molecules(molecule: Molecule, region: str, atoms_per_molecule: int) -> list[list[int]]:\n", " atom_indices = [i for i, atom in enumerate(molecule, start=1) if atom_region(atom) == region]\n", " if not atom_indices:\n", " raise ValueError(f\"No atoms found for region '{region}'\")\n", " if len(atom_indices) % atoms_per_molecule:\n", " raise ValueError(\n", " f\"Region '{region}' has {len(atom_indices)} atoms, which is not divisible by \"\n", " f\"{atoms_per_molecule} atoms per molecule\"\n", " )\n", " return [\n", " atom_indices[i : i + atoms_per_molecule]\n", " for i in range(0, len(atom_indices), atoms_per_molecule)\n", " ]\n", "\n", "\n", "def lattice_matrix(molecule: Molecule) -> np.ndarray:\n", " return np.asarray(molecule.lattice, dtype=float)\n", "\n", "\n", "def minimum_image_vectors(vectors: np.ndarray, lattice: np.ndarray) -> np.ndarray:\n", " fractional = vectors @ np.linalg.inv(lattice)\n", " return vectors - np.round(fractional) @ lattice\n", "\n", "\n", "def center_of_mass(molecule: Molecule, atom_indices: list[int]) -> np.ndarray:\n", " molecule_atoms = list(molecule)\n", " atoms = [molecule_atoms[i - 1] for i in atom_indices]\n", " masses = np.asarray([atom.mass for atom in atoms], dtype=float)\n", " coords = np.asarray([atom.coords for atom in atoms], dtype=float)\n", " lattice = lattice_matrix(molecule)\n", " unwrapped = coords[0] + minimum_image_vectors(coords - coords[0], lattice)\n", " com = np.average(unwrapped, axis=0, weights=masses)\n", " return com - np.floor(com @ np.linalg.inv(lattice)) @ lattice\n", "\n", "\n", "def pair_distances(points_a: np.ndarray, points_b: np.ndarray, lattice: np.ndarray, same_region: bool) -> np.ndarray:\n", " if same_region:\n", " vectors = [\n", " points_a[j] - points_a[i]\n", " for i in range(len(points_a))\n", " for j in range(i + 1, len(points_a))\n", " ]\n", " else:\n", " vectors = [point_b - point_a for point_a in points_a for point_b in points_b]\n", " if not vectors:\n", " return np.empty(0, dtype=float)\n", " return np.linalg.norm(minimum_image_vectors(np.asarray(vectors), lattice), axis=1)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "1f35181e", "metadata": {}, "outputs": [], "source": [ "def calculate_rdfs(\n", " job: AMSJob,\n", " pairs: tuple[tuple[str, str], ...],\n", " atoms_per_molecule: dict[str, int],\n", " rmax: float,\n", " dr: float,\n", " start_fraction: float,\n", ") -> tuple[np.ndarray, dict[str, np.ndarray]]:\n", " n_frames = job.results.readrkf(\"History\", \"nEntries\")\n", " start_frame = max(1, min(n_frames, int(n_frames * start_fraction) + 1))\n", " first_molecule = job.results.get_history_molecule(start_frame)\n", " max_r = min(np.linalg.norm(vector) for vector in lattice_matrix(first_molecule)) / 2.0\n", " if rmax > max_r:\n", " warnings.warn(f\"Capping rmax from {rmax:.3f} to {max_r:.3f} A\", RuntimeWarning)\n", " rmax = max_r\n", "\n", " bins = np.arange(0.0, rmax + dr, dr)\n", " histograms = {f\"{a}-{b}\": np.zeros(len(bins) - 1, dtype=float) for a, b in pairs}\n", " volumes: list[float] = []\n", " molecule_counts: dict[str, int] | None = None\n", "\n", " print(f\"Using RDF frames {start_frame}-{n_frames} of {n_frames}\")\n", " for frame in range(start_frame, n_frames + 1):\n", " molecule = job.results.get_history_molecule(frame)\n", " centers = {\n", " region: np.asarray(\n", " [\n", " center_of_mass(molecule, group)\n", " for group in region_molecules(molecule, region, atoms_per_molecule[region])\n", " ]\n", " )\n", " for region in atoms_per_molecule\n", " }\n", "\n", " counts = {region: len(points) for region, points in centers.items()}\n", " if molecule_counts is None:\n", " molecule_counts = counts\n", " print(\"Molecule counts:\", molecule_counts)\n", " elif counts != molecule_counts:\n", " raise ValueError(f\"Frame {frame} has a different number of molecules\")\n", "\n", " lattice = lattice_matrix(molecule)\n", " for region_a, region_b in pairs:\n", " distances = pair_distances(centers[region_a], centers[region_b], lattice, region_a == region_b)\n", " histograms[f\"{region_a}-{region_b}\"] += np.histogram(distances, bins=bins)[0]\n", " volumes.append(molecule.unit_cell_volume())\n", "\n", " if molecule_counts is None or not volumes:\n", " raise ValueError(\"No trajectory frames were used\")\n", "\n", " radii = 0.5 * (bins[:-1] + bins[1:])\n", " shell_volumes = 4.0 * np.pi / 3.0 * (bins[1:] ** 3 - bins[:-1] ** 3)\n", " volume = float(np.mean(volumes))\n", "\n", " rdf_data: dict[str, np.ndarray] = {}\n", " for region_a, region_b in pairs:\n", " n_a = molecule_counts[region_a]\n", " n_b = molecule_counts[region_b]\n", " pair_count = n_a * (n_a - 1) / 2.0 if region_a == region_b else n_a * n_b\n", " normalization = len(volumes) * shell_volumes * pair_count / volume\n", " rdf_data[f\"{region_a}-{region_b}\"] = histograms[f\"{region_a}-{region_b}\"] / normalization\n", "\n", " return radii, rdf_data\n" ] }, { "cell_type": "code", "execution_count": null, "id": "530a44c8", "metadata": {}, "outputs": [], "source": [ "atoms_per_molecule = {\n", " \"polymer\": len(Molecule(SYSTEM[\"polymer\"])),\n", " \"resin\": len(Molecule(SYSTEM[\"resin\"])),\n", "}\n", "radii, rdf_data = calculate_rdfs(\n", " production_job,\n", " RDF_PAIRS,\n", " atoms_per_molecule,\n", " RDF_RMAX,\n", " RDF_DR,\n", " ANALYSIS_START_FRACTION,\n", ")\n", "\n", "fig, ax = plt.subplots(figsize=(5.2, 3.4))\n", "for label, values in rdf_data.items():\n", " ax.plot(radii, values, label=label)\n", "ax.set_xlabel(\"r (A)\")\n", "ax.set_ylabel(\"g(r)\")\n", "ax.legend(frameon=False)\n", "fig.tight_layout()\n", "fig.savefig(\"RDF_com_regions.png\")\n", "\n", "rdf_table = np.column_stack([radii, *rdf_data.values()])\n", "np.savetxt(\"RDF_com_regions.txt\", rdf_table, header=\"r_A \" + \" \".join(rdf_data))\n" ] }, { "cell_type": "markdown", "id": "06ad0baf", "metadata": {}, "source": [ "## Resin density map\n", "\n", "The 2D density map bins resin atoms in the x-y plane over the second half of the production trajectory. This gives a simple view of resin-rich domains in the pseudo-2D cell.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "535ada89", "metadata": {}, "outputs": [], "source": [ "def region_indices(molecule: Molecule, region: str) -> list[int]:\n", " return [i for i, atom in enumerate(molecule, start=1) if atom_region(atom) == region]\n", "\n", "\n", "def density_map_xy(\n", " job: AMSJob,\n", " region: str = \"resin\",\n", " bins: int = DENSITY_MAP_BINS,\n", " start_fraction: float = ANALYSIS_START_FRACTION,\n", ") -> tuple[np.ndarray, np.ndarray, np.ndarray]:\n", " n_frames = job.results.readrkf(\"History\", \"nEntries\")\n", " start_frame = max(1, min(n_frames, int(n_frames * start_fraction) + 1))\n", "\n", " final_molecule = job.results.get_main_molecule()\n", " lx = final_molecule.lattice[0][0]\n", " ly = final_molecule.lattice[1][1]\n", " histogram = np.zeros((bins, bins), dtype=float)\n", "\n", " print(f\"Using density-map frames {start_frame}-{n_frames} of {n_frames}\")\n", " for used_frames, frame in enumerate(range(start_frame, n_frames + 1), start=1):\n", " molecule = job.results.get_history_molecule(frame)\n", " atoms = list(molecule)\n", " indices = region_indices(molecule, region)\n", " if not indices:\n", " raise ValueError(f\"No atoms found for region '{region}'\")\n", "\n", " coords = np.asarray([atoms[i - 1].coords[:2] for i in indices], dtype=float)\n", " coords[:, 0] %= lx\n", " coords[:, 1] %= ly\n", " histogram += np.histogram2d(coords[:, 0], coords[:, 1], bins=bins, range=[[0.0, lx], [0.0, ly]])[0]\n", "\n", " dx_nm = (lx / bins) / 10.0\n", " dy_nm = (ly / bins) / 10.0\n", " density = histogram / used_frames / (dx_nm * dy_nm)\n", " x_edges = np.linspace(0.0, lx / 10.0, bins + 1)\n", " y_edges = np.linspace(0.0, ly / 10.0, bins + 1)\n", " return x_edges, y_edges, density.T\n" ] }, { "cell_type": "code", "execution_count": null, "id": "fc2e7dd6", "metadata": {}, "outputs": [], "source": [ "x_edges, y_edges, resin_density = density_map_xy(production_job)\n", "np.savetxt(\"DENSMAP_resin.txt\", resin_density, header=\"resin atoms / nm^2 / frame\")\n", "\n", "fig, ax = plt.subplots(figsize=(5.0, 4.0))\n", "mesh = ax.pcolormesh(x_edges, y_edges, resin_density, shading=\"auto\")\n", "fig.colorbar(mesh, ax=ax, label=\"resin atoms / nm2 / frame\")\n", "ax.set_xlabel(\"x (nm)\")\n", "ax.set_ylabel(\"y (nm)\")\n", "ax.set_title(\"Resin density map\")\n", "fig.tight_layout()\n", "fig.savefig(\"DENSMAP_resin.png\")\n" ] }, { "cell_type": "markdown", "id": "dbbb09ed", "metadata": {}, "source": [ "## Summary\n", "\n", "You built a small PI/H-DCPD blend, compressed it into a pseudo-2D bulk cell, annealed it with ReaxFF, and inspected the miscibility of the resin with center-of-mass RDFs and a resin density map.\n", "\n", "For better RDF statistics, run a longer production trajectory, save frames more frequently, use independent replicas, or increase the system size.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "8ea52b4a", "metadata": {}, "outputs": [], "source": [ "finish()\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.16" } }, "nbformat": 4, "nbformat_minor": 5 }