Metadata-Version: 2.4
Name: pymdkit
Version: 1.6.4
Summary: A unified command-line toolkit for atomistic / MD structure workflows.
Author-email: Yueda Wang <ydwang0608@ustc.edu.cn>
License-Expression: GPL-3.0-or-later
Keywords: materials-science,molecular-dynamics,vasp,gpumd,ase,pymatgen
Classifier: Programming Language :: Python :: 3
Classifier: Operating System :: OS Independent
Classifier: Intended Audience :: Science/Research
Classifier: Topic :: Scientific/Engineering :: Chemistry
Classifier: Topic :: Scientific/Engineering :: Physics
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy
Requires-Dist: ase
Requires-Dist: pymatgen
Requires-Dist: mp-api
Requires-Dist: py4vasp-core
Requires-Dist: dpdata
Requires-Dist: matplotlib
Requires-Dist: umap-learn
Provides-Extra: gpu-cuda12
Requires-Dist: cupy-cuda12x; extra == "gpu-cuda12"
Dynamic: license-file

# pymdkit

A single command-line tool that bundles a collection of atomistic / molecular-dynamics
structure scripts behind one executable: `pmk`. Instead of copying individual
scripts into each working folder and running `python some_script.py`, you install
`pymdkit` once and call any tool from anywhere as `pmk <command> [options]`.

Every command exposes named `-flags` (no positional guessing), and each underlying
script is still runnable on its own.

## Install

Create a clean conda environment, activate it, then install `pymdkit` with pip:

```bash
conda create -n pymdkit python=3.10
conda activate pymdkit
pip install pymdkit
```

This installs the `pmk` command into the active conda environment, together
with its direct dependencies, including ASE, pymatgen, py4vasp-core, dpdata, matplotlib, and UMAP. Some of these packages may install their own transitive dependencies.

GPU-accelerated RMSD screening is optional. For a CUDA 12.x cluster (including
the supplied `module load cuda/12.6` setup), install the matching CuPy extra:

```bash
python -m pip install "pymdkit[gpu-cuda12]"
```

Verify:

```bash
pmk -version
pmk -help                   # lists every command
pmk <command> -help         # shows that command's flags
```

## Global configuration and logs

On first use, PMK creates `.pmkrc.yaml` in the active Python environment root.
For example, an executable at `~/anaconda3/envs/pymdkit/bin/pmk` uses
`~/anaconda3/envs/pymdkit/.pmkrc.yaml`. View the active file and settings with:

```bash
pmk config
pmk config silence -on
pmk config silence -off
```

With silence enabled, successful commands do not print to the terminal. Errors
that prevent an operation are still shown. Command details are never discarded:
every `pmk <command> ...` invocation appends its exact command, timestamp, exit
status, and screen output to `<command>.txt` in the directory where PMK was
invoked. Repeated runs remain in chronological order in the same file. This log
always stays in the invocation directory even when `-of` selects an output
folder. Existing command reports that use the same filename, such as `rmsd.txt`
and `ehull.txt`, are preserved as the report section of that run.

## Commands

Commands that transform structures accept either a single file (`-i`/`-o`) or a
whole folder (`-if`/`-of`); commands that analyse VASP runs scan the current
directory for job sub-folders automatically.

| Command | What it does |
|---|---|
| `config` | View or edit global PMK settings in `.pmkrc.yaml` |
| `add-config-type` | Add Config_type to one XYZ structure or every trajectory frame |
| `gpumd-group` | Tag atoms with a GPUMD group index by element order |
| `gpumd-relax` | Write GPUMD energy-minimization jobs for a file, folder, or trajectory |
| `ewald` | Compute CIF electrostatic energy with pymatgen EwaldSummation |
| `ehull` | Auto-detect VASP job folders and compute E_hull vs Materials Project |
| `vasp-fe` | List VASP job final energies from lowest to highest |
| `gather-fs` | Recursively gather final structures from converged VASP or GPUMD jobs |
| `gpumd-thermo` | Export and plot GPUMD thermo data and detect sustained equilibrium changes |
| `gpumd-trjcat` | Merge normal or all GPUMD trajectories and extract the first abnormal region |
| `pca` | Reduce descriptor data to two dimensions; optional FPS sampling |
| `trj-extract` | Extract a GPUMD XYZ trajectory by time range |
| `trj-sparse` | Deterministically retain a requested fraction of XYZ trajectory frames |
| `umap` | Reduce descriptor data to two dimensions with UMAP; optional FPS sampling |
| `gpumd-msd` | Recursively extract directional SDC and Nernst-Einstein conductivity from GPUMD jobs |
| `nep-rmse` | Compute NEP energy/force/stress RMSE with terminal plots and optional candidate selection |
| `perturb` | Generate perturbed structures with dpdata |
| `vasp2xyz` | Collect SCF-converged VASP job folders (any name) into one extxyz file |
| `submit-vasp` | Submit/resubmit VASP job folders while limiting active queue jobs |
| `rmsd` | Compare structures by RMSD and optionally remove higher-energy duplicates |
| `chemsys-entry` | Download stable Materials Project structures for a chemical system |
| `convert` | Convert structures/trajectories/folders, refine CIF symmetry, or extract CSV structures |
| `realloc-occ` | Generate charge-balanced CIF occupancy allocations for a target formula |
| `substitute` | Randomly substitute or remove selected atoms/sites from a structure |
| `supercell` | Build a supercell by matrix, maximum length, or nearest target length; optional GPUMD setup |

| `vasp-relax` | Write VASP relaxation inputs for a structure (or folder); INCAR tags overridable |
| `vasp-static` | Write VASP static / single-point inputs for a structure (or folder) |

## Examples

```bash
pmk add-config-type -i single-structure.xyz -n P-3m1
pmk add-config-type -it trajectory.xyz -n P-3m1
pmk gpumd-group -i opted.cif -elements Li Y Cl -o model.xyz
pmk gpumd-group -if cifs/ -elements Li Y Cl -of cifs-grouped/   # whole folder
pmk gpumd-group -elements Li Y Cl                              # scan subfolders, tag each model.xyz in place
pmk convert -i NaCl.vasp -o NaCl.cif                            # CIF defaults to P1
pmk convert -i NaCl.vasp -o NaCl-symm.cif -symm                 # detect/import symmetry
pmk convert -if vasp-opted -of cif-opted -oe cif                # folder CIFs default to P1
pmk convert -it train.xyz -ot train.extxyz                      # convert a complete trajectory
pmk convert -i model.cif -o matrix.cif -supercell-matrix 1x1x2
pmk convert -i model.cif -o capped.cif -supercell-max-abc 40
pmk convert -i model.cif -o close.cif -supercell-close-abc 20
pmk convert -if Li4YCl6-disorder -of Li4YCl6-order -d2o -supercell-matrix 1x1x1
pmk supercell -i opted.vasp -o matrix.vasp -matrix 1x1x2
pmk supercell -i opted.xyz -o grouped.xyz -matrix 2x2x2 -gpumd-group Li Y Cl
pmk supercell -i opted.vasp -o capped.vasp -max-abc 40
pmk supercell -i opted.vasp -o close.vasp -close-abc 20       # nearest target: e.g. 21,20,20 A
pmk supercell -if vasp-opted -max-abc 20 -individual          # per-structure ./<name>/<name>.<ext>
pmk supercell -if extxyz-opted -max-abc 24 -individual -temp 500 600 -md-if input-files -gpumd-group Li Y Cl
                                                                  # GPUMD: ./<name>/model.xyz (grouped) + ./<name>/<T>/ jobs
pmk vasp-relax  -i opted.vasp -supercell-max-abc 20           # expand cell, then write relax inputs
pmk vasp-relax  -if optimal_occupancy                          # one ./<name>/ job folder per structure
pmk gpumd-relax -if example -nep nep89_20250409.txt            # model.xyz + run.in + NEP per job folder
pmk gpumd-relax -it train.xyz -nep nep.txt -in run.in         # one ./frame_N/ GPUMD job per frame
pmk vasp-static -if cifs/ -supercell-close-abc 20 -custom-setting my_incar.txt
pmk vasp-static -it traj.xyz                                    # one ./frame_N/ job per trajectory frame
pmk gpumd-msd                                   # recursively scans GPUMD jobs -> per-job msd/ + sigma-sdc-summary.txt
pmk gpumd-msd -diffuse_ion Li -ion_charge 1     # choose the mobile ion and signed valence
pmk chemsys-entry -s Li La Ta Cl                 # MP stable entries -> Li-La-Ta-Cl-stable-entries/
pmk ehull -mp-api-key $MP_API_KEY              # scans ./ for VASP jobs -> ehull.txt
pmk ehull -local Li-La-Ta-Cl-stable-entries-opted
pmk vasp-fe                                # scans ./ for VASP jobs -> final-energy.txt with convergence status
pmk gather-fs -job vasp -fs-name CONTCAR -of vasp-opted
pmk gather-fs -job vasp -fs-name CONTCAR -of vasp-opted -ehull 0.028
pmk gather-fs -job gpumd -fs-name relaxed.xyz -of gpumd-opted
pmk vasp2xyz                                    # scans ./ for VASP output folders -> scf-converged.xyz
pmk vasp2xyz -position-only                    # write positions only, without energy/forces/stress
pmk submit-vasp -subscript sub_vasp -queue slurm -max-job-num 30
nohup pmk submit-vasp -subscript sub_vasp -queue slurm -max-job-num 30 > submit-vasp.log 2>&1 &
pmk realloc-occ -i 66-Li4YCl6-5.cif -of chg-balance -formula Li3YCl6
pmk realloc-occ -if cifs -of chg-balance -formula Li3YCl6
pmk realloc-occ -i 66-Li4YCl6-5.cif -of chg-balance -formula Li3YCl6 -occ 0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1
pmk realloc-occ -if cifs -of chg-balance -formula Li3YCl6 -occ 0.25 0.5 0.75 1 -max-supercell 2x2x2
pmk substitute -i Li3YCl6.cif -se Li -sn 3 -we Na -wn 3 -on 100
pmk substitute -i Li3YCl6.cif -se Li -sn 3 -we none -on 100
pmk substitute -i Li96Ta6La11Cl72.cif -se Li1 Li2 -sn 20 67 -we none -on 100
pmk substitute -i Li96Ta6La11Cl72.cif -se Li2 -we none -ref La Ta -d 1.01 1.02
pmk substitute -i Li96Ta6La11Cl72.cif -se Li2 -we none -ref La Ta
pmk ewald -i Li3YCl6.cif
pmk ewald -if Li3YCl6-all
pmk nep-rmse                                    # writes energy/force/stress train txt files, rmse_value.txt, and terminal plots
pmk nep-rmse -select -xyz train.xyz   # interactive candidate selection; writes candidate.xyz and accurate.xyz
pmk perturb -i example.xyz -atom 0.2 -lattice 0.03 -n 100 -o example-perturb-atom-0.2-lattice-0.03.xyz
pmk pca -i descriptor.out -o pca-descriptor.txt -it train.xyz  # also split descriptors by Config_type
pmk pca -i descriptor.out -o pca-descriptor.txt -fps 0.01 -it train.xyz
pmk umap -i descriptor.out -o umap-descriptor.txt -fps 0.01 -it train.xyz
pmk gpumd-thermo                              # current job, or recursively scan job folders
pmk gpumd-thermo -t thermo.out -in run.in
pmk gpumd-trjcat                            # stable frames -> train.xyz; abnormal regions -> abnormal.xyz
pmk gpumd-trjcat -i traj.xyz -o train.xyz -ao abnormal.xyz
pmk gpumd-trjcat -all                       # merge every readable trajectory into train.xyz
pmk trj-extract -it traj.xyz -in run.in -b 100 -e 200 -o traj-100ps-200ps.xyz
pmk trj-sparse -it example.xyz -r 0.5 -o sparse-example.xyz

pmk rmsd -i structure-1.vasp structure-2.vasp            # one pair -> rmsd.txt
pmk rmsd -if vasp-opted                                  # every pair -> rmsd.txt
pmk rmsd -if example -rm-duplicate -of unique-example   # RMSD <= 0.1: keep lowest energy
pmk rmsd -if example -rm-duplicate -of unique-example -gpu
pmk convert -i opted.cif -o opted-symm.cif -symm -add-oxidation yes
pmk convert -if my_cifs/ -of my_cifs-symm -oe cif -symm -symprec 0.1 -add-oxidation no
pmk convert -csv output_Li2YCl3_struct.csv -of Li2YCl3 -symm
```

VASP input commands (`vasp-relax`, `vasp-static`) always produce **individual**
jobs (one structure per folder): `-i` writes into the current dir (or `-o`),
`-if` creates one `./<name>/` folder per structure, and `-it` creates one
`./frame_N/` folder per trajectory frame - all directly in the current path.

Both commands accept exactly one of `-supercell-matrix`,
`-supercell-max-abc`, or `-supercell-close-abc` for single files, folders, and
trajectories. The matrix mode accepts `1x1x2`, three diagonal factors, or nine
integer entries. Maximum mode uses `max(1, floor(L / cell_length))`; close mode
uses the positive integer repeat nearest `L / cell_length`.

They start from sensible default INCAR settings; override them by passing a
settings file with `-custom-setting FILE`. The file may be a Python-dict block
or `KEY = VALUE` lines (a `None`/blank value clears a tag):

```text
custom_settings = {
    "ENCUT": "600.0",
    "ISIF": "3",
    "MAGMOM": None
}
```

`vasp-static -it traj.xyz` (also available on `vasp-relax`) reads a
multi-structure trajectory and writes one job sub-folder per frame
(`frame_1/`, `frame_2/`, ..., prefix configurable via `-frame-prefix`). Each
folder also keeps a `frame_N.xyz`, so `Config_type` survives for a later
`vasp2xyz`.
`gpumd-relax` accepts the same `-i`, `-if`, and `-it` structure modes. `-nep`
is required and the selected potential is copied into every job folder together
with `model.xyz` and `run.in`. `-in FILE` supplies a custom GPUMD input; its
first `potential` filename is synchronized to the copied NEP. Without `-in`,
the generated job-local input is:

```text
potential       <NEP filename>

minimize        fire 2.0e-2 1000000

ensemble nve
dump_xyz    -1 0 1 relaxed.xyz force
run 1
```

`gpumd-msd` recursively discovers job folders at any depth below `-entry-path`
(the current directory by default). A job must contain `msd.out`, `model.xyz`,
`run.in`, and `thermo.out`; the common accidental name `thermos.out` is also
accepted. Following the [official GPUMD `msd.out` format](https://gpumd.org/gpumd/output_files/msd_out.html), each group block is read as three MSD columns
followed by three directional SDC columns. PMK takes the arithmetic mean over
the last 50% of the SDC correlation-time range, converts A^2/ps to cm^2/s by
multiplying by `1e-4`, and reports `D_total = (D_x + D_y + D_z) / 3`.
Each job receives `msd/group_INDEX_ELEMENT_TEMP_msd.txt`, a matching `_sdc.txt`,
and a `_results.txt`. The root receives a compact `sigma-sdc-summary.txt` with
one fixed-width, right-aligned row per job folder. Conductivities use two
decimal places; diffusion coefficients use compact scientific notation with
two decimal places. If the
mobile ion occupies multiple groups, conductivities are summed and diffusion
coefficients are particle-count-weighted.

Ionic conductivity uses the Nernst-Einstein relation with the mean temperature
of the final SDC-producing MD run and the official 18-column
[`thermo.out` format](https://gpumd.org/gpumd/output_files/thermo_out.html):
for NPT ensembles (`npt_scr`, `npt_mttk`, and other `npt_*` forms), PMK uses
the mean triclinic-box volume over that run; for NVT ensembles (`nvt_bdp`,
`nvt_ber`, and other `nvt_*` forms), it uses the final fixed-cell volume.
`V(A^3) * 1e-24 -> V(cm^3)`, exact SI values for elementary charge and the
Boltzmann constant, and a final `* 1000` conversion from S/cm to mS/cm.

`realloc-occ -i INPUT.cif -of OUTPUT -formula FORMULA` computes atom counts as
Wyckoff multiplicity times occupancy. Because source `_cell_formula_units_Z`
belongs to the source formula, PMK infers a new target Z from the available site
capacities and current atom counts, then writes that target Z to each output CIF.
It enumerates finite site-wise
occupancy allocations, retains only charge-balanced candidates using the CIF
oxidation-number loop, and preserves the original symmetry, coordinates,
Wyckoff labels, and CIF layout. A missing or all-zero oxidation assignment is
first inferred from pymatgen's highest-ranked charge-balanced composition
guess. If a defective/non-stoichiometric composition has no charge-balanced
guess, PMK uses each element's first conventional pymatgen oxidation state.
Every input receives `OUTPUT/INPUT_STEM/`; for
the supplied Li4YCl6 example, the two outputs are `...-d1.cif` with `Li1=0.50`
and `...-d2.cif` with `Li2=0.50`. Folder mode uses `-if` with the same `-of` and
`-formula` options. Add `-occ` followed by the only permitted occupancy values;
input grids support three decimal places, including `0.125`, `0.375`, `0.625`,
and `0.875`. Generated CIF occupancy fields use two decimal places.
`-max-supercell AxBxC` searches bounded diagonal supercells, smallest volume
first, until every multiplicity times occupancy times supercell determinant is
an integer. Equal-volume choices favor the most isotropic resulting cell.
All selections, site atom counts, oxidation states, charges, failures, the
allowed occupancy grid, and the final summary are recorded in
`realloc-occ.txt` in the invocation directory for both single-file and folder modes.

`gather-fs` recursively scans all job folders below the current directory. For
VASP it follows the global `vaspout.h5 > vasprun.xml > OUTCAR` priority and
requires full convergence before copying the selected `-fs-name` as `.vasp`. The gathered VASP file's first line is replaced with `energy=<final_energy> eV` for later screening.
For GPUMD it requires `gpumd.out` to report a force tolerance and a final
`f_max` no greater than that tolerance before copying the final structure.
`-fs-name` defaults to `CONTCAR` for VASP and `relaxed.xyz` for GPUMD. The
`-ehull` and `-ehull-file` filters are available only with `-job vasp`.
`rmsd -i STRUCTURE_1 STRUCTURE_2` compares one pair, while `rmsd -if FOLDER`
records every unique pair in `rmsd.txt`. With `-rm-duplicate -of OUTPUT`, files
are processed from lowest to highest energy; a structure is removed only when
its RMSD to an already retained lower-energy representative is no greater than `0.1`. Duplicate removal supports `.vasp`,
`.xyz`, and `.extxyz`: VASP energy is read from the first line, and extended-XYZ
energy from the second-line `energy=` tag. The report records every RMSD,
duplicate decision, energy, group, and retained representative.

`rmsd -if FOLDER -rm-duplicate -of OUTPUT -gpu` is intended for a GPU batch
job such as the supplied Slurm `sub_pmk`. CuPy computes permutation-invariant,
periodic distance fingerprints in CUDA batches and skips clearly dissimilar
pairs; pymatgen `StructureMatcher` remains the final authority for every
candidate RMSD and duplicate decision. `rmsd.txt` labels pair rows as `exact`,
`gpu-filtered`, or `incompatible` and records the CUDA device and pair counts.
The `-gpu` option requires both `-if` and `-rm-duplicate`.

For the supplied CUDA 12.6 Slurm script, submit the GPU workflow with:

```bash
sbatch sub_pmk
```

`convert -symm` uses pymatgen's spglib-backed `SpacegroupAnalyzer` for authoritative
space-group detection and conventional-cell refinement. Its default Cartesian
symmetry tolerance is `0.1 Angstrom` (`-symprec`), with a `5 degree` angle
tolerance (`-angle-tolerance`). CIF occupancies are retained by the pymatgen
structure and written directly rather than reconstructed after symmetry finding.
With `-symm`, CSV mode reads serialized pymatgen `Structure` dictionaries from a required `cif`
column. `-csv INPUT.csv -of NAME -symm` writes `INDEX-NAME-SPACEGROUP.cif`, for
example `1-Li2YCl3-66.cif`. Every run also writes `symmetry.txt` with one row
per successfully generated CIF: filename, crystal system, space-group symbol,
and space-group number. Crystal-system and space-group summaries are sorted
from highest to lowest count, with alphabetical ordering for ties, followed by
the total number of structures. Single-file mode places the report beside the
output CIF; folder and CSV modes place it inside `-of`. Folder mode keeps each
input stem unchanged, writing `NAME.cif` rather than `NAME-symm.cif`.

Symmetry-refined CIFs use `pymatgen.core.Lattice` and the international
monoclinic convention, where beta is the non-right angle. The lattice matrix is
always interpreted row by row as vectors a, b, and c. For standardized cells,
the `_symmetry_equiv_pos_as_xyz` loop is sourced from pymatgen's `SpaceGroup`
database instead of serializing noisy numerical translations returned by
structure-specific symmetry detection.

Every CIF exported by `convert`, `substitute`, `chemsys-entry`,
`supercell`, or `gpumd-group` is written through the same CIF helper. `convert`
writes P1 by default and performs symmetry detection only when `-symm` is
present; the other commands retain their command-specific symmetry behavior.
Each CIF places the IUCr-defined `_space_group_crystal_system` item immediately
below its `data_...` header and contains both
`_atom_site_site_symmetry_multiplicity` and `_atom_site_Wyckoff_symbol` in the
atom-site loop. Multiplicity-letter
assignments use the conventional
International Tables setting from pymatgen/spglib and correspond to the
[Bilbao Crystallographic Server WYCKPOS tables](https://cryst.ehu.es/cryst/get_wp.html);
the provenance URL is also recorded inside each generated CIF.

For structural conversion, supercell creation, and GPUMD grouping, PMK reads
CIF input with pymatgen rather than ASE. ASE interprets
`_space_group_crystal_system` as a setting hint only for rhombohedral groups
and emits a potentially misleading warning for other crystal systems. Using
pymatgen avoids that ambiguity while PMK continues to use ASE for non-CIF
formats. Pymatgen intentionally snaps fractional coordinates within its
default `1e-4` tolerance to ideal values before applying symmetry. PMK keeps
that corrective behavior but hides only its corresponding finite-precision
notice; composition, occupancy, malformed-CIF, and all other warnings remain
visible. See the
[pymatgen CifParser documentation](https://pymatgen.org/pymatgen.io.html#pymatgen.io.cif.CifParser).
A partially occupied CIF must first be ordered with `convert -d2o` before it
can be represented as an explicit atom list.

CIF `_atom_type_symbol` and `_atom_site_label` values use per-element
inequivalent-site identifiers such as `Li1`, `Y1`, `Cl1`, and `Cl2`.
The final atom-site loop starts with `_atom_site_label`, followed by
`_atom_site_type_symbol`, whose values are plain elements such as `Li`, `Y`,
and `Cl`. Oxidation numbers remain in the separate
`_atom_type_oxidation_number` column. When a generated multi-element CIF has
no oxidation states or only zero-valued states, PMK first tries pymatgen's
highest-ranked charge-balanced assignment. If none exists, it uses the first
conventional oxidation state of each element. This writes `Li1.0`, `Y3.0`, and
`Cl-1.0` even for non-neutral defective compositions such as `Li2YCl3`.
Explicit nonzero oxidation states are preserved. PMK does not add an
`_audit_creation_method` block or a generated-by header comment.

Newly serialized floating-point values use eight digits after the decimal
point across PMK structure, trajectory, data, report, and screen output. CIF
`_atom_site_occupancy` values use two digits and
`_atom_type_oxidation_number` values use one digit. In extended XYZ, zero components
inside the second-line `Lattice="..."` value remain `0.0`; other generated
floating-point values use eight digits. Commands that select or concatenate
existing trajectory frames continue copying those frame bytes exactly.
`convert` writes `.vasp` structures with fractional (`Direct`) coordinates.

`add-config-type` updates its input atomically in place. With `-i`, the file must contain exactly one XYZ structure; use `-it` for a trajectory. Existing `Config_type` values are replaced, missing values are added, and all other extended-XYZ bytes remain unchanged.

`nep-rmse -select`, PCA FPS, and UMAP FPS copy selected extended-XYZ frame blocks directly from the input trajectory. No parser rewrites energy, stress, forces, positions, precision, or extra metadata.

When `pca` or `umap` receives `-it train.xyz`, it reads `Config_type` from every frame and writes one additional descriptor table per value, such as `pca-descriptor-P-3m1.txt` and `pca-descriptor-Pnma.txt`. With `-fps`, it also writes matching files such as `fps-0.01-pca-descriptor-P-3m1.txt` and `fps-0.01-train-P-3m1.xyz`. These files partition the single global FPS selection; FPS is not rerun independently for each `Config_type`. Frames without `Config_type` remain in the main outputs but do not receive a type-specific file.

For folder conversion, `convert -if INPUT -o OUTPUT` and
`convert -if INPUT -of OUTPUT` are equivalent and write `.xyz` by default. Set
another target extension with `-oe`, for example `-oe cif` or `-oe vasp`.
Single-file and trajectory modes infer the target format from `-o` and `-ot`.
Because Wyckoff metadata describes one crystal structure, each `.cif` output
file accepts exactly one frame.

For ordinary conversion, the three supercell options are mutually exclusive and
are resolved separately for every input structure. For example, a cell with
`a=3 A`, `b=4 A`, and `c=5 A` plus `-supercell-close-abc 20` uses `7x5x4`,
giving approximately `21 A`, `20 A`, and `20 A`. The numeric value is the target:
`-supercell-close-abc 40` targets approximately `40 A`.

`convert -d2o` orders partially occupied CIF structures.
`-supercell-matrix` defaults internally to `1x1x1` in this mode; three
diagonal factors or nine integer matrix entries are accepted. Because Ewald
candidates must be compared using one fixed cell, `-supercell-max-abc` and
`-supercell-close-abc` are not accepted with `-d2o`. Folder ordering scans
recursively, groups `BASE-dN.cif` files by `BASE`, and writes only the
lowest-total-Ewald candidate into the flat output folder as `BASE-oN.cif`,
where `N` remains tied to its source `-dN` file. Pymatgen
`EnumerateStructureTransformation` is used when enumlib is available; the
expected missing-enumlib runtime error quietly switches to
`OrderDisorderedStructureTransformation`. Every output CIF records total Ewald
energy on line two as `# E_e = VALUE eV`, and `d2o.txt` records all candidates,
methods, energies, failures, and selections. Add `-symm` to import symmetry into
the ordered CIF; otherwise it is P1. Without `-d2o`, the
`-supercell-matrix`, `-supercell-max-abc`, and `-supercell-close-abc` options
build an ordinary structural supercell and never enable ordering automatically.

After `supercell -temp` creates all requested jobs, it checks generated folders
that contain `run.in`, `model.xyz`, and `nep.txt`. Existing `ensemble npt_scr`
lines retain their temperature parameters but are rewritten with three
pressure/modulus components for 90-degree orthorhombic cells and all six
components for non-orthorhombic cells, following the
[official GPUMD ensemble syntax](https://gpumd.org/gpumd/input_parameters/ensemble_standard.html).
Other ensemble types are not modified. Use `-gpumd-group Li Y Cl` to embed the
standalone `gpumd-group` tagging behavior in generated `model.xyz` files. The
same option also tags ordinary `.xyz` or `.extxyz` supercell output; formats
that cannot preserve GPUMD group arrays are rejected before output is written.
`gpumd-thermo` writes `temperature.txt`, `potential-energy.txt`, `pressure.txt`, `lattice-parameters.txt`, `volume.txt`, `lattice-angles.txt`, and a headless `thermo.png` in every detected GPUMD job folder. For an abnormal run it also writes `abnormal-time.txt`, with one `BEGIN ps to END ps` interval per line; a normal rerun removes a stale file. Potential energy is the primary equilibrium signal. Detection requires a persistent change with a meaningful magnitude relative to the robust natural energy fluctuations, so ordinary correlated thermal noise and isolated spikes are not classified as abnormal. An abnormal classification always requires a mapped potential-energy change. When at least two lattice-parameter or volume signals show the same mapped transition, they may enable a more sensitive noise-relative energy test; this avoids a fixed eV or percentage-of-total-energy threshold that would depend on system size. Cell signals can corroborate and refine the energy interval, but can never classify a run by themselves. Temperature and all six pressure components remain diagnostics only and cannot classify an otherwise stationary run as abnormal by themselves. A lone cell-axis drift, an isolated spike, or an empty anomaly mask never triggers a whole-run fallback. Thermo plots use Arial when available and Matplotlib's built-in qualitative `Dark2` color cycle. They plot all six pressure components, use frameless upper-right legends with added vertical headroom, preserve correct Å units, and show abnormal regions as shaded legend entries on potential energy, lattice parameters, and volume. Potential energy and volume are divided by 1000 and labeled ×10³ only when their absolute plotted values reach 1000. The analysis uses the last `run` segment, or the last 80% when `run.in` has no `run` record; inspect the plot before making a final scientific judgment.

GPUMD box matrices are interpreted row by row by
`pymatgen.core.Lattice`, which supplies a, b, c, alpha, beta, gamma, and volume.

The thermo parser follows the [official GPUMD thermo.out format](https://gpumd.org/gpumd/output_files/thermo_out.html): 18 columns are required, `pressure.txt` exports `Pxx Pyy Pzz Pyz Pxz Pxy`, box vectors are interpreted as a full 3×3 matrix, `time_step` propagates between runs, and `dump_thermo` does not.

`gpumd-trjcat` recursively scans job subfolders for `traj.xyz`, `thermo.out`, and `run.in`. By default, normal trajectories are merged in deterministic folder order into `train.xyz` (or `-o`), while only frames from the earliest mapped abnormal interval of each abnormal job are written to `abnormal.xyz` (or `-ao`). With `-all`, every readable full trajectory is merged into `train.xyz` regardless of normal, abnormal, missing-thermo, or failed-thermo status; an analyzable abnormal job still contributes its earliest abnormal interval to `abnormal.xyz`. If no mapped abnormal frames are found, `abnormal.xyz` is not created and a stale file at the selected `-ao` path is removed. Frame times are obtained from `time_step` and `dump_xyz` in that job's `run.in`. In default mode, missing, unreadable, or unmappable jobs are skipped. When `<job>/model.xyz` contains `Config_type`, that value is applied to every exported frame while all other extended-XYZ bytes remain unchanged.

`trj-sparse` keeps `max(1, floor(frame_count × ratio))` frames without rewriting their extended-XYZ content. Selection is deterministic and starts with frame zero; for 100 frames and `-r 0.5`, it writes source indices `0, 2, ..., 98`.

Each command's full flag list is in `pmk <command> -help`.

Development rule: when an established dependency already provides the required
domain operation, PMK uses that supported API rather than maintaining a custom
implementation.

For long VASP batch submission, run `submit-vasp` with `nohup` and `&` if you
want it to keep sleeping, checking the queue, and submitting new jobs after you
exit the terminal:

```bash
nohup pmk submit-vasp -subscript sub_vasp -queue slurm -max-job-num 30 > submit-vasp.log 2>&1 &
```

The command itself controls the loop: it submits until the active queue reaches
`-max-job-num`, sleeps when the queue is full, checks again, and continues until
all needed jobs are submitted. `nohup ... &` is what makes that loop continue in
the background after logout.

`substitute -ref` removes selected sites near reference sites and writes one
`<input-stem>_substitute.cif` in the current path. If `-d` is omitted, each
cutoff is `0.7 * (selected covalent radius + reference covalent radius)` using
the covalent radii from Cordero et al., Dalton Trans., 2008, 2832-2838.

VASP-output readers (`vasp2xyz`, `ehull`, `vasp-fe`, and other future VASP-output
commands) use the global priority `vaspout.h5 > vasprun.xml > OUTCAR`.

`ehull` auto-detects every sub-folder of the current path that contains a
supported VASP output, groups them by chemical system (elements ordered by
electronegativity, e.g. `Li-Y-Cl`), and builds/reuses one `mp_cache_<system>.json`
per system - so a pure Li-Y-Cl batch yields a single `mp_cache_Li-Y-Cl.json`, while
a mixed Li-Y-Cl + La-O batch yields both `mp_cache_Li-Y-Cl.json` and
`mp_cache_La-O.json`. (Formation energy is reported alongside E_hull in
`ehull.txt`.)

## Layout

```
pymdkit/
|-- pyproject.toml              # package metadata + the `pmk` entry point
|-- README.md
`-- src/pymdkit/
    |-- pymdkit_main.py         # dispatcher: discovers and runs commands
    `-- commands/               # one module per command
        |-- _fileio.py          # shared -i/-o/-if/-of helper (not a command)
        |-- _gpumd.py           # shared GPUMD run.in and thermo analysis
        |-- _geometry.py        # shared minimum-image geometry helper
        |-- _vaspset.py         # shared VASP input-set helper (not a command)
        |-- gpumd_group.py
        |-- gpumd_msd.py
        |-- gpumd_relax.py
        |-- gpumd_thermo.py
        |-- gpumd_trjcat.py
        |-- compute_ehull.py
        |-- compute_rmsd.py
        |-- ewald.py
        |-- vasp_fe.py
        |-- perturb.py
        |-- nep_rmse.py
        |-- add_config_type.py
        |-- _cifio.py
        |-- chemsys_entry.py
        |-- submit_vasp.py
        |-- convert.py
        |-- realloc_occ.py
        |-- substitute.py
        |-- supercell.py
        |-- vasp2xyz.py
        |-- vasp_relax.py
        |-- vasp_static.py
        |-- ...
        `-- _symmetry.py       # private symmetry conversion/report helper
```

Modules whose name starts with `_` are shared helpers and are skipped by the
dispatcher, so they never appear as commands.

## Adding a new tool later

Drop a module in `src/pymdkit/commands/` that defines four things:

```python
COMMAND = "my-tool"                 # the subcommand name you'll type
HELP = "One-line description."

def add_arguments(parser):          # register flags
    parser.add_argument("-input", required=True)

def run(args):                      # do the work; return an exit code (0 = ok)
    ...
    return 0

if __name__ == "__main__":          # keeps the script runnable on its own
    import argparse
    _p = argparse.ArgumentParser(description=__doc__)
    add_arguments(_p)
    raise SystemExit(run(_p.parse_args()))
```

It will appear in `pmk -help` automatically - no central registration needed.
Put heavy imports (pymatgen, ase, ...) inside `run()` where practical; the dispatcher
reads each command's name and help without importing it, so `pmk -help` stays
fast and a missing optional dependency only affects the one command that needs it.

## Running a script standalone

Every command module still works directly, which is handy for debugging:

```bash
python src/pymdkit/commands/supercell.py -i in.cif -close-abc 20 -o sc.vasp
```
