[Bug]: empty index array `framework_indices` in `framework_disp = self.disp[framework_indices]`

[KylinGuo]

Email (Optional)

No response

Version

2025.11.15

Which OS(es) are you using?

• MacOS
• Windows
• Linux

What happened?

I attempted to plot the RMSD of Si using the DiffusionAnalyzer class. Physically, I expect the RMSD to increase roughly linearly with simulation time. However, the code raises the following error:

ValueError: zero-size array to reduction operation maximum which has no identity

Cause of the Error

The error occurs because framework_indices ends up being an empty list in the following code block (excerpt from DiffusionAnalyzer):

        indices: list = []
        framework_indices: list[int] = []  # improved here
        for i, site in enumerate(structure):
            if site.specie.symbol == specie:
                indices.append(i)
            else:
                framework_indices.append(i)
        if self.disp.shape[1] < 2:
            self.diffusivity = 0.0
            self.conductivity = 0.0
            self.diffusivity_components = np.array([0.0, 0.0, 0.0])
            self.conductivity_components = np.array([0.0, 0.0, 0.0])
            self.max_framework_displacement = 0
        else:
            framework_disp = self.disp[framework_indices]  # bug here
            drift = np.average(framework_disp, axis=0)[None, :, :]

and

            # Drift and displacement information.
            self.drift = drift
            self.corrected_displacements = dc
            self.max_ion_displacements = np.max(np.sum(dc**2, axis=-1) ** 0.5, axis=1)
            self.max_framework_displacement = np.max(self.max_ion_displacements[framework_indices])  # bug here
            self.msd = msd
            self.mscd = mscd
            self.haven_ratio = self.diffusivity / self.chg_diffusivity
            self.sq_disp_ions = sq_disp_ions
            self.msd_components = msd_components
            self.dt = dt
            self.indices = indices
            self.framework_indices = framework_indices

See source code here and here

XDATCAR file used to reproduce this bug:
XDATCAR.tar.gz

Why this happens

If the structure contains only the target species (Si in my case), then the framework_indices list is empty. As a result:

• framework_disp becomes a zero-length array.
• np.average() is called on this empty array.
• NumPy raises ValueError: zero-size array to reduction operation maximum which has no identity.

This makes DiffusionAnalyzer unusable for single-species systems or systems where all atoms are considered diffusing species.

Expected Behavior

DiffusionAnalyzer should:

• handle single-species systems gracefully, or
• provide a clear error message explaining that at least one “framework” atom is required for drift correction, or
• allow an option to disable drift correction when no framework atoms exist.

Suggested Fix

Before computing the drift, add a check such as:

# --- Fixed: guard against empty framework_indices ---
if len(framework_indices) == 0:
    # No framework atoms -> no drift correction from framework.
    # Make a zero drift array with the same time/coord dimensions as expected:
    _, nsteps, ndim = self.disp.shape
    framework_disp = np.zeros((0, nsteps, ndim), dtype=self.disp.dtype)
    drift = np.zeros((1, nsteps, ndim), dtype=self.disp.dtype)
else:
    framework_disp = self.disp[framework_indices]
    drift = np.average(framework_disp, axis=0)[None, :, :]

# --- Fixed: handle empty framework_indices for max retrieval ---
if len(framework_indices) == 0:
    self.max_framework_displacement = 0.0
else:
    self.max_framework_displacement = np.max(self.max_ion_displacements[framework_indices])

Code snippet

import numpy as np
from pymatgen.io.vasp import Xdatcar
from pymatgen.analysis.diffusion.analyzer import DiffusionAnalyzer as DiffusionAnalyzerPmg
import matplotlib.pyplot as plt


def plot_msd():
    xdatcar_file = "XDATCAR.tar.gz"
    xdatcar = Xdatcar(filename=xdatcar_file)

    specie = 'Si'
    diffusion_analyzer_pmg = DiffusionAnalyzerPmg.from_structures(
        structures=xdatcar.structures,
        specie=specie,
        temperature=2000,
        time_step=3,  # POTIM = 3.0
        step_skip=1,  # NBLOCK = 1
        smoothed=False
    )
    x = diffusion_analyzer_pmg.dt
    y = np.sqrt(diffusion_analyzer_pmg.msd)

    fig, ax = plt.subplots()
    ax.plot(x, y, label=specie)
    ax.legend()
    ax.set_xlabel(f'Time (fs)')
    ax.set_ylabel(f'RMSD (Angstrom)')
    ax.set_xlim(left=0, right=None)
    ax.set_ylim(bottom=0, top=None)
    plt.show()


def main():
    plot_msd()


if __name__ == '__main__':
    main()

Log output

Traceback (most recent call last):
  File "/home/USER/PycharmProjects/OVITO-Project/source/demo.py", line 38, in <module>
    main()
    ~~~~^^
  File "/home/USER/PycharmProjects/OVITO-Project/source/demo.py", line 34, in main
    plot_msd()
    ~~~~~~~~^^
  File "/home/USER/PycharmProjects/OVITO-Project/source/demo.py", line 12, in plot_msd
    diffusion_analyzer_pmg = DiffusionAnalyzerPmg.from_structures(
        structures=xdatcar.structures,
    ...<4 lines>...
        smoothed=False
    )
  File "/home/USER/miniconda3/envs/ovito_env/lib/python3.13/site-packages/pymatgen/analysis/diffusion/analyzer.py", line 657, in from_structures
    return cls(
        structure,
    ...<6 lines>...
        **kwargs,
    )
  File "/home/USER/miniconda3/envs/ovito_env/lib/python3.13/site-packages/pymatgen/analysis/diffusion/analyzer.py", line 377, in __init__
    self.max_framework_displacement = np.max(self.max_ion_displacements[framework_indices])
                                      ~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/USER/miniconda3/envs/ovito_env/lib/python3.13/site-packages/numpy/_core/fromnumeric.py", line 3164, in max
    return _wrapreduction(a, np.maximum, 'max', axis, None, out,
                          keepdims=keepdims, initial=initial, where=where)
  File "/home/USER/miniconda3/envs/ovito_env/lib/python3.13/site-packages/numpy/_core/fromnumeric.py", line 86, in _wrapreduction
    return ufunc.reduce(obj, axis, dtype, out, **passkwargs)
           ~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
ValueError: zero-size array to reduction operation maximum which has no identity

Process finished with exit code 1

Code of Conduct

• I agree to follow this project's Code of Conduct