Use the ams2026 skill! 

# Glycine diffusion with ReaxFF at 380 K

## Summary

Build a periodic cubic system containing one neutral glycine molecule
(non-zwitterionic) and 40 water molecules at 1.0 g/cm3. 

Use ReaxFF Glycine.ff for iterative NVT equilibration, followed by a separate
NVT production trajectory initialized from the final equilibrated coordinates
and velocities. Report the glycine center-of-mass diffusion coefficient from
its MSD.

Run all ReaxFF simulations with nproc=1 to disable MPI (it will automatically
use OpenMP).

## Implementation

Be verbose with `plams.log("message")` lines, especially when testing
equilibration and production convergence (criteria vs measured values,
failure/success, ...).

Make the equilibration and production scripts restartable from themselves, in
case there are errors in the scripts, so that any completed MD simulations can
be reused.

### Script 1: initial system

- Pack the molecules, retain explicit glycine and water regions, remove bad
  contacts with a ReaxFF geometry optimization (Convergence.Quality basic).

### Script 2: equilibration

- Timestep 0.5 fs, sample every 100 fs
- Equilibrate at 380 K with a Berendsen thermostat and 100 fs coupling time.
  Extend in 5 ps increments with Restart and CopyRestartTrajectory Yes, up to
  100 ps.
- After at least 10 ps, examine the trailing half of the accumulated
  trajectory. Divide it into five blocks and fit block-mean potential energy
  against time. Accept equilibration after two consecutive checks where the
  fitted change across the window is at most 0.5 trajectory standard
  deviations. Also require the mean temperature to be within 2% of 380 K.

### Script 3: production

- Timestep 0.5 fs, sample every 100 fs
- Start a production NVT job from the final equilibrated structure and
  velocities without copying equilibration frames. Use an NHC thermostat with
  100 fs coupling.
- Run at least 100 ps of production, then extend in 50 ps increments using
  CopyRestartTrajectory, with a 1 ns cap.
- Use the `UseMolecularCentersOfMass` input block/option to trajectory analysis
  to get center-of-mass diffusion. 
- Calculate the MSD in a 20 ps window starting the fit at 5 ps.
- Declare convergence when the last three cumulative diffusion estimates each
  lie within 15% of their mean. If the 1 ns cap is reached, report the best
  estimate as unconverged.
- The MSD calculation center-of-mass fails on glycine fragmentation or atom
  exchange. If this happens, do not equilibrate further but report the final
  results and the fragmentation status.

The production simulation may take a long time to finish. Once everything seems
to be working correctly, you can monitor/poll the simulation less frequently.

### Reporting script

- Produce tables of convergence history, potential-energy and temperature
  traces, MSD and fitted-region figures, and an image of the packed system.
- Include the final diffusion coefficient in m2/s (exponential notation), the
  fit interval, trajectory length, convergence status, observed temperature and
  pressure, and the complete AMS inputs.
- Note the finite-size limitation from using only 40 water molecules and the
  fixed-density NVT state.

## Assumptions

- “Diffusion coefficient” means translational diffusion of the glycine
  molecular center of mass.
- The initial glycine form is neutral NH2-CH2-COOH.
- No finite-size correction or independent replicate trajectories are
  requested.