Resin Blend from Molecular Dynamics

ReaxFF Molecular Dynamics Polymers RDF Packmol

This tutorial will teach you how to:

• build a small polyisoprene/H-DCPD resin blend based on Ref. [1];

• equilibrate a pseudo-2D bulk polymer cell with UFF and ReaxFF;

• monitor the density during annealing;

• analyze the miscibility of the resin.

Introduction

Rissanou et al. [1] investigated hydrogenated dicyclopentadiene (H-DCPD) resin blended with polyisoprene (PI) and SBR. For the PI system at high resin loading, the paper used 100 PI chains with 232 H-DCPD resin molecules. This tutorial scales that composition down to 20 PI chains and 46 H-DCPD molecules, which gives a system of about 10,000 atoms.

The goal is not to reproduce the statistical quality of the production paper. Instead, the tutorial reproduces the workflow: packing the blend, compressing the low-density starting structure, annealing the material, and analyzing resin-rich structure.

You can run this tutorial with Python from the companion notebook.

Downloads

Use these structures as inputs for the tutorial:

• PI_30mer.xyz

• HDCPD_trimer.xyz

Note

Refer to Building Polymers to build the PI chain. You can define the cis-1,4-isoprene molecule, define the anchors, and generate a 30-mer.

Part A: Build the Blend

Create the packed PI/H-DCPD structure

1. Start AMSinput: SCM → New Input

2. Switch to the ForceField panel

3. Choose Builders → Packed Molecules…

4. Set Density to 0.05 g/cm ³

5. Tick Use number of molecules

6. In the first SMILES or file row, use the folder button to select PI_30mer.xyz

7. Set N mols for PI to 20 and Region to polymer

8. Click the add-row button

9. In the second row, use the folder button to select HDCPD_trimer.xyz

10. Set N mols for H-DCPD to 46 and Region to resin

11. Click Generate molecules

12. Close the Packed Molecules builder

Note

The paper used 100 PI chains and 232 resin molecules. The 20:46 system used here preserves the same approximate composition while keeping the tutorial small enough for a local workstation.

Part B: Equilibrate the Polymer Bulk

Run a UFF geometry relaxation

1. In AMSinput, keep the ForceField engine

2. Set Task to Geometry Optimization

3. Save the job as UFF_InitialRelax

4. Run the calculation.

Pre-shrink to a pseudo-2D bulk cell

We will intentionally create a pseudo-2D periodic bulk cell with a' = b' = λc' (λ > 1). This makes the later resin-density map easier to interpret while still using 3D periodic boundary conditions.

The target volume V as a function of the initial volume, initial density, and target density is therefore:

\[V = \rho_0 V_0 / \rho\]

We can then deduce the target cell dimensions as a function of the c-scaling factor λ:

\[\begin{split}&a' = b' = \left(V\lambda\right)^\left(1/3\right) \\ &c' = a' / \lambda\end{split}\]

We will use λ = 1.5.

Tip

If you open UFF_InitialRelax with AMSinput, you can easily read the cell volume and density from Properties → Properties (Including Estimated).

1. Re-open UFF_InitialRelax and Update Molecule

2. Keep the ForceField engine and set Task to Molecular Dynamics

3. In Model → MD, set Number of steps to 20000

4. Keep Time step at 0.25 fs

5. Set Sampling frequency to 1000

6. Set Initial velocities temperature to 5 K

7. Add a Berendsen Thermostat, set the Temperature to 5 K, and the Damping Constant to 10 fs

8. Switch to Model → MD → Deformation and set the target lattice with diagonal values 53 53 35

9. Save as UFF_Shrink

10. Run the calculation

The calculation takes about 20 minutes on an Apple M2 with 16 GB of memory.

Run a ReaxFF relaxation

1. Open the final structure from UFF_Shrink in AMSinput

2. Switch to the ReaxFF engine

3. Set Force field to dispersion/CHONSSi-lg.ff

4. Set Task to Geometry Optimization

5. Save as ReaxFF_InitialRelax

6. Run the calculation

Note

If the geometry optimization reaches a plateau with small gradients, you can request an early stop in AMSjobs by selecting Jobs → Request Early Stop. Tight convergence is not required here because an MD simulation will be run immediately afterward. However, do not skip this step: the forces should still be minimized before starting the MD simulation to help prevent the molecule from decomposing.

Part C: Anneal and Produce a Trajectory

Run ReaxFF annealing at 413 K

1. Open the final structure from ReaxFF_InitialRelax in AMSinput

2. Keep the ReaxFF engine and dispersion/CHONSSi-lg.ff force field

3. Set Task to Molecular Dynamics

4. Set Number of steps to 25000

5. Set Sampling frequency to 500

6. Set Initial temperature to 413 K

7. Set Thermostat to Berendsen, Temperature to 413 K, and Tau to 100 fs

8. Set Barostat to Berendsen, Pressure to 101325 Pa, and Tau to 500 fs

9. Make sure deformation is disabled in Model → MD → Deformation

10. Save as ReaxFF_Anneal_1 and run

Repeat the procedure from the final structure (update the molecule and save) until the last-third average density changes by less than about 0.5 percent.

Tip

To get the density, open the job in AMSmovie, Graph → Add Graph, plot the density with MD Properties → Density, and save it to a file with Graph → Save as XY. Then you can get the last-third average from a terminal (from the same folder as the data file):

awk 'NR==FNR && FNR>1 {total++; next} FNR>1 {i++; if (i >= int(2*total/3)+1) {sum+=$2; n++}} END {print sum/n}' rho_1.xy rho_1.xy

Table 12 Summary of the densities (in g/cm ³) obtained from the successive annealing procedures.

Annealing Step	Initial Density	Final Density	Last Third Average Density	% Change
1	0.8988	0.7290	0.7168	_
2	0.7290	0.7962	0.7917	9.46
3	0.7962	0.8377	0.8362	5.32
4	0.8377	0.8595	0.8580	2.54
5	0.8595	0.8744	0.8704	1.42
6	0.8744	0.8893	0.8871	1.88
7	0.8894	0.8932	0.8920	0.55

Run a production trajectory

1. Open the final annealed structure in AMSinput

2. Keep ReaxFF with dispersion/CHONSSi-lg.ff

3. Set Task to Molecular Dynamics

4. Set Number of steps to 50000

5. Set Sampling frequency to 100

6. Set Initial velocities temperature to 413 K

7. Use a Berendsen thermostat at 413 K with Tau 100 fs

8. Remove the barostat

9. Save as ReaxFF_Production and run

Part D: Analyze the Trajectory

Radial distribution function

To understand the miscibility of the resin in the polymer matrix, we compute and report the resin-resin and resin-polymer center-of-mass RDFs. You can use this PLAMS script and simply modify the path so the trajectory points to the ams.rkf of the final production run.

Density map

The miscibility can also be appreciated by plotting the 2D density map of the resin atoms. You can use this PLAMS script to generate and plot the density map.

Note

To get better statistics, you should run a longer production trajectory, increase the sampling frequency, or perform the simulations with a larger system size.

Summary

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.

References

[1] (1,2)

1. 14. Rissanou, R. Ghanta, M. Doxastakis, V. Harmandaris, Investigating the Role of Miscibility in Hydrogenated Dicyclopentadiene Resin/Polymer Blends: A Molecular Dynamics Study. Polymers, 18, 594 (2026)