The density was right to three decimal places. The stiffness was wrong by a factor of fifteen.

The cell

Fifty polyethylene chains of twenty carbons each, packed into a periodic cube and relaxed. This is the actual structure the elastic constants were computed on — 3,100 atoms, exported from the simulation and rendered here in your browser. Chains extend past the box edge because coordinates are unwrapped: the cell is periodic, so a chain leaving one face re-enters the opposite one.

What it got right

The equilibrated cell settled at a density of 0.854 g/cm³. Measured amorphous polyethylene sits at roughly 0.85 g/cm³. That is not a lucky coincidence — it is the kind of property atomistic simulation is genuinely good at, because density is an equilibrium packing property. It depends on how molecules sit next to each other, the forcefield was parameterised against exactly that, and it converges in tens of picoseconds.

Worth holding onto, because it is the reason the rest of the result is seductive. A simulation that nails one number invites you to trust the next one.

What it got wrong

The full stiffness tensor, from a constant-strain calculation at ten strain magnitudes up to 2.5%:

Reducing it: Hill-averaged bulk modulus 16.0 GPa, shear modulus 9.4 GPa, and directional Young's moduli of 13.9, 25.3 and 23.5 GPaalong X, Y and Z.

Look at those three numbers again. This is an amorphous cell — no crystallinity, no orientation, no processing history. It should be isotropic, and the three should agree. They differ by 80%. The anisotropy is not physics; it is the fingerprint of a single configuration that has not relaxed.

Why — three reasons, all visible in the data

1. The cell was never equilibrated, and the temperature hid it

The usual equilibration check is temperature, and temperature looks perfect. It starts wild — a standard deviation of 70 K over the first 2 ps, swinging between 199 K and 396 K as the thermostat fights random initial velocities — and settles to 297.8 ± 4.3 K over the final 100 ps. Textbook.

Now the potential energy over the same run. It is still falling when the simulation stops. Between the 50–100 ps window and the 100–200 ps window the mean drops a further 67.8 kcal/mol, about 3.5%, with no plateau anywhere in sight.

2. The chains are far too short to entangle

C₂₀H₄₂ is 282 g/mol. The entanglement molecular weight of polyethylene is roughly 1,000–1,400 g/mol. These chains are four to five times below it, which means the cell cannot contain an entanglement network — not because it was not equilibrated long enough, but because the molecules are not long enough for one to exist.

Entanglement is most of what governs the mechanical response of real polyethylene. Simulating twenty-carbon chains and reporting a modulus is closer to measuring a wax than a thermoplastic.

3. Constant-strain elastic constants measure the unrelaxed response

The method applies a strain, geometry-optimises the atoms into the deformed box, and reads the stress. There is no time axis in that procedure at all — it is effectively an infinite strain rate with no segmental relaxation permitted.

Polyethylene at room temperature is roughly 140 °C above its glass transition. It is viscoelastic, and its stiffness depends entirely on how long you give it to relax. Measuring it with relaxation switched off returns the glassy, instantaneous modulus — which is genuinely much higher than the engineering modulus, and is not the number anyone specifies against.

What I would change

• Chains above the entanglement length — C₂₀₀ or longer at atomistic resolution, or a coarse-grained bead model if the timescales make atomistic impractical, which they usually do.
• Equilibrate against potential energy, not temperature, under NPT, and run until the block means genuinely plateau. On this system that is nanoseconds, not 200 picoseconds.
• Five to ten independent cells, averaged, with the scatter reported. One cell gives you a number with no error bar, and the 80% anisotropy above shows what that number is worth.
• A finite-rate deformation protocol — or a stress-relaxation or oscillatory calculation — so the result carries a strain rate and can be compared with a real DMA measurement.
• Anchor to one measured value before trusting any prediction.The density agreement here was the anchor that existed; there was no mechanical anchor, and that is precisely where it went wrong.

The point

I publish this one because the failure is more useful than a success would have been. It is a clean demonstration of something I tell clients constantly: simulation output is a hypothesis, not a measurement, and a number that arrives with no uncertainty and no validation against something you actually measured is not a result.

It is also the reason the first thing I do on any engagement is ask what has been measured, not what has been modelled. Models are cheap. Anchors are not.