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, h5py, dpdata, matplotlib, and UMAP. Some packages also install 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:

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

Verify:

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:

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 a single file (-i/-o), a folder (-if/-of), or a trajectory where applicable. PMK delegates ordinary format recognition to ASE, while CIF uses the specialized pymatgen path needed to preserve occupancy, oxidation, symmetry, and Wyckoff metadata.

Examples

pmk convert -i NaCl.vasp -o NaCl.xyz                    # always extended XYZ
pmk convert -i NaCl.vasp -o NaCl.xyz -config-type Fm-3m
pmk convert -if structures -of structures-pdb -oe pdb
pmk convert -it trajectory.traj -ot trajectory.xyz      # ASE trajectory conversion
pmk convert -i NaCl.vasp -o NaCl-symm.cif -symm
pmk convert -if cifs -of ordered -d2o -supercell-matrix 1x1x1 -ncore 64

pmk supercell -i model.vasp -o model-sc.vasp -matrix 1x1x2
pmk supercell -if structures -of capped -max-abc 40
pmk supercell -if structures -of close -close-abc 50

pmk gpumd-job -if Li3YCl6-all -supercell-close-abc 50 -nep nep.txt -in run.in -temp 300 320 340 360 -group Li Y Cl
pmk gpumd-job -it train.xyz -supercell-matrix 2x2x2 -nep nep.txt -group Li Y Cl

pmk vasp-job -if Li3YCl6-all -incar INCAR
pmk vasp-job -it train.xyz -incar INCAR -supercell-max-abc 40
pmk vasp-job -submit -subscript sub_vasp -queue slurm -max-job-num 30
pmk vasp-job -status -queue slurm
nohup pmk vasp-job -submit -subscript sub_vasp -queue slurm -max-job-num 30 > vasp-job.log 2>&1 &

pmk trjconv -i example.xyz -b 10 -e 20 -o frames-10-20.xyz
pmk trjconv -i example.xyz -dt 10 -o every-10-frames.xyz
pmk trjconv -i example.xyz -b 10 -e 20 -timestep 1 -o time-10-20ps.xyz
pmk trjconv -i example.xyz -dt 10 -timestep 1 -o every-10ps.xyz

pmk rd -pca -i descriptor.out -o pca-descriptor.txt -it train.xyz
pmk rd -pca -i descriptor.out -fps 0.01 -it train.xyz
pmk rd -umap -i descriptor.out -o umap-descriptor.txt -fps 0.01 -it train.xyz

pmk mp-chemsys -s Li La Ta Cl
pmk gpumd-msd
pmk gpumd-thermo
pmk vasp-out -energy -force -stress
pmk vasp-out -born-effective-charge -dielectric-tensor
pmk vasp-out -all

pmk gpumd-trjcat -all
pmk ehull -local Li-La-Ta-Cl-stable-entries-opted
pmk ehull -ref-form-e -1.1
pmk gather-fs -job vasp -fs-name CONTCAR -of vasp-opted
pmk vasp2xyz -position-only
pmk ewald -if cifs -ncore 64
pmk realloc-occ -if cifs -of disorder -formula Li3YCl6 -occ 0.25 0.5 0.75 1 -ncore 64
pmk rmsd -if vasp-opted -rm-duplicate -of unique -gpu -ncore 64

File I/O and XYZ

convert uses ASE read, iread/multi-frame indexing, format detection, and write for ordinary formats supported by the installed ASE version. This covers single structures, files in a folder, and trajectory conversion whenever the selected format supports multiple configurations. Run ase info --formats to inspect the exact formats available in that environment. Specialized CIF reading/writing remains pymatgen-based to preserve PMK's disorder and symmetry metadata.

For CIFs with explicit symmetry multiplicities, PMK validates every labeled Wyckoff orbit and the occupied composition after pymatgen expands the asymmetric unit. If finite coordinate noise splits a special position, PMK retries with the smallest bounded site tolerance that restores the declared multiplicities. It never accepts an expansion with different stoichiometry; an irreconcilable CIF stops with an explicit error instead of continuing with the wrong atoms.

Every newly exported .xyz file is written as extended XYZ (extxyz), including single structures and trajectories. Commands that select existing trajectory frames preserve already-extended source blocks byte-for-byte; if the source is plain XYZ, selected frames are upgraded to extxyz. convert -config-type NAME adds or replaces Config_type in every output XYZ frame and supports -i, -if, and -it conversion modes.

Trajectory conversion

Without -timestep, trjconv -b/-e uses one-based, inclusive structure numbers. -dt N retains the first frame and then every Nth frame. With -timestep, every frame must contain a numeric Time= value and physical time is Time * timestep in ps; -b, -e, and -dt are then interpreted in ps. The output format is inferred from -o, and .xyz output is always extxyz.

Dimension reduction

rd requires exactly one of -pca or -umap. Descriptor rows remain aligned with -it trajectory frames. Config_type-specific descriptor files and FPS trajectory/descriptor outputs retain the earlier PCA/UMAP behavior; FPS is run once globally and the selection is then partitioned by Config_type.

GPUMD jobs

gpumd-job accepts exactly one of -i, -if, or -it and optionally one of -supercell-matrix, -supercell-max-abc, or -supercell-close-abc. -nep and -in default to nep.txt and run.in in the invocation directory. Each job receives model.xyz as extxyz, the potential as nep.txt, and the input as run.in; the potential line is synchronized to nep.txt. With -temp, PMK creates <structure>/<T>/, replaces every target_t token with T, applies -group to all model files, and normalizes npt_scr pressure components for the actual orthorhombic or triclinic cell.

The standalone supercell command is now structure-only. Its removed -individual, -temp, -md-if, and -gpumd-group workflows belong to gpumd-job. The standalone gpumd-group command remains available for tagging existing files.

E_hull reference

ehull -ref-form-e EV_PER_ATOM keeps the normal Materials Project or -local formation-energy calculation, but replaces the phase-diagram hull distance with E_hull = E_form - EV_PER_ATOM. All scanned structures must have the same reduced formula because one scalar reference applies to one composition. The exact difference is retained, including a negative value when a structure lies below the supplied reference.

VASP jobs

vasp-job -i/-if/-it -incar INCAR creates one job folder per structure. PMK uses pymatgen to generate POSCAR, POTCAR, and KPOINTS, then preserves the supplied INCAR file exactly in each job. All three supercell modes are available during generation. vasp-job -submit retains the previous convergence-aware Slurm/PBS queue limiter: converged jobs are skipped, unconverged relaxations restart from CONTCAR, and the process sleeps while the active queue is full. Use nohup ... & when the submit/check loop must survive logout.

vasp-job -status recursively reports each prepared job and summary counts for total, finished, converged, failed, and running. Submitted scheduler IDs are stored in each job folder so running jobs can be identified with squeue or qstat. For vaspout.h5, status uses h5py to read only the final energy, final positions, required INCAR scalars, and the final force slice when checking a relaxation. It does not construct a py4vasp Calculation or load complete ionic trajectories. VASP output parsing still follows vaspout.h5 > vasprun.xml > OUTCAR.

VASP quantities

vasp-out scans VASP output files in the current path and every descendant folder. Select any combination of quantity flags, or use -all to attempt all 39 quantities listed by py4vasp Calculation, including namespaced quantities such as -electron-phonon-transport and -phonon-mode. Underscore/dotted spellings from the API are also accepted; see pmk vasp-out -h for the exact list.

Reader priority is applied independently for every requested quantity: vaspout.h5, then vasprun.xml with pymatgen, then OUTCAR with pymatgen. For vaspout.h5, final energy is read directly with h5py as one scalar hyperslab; other quantities use py4vasp. The fast path follows VASP's existing dataset and chunk layout and never rewrites or rechunks the file. If a preferred file exists but does not contain that quantity, PMK continues to the next reader. Energy always means the final energy(sigma->0) value (e_0_energy in pymatgen Vasprun data). Older or nonstandard HDF5 layouts fall back to py4vasp.

Results are stored as one text report per quantity under vasp-out/, such as energy.txt, force.txt, and stress.txt. Numeric values use eight decimal places. Forces include element-number labels; stresses are reported in both GPa and eV/ų using 1 eV/ų = 160.2176621 GPa. Multiple calculations are identified by relative folder path. -all skips unavailable datasets per calculation instead of aborting the scan.

Other behavior

convert -d2o first symmetrizes the disordered structure so each distinct Wyckoff orbit keeps its own occupancy constraint. It then uses pymatgen's branch-and-bound Ewald ordering to find the lowest-energy ordered structure for the fixed supercell matrix. This avoids the previous greedy solution and prevents atoms from being redistributed between unrelated but chemically identical orbits.

There is deliberately no convert -d2o -gpu flag: symmetry grouping and the irregular branch-and-bound search are CPU/control-flow work, so transferring small candidate sets to CUDA would add overhead without accelerating the search.

Multicore CPU

Independent CPU-heavy inputs are processed concurrently by convert -d2o, ewald -if, realloc-occ -if, and exact folder-pair comparisons in rmsd. Use -ncore N to choose the number of worker processes, or -ncore 1 to disable multiprocessing. Without -ncore, PMK uses the single-node CPU allocation reported by Slurm, PBS, or SGE and never creates more workers than there are independent tasks. Outside a scheduler, automatic mode is capped at four workers to avoid unexpectedly occupying an entire workstation. Numerical libraries are limited to one native thread in each worker to prevent nested thread oversubscription.

For one PMK process using 64 CPU cores on one Slurm node:

#!/bin/bash
#SBATCH -J Nwyd-pmk
#SBATCH -p normal
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=64

source ~/softwares-env/anaconda.sh
conda activate pymdkit

pmk convert -if Li4YCl6-disorder/ -of Li4YCl6-disorder-order -d2o

The older single-node form -N 1 with -n 64 is also detected through SLURM_NTASKS. This is local multiprocessing, not distributed multi-node execution, so keep the job on one node. The project audit intentionally leaves HDF5/VASP scans, GPUMD plotting, trajectory copying, scheduler submission, network/API commands, and NumPy/UMAP reductions outside this process pool: those paths are I/O-bound, own their native parallelism, or write shared state. realloc-occ.txt and every generic <command>.txt log stay in the invocation directory. Silence mode suppresses successful command output, but pmk <command> -h always displays help and operation-blocking errors are always visible.

Each command's full flag list is available with pmk <command> -h. Automatic scan modes include matching files directly in the current path as well as files in folders, subfolders, and deeper descendants. This rule applies to shared structure-folder processing, VASP calculations, GPUMD outputs, and trajectory gathering; explicit -i and -if modes remain available.

Layout

pymdkit/
|-- pyproject.toml              # package metadata and `pmk` entry point
|-- README.md
`-- src/pymdkit/
    |-- pymdkit_main.py         # command discovery and dispatch
    |-- _config.py              # environment-wide .pmkrc.yaml settings
    |-- _logging.py             # append-only command reports
    `-- commands/
        |-- convert.py          # ASE/pymatgen conversion and CIF workflows
        |-- trjconv.py          # trajectory extraction and stride
        |-- rd.py               # PCA/UMAP dimensional reduction
        |-- gpumd_job.py        # GPUMD job generation
        |-- vasp_job.py         # VASP generation, submission, and status
        |-- chemsys_entry.py    # `mp-chemsys` implementation
        |-- ...                 # remaining public commands
        `-- _*.py               # private shared helpers

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 -close-abc 20 -o sc.vasp