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:

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.

Verify:

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

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.

Examples

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 opted.vasp -o opted.xyz                          # convert one structure
pmk convert -if vasp-opted -of cif-opted -oe cif                # folder conversion target format
pmk convert -it train.xyz -ot train.extxyz                      # convert a complete trajectory
pmk supercell -i opted.vasp -o sc.vasp -max-abc 20            # cell lengths <= 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 -add-groups Li Y Cl
                                                                  # GPUMD: ./<name>/model.xyz (grouped) + ./<name>/<T>/ jobs
pmk vasp-relax  -i opted.vasp                                  # relax inputs in current dir
pmk vasp-relax  -if optimal_occupancy                          # one ./<name>/ job folder per structure
pmk vasp-static -if cifs/ -custom-setting my_incar.txt         # static inputs, custom INCAR
pmk vasp-static -it traj.xyz                                    # one ./frame_N/ job per trajectory frame
pmk msd                                         # scans <structure>/<temp>/ -> per-job msd/ + msd_summary.txt
pmk msd -diffuse_ion Li -ion_charge 1         # choose the mobile ion for conductivity
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-contcar -of vasp-opted               # CONTCARs -> vasp-opted/<folder>.vasp
pmk gather-contcar -of vasp-opted -ehull 0.028  # only structures with E_hull < 0.028 eV/atom
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 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 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 a.cif b.cif                            # RMSD of two files -> rmsd.txt
pmk rmsd vasp-opted/                            # all pairs in a folder -> rmsd.txt
pmk symmetrize -i opted.cif -symprec 0.01 -add_oxidation yes -o opted-symm.cif
pmk symmetrize -if my_cifs/ -symprec 0.01 -add_oxidation no -of my_cifs-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.

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):

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.

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 -of OUTPUT writes .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.

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. 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. When multiple lattice or volume measures independently show a transition, a more sensitive noise-relative potential-energy test is enabled; this avoids a fixed eV or percentage-of-total-energy threshold that would depend on system size. Lattice parameters and volume can corroborate an energy transition. Temperature and all six pressure components remain diagnostics only and cannot classify an otherwise stationary run as abnormal by themselves. Thermo plots use Arial when available and the Morandi color cycle. They plot all six pressure components, use frameless upper-right legends with added vertical headroom, preserve correct Å units, and shade abnormal regions 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.

The thermo parser follows the official GPUMD thermo.out format: 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. Stable trajectories are merged in deterministic folder order into train.xyz (or -o). For an abnormal job, frames from the structurally corroborated potential-energy region are written to abnormal.xyz (or -ao); if no energy region can be resolved, the dominant structural transition is used as a fallback so an abnormal job cannot produce an empty abnormal output. Frame times are obtained from time_step and dump_xyz in that job's run.in. 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.

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:

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_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
        |-- chemsys_entry.py
        |-- submit_vasp.py
        |-- convert.py
        |-- substitute.py
        |-- supercell.py
        |-- vasp2xyz.py
        |-- vasp_relax.py
        |-- vasp_static.py
        |-- ...
        `-- symmetrize.py

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:

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:

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