Regression: SPC/E water#

The tests in tests/ check the analysis on synthetic data with a known viscosity. They cannot catch problems in the MD itself, such as the wrong constraints, or in the chain from the LAMMPS step through md.h5md to the Viscosity and Diffusivity steps. This regression does, by computing the viscosity and diffusion coefficient of rigid SPC/E water and comparing them with published values. It takes too long for CI, so it is run by hand: after changes to this step, to Diffusivity, to seamm_md_analysis, or to how the LAMMPS step samples the pressure tensor or applies constraints.

The files are in regression/spce/:

spce.yaml

The flowchart, as a specification for seamm-flowchart build.

spce.flow

The flowchart built from it, ready to run.

check.py

Reads the results from job.out and checks them against the expected ranges.

The calculation#

  • 512 SPC/E waters (nacl_water.frc), rigid, held by SHAKE with a tolerance of 1e-6, in a cubic cell at 0.997 g/mL (24.86 Å), packed by Packmol.

  • No minimization: the hydrogens of SPC/E have no Lennard-Jones repulsion, and without constraints the minimizer can collapse them onto oxygens.

  • LAMMPS: PPPM (1e-5), a 10 Å cutoff with the long-range dispersion correction, Nosé– Hoover NVT at 298.15 K, 2 fs timestep. 100 ps of equilibration, then 3 ns of production sampling the pressure tensor every 4 fs and the centres of mass every 1 ps.

  • Viscosity and then Diffusivity (mean square displacement of the molecules, with the Yeh–Hummer correction from the Viscosity step) analyze the production run in 5 blocks of 600 ps.

It takes about 35 minutes on 4 cores of an Apple M5 laptop, and 50 on an M3 Pro.

Running it#

From an empty directory:

run_flowchart /path/to/viscosity_step/regression/spce/spce.flow
python /path/to/viscosity_step/regression/spce/check.py .

check.py takes one or more job directories (or job.out files); with several it also checks their average, which reduces the noise by the square root of the number of runs.

Expected values#

The ranges are the references widened by the 95% error of a single 3 ns run, about 12% for the viscosity:

Quantity

Accepted

References

Viscosity, Green–Kubo average (mPa s)

0.62–0.82

0.72 ± 0.02 (our LAMMPS baseline, 3 × 1.7 ns, 1 fs); 0.68 ± 0.02 (Tazi et al., J. Phys.: Condens. Matter 24, 284117 (2012))

D0, corrected (10-9 m2/s)

2.75–3.05

2.86 (baseline; 2.51 ± 0.05 before the correction); 2.97 ± 0.05 (Tazi et al.)

The baseline was computed outside SEAMM, directly with LAMMPS, at the NPT density of the model (0.9995 g/mL) with a 1 fs timestep. The small differences in density and timestep do not matter at this precision.

Results#

Date

Versions, machine

η, GK (mPa s)

η, Helfand

D0

2026-10-15

lammps-step 2026.10.15.4, viscosity-step 2026.10.15; Apple M5

0.716 ± 0.085

0.72 ± 0.12

2.91 ± 0.08

2026-10-15

the same; Apple M3 Pro, a second run

0.725 ± 0.036

0.698 ± 0.087

2.80 ± 0.08

2026-10-15

the average of the two runs

0.720

2.86

The five shear components gave 0.74, 0.715, 0.76, 0.635 and 0.73 mPa s in the first run and 0.79, 0.77, 0.709, 0.789 and 0.61 in the second; a spread of this size is normal for runs of this length. The average of the two runs reproduces the baseline, 0.72 mPa s and 2.86 × 10-9 m2/s.

History#

The first run of this regression, with lammps-step 2026.10.15.3, failed: η came out 1.42 mPa s, and D in y and z was up to 60 times too large. Rigid water was then held by LAMMPS’s RATTLE, which for this system gave a pressure about 850 atm too high, Pxx about 150 atm above Pyy and Pzz for the whole run, and a center of mass that drifted ever faster. Because the mean of the shear stress is not subtracted (it is zero for a fluid), the persistent (Pxx - Pyy)/2 integrated into a viscosity of 6.3 for that component. lammps-step 2026.10.15.4 uses SHAKE instead, and the regression passes.