This repository provides Python scripts for analyzing the spatial distributions of Li⁺ and I⁻ around PS₄³⁻ tetrahedra in Li₂S–P₂S₅–LiI glass trajectories.

The workflow aligns individual PS₄³⁻ units into a common reference frame, accumulates the positions of nearby Li⁺ ions and I⁻ ions, and exports three-dimensional spatial-density fields in Gaussian cube format.

The generated cube files can be visualized as isosurfaces using Mayavi.

Workflow

LAMMPS trajectory
      |
      v
dump2cube.py
      |
      +-- PS₄_Li.cube : Li⁺ spatial-density field
      |
      +-- PS₄_I.cube  : I⁻ spatial-density field
      |
      v
cube2mayavi.py
      |
      v
3D isosurface visualization

Files

File	Description
`dump2cube.py`	Reads a LAMMPS trajectory, aligns PS₄³⁻ units, and generates Li⁺/I⁻ cube files.
`cube2mayavi.py`	Reads a Gaussian cube file and visualizes the density field using Mayavi.
`050Li3PS4-050LiI.lammpstrj`	Example LAMMPS trajectory for a Li₂S–P₂S₅–LiI glass.
`050Li3PS4-050LiI_PS4_I.cube`	Example cube file containing the I⁻ spatial-density field.

Requirements

Python packages

The scripts require Python 3 and the following packages:

NumPy
Mayavi
VTK
Traits
PyQt5 or PySide6

Mayavi is most easily installed through conda-forge.

conda create -n ps4-density python=3.10
conda activate ps4-density
conda install -c conda-forge numpy mayavi pyqt

If NumPy is not installed, it can also be installed with pip.

pip install numpy

Input Trajectory Format

dump2cube.py expects a LAMMPS trajectory containing the following atom columns:

ITEM: ATOMS id type element xu yu zu mol

The trajectory must contain the following elements:

The script assumes that P and S atoms belonging to the same PS₄³⁻ unit share the same molecule ID (mol).

The current workflow has been tested primarily for orthorhombic simulation cells. Periodic boundary conditions are applied when calculating relative atomic positions.

Analysis Method

For every frame in the trajectory, the following procedure is performed.

Each P atom is selected as the center of a reference PS₄³⁻ tetrahedron.
The four nearest I⁻ ions around the selected P atom are identified.
The nearest I⁻ ion is used to define the orientation of the PS₄³⁻ unit.
The four S atoms belonging to the PS₄³⁻ unit are rotated into a common reference coordinate system.
Li⁺ and I⁻ ions within the specified cutoff distance from P are collected.
The accumulated Li⁺ and I⁻ positions are converted into three-dimensional histograms.
The histograms are exported in Gaussian cube format.

Definition of the reference orientation

The PS₄³⁻ tetrahedron is aligned using the nearest I⁻ ion.

The S atom farthest from the nearest I⁻ ion is aligned with the z axis.
A second S atom is used to define the rotation around the z axis.
The P atom is placed at the origin.

This alignment enables the spatial distributions of Li⁺ and I⁻ around many PS₄³⁻ units and trajectory frames to be accumulated in a common coordinate system.

Usage

Generate Cube Files

python dump2cube.py \
  -i 050Li3PS4-050LiI.lammpstrj \
  -m 160 160 160 \
  -cut 8

Arguments:

Argument	Description
`-i`, `--trjfile`	Input LAMMPS trajectory file.
`-m`, `--mesh`	Number of grid points along the x, y, and z directions.
`-cut`, `--cutoff`	Cutoff distance around the reference P atom in Å.

For the example above, the expected output files are:

050Li3PS4-050LiI_PS4_Li.cube
050Li3PS4-050LiI_PS4_I.cube

Visualize the I⁻ Spatial-Density Field

python cube2mayavi.py \
  -i 050Li3PS4-050LiI_PS4_I.cube \
  -atom I \
  -iso 1.26483e-09

Visualize the Li⁺ Spatial-Density Field

python cube2mayavi.py \
  -i 050Li3PS4-050LiI_PS4_Li.cube \
  -atom Li \
  -iso 1.0e-09

The appropriate isovalue depends on the number of trajectory frames, the number of PS₄³⁻ units, the mesh size, and the cutoff distance. Therefore, the -iso value should be adjusted for each dataset.

Output

Each generated cube file contains:

A reference PS₄³⁻ tetrahedron
- One P atom located at the origin
- Four S atoms in the aligned coordinate system
A three-dimensional spatial-density field for Li⁺ or I⁻

The cube files can also be viewed with software supporting Gaussian cube format, including:

Mayavi
VMD
ParaView
PyMOL

Visualization Settings

The default colors in cube2mayavi.py are:

Object	Color
P atom	Purple
S atom	Yellow
P–S bond	Gray
Li⁺ density	Blue
I⁻ density	Red

For Li⁺ visualization, the script separates the density field into an inner region and an outer region:

Inner region: within 5 Å of the P atom
Outer region: between 5 Å and 9 Å from the P atom

Example Interpretation

The I⁻ cube file can be used to examine whether I⁻ ions are distributed uniformly around PS₄³⁻ tetrahedra or occupy preferred directions.

An anisotropic I⁻ density distribution indicates that I⁻ ions form characteristic local geometries relative to PS₄³⁻ units rather than being randomly distributed in the glass network.

Notes and Limitations

The input trajectory must include id type element xu yu zu mol columns.
The P and S atoms of a PS₄³⁻ tetrahedron must have the same molecule ID.
The cube field is generated from accumulated histogram counts and is intended for relative spatial-distribution analysis.
The output file name is generated from the input trajectory name. The current implementation expects a filename containing LiI.
The visualization script opens an interactive Mayavi window. A GUI-capable Python environment is required.
Before use, confirm that the trajectory format and simulation-cell definition are compatible with the parser implemented in dump2cube.py.

License

This repository is intended for academic and research use.

README.org Unescape Escape

PS₄³⁻-Centered Li⁺/I⁻ Spatial Density Analysis in Li₂S–P₂S₅–LiI Glasses

Overview