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.

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

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 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 -if Li4YCl6-disorder -of Li4YCl6-order -d2o -supercell 1x1x1
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 -group Li Y Cl
                                                                  # GPUMD: ./<name>/model.xyz (grouped) + ./<name>/<T>/ jobs
pmk vasp-relax  -i opted.vasp -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/ -max-abc 20 -custom-setting my_incar.txt # expand cell, static inputs
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 symmetrize -i opted.cif -add_oxidation yes -o opted-symm.cif
pmk symmetrize -if my_cifs/ -symprec 0.1 -add_oxidation no -of my_cifs-symm
pmk symmetrize -csv output_Li2YCl3_struct.csv -of Li2YCl3

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 -max-abc L for single files, folders, and trajectories. Before writing each VASP job, each lattice vector is repeated max(1, floor(L / cell_length)) times, matching the supercell command.

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

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, 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: 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, and failures are recorded in OUTPUT/realloc-occ.txt.

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:

sbatch sub_pmk

symmetrize 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. CSV mode reads serialized pymatgen Structure dictionaries from a required cif column. -csv INPUT.csv -of NAME 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.

Every CIF exported by symmetrize, 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; the provenance URL is also recorded inside each generated CIF.

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.

convert -d2o orders partially occupied CIF structures. -supercell defaults to 1x1x1; three diagonal factors or nine integer matrix entries are also accepted. 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, -supercell only builds an ordinary structural supercell and never enables 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. Other ensemble types are not modified. Use -group Li Y Cl to add GPUMD groups.

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.

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. 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.

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_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
        |-- ...
        `-- 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