2026-10-14 MD analysis, phase 1: the shared library#

The campaign’s design is in the SEAMM developer guide: MD analysis: H5MD output, a shared library and analysis steps. This package is its phase 1.

Phase 1#

The code that diffusivity_step and thermal_conductivity_step had each copied and let drift is gathered here, and both steps import it, with no change for users.

Here

Was

stats.fmt_err

fmt_err in both steps. They differed in the decimals used when the error is unusable (3 in Diffusivity, 2 in Thermal Conductivity): now the fallback_decimals argument.

correlate.ccovf, acf_err, ccf_err, moving_average, tensor_correlation_functions

thermal_conductivity_step.analysis: statsmodels’ ccovf, acf_err, ccf_err, moving_average, create_correlation_functions. ccovf uses the same definition and scipy’s FFT correlation, so statsmodels is no longer needed.

transport.collective_helfand_moments

create_helfand_moments in Thermal Conductivity.

transport.species_helfand_moments, msd_by_species

create_helfand_moments and compute_msd in Diffusivity.

transport.fit_green_kubo_integral, fit_helfand_moment, helfand_slope

The same in Thermal Conductivity (get_helfand_slope).

transport.fit_slope

The same in Diffusivity.

transport.yeh_hummer_correction

Inline in Diffusivity (twice).

plots

The plot functions of both steps, now one helper (add_component_traces) with thin wrappers.

legacy.read_trajectory, read_vector_trajectory, read_dump_trajectory

The pandas reader in Thermal Conductivity; the two readers in Diffusivity.

Changes in behaviour#

  • Speed. The Helfand moments and MSDs are computed with FFTs. The old code looped over every time origin in Python, which took O(n m) time. For a heat flux of a million points with m = 10,000 that is about 2 minutes, against 0.4 s now (this Mac). For the Helfand moments of 512 molecules over 1,001 frames it was 3.3 s, against 0.09 s now. The results are the same to rounding: the tests check this against the original loops, and on real LAMMPS output from earlier jobs they agree to about 1e-15 relative.

  • Precision. A flux with a non-zero mean has a running integral that drifts linearly. The straight line is removed before the FFTs and its contribution added back exactly, so the rounding error does not grow with the length of the series.

  • Latent bugs fixed:

    • fit_slope with weights passed sigma[i:] with y[i:j], a shape error whenever the fit ended before the last point. No caller used weights.

    • fit_green_kubo_integral failed without sigma, with a NameError on very short series, and when the sampling was coarser than 0.4 ps (a zero-width moving average; this is what ended the 2025-11-19 argon demo job, Job_003234).

    • read_dump_trajectory failed on a dump with a single timed frame.

    • The MSD trace had units Å instead of Å^2 in its metadata.

  • Defects corrected (Paul, 2026-10-14: fix the older problems in phase 1). These change results:

    • fit_green_kubo_integral accepted any second exponential whose plateau was positive and less than twice the first. On an OU test series it kept one with a negative amplitude and a time constant of 46 ps, beyond the 15 ps of data, and gave 0.018 instead of about 0.5. The second exponential is now kept only if both amplitudes are positive, the two time constants differ by more than 10 %, and three times the slower one lies within the fitted data.

    • Its error was the sum of the errors of the two amplitudes. When the two time constants are nearly equal the amplitudes are individually undetermined and strongly anticorrelated, so the error came out as large as 10^6 for a value of 0.5. The error now includes their covariance.

    • All the fits treated the errors as weights even when they were zero, as they are for a single run (the spread over runs). Thermal Conductivity’s Helfand-moment fit then returned its starting guess, 1.00. A fit with missing, zero or non-finite errors is now unweighted.

    • helfand_slope stored the central difference over points k and k+2 at k+2 instead of k+1, which shifted the curve by one sample, and its error was the error of the moment, not divided by 2 dt. Both are corrected; the slope at the origin is zero, as the running Green–Kubo integral is.

    • The fourth column of species_helfand_moments, the total, was the moment of the summed velocity components, (vx + vy + vz)^2. Its cross terms average to zero only for an isotropic system, and add noise in any case. It is now the sum of the three moments, as for the MSD.

Defects found and corrected in the steps at the same time:

  • LAMMPS Heat Flux sub-step (lammps_step). To remove each atom’s average stress, the callback passed the average stress, negated, as a virial to fix external. compute stress/atom multiplies a virial by -nktv2p (-68,568 in real units), so the stress was changed by +68,568 times its average. For argon this made the heat flux about 10^4 times too large and the thermal conductivity about 4e7 W/m/K; it is now 0.1 W/m/K, as measured. For molecular systems the correction reached neither stress compute, and the convective term was counted twice (both heat/flux computes include it).

  • Diffusivity. The per-direction coefficients Dx, Dy and Dz used the 1/6 (MSD) or 1/3 (Helfand) of the total, so they were a third of the diffusion coefficient, and they got a third of the Yeh–Hummer correction. They are now D_alpha = slope / 2 (MSD) or slope (Helfand moments), and equal D for an isotropic system.

Phase 2: H5MD#

Done on 2026-10-15, following the design (“The data contract” and “Non-periodic systems” in the plan).

  • The library gained h5md (H5MDWriter), io (open_md(), MDData, TimeSeries) and units (H5MD unit strings and pint), and streaming readers for LAMMPS dumps and vector trajectories. The layout is specified in The H5MD file.

  • The LAMMPS step writes md.h5md at the end of each NVE, NVT, NPT and Heat Flux sub-step (lammps_step/md_h5md.py) and appends the run to _md_runs. The native files are kept.

  • Diffusivity and Thermal Conductivity read md.h5md and fall back to the text files of older runs. Both stop before the MD for a non-periodic system.

Decisions taken while implementing:

  • A fixed box is written as time-dependent edges, repeated for each frame of the positions. MDAnalysis reads box/edges/value[frame] and fails on the fixed form that H5MD also allows.

  • When both atomic positions and a trajectory are sampled, the atomic positions (and velocities) go into particles/all, and the trajectory into its own group, particles/trajectory. MDAnalysis reads the first group alphabetically, “all”.

  • The off-diagonal pressure sampled for viscosity (Pxy Pxz Pyz) is written as pressure_offdiagonal (attribute components = “xy xz yz”) until phase 4 replaces it by the full pressure_tensor.

  • Thermal Conductivity averages the temperature over the same samples as the heat flux: the whole run. Before, the temperature came from state.json, which drops the equilibration period detected by pymbar and keeps only uncorrelated samples (494 of 5,001 in an argon test), while the heat flux was used whole. The two were thus taken over different parts of the run. For that argon NVE run the temperature changed by 0.03 K, or 0.05% in the thermal conductivity. Dropping an equilibration period, if wanted later, must cut the heat flux and the temperature at the same step.

Defects found on the way:

  • LAMMPS, non-periodic systems: sampling the trajectory of a molecule failed (“Variable name sxx … does not exist”), since the stress columns were added regardless of periodicity. The pressure was reported though meaningless.

  • Diffusivity: without centre-of-mass data it used the atomic positions or velocities with the molecules’ indices, which is wrong unless every molecule is one atom. It now checks the number of particles.

  • Thermal Conductivity: a missing state.json gave a NameError; the sampling interval was taken as fs even when the file gave ps.

Phase 3: requirements and analysing earlier MD#

Done on 2026-10-15.

  • The library gained requirements (Need, Requirements, Provides, check), select_md_runs, and step_support, the code the analysis steps share. The protocol is described in MD steps and analysis steps.

  • seamm: builder.validate calls each step’s check_flowchart(preceding); the editor validates before submitting. seamm_exec: run_flowchart validates before running any step, except in an iteration of a parallel loop, which its parent job checked. A failing check is logged, never reported, so a bug in a step cannot stop every job.

  • LAMMPS: md_provides and request_sampling on NVE, NVT, NPT and Heat Flux (md_provides.py); before running, the LAMMPS step calls check_configuration on the steps after it up to the next LAMMPS step, following the flowchart’s edges without calling next(), which a Loop may override.

  • Diffusivity and Thermal Conductivity: the md source, md files and number of blocks parameters, both modes, automatic sampling in mode A, and the dialogs hide what does not apply, including the sub-flowchart tab in mode B.

Decisions taken while implementing:

  • Mode B reads the last run that samples the first need, not simply the last run, both when checking the flowchart and when running. A LAMMPS step with an NVT run sampling centres of mass followed by a Heat Flux run then feeds both Diffusivity (the NVT run) and Thermal Conductivity (the Heat Flux run).

  • A single run is divided into blocks (5 by default) analyzed as independent runs for the errors. Thermal Conductivity requires at least 2000 heat-flux samples per block; Diffusivity at least 10 frames.

  • Sources: “the last MD run”, “the last MD run of each loop iteration” (runs by the same step, i.e. the iterations of a Loop), and files. “All runs of the last MD step”, for production split over several sub-steps, waits for the segments of phase 8.

  • Diffusivity’s requirements: positions of the molecules (or atoms, for monatomic systems) at least every 1 ps, turned on at 100 fs; velocities at least every 20 fs, turned on at 10 fs. Thermal Conductivity: the heat flux every two steps. Mode A turns on what is missing in the last MD sub-step of the sub-flowchart and says so.

  • Thermal Conductivity in mode B uses the Green–Kubo approach only (the only one implemented) at the temperature of the MD; the approach and temperature parameters apply only when it runs its own MD.

Found on the way: when the correlation window is shorter than the start of the Helfand-moment fit (at least 1 ps), Thermal Conductivity silently dropped the Helfand Moments rows; it now says why.

Phase 4: the viscosity#

Done on 2026-10-15.

  • The library gained shear_stress_components, autocorrelation_functions, component_helfand_moments, publish_viscosity / preceding_viscosity (the workspace variable _viscosity) and plots.plot_components.

  • LAMMPS: the NVE, NVT and NPT sub-steps sample the full pressure tensor (“Sample the pressure tensor”) instead of its three off-diagonal components. The parameters shear stress* were renamed pressure tensor*; NVE_Parameters translates the old names when a flowchart is read (format 2.0 and 3.0), as NPT already did for “keep orthorhombic”. The sampled tensor is written to md.h5md as pressure_tensor, replacing the coarser one from the thermodynamic state, so pressure_offdiagonal is no longer written.

  • Viscosity (new, viscosity_step): Green-Kubo and Helfand moments of the five shear components and their average, in both modes; needs pressure_tensor at least every 5 fs and a periodic system. It publishes the Green-Kubo average for later steps and stores it as viscosity#Green-Kubo#<model>.

  • Diffusivity: “Viscosity from” chooses the preceding Viscosity step (the default) or a given value for the Yeh-Hummer correction, which is now on by default and skipped with a note when there is no viscosity. Older flowcharts, which always give the viscosity, keep it.

Decisions taken while implementing:

  • The mean of the shear stress is not subtracted. For a fluid it is zero, and subtracting the sample mean of a block of length T lowers the correlation function by about 2 sigma^2 tau / T at every lag. With synthetic stresses of known viscosity (0.9 mPa s, 200 ps blocks) the Green-Kubo fit averaged 0.837 with the mean subtracted and 0.878 without, over eight seeds.

  • The average is fitted, not averaged. The Green-Kubo and Helfand values of the average are fits to the curve averaged over the five components, which is less noisy than any one of them; its error comes from the spread of the fits to each run or block.

  • The Viscosity step resets its data at the start of each run, so it gives the right answer inside a loop. Thermal Conductivity does not, which is worth fixing.