Molecular dynamics · LAMMPS · OPLS-AA

Polyethylene and polypropylene have almost the same solubility parameter. They still don't mix.

All-atom simulations of a PE melt, a PP melt and a 50/50 blend, built from SMILES and run to 500 ps of production dynamics. The measured solubility parameters differ by less than half a percent — which, on the classical argument, would make them freely miscible. Every recycler knows they are not. The gap between those two statements is the whole problem with predicting blend behaviour.

Force field: OPLS-AA, assigned with Foyer
Engine: LAMMPS, real units, atom_style full
Conditions: 400 K, 1 atm, PPPM electrostatics
Protocol: 400 ps NPT → 100 ps volume average → 500 ps NVT
Observable: Cohesive energy density from intermolecular energy
Frames analysed: 501 per system

The blend

Thirty polyethylene chains and thirty polypropylene chains in a periodic cube, after 400 ps of constant-pressure equilibration and 500 ps of production dynamics at 400 K. This is the final configuration, rendered in your browser from the simulation output. Use the buttons to show one component at a time — the question you are looking at is whether the two have started to separate.

Loading structure…

PE/PP 50/50 blend, 2,370 atoms, 29.67 Å cube. Carbon backbones only — hydrogens are simulated but hidden, because at this density they obscure the chain paths. Isolating one component is the fastest way to see whether domains have formed.

A caution on reading it: at 5 repeat units per chain and a 30 Å box, this system is far too small and far too short-chained to develop the domain structure a real blend forms. Absence of visible separation here is not evidence of miscibility — it is a consequence of the model. More on that below.

The two pure melts

Each component is simulated separately under identical conditions, because the cohesive energy density of the blend only means something next to the two references.

Loading structure…

Pure PE melt — 60 chains, 1,920 atoms, 28.08 Å.

Loading structure…

Pure PP melt — 60 chains, 2,820 atoms, 31.15 Å.

What was measured

Cohesive energy density is the energy needed to pull a condensed phase apart into isolated chains, per unit volume. It is obtained here by computing the intermolecular non-bonded energy of every stored frame — the energy of the periodic melt minus the energy of the same chains, in the same conformations, placed in isolation. Its square root is the Hildebrand solubility parameter.

System	Chains	Atoms	Box (Å)	ρ (g cm⁻³)	CED (J cm⁻³)	δ (J cm⁻³)^½
PE melt	60 × C₁₀H₂₂	1,920	28.08	0.640	188.20 ± 0.15	13.72
PP melt	60 × C₁₅H₃₂	2,820	31.15	0.700	189.74 ± 0.08	13.78
50/50 blend	30 + 30	2,370	29.67	0.668	190.08 ± 0.13	13.79

Cohesive energy density during NVT production

J cm⁻³ · 10 ps block means

Read the vertical axis before anything else. Individual 10 ps block means scatter by one to two J cm⁻³, while the three systems are separated by less than that — the standard error on each run mean is 0.08–0.15 J cm⁻³ only because 500 ps averages a great deal of that scatter away. The quantity that determines whether these two polymers mix is smaller still, which is the honest reason it is not quoted on this page.

The melt densities, 0.640 and 0.700 g cm⁻³ for PE and PP at 400 K, are where they should be for short-chain polyolefin melts, and the blend volume is within 0.2% of the sum of the pure-component volumes. The energy bookkeeping closes.

Why the solubility parameters do not settle it

δ_PE = 13.72 and δ_PP = 13.78 (J cm⁻³)^½. On the regular-solution argument, the enthalpy of mixing goes as (δ_A − δ_B)², and a difference of 0.06 gives a number indistinguishable from zero. Take that at face value and PE and PP should be as miscible as two batches of the same grade.

They are not, and the reason is not enthalpy. It is that the entropy which normally drives mixing disappears in a polymer. In Flory–Huggins the combinatorial entropy per lattice site carries a factor 1/N, so for a symmetric blend the critical interaction parameter is

χ_c = 2 / N

For the five-unit oligomers simulated here, χ_c = 0.4 — an enormous tolerance. Decane and pentadecane really do mix. For a commercial polyolefin with N in the thousands, χ_c falls to around 0.002. A χ small enough to be invisible in a solubility-parameter comparison is then more than enough to drive full phase separation.

That is the practical lesson for anyone assessing a contaminated recyclate:the two polymers being chemically similar is not reassurance. PE and PP are about as similar as two different polymers can be, and a PP part with PE in it still fails differently from a clean one.

What this run does not yet tell you

The obvious next step is to extract χ itself from these three cohesive energy densities. I am not publishing that number yet, and the reason is worth stating plainly.

The signal is about half a percent of the measurement. The mixing energy that sets χ is a difference of order 1 J cm⁻³ between quantities of order 190 J cm⁻³. Getting χ to the precision that matters for a polymer — around 0.001 — means resolving the CEDs to roughly 0.01 J cm⁻³.
One replica of 500 ps does not reach that. Splitting the production run in half gives two answers that disagree by more than the quoted statistical error, which is the clearest possible sign that the frames are more correlated than the block analysis assumes and the blend is still relaxing.
Packing can masquerade as chemistry. The blend equilibrated 0.2% denser than the ideal mixture. Since CED scales with density, a systematic error of that size accounts for a good fraction of any apparent mixing energy before any chemistry is invoked.
Five-mers are not polymers. The result is an oligomer χ. Chain ends are a large fraction of the system, and PP tacticity is not controlled by the builder.

What would fix it: three to five independent replicas per system with different packing and velocity seeds, production runs extended into nanoseconds, and longer chains. That is queued. When the two halves of a run agree and the replicas agree with each other, the number goes up here with its error bar.

Publishing the pipeline and the converged quantities, while holding back the one that is not converged, is not modesty. It is the difference between a measurement and a guess with decimal places.

Method

Chain building — repeat units defined as SMILES with placeholder atoms, joined into five-unit oligomers with RDKit and PySoftK, MMFF-relaxed and end-capped, giving C₁₀H₂₂ (PE) and a branched C₁₅H₃₂ (PP).
Typing — OPLS-AA applied through Foyer's SMARTS rules; parameters carried in a ParmEd structure and written directly to LAMMPS data files, with Ryckaert–Bellemans torsions converted to the OPLS Fourier form.
Packing — Packmol at 0.5 g cm⁻³ with a 2 Å minimum separation, deliberately below melt density so packing always converges and NPT does the compression.
Dynamics — LAMMPS, 1 fs timestep, Nosé–Hoover thermostat and barostat, lj/cut/coul/long at 10 Å with analytic tail correction, PPPM to 10⁻⁴. Minimisation, then 400 ps NPT, then 100 ps of volume averaging, then the box rescaled to the mean volume, then 500 ps NVT with frames every 1 ps.
Energy decomposition — each production frame re-evaluated twice with the LAMMPS rerun command: once as the periodic melt, once with the same chains made whole and isolated. The difference is the intermolecular energy, and bonded terms cancel exactly. This two-pass scheme avoids the intramolecular contamination that a cutoff Coulomb sum with neighbour exclusions would introduce.
Analysis — per-frame CED, block averaging over five contiguous blocks for the standard error, and error propagation into the derived quantities.

Dealing with mixed polyolefin streams? The practical question is rarely what χ is — it is how much of the wrong polymer is in your feed and what it costs you. That is measurable, and it is the subject of the companion article.

Get in touch