Skip to content

🧮 Calculator Scripts

The calculator module contains scripts for computing properties from MD trajectories, NEP models, and structure files.

Script Location: Scripts/calculators/

Overview

The calculator module can be read as three groups:

  • Trajectory properties: compute time-dependent quantities such as MSD and ionic conductivity from GPUMD or extxyz trajectories, and calculate X-ray diffraction from extxyz trajectories;
  • Phonon properties: calculate phonon force constants and band structures from a primitive cell, a NEP model, and a QPOINTS path;
  • NEP-assisted calculations: use a NEP model to predict energies/forces/stresses, extract descriptors, calculate DOAS, run NEB, or minimize structures;
  • Polar-material analysis: build neighbor lists and calculate local displacement, averaged structures, octahedral tilt, and local polarization for perovskite or polar systems.

If you are not sure about the required arguments, start from the interactive menu to see the prompt. If the arguments are already clear, use gpumdkit.sh -calc ... directly for calculators that provide a CLI; XRD and phonon calculations are intentionally interactive only.

Task Command Main Input
Ionic conductivity gpumdkit.sh -calc ionic-cond <element> <charge> msd.out, thermo.out, model.xyz
MSD from trajectory gpumdkit.sh -calc msd <trajectory.xyz> <element> <dt_fs> extxyz trajectory
X-ray diffraction 4 → 413 (interactive only) extxyz trajectory with Lattice/pbc
Phonon band structure 4 → 414 (interactive only) PRIMCELL.vasp, nep.txt, QPOINTS
NEP prediction gpumdkit.sh -calc nep <input.xyz> <output.xyz> <nep.txt> extxyz + NEP model
NEP descriptors gpumdkit.sh -calc des <input.xyz> <output.npy> <nep.txt> <element> extxyz + NEP model
DOAS gpumdkit.sh -calc doas <input.xyz> <nep.txt> <output.txt> extxyz + NEP model
NEB gpumdkit.sh -calc neb <initial.xyz> <final.xyz> <n_images> <nep.txt> initial/final structures
Minimization gpumdkit.sh -calc minimize <structure> <nep.txt> [fmax] [max_steps] structure + NEP model
Neighbor list gpumdkit.sh -calc nlist ... see Polar Material Analysis
Displacement gpumdkit.sh -calc disp ... see Polar Material Analysis
Average structure gpumdkit.sh -calc avg-struct ... see Polar Material Analysis
Octahedral tilt gpumdkit.sh -calc oct-tilt ... see Polar Material Analysis
ABO3 polarization gpumdkit.sh -calc pol-abo3 ... see Polar Material Analysis

For a full command list:

gpumdkit.sh -calc -h

The command-line help table looks like:

+-------------------------------------------------------------------------------------------------------+
|                                      CALCULATOR TOOLS                                                 |
+-------------------------------------------------------------------------------------------------------+
| Usage: gpumdkit.sh -calc <type> [args...]                                                             |
+-------------------------------------------------------------------------------------------------------+
| ionic-cond <element> <charge>                 Calculate ionic conductivity from MSD data              |
| nep <input.xyz> <output.xyz> <nep_model>      Calculate energy/force/virial with a NEP model          |
| des <input.xyz> <output.npy> <nep_model> <el> Calculate NEP descriptors for one element               |
| doas <input.xyz> <nep_model> <output.txt>     Calculate density of atomistic states                   |
| neb <initial.xyz> <final.xyz> <n_images> <nep> Run NEB calculation with a NEP model                   |
| minimize <structure> <nep_model> [fmax] [n]   Minimize a structure with a NEP model                   |
| msd <trajectory.xyz> <element> <dt_fs> [n]    Calculate MSD from an extxyz trajectory                 |
| nlist [script args...]                        Build neighbor lists                                    |
| disp [script args...]                         Calculate displacement from trajectory                  |
| avg-struct [script args...]                   Calculate averaged structure                            |
| oct-tilt [script args...]                     Calculate octahedral tilt                               |
| pol-abo3 [script args...]                     Calculate local polarization for ABO3                   |
+-------------------------------------------------------------------------------------------------------+

In interactive mode, choose 4) Calculators. The menu is:

+----------------------------------------------------------+
|                     CALCULATOR TOOLS                     |
+----------------------------------------------------------+
| 401) Calc ionic conductivity                             |
| 402) Calc properties by nep                              |
| 403) Calc descriptors of specific elements               |
| 404) Calc density of atomistic states (DOAS)             |
| 405) Calc nudged elastic band (NEB) by nep               |
| 406) Build neighbor list                                 |
| 407) Calc displacement from trajectory                   |
| 408) Calc averaged structure                             |
| 409) Calc octahedral tilt                                |
| 410) Calc polarization for ABO3                          |
| 411) Minimize structure by nep                           |
| 412) Calc mean square displacement (MSD) from trajectory |
| 413) Calc XRD from extxyz trajectory                     |
| 414) Calc phonon band structure                          |
+----------------------------------------------------------+
| 000) Return to the main menu                             |
+----------------------------------------------------------+
Input the function number:

Ionic Conductivity

calc_ion_conductivity.py calculates ionic diffusivity and conductivity from msd.out.

Required and Optional Files

File Role
msd.out Required MSD data
thermo.out Optional, used for automatic temperature detection
model.xyz Optional, used for volume and ion-count detection
run.in Optional, used to detect replication

Usage

gpumdkit.sh -calc ionic-cond Li 1
Argument What it identifies Check before running
element the mobile species whose MSD is analysed The symbol must match the species represented in msd.out and model.xyz.
charge the charge magnitude used in the Nernst–Einstein conversion Supply the charge definition appropriate to your system. The implementation uses its square, so changing only the sign does not change the reported conductivity.

The example Li 1 is only an argument-format example; it is not a default choice for every system. The current implementation fits the middle 40–80% of the available MSD points. Inspect the MSD curve and its convergence before using the resulting diffusion coefficient or conductivity in an interpretation.

From interactive mode, choose 401. You will see:

>-------------------------------------------------<
| This function calls the script in calculators   |
| Script: calc_ion_conductivity.py                |
| Developer: Zihan YAN (yanzihan@westlake.edu.cn) |
>-------------------------------------------------<
Input <element> <charge> (eg. Li 1)
------------>>

If automatic files are missing, the script will ask for temperature, volume, and ion count interactively.

Interactive prompts in manual mode look like:

Files 'thermo.out' and 'model.xyz' are not found.
Please provide the following values:
--------------------------->
Enter average temperature (in K):
Enter system volume (in A^3):
Enter number of ions:

Example Output

Diffusivity (D):
  D_x: 4.153e-07 cm^2/s
  D_y: 4.174e-07 cm^2/s
  D_z: 2.610e-07 cm^2/s
  D_total: 3.646e-07 cm^2/s
------------------------------
Ionic Conductivity:
  Sigma_x: 2.576e-02 mS/cm
  Sigma_y: 2.589e-02 mS/cm
  Sigma_z: 1.619e-02 mS/cm
  Sigma_total: 2.261e-02 mS/cm

Mean Square Displacement

calc_msd.py computes MSD directly from an extxyz trajectory.

gpumdkit.sh -calc msd dump.xyz Li 10
gpumdkit.sh -calc msd dump.xyz Li 10 5000

From interactive mode, choose 412. You will see:

Input <extxyz_file> <element_symbol> <dt_fs> [max_corr_steps]
  Optional argument: max_corr_steps (default: frame number)
Example: dump.xyz Li 10
------------>>

Arguments:

Argument Meaning
dump.xyz input trajectory
Li target mobile species
10 time interval between frames, in fs
5000 optional maximum correlation steps

max_corr_steps limits the largest time lag included in the correlation. If it is omitted, the script uses the number of available frames. It changes the calculation window, so select it from the length and sampling of your own trajectory rather than copying a value from another system.

Output:

  • msd.out

The beginning of msd.out is a text table with time and MSD columns. After generating it, use the plot commands below.

You can then plot:

gpumdkit.sh -plt msd
gpumdkit.sh -plt sdc
MSD plot SDC plot

X-ray Diffraction (XRD)

calc_xrd.py calculates and averages LAMMPS-compatible XRD intensities from an extended XYZ trajectory. It is intentionally available through the interactive menu only; no CLI flag is provided.

Choose 4) Calculators, then 413) Calc XRD from extxyz trajectory. The Python page asks for the following values:

Input extended XYZ trajectory
Output XRD file
X-ray wavelength (Angstrom)
2theta range (degrees; min max)
Number of bins in this 2theta range
Elements to include (all or comma-separated) [all]
CPU workers (0 means automatic) [0]

The bins always cover the selected 2theta interval. For example, a range of 10 60 with 500 bins produces 500 equal-width bins from 10 to 60 degrees; there is no separate output-bin range.

Use all to include every atom, or enter element symbols such as Li or Li,Cl. Element matching is case-insensitive. The default calculation uses the input Lattice and pbc, standard LAMMPS scattering factors, and the Lorentz-polarization factor. Orthogonal cells are supported; triclinic cells are rejected because this calculator follows the current LAMMPS-compatible mesh convention.

The output contains metadata headers followed by four columns:

Column Meaning
Bin one-based histogram-bin index
Coord center of the selected 2theta bin, in degrees
Count averaged XRD intensity
Count/Total intensity normalized by the total selected-range intensity

For faster runs, the script reads the trajectory once, reuses the reciprocal mesh when the cell is unchanged, groups atoms by element, and supports ordered threaded workers. These optimizations preserve the selected atoms, scattering formula, and 2theta bin definition.

Phonon Band Structure

calc_phonon.py calculates force constants with a NEP model and evaluates the phonon band structure along the path defined in QPOINTS. It is available through the interactive menu only:

gpumdkit.sh -> 4) Calculators -> 414) Calc phonon band structure

The Python page asks for the following values:

Primitive cell structure [PRIMCELL.vasp]
NEP model [nep.txt]
QPOINTS file [QPOINTS]
Supercell dimensions [1 1 1]
Displacement distance [0.015]
Output phonon file [phonon_NEP.dat]

The QPOINTS file must use line-mode endpoint pairs. The number of points per segment is read from its second line, and the resulting output rows are later used by plt_phonon.py and plt_phonon_comp.py. Install phonopy before running this calculator:

pip install phonopy

NEP Property Prediction

calc_properties_with_nep.py calculates energy, force, and stress for structures using a NEP model.

Dependency: calorine

pip install calorine

gpumdkit.sh -calc nep structures.xyz predictions.xyz nep.txt

Use this function to run predictions with a trained NEP model. For best results, validate model quality on your target structures before relying on the output.

Tip: Before prediction, you may want to clean the extxyz metadata:

gpumdkit.sh -clean_xyz train.xyz clean_train.xyz

NEP Descriptors

calc_descriptors.py extracts NEP descriptors for a selected element.

gpumdkit.sh -calc des train.xyz descriptors.npy nep.txt Li

Use cases:

  • visualize chemical environments with PCA/UMAP;
  • compare training and candidate structures;
  • inspect whether new data expands descriptor space.

Plot descriptors with:

gpumdkit.sh -plt des pca
gpumdkit.sh -plt des umap
Descriptor UMAP

Density of Atomistic States

calc_doas.py calculates density of atomistic states (DOAS), following the idea proposed by Wang et al..

gpumdkit.sh -calc doas structures.xyz nep.txt doas.out
gpumdkit.sh -plt doas doas.out Li

The script:

  1. reads all structures;
  2. relaxes each structure with a NEP calculator;
  3. extracts per-atom energies;
  4. groups atomic energies by element;
  5. writes the grouped values to the output file.

For very large systems, running minimization and atomistic-energy extraction directly in GPUMD can be more efficient.

Density of atomistic states

NEB with a NEP Model

gpumdkit.sh -calc neb initial.xyz final.xyz 9 nep.txt

This runs a NEB calculation with 9 intermediate images. During execution, the script asks how atoms should be fixed:

  • none: no atoms fixed;
  • index: fix atoms by index;
  • element: fix all atoms of one element;
  • position: fix atoms inside a coordinate range.

Structure Minimization

calc_minimize.py minimizes a structure using a NEP model through calorine.

Dependency: calorine

pip install calorine

gpumdkit.sh -calc minimize POSCAR nep.txt 0.01 1000

Arguments:

Argument Meaning
POSCAR input structure, POSCAR/CONTCAR or extxyz
nep.txt NEP model
0.01 optional force convergence threshold in eV/Ang
1000 optional maximum optimization steps

Output:

  • minimize.xyz
  • minimize.log

RDF Calculation with OVITO

rdf_calculator_ovito.py calculates radial distribution function using OVITO's analysis tools.

Dependency: OVITO

pip install ovito

Input file: Structure file (single frame or trajectory)

python Scripts/calculators/rdf_calculator_ovito.py trajectory.xyz 6.0 400

Parameters:

Argument Meaning
trajectory.xyz Input structure file
6.0 Maximum distance for RDF calculation (Å)
400 Number of histogram bins

Visualization:

gpumdkit.sh -plt rdf

Note: It is recommended to use the compute_rdf command directly in GPUMD when possible.

Ferroelectric and Polar Material Tools

Options 406–410 are for perovskite and polar-material analysis. They are commonly used to extract local structural information from MD trajectories before analyzing phase transitions, domain patterns, or polarization textures.

This page only gives a quick index for these scripts because they are usually used together. For full workflows and argument details, see Polar Material Analysis.

Menu CLI subcommand Purpose Details
406 gpumdkit.sh -calc nlist ... Build neighbor lists between selected center and neighbor atoms Polar Material Analysis
407 gpumdkit.sh -calc disp ... Calculate local displacement from a trajectory and neighbor list Polar Material Analysis
408 gpumdkit.sh -calc avg-struct ... Calculate an averaged structure from a trajectory Polar Material Analysis
409 gpumdkit.sh -calc oct-tilt ... Calculate octahedral tilt angles Polar Material Analysis
410 gpumdkit.sh -calc pol-abo3 ... Estimate local ABO3 polarization from Born effective charges Polar Material Analysis

Some scripts require ferrodispcalc:

pip3 install git+https://github.com/MoseyQAQ/ferrodispcalc.git

A typical sequence is to build neighbor lists with nlist, then reuse those lists for displacement, tilt, or polarization calculations. For averaged structures, you can start directly from the trajectory and use the length tolerance to control which frames are included.