AIMD / Trajectory Analysis

The AIMD / trajectory analysis module provides post-processing tools for VASP molecular-dynamics trajectories.

It is available from:

VQkit
 -> Post-processing
 -> [6] AIMD / trajectory analysis

The current implementation analyzes VASP trajectory data from:

XDATCAR
OSZICAR
OUTCAR
INCAR

The main menu contains:

[1] Trajectory summary
[2] Atomic position coordinates vs time (XDATCAR)
[3] Energy vs time
[4] Volume / lattice parameters vs time
[5] Bond distance vs time
[6] Mean-square displacement (MSD)
[7] Diffusion coefficient
[8] Radial distribution function (RDF)
[9] Coordination number vs time

[0] Back
[99] Exit VQkit

This page is the detailed post-processing reference. A shorter hands-on AIMD tutorial may be maintained separately under Tutorials.

General AIMD requirements

The main trajectory source is:

XDATCAR

Depending on the selected analysis, VQkit may also use:

OSZICAR
OUTCAR
INCAR

The working directory should contain the files from the same AIMD simulation.

Typical VASP AIMD settings include:

IBRION = 0
NSW    = ...
POTIM  = ...

The exact thermostat, ensemble, temperature-control, and electronic settings depend on the physical simulation and are not inferred by the post-processing module.

Time-step detection

VQkit attempts to determine the molecular-dynamics time step automatically.

The most important quantity is:

POTIM

in femtoseconds.

If the time step is detected, VQkit displays a prompt such as:

MD time step POTIM in fs [1, detected from INCAR]:

The detected value can be accepted or overridden.

If VQkit cannot determine the time step automatically, it asks:

MD time step POTIM in fs:

A positive value must be supplied.

Time convention

For XDATCAR-based trajectory analysis, VQkit uses:

\[t_i = N_i \Delta t\]

where:

  • \(N_i\) is the XDATCAR configuration number;

  • \(\Delta t\) is POTIM in femtoseconds.

Therefore:

time_fs = XDATCAR configuration number * POTIM

This convention is written explicitly in the generated summaries.

Variable-cell / NPT trajectories

For XDATCAR-based analyses, VQkit asks:

Variable-cell / NPT XDATCAR? [y/N]:

Choose:

y

when the trajectory contains frame-dependent lattice information.

For ordinary fixed-cell AIMD, press Enter or answer:

N

VQkit then uses the appropriate lattice for each frame.

Common XDATCAR prompt

Most trajectory analyses begin with:

XDATCAR file in current directory [XDATCAR]:

Press Enter when the standard filename is used.

Otherwise, enter another XDATCAR filename in the current directory.

Coordinate handling

VQkit converts trajectory coordinates into a common representation.

For each frame it can obtain:

  • fractional coordinates;

  • wrapped Cartesian coordinates;

  • PBC-unwrapped Cartesian coordinates.

The frame lattice is converted to Angstrom.

For VASP Direct coordinates:

\[\mathbf{r} = \mathbf{s}\mathbf{A}\]

where:

  • \(\mathbf{s}\) is the fractional coordinate;

  • \(\mathbf{A}\) is the lattice matrix.

Periodic unwrapping

For trajectory-dependent displacement analysis, VQkit removes artificial jumps caused by periodic boundary wrapping.

For a fractional displacement:

\[\Delta\mathbf{s} = \mathbf{s}_{i} - \mathbf{s}_{i-1}\]

VQkit applies the nearest periodic displacement:

\[\Delta\mathbf{s} \leftarrow \Delta\mathbf{s} - \mathrm{round}(\Delta\mathbf{s})\]

This allows atoms crossing a cell boundary to be tracked continuously.

[1] Trajectory summary

Choose:

[1] Trajectory summary

VQkit reads the XDATCAR and reports basic trajectory information.

The summary includes:

  • XDATCAR path;

  • number of frames;

  • number of atoms;

  • species and counts;

  • MD time step;

  • time-step source;

  • first configuration number;

  • last configuration number;

  • final trajectory time;

  • variable-cell status.

For variable-cell trajectories, it also examines the volume evolution and reports minimum, maximum, and mean volume.

Output directory

VQkit writes:

VQkit_AIMD_summary/

and inside it:

VQkit_AIMD_trajectory_summary.txt

Representative summary

A summary may contain:

VQkit AIMD trajectory summary
XDATCAR: /path/to/XDATCAR
Frames: 5000
Atoms: 64
Species: Li32 S32
Time step: 1 fs (INCAR)
First configuration: 1
Last configuration: 5000
Last-frame time: 5000 fs
Variable-cell mode: no

[2] Atomic position coordinates vs time

Choose:

[2] Atomic position coordinates vs time (XDATCAR)

This workflow extracts the trajectory of selected atoms.

Atom selection

VQkit asks for 1-based atom indices.

Examples:

1

1,5,10-12

The selection parser accepts:

  • individual indices;

  • comma-separated indices;

  • inclusive index ranges.

Every requested index must lie within:

1 .. number_of_atoms

Coordinate output

For each selected atom and frame, VQkit exports:

  • configuration number;

  • time;

  • global atom index;

  • element symbol;

  • element-local index;

  • wrapped fractional coordinates;

  • wrapped Cartesian coordinates;

  • PBC-unwrapped Cartesian coordinates.

The columns include:

configuration
time_fs
atom_index
element
element_index
frac_x
frac_y
frac_z
cart_x_A
cart_y_A
cart_z_A
unwrap_x_A
unwrap_y_A
unwrap_z_A

Output directory

Default:

VQkit_AIMD_positions/

Generated files include:

VQkit_XDATCAR_position_vs_time.csv
VQkit_XDATCAR_position_vs_time.dat
VQkit_XDATCAR_position_vs_time_summary.txt

For each selected atom, VQkit also generates position-vs-time plots.

Typical plot names are:

VQkit_AIMD_atom_1_position_vs_time.png
VQkit_AIMD_atom_1_position_vs_time.pdf

The plotted x, y, and z trajectories use the unwrapped Cartesian coordinates.

Why unwrapped coordinates matter

Suppose an atom crosses from:

fractional x = 0.98

to:

fractional x = 0.02

The wrapped trajectory appears to jump by approximately one full cell length.

The unwrapped trajectory instead follows the physically short crossing.

This is essential for diffusion and migration analysis.

[3] Energy vs time

Choose:

[3] Energy vs time

VQkit first tries to read AIMD records from:

OSZICAR

If suitable AIMD records are unavailable, it falls back to:

OUTCAR

OSZICAR parsing

For MD-style OSZICAR lines, VQkit can extract quantities including:

T
E
F
E0
EK

where available.

The exported columns are:

step
time_fs
temperature_K
E_eV
F_eV
E0_eV
EK_eV

OUTCAR fallback

If no suitable OSZICAR AIMD records are found, VQkit extracts:

free energy TOTEN

from OUTCAR.

Output directory

VQkit writes:

VQkit_AIMD_energy/

Typical outputs are:

VQkit_AIMD_energy_vs_time.csv
VQkit_AIMD_energy_vs_time.dat
VQkit_AIMD_energy_vs_time_summary.txt
VQkit_AIMD_energy_vs_time.png
VQkit_AIMD_energy_vs_time.pdf

The plot includes the available energy series.

Scientific interpretation

Energy fluctuations are expected in AIMD.

Their interpretation depends on the ensemble and thermostat/barostat.

Do not assess AIMD equilibration from energy alone.

Temperature, structure, volume, and relevant dynamical observables should also be examined.

[4] Volume / lattice parameters vs time

Choose:

[4] Volume / lattice parameters vs time

For each XDATCAR frame, VQkit calculates:

  • lattice lengths a, b, c;

  • angles alpha, beta, gamma;

  • cell volume.

The cell volume is:

\[V = |\det(\mathbf{A})|\]

The lattice lengths are reported in Angstrom and the angles in degrees.

Output directory

VQkit writes:

VQkit_AIMD_cell/

The exported table contains:

configuration
time_fs
a_A
b_A
c_A
alpha_deg
beta_deg
gamma_deg
volume_A3

Plots are generated for lattice evolution and volume evolution.

This analysis is especially useful for:

  • NPT trajectories;

  • variable-cell MD;

  • thermal expansion;

  • detecting cell instability.

For a fixed-cell trajectory, the lattice/volume curves should remain constant.

[5] Bond distance vs time

Choose:

[5] Bond distance vs time

VQkit asks for:

First atom index
Second atom index

The bond distance is calculated for every frame using periodic minimum-image geometry.

For fractional positions \(\mathbf{s}_i\) and \(\mathbf{s}_j\):

\[\Delta\mathbf{s} = \mathbf{s}_j - \mathbf{s}_i\]

with periodic wrapping:

\[\Delta\mathbf{s} \leftarrow \Delta\mathbf{s} - \mathrm{round}(\Delta\mathbf{s})\]

The Cartesian displacement is:

\[\Delta\mathbf{r} = \Delta\mathbf{s}\mathbf{A}\]

and the bond distance is:

\[d = |\Delta\mathbf{r}|\]

This avoids artificial distance jumps when one atom crosses a periodic boundary.

Output directory

The default output directory is:

VQkit_AIMD_bond/

The module writes numerical data and PNG/PDF distance-vs-time plots.

Bond-distance analysis is useful for:

  • bond formation;

  • bond breaking;

  • adsorption/desorption events;

  • coordination changes;

  • reaction-coordinate monitoring.

[6] Mean-square displacement (MSD)

Choose:

[6] Mean-square displacement (MSD)

Atom selection

The current AIMD group-selection syntax accepts:

all

an element symbol:

Li

or atom indices:

1,5,10-12

Multiple element symbols can also be supplied where supported by the selection parser.

MSD definition

For selected atoms:

\[\mathrm{MSD}(t) = \frac{1}{N} \sum_{i=1}^{N} |\mathbf{r}_i(t)-\mathbf{r}_i(0)|^2\]

VQkit uses:

  • the first XDATCAR frame as the reference;

  • PBC-unwrapped Cartesian trajectories.

This avoids artificial jumps at periodic boundaries.

Output directory

VQkit writes:

VQkit_AIMD_MSD/

Typical files are:

VQkit_AIMD_MSD.csv
VQkit_AIMD_MSD.dat
VQkit_AIMD_MSD_summary.txt
VQkit_AIMD_MSD.png
VQkit_AIMD_MSD.pdf

The data columns are:

configuration
time_fs
msd_A2

Variable-cell caution

For a variable-cell trajectory, the Cartesian MSD also contains displacement associated with cell deformation.

VQkit explicitly records this caution in the summary.

For strongly deforming NPT cells, interpret Cartesian MSD carefully.

[7] Diffusion coefficient

Choose:

[7] Diffusion coefficient

The diffusion analysis first calculates the MSD for the selected atoms.

Einstein relation

VQkit uses:

\[D = \frac{1}{2d} \frac{d\,\mathrm{MSD}(t)}{dt}\]

where:

  • \(D\) is the diffusion coefficient;

  • \(d\) is the diffusion dimensionality;

  • the derivative is obtained from a linear fit to the chosen part of the MSD curve.

Dimensionality

VQkit asks:

Diffusion dimensionality [3]:

Allowed values are:

1
2
3

Examples:

  • 3 for bulk 3D diffusion;

  • 2 for in-plane diffusion;

  • 1 for essentially one-dimensional transport.

Fit range

VQkit asks:

Fraction of trajectory to discard before linear fit [0.5]:

The default:

0.5

means the first 50% of the trajectory is excluded from the linear fit.

The fit is then applied to the remaining MSD points.

This is only a user-controlled heuristic.

VQkit does not automatically identify the physically diffusive regime.

Diffusion coefficient

If the fitted MSD slope is:

\[m = \frac{d\,\mathrm{MSD}}{dt}\]

then:

\[D = \frac{m}{2d}\]

VQkit first obtains:

D in Angstrom^2/fs

and converts it to:

cm^2/s

using:

\[1\ \mathrm{\AA^2/fs} = 0.1\ \mathrm{cm^2/s}\]

Therefore:

\[D_{\mathrm{cm^2/s}} = 0.1 D_{\mathrm{\AA^2/fs}}\]

Fit quality

VQkit also reports:

R^2

for the linear fit.

A high R^2 alone does not guarantee physically valid diffusion.

The selected region must also correspond to a genuine linear diffusive regime.

Output directory

VQkit writes:

VQkit_AIMD_diffusion/

Typical outputs are:

VQkit_AIMD_diffusion_fit.csv
VQkit_AIMD_diffusion_fit.dat
VQkit_AIMD_diffusion_summary.txt
VQkit_AIMD_diffusion_fit.png
VQkit_AIMD_diffusion_fit.pdf

The summary reports:

  • selected dimensionality;

  • fit-start fraction;

  • MSD slope;

  • \(R^2\);

  • diffusion coefficient in Angstrom²/fs;

  • diffusion coefficient in cm²/s.

Diffusion-analysis caution

A reliable diffusion coefficient normally requires:

  • equilibrated dynamics;

  • sufficiently long trajectory;

  • statistically meaningful number of diffusing atoms;

  • linear MSD regime;

  • appropriate dimensionality;

  • sufficient sampling of hopping events.

A short trajectory can give a numerically fitted D that is not physically meaningful.

[8] Radial distribution function (RDF)

Choose:

[8] Radial distribution function (RDF)

VQkit asks for:

Group A
Group B

Selections can be:

  • all atoms;

  • element symbol(s);

  • atom indices.

Examples:

Li

S

1-16

For unambiguous normalization, the two groups must be either:

  • identical;

  • disjoint.

Partially overlapping but non-identical groups are rejected.

RDF definition

VQkit bins periodic minimum-image pair distances and normalizes them to the ideal-shell population.

The shell volume between radii \(r_1\) and \(r_2\) is:

\[V_{\mathrm{shell}} = \frac{4\pi}{3} \left( r_2^3-r_1^3 \right)\]

The normalization is evaluated using the instantaneous frame volume.

This is important for variable-cell trajectories.

Maximum radius

VQkit asks:

Maximum r in Angstrom [automatic]:

If left blank, VQkit chooses a conservative maximum based on approximately half of the shortest lattice-vector length across the trajectory.

This reduces problems from sampling beyond the minimum-image range.

Bin width

The default is:

0.05 Angstrom

A smaller bin width gives finer resolution but noisier statistics.

A larger bin width smooths the RDF.

Output directory

VQkit writes:

VQkit_AIMD_RDF/

Typical files include:

VQkit_AIMD_RDF.csv
VQkit_AIMD_RDF.dat
VQkit_AIMD_RDF_summary.txt
VQkit_AIMD_RDF.png
VQkit_AIMD_RDF.pdf

The data columns are:

r_A
g_r
pair_counts

RDF interpretation

Typical interpretation:

  • first strong peak -> nearest-neighbor shell;

  • peak position -> characteristic pair distance;

  • peak width -> structural disorder / thermal broadening;

  • long-range oscillations -> persistent structural ordering.

For liquids, \(g(r)\) should usually approach approximately 1 at large distance when sampling and normalization are adequate.

[9] Coordination number vs time

Choose:

[9] Coordination number vs time

VQkit asks for:

Center atoms
Neighbor atoms
Coordination cutoff in Angstrom [3.0]

Both groups can be selected using:

  • all;

  • element symbols;

  • atom indices.

For each center atom and frame, VQkit counts neighboring atoms satisfying:

\[d_{ij} \le d_{\mathrm{cutoff}}\]

Self-pairs are excluded.

Periodic minimum-image distances are used.

Frame statistics

For each trajectory frame, VQkit reports:

  • mean coordination number;

  • minimum coordination number;

  • maximum coordination number.

Output directory

VQkit writes:

VQkit_AIMD_coordination/

Typical files include:

VQkit_AIMD_coordination_vs_time.csv
VQkit_AIMD_coordination_vs_time.dat
VQkit_AIMD_coordination_summary.txt
VQkit_AIMD_coordination_vs_time.png
VQkit_AIMD_coordination_vs_time.pdf

The data columns are:

configuration
time_fs
mean_coordination
min_coordination
max_coordination

Coordination cutoff

The default:

3.0 Angstrom

is only a starting value.

The cutoff should normally be chosen from:

  • a known bonding criterion;

  • a minimum after the first RDF peak;

  • chemical knowledge of the system.

Using an arbitrary cutoff can give misleading coordination numbers.

Atom and group selection

Different AIMD analyses use slightly different selection syntax.

Position analysis

Atomic position vs time uses explicit 1-based atom indices:

1
1,4,8-10

Group analyses

MSD, diffusion, RDF, and coordination analyses accept:

all

or element symbol(s):

Li

or atom indices:

1,5,10-12

This allows species-specific transport and structural analysis.

Example: Li-ion diffusion

Suppose an AIMD simulation contains Li atoms and a host lattice.

Step 1

Run:

Post-processing
AIMD / trajectory analysis
Mean-square displacement

Step 2

For the atom selection, enter:

Li

Step 3

Inspect:

VQkit_AIMD_MSD.dat
VQkit_AIMD_MSD.png

Step 4

Run:

Diffusion coefficient

Step 5

Again select:

Li

Step 6

Choose the correct diffusion dimensionality.

For bulk:

3

For a layered material with in-plane diffusion:

2

Step 7

Inspect the MSD region used for the fit.

Step 8

Record:

D
R^2

from:

VQkit_AIMD_diffusion_summary.txt

Example: bond-breaking trajectory

For a reaction involving atoms 12 and 18:

  1. Choose Bond distance vs time.

  2. Enter atom 12.

  3. Enter atom 18.

  4. Confirm POTIM.

  5. Plot the bond distance through the trajectory.

A rapid increase in distance may indicate bond cleavage.

Combine this analysis with:

  • energy vs time;

  • coordination number;

  • structural snapshots.

Example: liquid structure

For a liquid or solvated system:

  1. Run RDF for the pair of interest.

  2. Choose Group A and Group B by element.

  3. Use automatic r_max initially.

  4. Start with 0.05 Angstrom bin width.

  5. Identify the first RDF minimum.

  6. Use that minimum as a physically motivated coordination-number cutoff.

This provides a consistent link between RDF and coordination analysis.

Output directory summary

Trajectory summary

VQkit_AIMD_summary/

Atomic positions

VQkit_AIMD_positions/

Energy

VQkit_AIMD_energy/

Cell / volume

VQkit_AIMD_cell/

Bond distance

VQkit_AIMD_bond/

MSD

VQkit_AIMD_MSD/

Diffusion

VQkit_AIMD_diffusion/

RDF

VQkit_AIMD_RDF/

Coordination

VQkit_AIMD_coordination/

Plot generation

When Matplotlib is available, VQkit generates PNG and PDF plots for the relevant AIMD analyses.

The numerical .csv and .dat files are retained as the primary plot-ready outputs.

These can also be analyzed using:

  • Origin;

  • gnuplot;

  • Python;

  • spreadsheet software.

Important scientific cautions

Equilibration

Do not mix non-equilibrated trajectory regions with production sampling without justification.

Temperature stability

Check the temperature evolution before interpreting dynamic properties.

Trajectory length

Diffusion and RDF statistics may require substantially longer trajectories than simple structural inspection.

Finite-size effects

Small simulation cells can distort:

  • diffusion coefficients;

  • RDFs;

  • collective motion;

  • coordination statistics.

Time-step choice

AIMD dynamics depend on a physically appropriate integration time step.

Post-processing cannot correct an unstable or excessively large POTIM.

Variable-cell effects

NPT trajectories require extra care because the simulation cell itself evolves.

In particular, Cartesian MSD can include motion associated with cell deformation.

Diffusion regime

VQkit performs a linear fit over the user-selected portion of MSD.

It does not automatically determine whether that region is truly diffusive.

Always inspect the MSD curve.