#!/usr/bin/env python # coding: utf-8 # ## Overview # # This notebook shows a PLAMS workflow for an accelerated thermo-oxidative degradation study with AMS/ReaxFF. The model system is a poly(acrylate) proxy oligomer mixed with molecular oxygen. The chemistry we want to observe is the onset of oxidation and fragmentation in an oxygen-rich polymer environment: oxygen should gradually be consumed, oxidized products such as CO2 can appear, and the original large organic material can break into smaller oxygenated and hydrocarbon fragments. # The workflow is intentionally qualitative rather than quantitative. The temperatures are chosen high enough that reactive events become visible on short MD timescales, so the purpose is not to reproduce an experimental degradation rate directly, but to show how to set up and analyze this type of reactive simulation. # # The main properties followed in the analysis are the amount of `unreacted_oxygen` (`O2`), the amount of `co2`, the populations of `small_neutral`, `small_oxygenated`, `medium_oxygenated_fragment`, `oxygenated_fragment`, `parent_or_large_fragment`, and `hydrocarbon_fragment`, and a `scission_proxy` that acts as a qualitative indicator for fragmentation. These quantities are important because they summarize, in a compact way, how much oxidant is left, how far oxidation has progressed, and how strongly the original polymer-like material is being broken down. # # In this notebook the family labels have the following practical meaning: # # - `unreacted_oxygen`: molecules identified as `O2`. # - `co2`: molecules identified as `CO2`. # - `small_neutral`: very small neutral molecules tracked explicitly here as `CO`, `H2O`, or `H2`. # - `small_oxygenated`: carbon-and-oxygen-containing fragments with at most 4 carbon atoms. # - `medium_oxygenated_fragment`: oxygen-containing organic fragments with more than 4 carbon atoms but fewer than 8. # - `oxygenated_fragment`: oxygen-containing fragments that do not fall into the previous carbon-based oxygenated groups. # - `parent_or_large_fragment`: large carbon-containing fragments with at least 8 carbon atoms. This group includes material still close to the starting oligomer size. # - `hydrocarbon_fragment`: carbon-containing fragments without oxygen that do not fall into the large-fragment group. # # The logic of the workflow is therefore: first build one representative polymer proxy oligomer, then place it in contact with explicit oxygen, and finally compare several temperatures. The oligomer is cleaned before packing so that the reactive box starts from a sensible structure. Oxygen is added explicitly because the goal is to study oxidation chemistry rather than pure thermal fragmentation. Running four temperatures lets us see whether the same qualitative chemistry becomes stronger or faster as the temperature increases. # # When you reach the time-dependent plots, you should be able to read them as a degradation story: a faster decrease of unreacted oxygen suggests more rapid oxidation, increasing CO2 suggests deeper oxidative breakdown, increasing oxygenated fragments points to oxidized intermediates, and a growing scission proxy suggests stronger fragmentation of the original material. # # So, let's start! # ## Imports and global settings import matplotlib.pyplot as plt import numpy as np import pandas as pd from scm.plams import AMSJob, JobRunner, Settings, config, from_smiles, packmol, plot_image_grid, view # ## Practical notes # # - ReaxFF is often efficient in serial for systems around this size, so each MD job uses `nproc = 1` and `OMP_NUM_THREADS = 1`. # - The different temperatures are still run concurrently through the PLAMS `JobRunner`. # - The trajectories are deliberately short and accelerated. Increase the number of MD steps for a longer study, or decrease the trajectory sampling frequency if you want finer reaction-event resolution. # ## Reusable helper functions # # The helper functions for this notebook live in the local file `thermo_oxidative_polyacrylate_reaxff_helpers.py`. # # They convert AMS molecule histories into tables, group species into the product families, and prepare representative molecules for display. # # You can continue the notebook without getting into the details of those functions. Only users who want to inspect or modify the implementation need to open that helper file. from thermo_oxidative_polyacrylate_reaxff_helpers import ( add_degradation_metrics, find_representative_molecule, molecule_counts_from_ams_results, prepare_molecule_for_display, species_family, summarize_product_families, ) # The product-family labels introduced in the overview are only analysis categories. They are useful for comparing temperature trends, but they should not be interpreted as a complete mechanistic description of every reaction event in the trajectory. # ## Build and clean the poly(acrylate) proxy oligomer # # We begin by building one poly(acrylate) proxy oligomer from SMILES. This oligomer is the organic material whose oxidative degradation we want to follow later in the reactive box. Before mixing it with oxygen, we run a short ReaxFF geometry optimization so the oligomer starts from a reasonable, clash-free structure. The goal here is not an exhaustive conformer search, but simply to create a sensible building block for the packed simulation cell. config.default_jobrunner = JobRunner(parallel=True, maxjobs=4) config.default_jobmanager.hashing = None polymer = from_smiles("CC(C(=O)OC)CC(C(=O)OCC)CC(C(=O)OCCCC)CC(C(=O)OC)CC(C(=O)OCC)CC(C(=O)OCCCC)") oxygen = from_smiles("O=O") print("Oligomer atoms:", len(polymer)) print("O2 atoms:", len(oxygen)) print("Target total atoms:", 8 * len(polymer) + 80 * len(oxygen)) go_settings = Settings() go_settings.input.ReaxFF.ForceField = "CHO.ff" go_settings.input.ams.Task = "GeometryOptimization" go_settings.runscript.nproc = 1 go_settings.runscript.preamble_lines = ["export OMP_NUM_THREADS=1"] go_job = AMSJob(name="00_clean_polyacrylate_oligomer", molecule=polymer, settings=go_settings) go_results = go_job.run(); if go_results.ok(): polymer_clean = go_results.get_main_molecule() print("Oligomer cleanup ok: True") else: polymer_clean = polymer print("Oligomer cleanup ok: False") print("Continuing with the SMILES-built oligomer.") view(polymer_clean, direction="along_pca3", picture_path="picture1.png") # ## Pack the polymer and oxygen box # # We now add explicit molecular oxygen and use Packmol to build the periodic simulation box. This step defines the reactive environment: the oligomers provide the organic material and O2 provides the oxidant. The chosen composition keeps the system close to the 1000-atom scale, which is still practical for short ReaxFF trajectories while being large enough to display competing oxidation and fragmentation events. box = packmol( molecules=[polymer_clean, oxygen], n_molecules=[8, 80], density=0.65, ) view(box, direction="tilt_z", picture_path="picture2.png") # ## Run short ReaxFF MD simulations at several temperatures # # We run the same reactive system at four temperatures so that we can compare how strongly the oxidation and fragmentation chemistry is activated as the temperature increases. Each simulation uses one core, while PLAMS runs the different temperatures concurrently through a parallel `JobRunner`. The trajectories are intentionally short so the notebook remains practical, but they are still long enough to reveal differences in oxygen consumption, fragment formation, and overall degradation activity. A fixed AMS random seed is also assigned so the initial velocity draw is reproducible when the notebook is rerun. # # The MD setup is deliberately simple. Random initial velocities are drawn at the target temperature, and a Berendsen thermostat is used to keep each trajectory close to that temperature during these short reactive runs. This is convenient here because the goal is comparative screening of oxidation activity at fixed temperature, not rigorous sampling of a particular statistical ensemble. The thermostat time constant `Tau = 100` fs gives fairly gentle coupling: it is short enough to control the temperature over a 12.5 ps run, but not so tight that every fast vibrational fluctuation is immediately damped away. jobs = {} for temperature in [1700, 1850, 2000, 2150]: settings = Settings() settings.input.ReaxFF.ForceField = "CHO.ff" settings.input.ams.Task = "MolecularDynamics" settings.input.ams.Properties.BondOrders = "Yes" settings.input.ams.RNGSeed = [10000 + temperature] md = settings.input.ams.MolecularDynamics md.NSteps = 50000 md.TimeStep = 0.25 md.InitialVelocities.Type = "Random" md.InitialVelocities.Temperature = float(temperature) md.Thermostat.Type = "Berendsen" md.Thermostat.Temperature = float(temperature) md.Thermostat.Tau = 100 md.Trajectory.SamplingFreq = 100 md.Trajectory.WriteMolecules = "Yes" md.Trajectory.WriteBonds = "Yes" md.Trajectory.WriteCharges = "Yes" settings.runscript.nproc = 1 settings.runscript.preamble_lines = ["export OMP_NUM_THREADS=1"] job = AMSJob( name=f"MD_T{temperature}K", molecule=box, settings=settings, ) jobs[temperature] = job for job in jobs.values(): job.run() for temperature, job in jobs.items(): if not job.results.ok(): raise RuntimeError(f"MD failed at {temperature} K") # ## Convert the trajectories into family summaries # # The next cell converts each trajectory into a compact table of family populations and then extracts a few summary metrics for `final_summary`. Here `final_*` means the value in the last stored frame of a trajectory, while `peak_*` means the largest value reached anywhere along that trajectory. # # This distinction matters because some species are transient. For example, `final_small_oxygenated` tells you how many short oxygen-containing fragments are still present at the end, whereas `peak_small_oxygenated` tells you whether the trajectory passed through a stage with stronger buildup of such intermediates even if they were later consumed. Likewise, `final_scission_proxy` is only the instantaneous fragmentation signal in the last frame, while `cumulative_scission_proxy` adds that signal over time and is therefore a more robust indicator of how much fragmentation activity occurred during the whole run. # # The `final_summary` table is important because it condenses the four short MD trajectories into a form that is easy to compare across temperature. It does not replace the time-dependent plots, but it makes it much easier to spot which temperatures leave more oxygen unreacted, produce more CO2, or generate stronger overall fragmentation. Notice that we also store one representative molecule for each family while scanning the trajectories, so those categories can be illustrated later with actual fragments from the MD history. all_family_rows = [] final_rows = [] representative_molecules = {} representative_formulas = {} time_step = 0.25 # ps sampling_freq = 100 fs_to_ps = 0.001 for temperature, job in jobs.items(): counts = molecule_counts_from_ams_results(job.results) counts["time_ps"] = ( (counts["frame"].astype(int) - 1) * time_step * sampling_freq ) * fs_to_ps families = summarize_product_families(counts) families = add_degradation_metrics(families) families["temperature_K"] = temperature all_family_rows.append(families) species_columns = [column for column in counts.columns if column not in {"frame", "time_ps"}] for species in species_columns: family = species_family(species) if family not in representative_molecules and counts[species].max() > 0: molecule = find_representative_molecule(job.results, species) if molecule is not None: representative_molecules[family] = molecule representative_formulas[family] = species final_row = families.iloc[-1] peak_row = families.max(numeric_only=True) final_rows.append( { "temperature_K": temperature, "time_ps": final_row["time_ps"], "final_unreacted_oxygen": final_row["unreacted_oxygen"], "final_co2": final_row["co2"], "peak_co2": peak_row["co2"], "final_small_oxygenated": final_row["small_oxygenated"], "peak_small_oxygenated": peak_row["small_oxygenated"], "final_medium_oxygenated": final_row["medium_oxygenated_fragment"], "peak_medium_oxygenated": peak_row["medium_oxygenated_fragment"], "final_scission_proxy": final_row["scission_proxy"], "cumulative_scission_proxy": final_row["cumulative_scission_proxy"], } ) family_long = pd.concat(all_family_rows, ignore_index=True) final_summary = pd.DataFrame(final_rows).sort_values("temperature_K") final_summary # The current run, the `final_summary` table already shows a clean temperature trend. The amount of unreacted `O2` decreases from 76 at 1700 K to 62 at 2150 K, while final `CO2` rises from 5 to 20 and the cumulative scission proxy rises from 16 to 52. The increase in `small_oxygenated` fragments is also strong, especially between 2000 and 2150 K. Taken together, this suggests a progression from relatively mild oxidation at 1700 K to much more extensive oxidative fragmentation at the highest temperatures. # ## Representative molecules for the product families # # It is useful to connect the family labels to actual fragments found in the trajectories. The panel below shows one representative molecule for each family that appeared in the current MD runs. # # These figures are helpful because they anchor the abstract labels to concrete structures. They let you see, for example, whether a family contains oligomer-like material, short oxygenated chain fragments, or small stable products such as `CO2`. This is still only an illustrative snapshot: each family can contain many different structures, especially at higher temperature, so the panel should be read as a visual guide rather than a complete mechanistic inventory. # # Some families may be absent in a short trajectory set, so only the families that were actually found are shown. family_images = { f"{family}\n{representative_formulas[family]}": view( prepare_molecule_for_display(representative_molecules[family]), width=220, height=180, guess_bonds=False ) for family in representative_molecules } plot_image_grid(family_images, rows=2, save_path="picture3.png"); # ## Plot the time-dependent trends # # Here we make a figure that shows how a few coarse observables evolve over time for each temperature: the remaining oxygen, the amount of `CO2`, the population of oxygenated fragments, and the cumulative scission proxy. A faster drop in unreacted oxygen suggests more rapid oxidation, increasing `CO2` points to deeper oxidative breakdown, rising oxygenated-fragment counts indicate the build-up of oxidized intermediates, and a growing scission proxy signals stronger fragmentation of the original polymer material. fig, axes = plt.subplots(2, 2, figsize=(12, 8), sharex=True) axes = axes.ravel() for temperature, group in family_long.groupby("temperature_K"): axes[0].plot(group["time_ps"], group["unreacted_oxygen"], label=f"{temperature} K") axes[1].plot(group["time_ps"], group["co2"], label=f"{temperature} K") axes[2].plot( group["time_ps"], group["small_oxygenated"] + group["medium_oxygenated_fragment"], label=f"{temperature} K", ) axes[3].plot(group["time_ps"], group["cumulative_scission_proxy"], label=f"{temperature} K") axes[0].set_title("Unreacted O2") axes[1].set_title("CO2") axes[2].set_title("Oxygenated fragments") axes[3].set_title("Cumulative scission proxy") for ax in axes: ax.set_xlabel("time / ps") ax.set_ylabel("count") ax.legend(fontsize=8) fig.tight_layout() axes; import matplotlib.pyplot as plt plt.gcf().savefig("picture4.png") # We can see a few interesting things from the figure above. Comparing when the curves in the different panels start to move gives us a qualitative timeline of the reaction sequence. If we look at the blue line in the top-left panel ("Unreacted O2" at 1700 K), it remains relatively flat at the beginning, indicating an induction period. Moving to the bottom-right panel ("Cumulative scission proxy"), notice how its stair-step lines begin to rise almost when the oxygen levels drop, illustrating how oxygen consumption triggers polymer fragmentation. Furthermore, if we check the CO2 panel (top-right), there is a time lag; these curves take much longer to take off compared to the immediate, early growth of the oxygenated fragments (bottom-left panel). This visual difference shows that deep mineralization into CO2 is a secondary process that requires the prior fragmentation of the polymer chain. # ## Compare the final metrics versus temperature # # A second plot compresses each trajectory into a few end-point or cumulative values. This makes it easier to compare the short accelerated runs at a glance and to see which temperatures lead to the strongest overall oxidation and fragmentation within the same simulation time. Two of the most useful quantities here are `peak_small_oxygenated` and `cumulative_scission_proxy`. A high `peak_small_oxygenated` value suggests that the trajectory passes through a stage rich in short oxidized fragments, which is often a sign that larger oxygenated pieces are being broken into smaller intermediates. A high `cumulative_scission_proxy` means that fragmentation activity is strong over the full trajectory, not just in the final frame. Together they help distinguish a system that briefly forms oxidized fragments from one that keeps fragmenting throughout the run. # # In the present figure, all three curves increase with temperature, so the short ReaxFF runs give a consistent phenomenological picture: hotter trajectories leave less unreacted oxygen, form more `CO2`, and generate more fragmentation. The strongest jump is between 2000 and 2150 K, where both `final_co2` and the scission proxy increase markedly. That makes the plot useful as a compact screening summary: if the goal is to identify a temperature window where oxidative degradation becomes much more aggressive, this panel does that immediately. fig, ax = plt.subplots(figsize=(7, 4.5)) ax.plot(final_summary["temperature_K"], final_summary["final_co2"], marker="o", label="final CO2") ax.plot( final_summary["temperature_K"], final_summary["cumulative_scission_proxy"], marker="o", label="cumulative scission proxy", ) ax.plot( final_summary["temperature_K"], final_summary["peak_small_oxygenated"], marker="o", label="peak small oxygenated", ) ax.set_xlabel("temperature / K") ax.set_ylabel("count") ax.set_title("Temperature effect in accelerated ReaxFF MD") ax.legend() fig.tight_layout() ax; ax.figure.savefig("picture5.png") # The diverging trends in this final graph perfectly illustrate how temperature alters competing mechanistic pathways. Following the orange line (cumulative scission proxy), we observe a steady, moderate slope up to 2000 K, followed by a drastic, sharp jump toward 2150 K. In stark contrast, if we look at the blue line (final CO2), its growth visibly slows down at those same extreme temperatures. This visual separation, where the orange curve shoots up while the blue curve plateaus, shows, that excessive heat shreds the polymer backbone much faster than the system can sequentially oxidize those newly formed fragments. # ## Fit an apparent Arrhenius law to the degradation proxy # # As a simple phenomenological model, we can treat the degradation signal as Arrhenius-like over this temperature window: # # $$k_{\mathrm{app}}(T) = A_{\mathrm{app}} e^{-E_{a,\mathrm{app}} / RT}$$ # # or, equivalently, # # $$\ln k_{\mathrm{app}} = \ln A_{\mathrm{app}} - \left(\frac{E_{a,\mathrm{app}}}{R}\right)\frac{1}{T}$$ # # Here `k_app` is not a fundamental rate constant. It is only the chosen trajectory-level proxy, taken as `cumulative_scission_proxy / time`. The fit therefore gives an *apparent* activation energy and an *apparent* preexponential factor for this accelerated ReaxFF observable. arr = final_summary.copy() metric = "cumulative_scission_proxy" arr["rate_proxy_per_ps"] = arr[metric] / arr["time_ps"].clip(lower=1e-12) arr = arr[arr["rate_proxy_per_ps"] > 0].copy() R = 8.31446261815324 # J/K/mol Gas constant J_to_kJ = 0.001 ps_to_s = 1.0e12 x = 1.0 / arr["temperature_K"].to_numpy(dtype=float) y = np.log(arr["rate_proxy_per_ps"].to_numpy(dtype=float)) slope, intercept = np.polyfit(x, y, 1) apparent_Ea_kJ_mol = -slope * R * J_to_kJ apparent_A_per_ps = np.exp(intercept) apparent_A_per_s = apparent_A_per_ps * ps_to_s arr["fit_ln_rate"] = slope * x + intercept print(f"Apparent Ea from {metric}: {apparent_Ea_kJ_mol:.1f} kJ/mol") print(f"Apparent A from {metric}: {apparent_A_per_ps:.2f} ps^-1 ({apparent_A_per_s:.2e} s^-1)") fig, ax = plt.subplots(figsize=(6.5, 4.2)) ax.scatter(x, y, label="rate proxy") order = np.argsort(x) ax.plot(x[order], arr["fit_ln_rate"].to_numpy()[order], label="linear fit") ax.set_xlabel("1 / T / K$^{-1}$") ax.set_ylabel("ln(proxy count ps$^{-1}$)") ax.set_title("Apparent Arrhenius fit for the scission proxy") ax.legend() fig.tight_layout() ax; ax.figure.savefig("picture6.png") # From the figure, we can see that the points are roughly linear. So, the chosen proxy behaves in a way that is at least consistent with thermally activated behavior over this temperature range. This graph shows how we can use super-fast, high-temperature simulations to predict how a material ages slowly in the real world. By measuring the slope of the solid blue line, we can calculate the activation energy (~75 kJ/mol), think of this as the minimum energy hurdle the reaction must jump over to happen. # ## Continue in AMSmovie and ChemTraYzer2 # # At this point the notebook has prepared trajectories and summary tables, but detailed reaction-event analysis is still best done interactively. A practical next step is to open one of the trajectories in AMSmovie and use the ChemTraYzer2 reaction detection tools there. # # Suggested workflow: # # 1. Open one of the printed `amsmovie` commands. # 2. In AMSmovie, choose **MD Properties -> Reaction Detection: ChemTraYzer2**. # 3. Inspect the detected reactions, populations, and individual reactive frames. # 4. If you need finer reaction-event resolution, rerun the MD with a smaller trajectory sampling frequency.