Getting Started

seamm_md_analysis is the shared library behind SEAMM’s analysis of molecular dynamics: the Diffusivity, Thermal Conductivity and Viscosity steps now, and Structure and the Dielectric constant later. The steps run the MD (or read an earlier run), and this library does the numerical work.

Installing

pip install seamm-md-analysis

It needs numpy, scipy and pymbar.

What is in it

seamm_md_analysis.stats

mean_and_stderr over replicas, block_average, statistical_inefficiency and detect_equilibration (pymbar, imported quietly when first used), mean_and_stderr_correlated, and fmt_err to print a value with its error.

seamm_md_analysis.correlate

ccovf, the cross-covariance of two series by FFT (the definition of statsmodels), the errors of correlation functions, and tensor_correlation_functions for the six components of a vector flux with their running integrals (the Green–Kubo input), and autocorrelation_functions for independent series such as the components of the shear stress.

seamm_md_analysis.transport

collective_helfand_moments for a flux such as the heat flux; shear_stress_components, the five independent shear components of a pressure tensor, and component_helfand_moments for each of them, for the viscosity; species_helfand_moments and msd_by_species for single-molecule diffusion; fit_green_kubo_integral, fit_helfand_moment, helfand_slope and fit_slope; and yeh_hummer_correction.

The Helfand moments and MSDs are computed with FFTs. They equal the direct sums over time origins to rounding, but take O(n log n) rather than O(n m) time.

seamm_md_analysis.electric

The static dielectric constant and the ionic conductivity from the trajectory and fixed charges: classify_molecules (neutral molecules, ions, and molecules with a non-integer charge), cell_dipole, dielectric_constant (the fluctuation formula for Ewald or PPPM electrostatics with conducting boundaries, with errors from blocks and the running estimate) and dipole_correlation; collective_msd and conductivity_from_msd (Einstein–Helfand, from the collective displacement of the ions’ charges) and nernst_einstein. For an electrolyte the dielectric constant uses the dipole of the neutral molecules: the ions carry a current, not a bounded dipole.

seamm_md_analysis.plots

Plots of correlation functions, Green–Kubo integrals, Helfand moments and MSDs, for seamm_util figures.

seamm_md_analysis.legacy

Readers for the files the LAMMPS step writes today: !MolSSI trajectory time series (read_trajectory), per-molecule vector_trajectory files (read_vector_trajectory) and dump files (read_dump_trajectory).

Example

The thermal conductivity from a heat flux sampled every dt:

import numpy as np
from seamm_md_analysis import (
    collective_helfand_moments,
    fit_green_kubo_integral,
    read_trajectory,
    tensor_correlation_functions,
)

metadata, columns = read_trajectory("heat_flux.trj")
J = np.stack((columns["Jx"], columns["Jy"], columns["Jz"]))
m = min(J.shape[1] // 10, 10000)

# Correlation functions and their integrals, per sampling interval
ccf, integral = tensor_correlation_functions(J, m=m)

# Helfand moments, per sampling interval squared
M = collective_helfand_moments(J, m=m)

Multiply by V / (k_B T²) and the sampling interval (Green–Kubo), or V dt² / (2 k_B T²) (Helfand), to get W/m/K.

The shear viscosity from a pressure tensor P (n × 3 × 3) works the same way, averaging the five components:

from seamm_md_analysis import (
    autocorrelation_functions,
    component_helfand_moments,
    shear_stress_components,
)

X = shear_stress_components(P)  # Pxy, Pxz, Pyz, (Pxx-Pyy)/2, (Pyy-Pzz)/2
acf, integral = autocorrelation_functions(X, m=m)
M = component_helfand_moments(X, m=m)

Multiply by V / (k_B T) and the sampling interval (Green–Kubo), or V dt² / (2 k_B T) (Helfand), to get the viscosity.