md-analysis-planner

Files (5)

• CHANGELOG.md
• SKILL.md
• evals.json · evals/
• md_analysis_checks.md · references/
• md_analysis_planner.py · scripts/

SKILL.md (first lines)

# MD Analysis Planner

## Goal

Choose the right MD trajectory analyses and prerequisites before writing post-processing code.

## Requirements

- Python 3.10+
- No external dependencies
- Works on Linux, macOS, and Windows

## Inputs to Gather

| Input | Description | Example |
|-------|-------------|---------|
| System | Material or molecular system | `oxide glass`, `liquid water` |
| Goals | Analysis goals | `rdf,diffusion,coordination` |
| Trajectory format | Dump, DCD, XYZ, H5MD, etc. | `LAMMPS dump` |
| Velocities | Whether velocities are stored | `true` |
| Stress | Whether stress/virial is stored | `true` |
| Unwrap needed | Whether atoms cross PBC | `true` |
| Timestep | fs per saved frame | `10` |

## Decision Guidance

- Use **RDF and coordination** for local structure.
- Use **MSD** for diffusion, but unwrap trajectories and verify diffusive regime.
- Use **VACF/VDOS** only when velocities or reliable finite-difference velocities exist.
- Use **stress-strain** only if stress/virial and deformation history are available.
- Always perform equilibration checks before fitting transport or thermodynamic properties.

## Script Outputs

`scripts/md_analysis_planner.py` emits these fields (JSON under `results`):

| Field | Description |
|-------|-------------|
| `analysis_plan` | One entry per goal: `goal`, `method`, and `status` |
| `required_data` | Sorted, de-duplicated data needed across all goals |
| `equilibration_checks` | Standard pre-fit equilibration checklist |
| `pbc_handling` | `unwrap_needed`, `minimum_action`, `format_note` |
| `warnings` | Safety-critical caveats and blockers |

### Status values

`status` is one of, from most to least severe: `blocked` > `needs time axis` >
`needs review` > `ready`. A more-severe status is never demoted by a less-severe
one. For example, a VACF/VDOS goal with no stored velocities reports `blocked`
even when the timestep is also missing (both warnings are still emitted).

### Default (non-JSON) output

Without `--json` the script prints the plan lines, a `Required data:` list, a
one-line `PBC:` note, and a `Warnings:` section (also mirrored to stderr) so the
safety-critical caveats are visible even when stdout is piped.

## Workflow

```bash
python3 skills/simulation-workflow/md-analysis-planner/scripts/md_analysis_planner.py \
  --system "oxide glass" \
  --goals rdf,coordination,bond-angle \
  --trajectory-format dump \
  --unwrap-needed \
  --timestep-fs 10 \
  --json
```

## Error Handling

If velocities, stress, or timestep information is missing, downgrade dependent
analyses and report warnings. The script exits with code `2` and a message on
stderr for invalid input (empty system, no goals, non-positive or non-finite
timestep, or inputs exceeding the size caps below).

## Limitations

This skill plans analysis and prerequisites; it does not parse large trajectories directly.

## Verification checklist

Do not trust trajectory-derived results until each applicable item below is
recorded against the planner's own output:

- [ ] Ran `scripts/md_analysis_planner.py` and confirmed no `analysis_plan` entry is `blocked` or `needs time axis`; record any `needs review` goal and how its custom analysis was resolved.
- [ ] Confirmed every item in the planner's `required_data` list is actually present in the trajectory (e.g. `velocities` for VACF/VDOS, `stress or virial` + `strain history` for stress-strain) before running the corresponding analysis.
- [ ] For diffusion/MSD: recorded the log-log MSD-vs-time slope and confirmed it is ~1 (diffusive regime) over the fit window, excluding ballistic/sub-diffusive transients, before quoting D = lim MSD/(2*d*t).
- [ ] For diffusion/MSD: assessed whether a hydrodynamic finite-size correction is applicable. For isotropic bulk fluids in 3D periodic cubic cells, record L, T and shear viscosity eta if applying Yeh-Hummer (D_0 = D_PBC + k_B*T*xi/(6*pi*eta*L), xi~=2.837). Report the measured D_PBC, the correction and its assumptions separately. Do not apply this formula automatically to solids, confined or anisotropic systems, non-cubic cells or ballistic trajectories.
- [ ] For any displacement-based analysis when `pbc_handling.unwrap_needed` is true: confirmed positions were unwrapped using `cell` + `image flags` before computing displacements (not raw wrapped coordinates).
- [ ] Worked through the planner's `equilibration_checks`: discarded the startup transient, confirmed temperature/pressure plateaus, and compared first-half vs second-half property estimates before any transport/thermodynamic fit.
- [ ] Reported an uncertainty from block averaging (or independent trajectories) for every quoted transport/thermodynamic property, not a single-window point value.

## Common pitfalls & rationalizations

| Tempting shortcut | Why it's wrong / what to do |
|-------------------|-----------------------------|
| "MSD vs time looks straight, so it's diffusive — just fit a line." | A straight-looking MSD can still include ballistic/sub-diffusive transients. Confirm the log-log slope is ~1 over the fit window first, then fit D = lim MSD/(2*d*t). |
| "D from the simulation is the diffusion coefficient." | Assess size dependence and report D_PBC. The Yeh-Hummer 1/L hydrodynamic correction applies under bulk-fluid, 3D periodic cubic-cell assumptions; justify those assumptions and the viscosity before separately reporting a corrected D_0. |
| "Positions are in the dump, so MSD is fine." | Wrapped coordinates make atoms jump across the box and corrupt displacements. When `unwrap_needed` is true, unwrap using `cell` + `image flags` before any displacement analysis. |
| "No velocities stored, but I can get VACF/VDOS from positions." | The planner marks VACF/VDOS `blocked` without velocities for a reason. A finite-difference velocity estimate must be explicitly justified (frame spacing, aliasing); otherwise the result is unreliable. |
| "The run finished, so I can fit transport properties on the whole trajectory." | Run completion is not equilibration. Discard the startup transient and verify temperature/pressure plateaus (the `equilibration_checks`) before fitting. |
| "One window gave a clean number, no need for error bars." | A single window hides correlation-driven variance. Use block averaging or independent trajectories so every transport/thermodynamic value carries an uncertainty. |

## Security

### Input Validation

Inputs are scalar CLI values and booleans only. `plan_md_analysis` validates and
bounds every field before use; any violation raises `ValueError`, which `main`
catches to print a message on stderr and exit with code `2`:

- `system` must be non-empty (after stripping) and at most 256 characters.
- At most 64 goals are allowed; each goal at most 256 characters.
- `trajectory_format` at most 256 characters.
- `timestep_fs`, if given, must be a positive, finite number (`math.isfinite`
  and `> 0`).

Goals are not allowlisted: an unrecognized goal is not rejected but is reported

From HeshamFS/materials-simulation-skills, Apache-2.0.

More from HeshamFS/materials-simulation-skills