atomipy: the atom library in Python

A modular Python toolbox for handling and analyzing molecular structures, particularly for mineral slabs with periodic boundary conditions, PBC. This toolbox is a light version of the MATLAB atom library and can in particular be used to generate molecular topology files for the MINFF forcefield with a streamlined Python interface. An online molecular system and topology builder running www.atomipy.io is now available. Test cases for hydrated montmorillonite using the general and tailored MINFF parameters (angle force constant 500 kJ/mol/rad²) can be found in the example cases of the atom Toolbox.

The package now supports generating GROMACS n2t (atom name to type) files for both MINFF and CLAYFF forcefields, enabling seamless integration with GROMACS utilities like gmx x2top for enhanced topology handling.

Contents

• Overview
• Two ways to use atomipy (web-module app vs. Python library)
• Project Structure
• Common Variables
  • Structure Containers
  • Atom dictionary fields
• Configuration & Performance
• Key Features
• Requirements
• Installation
• MINFF atom typing & topology generation
• XRD pattern simulation
• CIF to PDB helper script
• Bond valence analysis
• Ionic/crystal radii and bond estimates (Shannon radii)

Overview

This toolbox is designed to import, export, and analyze molecular structures with a focus on mineral slabs containing the elements Si, Al, Fe, Mg, Ca, Ti, Li, F, O, H. It handles periodic and triclinic simulation cells, and provides functions for calculating bonds, angles, and distances while taking periodic boundary conditions into account, and hence is ideal for generating molecular topology files for mineral bulk/slab systems that can be modelled using the MINFF forcefield. However it also has the capability to handle clay minerals, hydroxides, and oxyhydroxides using CLAYFF (Cygan, R.T.; Liang, J.J.; Kalinichev, A.G. Molecular Models of Hydroxide, Oxyhydroxide, and Clay Phases and the Development of a General Force Field. J. Phys. Chem. B 2004, 108, 1255-1266).

Built-in atom typing for MINFF and CLAYFF is a core feature: you can assign atom types, charges, and generate ITP/PSF/LAMMPS topologies directly from structures. For quick runs without local setup, use the molecular system and topology builder at www.atomipy.io.

The molecular structure information is stored in dictionaries where each atom has fields for coordinates, neighbors, bonds, angles, element type, and more.

Two ways to use atomipy

atomipy can be used in two complementary ways that share the same core:

1. Interactively, through the web-module app — the easiest way to get started. The visual node-graph builder at www.atomipy.io (open source: atomipy-web-module) lets you build, solvate, force-field, simulate and analyze systems by wiring nodes together — no Python required. Every node maps onto atomipy functions, and the app shows and exports the exact Python script it runs, so any workflow can be reproduced locally or on Google Colab. It covers:

• Import or build structures: presets, a bundled crystal library, parametric lattices (make_lattice/build_cell), and a dedicated Water-models library (SPC/E, TIP3P, TIP4P, TIP5P; hexagonal-ice, grid, and equilibrated boxes).
• Solvate, add ions, isomorphic substitution, supercell replication, box/PBC and coordinate tools.
• MINFF / CLAYFF force-fielding and topology export (GROMACS, LAMMPS, …).
• Simulation with OpenMM (CPU/GPU) or a local GROMACS engine (grompp + mdrun) — energy minimization and NVT/NPT, chainable in any order, with a fully editable .mdp.
• Trajectory analysis nodes: RDF g(r) + running coordination n(r), density profiles (x/y/z), MSD / self-diffusion, VACF / vibrational power spectrum, and hydrogen-bond analysis — plus interactive plotting, a 3D viewer, and an Inspector that reports the variables/files at any point in the graph.
• GROMACS on a Colab GPU — enable the engine inside the Colab app via the launcher's optional GROMACS cell.

2. As a Python library — full scripting control. Import atomipy and call its functions directly, as shown throughout this README. This is the way to automate, batch-process, and integrate atomipy into your own pipelines.

The web-module is the fastest way to get started and to discover the API (download the generated script to learn the calls); the library is the way to automate and extend it. Both run the same atomipy code.

Project Structure

atominpython/                # Repository root
├── atomipy/                 # Core package files
│   ├── data/                # Reference data (Shannon radii, BV params)
│   ├── ffparams/            # Force field parameters (GMINFF/TMINFF JSONs & ITPs)
│   ├── structures/          # Bundled library of mineral structures
│   ├── gromacs/             # Local GROMACS engine (detect gmx, .mdp, grompp/mdrun, trjconv, energy)
│   ├── __init__.py          # Package init, public API exports
│   ├── analysis.py          # High-performance analysis (RDF, CN, unwrap)
│   ├── bond_angle.py        # Topology analysis (bonds, angles, dihedrals)
│   ├── bond_valence.py      # Bond valence sum analysis (BVS, GII)
│   ├── build.py             # Structure manipulation, merge & substitutions
│   ├── cell_utils.py        # Box_dim ↔ Cell conversion utilities
│   ├── charge.py            # Charge assignment logic
│   ├── composition.py       # System composition / atom classification (water/ion/organic)
│   ├── dcd.py               # Zero-dependency pure-Python DCD trajectory reader
│   ├── diffraction.py       # XRD pattern calculation
│   ├── distances.py         # O(N) Sparse & O(N²) Direct distance dispatchers
│   ├── dummy_mineral.py     # Dummy FF (frozen framework) for non-MINFF/CLAYFF inorganics
│   ├── element.py           # Chemical element assignment
│   ├── forcefield.py        # MINFF/CLAYFF force field logic
│   ├── merge_top.py         # Merge mineral + organic/water/ion topologies into one .top
│   ├── move.py              # Translate, rotate, place, center, bend
│   ├── openmm_interface.py  # OpenMM topology/system loader
│   ├── oxidation.py         # Oxidation-state guessing (Shannon/ionic)
│   ├── replicate.py         # Supercell replication
│   ├── solvent.py           # Solvation & water molecule repairs
│   ├── transform.py         # Coordinate transformations (frac/cart/wrap)
│   ├── write_conf.py        # File exporters (PDB, GRO, XYZ, etc.)
│   ├── import_conf.py       # File importers (PDB, GRO, XYZ, CIF, etc.)
│   └── write_top.py         # Topology exporters (ITP, PSF, LMP)
├── scripts/                 # Example workflow scripts (run_*.py)
├── setup.py                 # Installation script
└── README.md                # This documentation

Common Variables

Understanding the core containers and fields used in atomipy makes it easier to follow the examples and extend the code. The tables below summarise the names, expected types, and roles of the variables you will encounter most often.

Structure Containers

Variable	Type	Description
atoms	list[dict]	Primary structure container; each dict holds coordinates, typing, charge, and connectivity for one atom.
Box	list[float] (len 3/6/9)	Canonical Cell dimensions accepted by most APIs; supports [lx, ly, lz], [a, b, c, alpha, beta, gamma], or the GROMACS-style [lx, ly, lz, 0, 0, xy, 0, xz, yz] in Angstrom and degrees.
Box_dim	list[float] (len 3/9)	Box dimensions returned by GRO import (1x3 orthogonal or 1x9 triclinic) and used for GRO/ITP/LAMMPS writers.
Cell	list[float] (len 6)	Unit-cell parameters [a, b, c, alpha, beta, gamma] in Angstrom and degrees.
Bond_index	numpy.ndarray or list	Bond records [i, j, distance] returned by bond_angle; indices are 0- or 1-based depending on export step.
Angle_index	numpy.ndarray or list	Angle records [i, j, k, angle, dx12, dy12, dz12, dx23, dy23, dz23] returned by bond_angle.
distances	numpy.ndarray (N, N)	Full pairwise distance matrix produced by dist_matrix.
dx, dy, dz	numpy.ndarray (N, N)	Cartesian component offsets that accompany distances.
cart_coords	numpy.ndarray (N, 3)	Cartesian coordinates created by transformation helpers such as fractional_to_cartesian.
frac_coords	numpy.ndarray (N, 3)	Fractional coordinates created by cartesian_to_fractional and wrapping utilities.
total_charge	float	Net system charge computed during forcefield or charge analysis.
atom_types	list[str]	Unique atom or forcefield type labels gathered from the atom dictionaries.

Atom dictionary fields

Key	Type	Description
index	int	1-based atom index assigned during import or when structures are combined.
molid	int	Molecule or residue identifier used to group atoms.
resname	str	Residue or component name (e.g. MMT, ION, SOL).
resid	int	Optional residue sequence number preserved from input files.
type	str	Working atom name used for matching templates and neighbour searches.
fftype	str	Forcefield-specific atom type assigned by minff, clayff, or related routines.
element	str	Chemical element symbol inferred by element.
x, y, z	float	Cartesian coordinates in Angstrom.
vx, vy, vz	float or None	Optional velocities in Angstrom per picosecond carried over from GRO files.
xfrac, yfrac, zfrac	float	Fractional coordinates stored when using the transformation helpers.
neigh	list[int]	Indices of bonded neighbours (0-based) populated by bond_angle.
bonds	list[tuple[int, float]]	(neighbour_index, distance_angstrom) pairs generated by bond_angle.
angles	list[tuple[tuple[int, int], float]]	((neighbour1, neighbour2), angle_deg) entries describing local angles.
cn	int	Coordination number saved when calculate_coordination=True.
charge	float	Partial charge in elementary charge units; initial imports may store strings before reassignment.
mass	float	Atomic mass in atomic mass units assigned via set_atomic_masses() or forcefield setup.
com	dict	Center of mass coordinates {'x': x, 'y': y, 'z': z} added by the com() function.
occupancy	float	PDB occupancy value (defaults to 1.0 when unavailable).
temp_factor	float	PDB temperature factor or B-factor (defaults to 0.0).
chain_id	str	PDB chain identifier used during export.
icode	str	PDB insertion code, if present.
is_nm	bool	Flag noting whether coordinates originated in nanometres before conversion to Angstrom.

Key Features

• Mixed Organic/Mineral Simulation Integration:
  • Parametrize organic molecules directly from SMILES using GAFF and CGenFF via openmmforcefields and ParmEd.
  • Mix organic structures seamlessly with CLAYFF/MINFF mineral slabs, handling combination rules and topological scaling automatically.
  • Dedicated Dockerized OpenFF worker microservice for SMIRNOFF-based (Sage, Parsley) parameterization.
• Robust Topology Merging & System Composition:
  • Merge multiple complex topology systems and linkages cleanly (merge_top.py).
  • Perform residue-level composition counting and system-wide mass calculations (composition.py).
• Support for multiple forcefields: MINFF and CLAYFF atom typing and parameter assignment
• Import/export PDB, Gromacs GRO, XYZ, CIF/mmCIF, PQR, POSCAR, SDF, and trajectories (multi-frame PDB/GRO plus GROMACS .xtc/.trr and .dcd import)
• CIF/mmCIF to PDB conversion helper (scripts/run_cif2pdb.py) with symmetry expansion
• Generation of GROMACS n2t (atom name to type) files for use with gmx x2top
• High-performance X-ray diffraction pattern simulation and plotting
• Distance matrices and neighbor lists with PBC corrections:
  • Extreme Scalability: Optimized $O(N)$ cell-lists for systems with >100,000 atoms.
  • Memory Efficiency: Sparse algorithms prevent RAM exhaustion on large-scale supercells.
  • Automatic Dispatch: Smart switching between Direct ($O(N^2)$) and Sparse ($O(N)$) methods.
• Local GROMACS engine (atomipy.gromacs): run grompp + mdrun directly from an atomipy system — energy minimization and NVT/NPT, MINFF/CLAYFF .mdp generation, define-aware force-field staging, and trajectory conversion. Used by the web-module's GROMACS Simulate node, including on a Colab GPU.
• Advanced Structural Analysis:
  • RDF/g(r) + running coordination n(r): radial distribution functions (single structure or ensemble-averaged over a trajectory) with the cumulative coordination number.
  • Coordination Numbers: Fast atom-by-atom environment analysis.
  • Density profiles (x/y/z): number / mass / charge density along a box axis.
  • MSD / self-diffusion: PBC-unwrapped, multi-time-origin mean-square displacement and the van Hove displacement distribution.
  • VACF / power spectrum: velocity autocorrelation (finite-difference), Green–Kubo diffusion and the vibrational density of states.
  • Hydrogen bonds: geometric donor–H···acceptor analysis (gmx-hbond convention) with per-molecule counts.
  • Unwrap: Restore molecular continuity for systems split by PBC.
• Progress tracking for computationally intensive calculations
• Consolidated charge module with support for formal, MINFF, and CLAYFF charge assignments
• Unified coordinate transformation system:
  • Cartesian-to-fractional and fractional-to-cartesian conversions
  • Orthogonal-to-triclinic and triclinic-to-orthogonal transformations
  • Direct transformation methods for crystallographic calculations
  • Coordinate wrapping for periodic boundary conditions
• Supercell replication with proper handling of triclinic cells
• Support for case-insensitive element matching in charge assignments
• Isomorphous substitution for creating defects in mineral structures:
  • Replace atom types in octahedral and/or tetrahedral sites (e.g., Al→Mg, Si→Al)
  • Ensure minimum separation distances between substituted atoms
  • Distribute substitutions evenly across the structure
  • Optional spatial limits for targeted substitution regions

Requirements

• NumPy (>=1.18.0)
• tqdm (>=4.45.0) - for progress bars
• Numba (>=0.50.0, optional) - for performance optimization via JIT compilation
• Matplotlib + SciPy (optional, required for XRD plotting and .mat export)
• GEMMI (>=0.7.0, required — core dependency; used for CIF/mmCIF import + symmetry expansion and the bundled inorganic crystal library)

Installation

If you're new to Python, follow these simple steps to get started with atomipy:

Step 1: Install Python

1. Download and install Python from python.org (version 3.7 or newer recommended)
2. During installation on Windows, make sure to check "Add Python to PATH"

Step 2: Install atomipy

Method 1: Install from GitHub (recommended)

1. Open a terminal or command prompt
2. Install Git if you don't have it already: git-scm.com
3. Clone the repository and install:

   git clone https://github.com/mholmboe/atomipy.git
   cd atomipy
   pip install -e .
   # Optional XRD extras (Matplotlib + SciPy):
   pip install -e ".[xrd]"
   # (gemmi for CIF/mmCIF is a core dependency, installed by default;
   #  the legacy ".[cif]" extra is now a no-op.)
   

Method 2: Manual Installation

1. Download this repository as a ZIP file (click the green "Code" button on GitHub and select "Download ZIP")
2. Extract the ZIP file to a folder on your computer
3. Open a terminal or command prompt
4. Navigate to the extracted folder:

   cd path/to/extracted/atomipy
   

5. Install the package and its dependencies:

   pip install -e .
   # Optional XRD extras (Matplotlib + SciPy):
   pip install -e ".[xrd]"
   # (gemmi for CIF/mmCIF is a core dependency, installed by default;
   #  the legacy ".[cif]" extra is now a no-op.)
   

Step 3: Verify Installation

To verify that atomipy is installed correctly, run this one-liner:

python - <<'PY'
import atomipy as ap
print("atomipy imported successfully, version:", ap.__version__)
PY

Configuration & Performance

atomipy is designed to handle systems ranging from small molecules to large mineral slabs with over 100,000 atoms. To balance speed and memory efficiency, it uses an $O(N)$ cell-list based configuration and dispatching system to control distance calculation thresholds.

The package features a central function, ap.distances.get_neighbor_list(), which acts as the decision-making "brain" for all distance-based calculations. Functions like bond_angle, solvate, and RDF all defer to this dispatcher.

Global Thresholds

By default, the dispatcher uses a size-based switching logic:

• Systems < 5,000 atoms: Use the Direct $O(N^2)$ method (faster for small systems due to NumPy vectorization).
• Systems ≥ 5,000 atoms: Use the Sparse Neighbor List $O(N)$ method (memory-efficient and scalable to hundreds of thousands of atoms).

Customizing Settings

You can easily customize this threshold or force a specific method at runtime:

import atomipy as ap

# Change the sparse method threshold globally
ap.config.SPARSE_THRESHOLD = 5000

# Or force a specific method for a single call
atoms, bonds, angles = ap.bond_angle(atoms, Box=cell, dm_method='sparse')

This ensures that even on memory-constrained systems, you can safely process large structures without risk of $O(N^2)$ memory bottlenecks.

Getting Started for Python Beginners

Here's a simple example to help you get started with atomipy:

# Create a file named my_first_atomipy.py with these contents:
import atomipy as ap

# Step 1: Load a structure file
print("Loading a GRO file...")
atoms, Box_dim = ap.import_gro("example.gro")  # Replace with your GRO file
print(f"Loaded {len(atoms)} atoms")

# Step 2: Assign elements based on atom names
print("Assigning elements to atoms...")
atoms = ap.element(atoms)

# Step 3: Calculate bonds and angles
print("Calculating bonds and angles...")
atoms, bonds, angles = ap.bond_angle(atoms, Box=Box_dim)

# Step 4: Save as a new file
print("Saving processed structure...")
ap.write_gro(atoms, Box=Box_dim, file_path="processed.gro")
print("Done!")

Run this script with:

python my_first_atomipy.py

Gromacs .n2t File Generation Example

Here's an example of generating a GROMACS n2t file for use with gmx x2top:

import atomipy as ap

# Load structure (auto-detect format)
auto_result = ap.import_auto("structure.gro")
if len(auto_result) == 2:
    atoms, box_like = auto_result
else:
    atoms, box_like = auto_result[0], auto_result[-1]

# Normalise Box (Cell -> Box_dim if needed)
if len(box_like) == 6:
    Box = ap.Cell2Box_dim(box_like)
else:
    Box = box_like

# Process with forcefield (optional)
atoms = ap.minff(atoms, Box)

# Generate n2t file (Box argument now precedes the optional output path)
n2t_path = ap.write_n2t(atoms, Box=Box, n2t_file="minff_atomtypes.n2t")
print(f"N2T file saved to: {n2t_path}")

The Box argument accepts orthogonal [lx, ly, lz] vectors, 1×6 Cell parameter lists, or the 1×9 Box_dim layout and will be normalised internally. This applies across writers (write_itp, write_psf, write_lmp, write_pdb, write_gro, write_xyz).

Dihedrals and 1–4 pair detection

Use bond_angle_dihedral to build bonds/angles and derive dihedrals plus 1–4 pairs in one call:

atoms, Box_dim = ap.import_gro("NMA_element.gro")
ap.element(atoms)
atoms, Bond_index, Angle_index, Dihedral_index, Pairlist = ap.bond_angle_dihedral(
    atoms, Box_dim, same_molecule_only=False, same_element_bonds=False
)
print(len(Bond_index), len(Angle_index), len(Dihedral_index), len(Pairlist))

Import existing .itp files

To read a GROMACS topology and inspect interaction sections:

from atomipy import import_itp_topology
itp = import_itp_topology("NMA.itp")
print("Pairs:", len(itp.get("pairs", {}).get("ai", [])))
print("Dihedrals:", len(itp.get("dihedrals", {}).get("ai", [])))

Alternatively, you can use the included MINFF helper script:

python scripts/run_minff2n2t.py structure.gro --output minff_atomtypes.n2t

For the simplest possible conversion you can call the structure-only helper:

python scripts/run_struct2n2t.py structure.gro

This script auto-detects the input format, forwards the Box dimensions, and writes <structure>.n2t.

MINFF atom typing & topology generation

Assign MINFF atom types/charges and write topology files directly from a structure:

import atomipy as ap

atoms, cell = ap.import_pdb("Kaolinite_GII_0.0487.pdb")
typed = ap.minff(atoms, Box=cell, log=True)

# Write out topologies
ap.write_itp(typed, Box=cell, file_path="minff_Kao.itp")
ap.write_psf(typed, Box=cell, file_path="minff_Kao.psf")
ap.write_lmp(typed, Box=cell, file_path="minff_Kao.data")

If your input is GRO, pass the returned Box_dim to minff/writers (they accept 1×3 Box_dim, 1×6 Cell, or 1×9 triclinic boxes interchangeably):

atoms, Box_dim = ap.import_gro("structure.gro")
typed = ap.minff(atoms, Box=Box_dim)
ap.write_itp(typed, Box=Box_dim, file_path="structure.itp")

Convenient helpers:

• scripts/run_create_itp_example.py shows a minimal end-to-end MINFF + ITP + typed PDB workflow.
• scripts/run_minff2n2t.py generates a MINFF-typed .n2t mapping for gmx x2top.

Force Field Parameters (JSON)

For greater flexibility and consistency (especially with LAMMPS), we now provide JSON parameter files for GMINFF and TMINFF in the atomipy/ffparams directory. You can load these using ap.load_forcefield:

import atomipy as ap

# Load simulation box and atoms
atoms, Box = ap.import_gro("system.gro")

# Load force field parameters (sigma, epsilon) from JSON
# This automatically handles unit conversion (Gromacs -> LAMMPS real by default)
ff = ap.load_forcefield(
    'GMINFF/gminff_all.json', 
    blocks=['GMINFF_k500', 'OPC3', 'OPC3_HFE_LM']
)

# Write LAMMPS data file with Pair Coeffs included
ap.write_lmp(atoms, Box, "system.data", forcefield=ff)

Available JSON files in ffparams:

• GMINFF: GMINFF/gminff_all.json (General MINFF)
• TMINFF: TMINFF/tminff_k*.json (Tailored MINFF for specific minerals)

See atomipy/ffparams/README.md for a full list of available blocks (minerals, ions, water models).

XRD pattern simulation

atomipy includes a fast X-ray diffraction module (atomipy.diffraction.xrd) that can turn a PDB/GRO/XYZ structure plus its box/cell into a calculated powder pattern with optional plotting and data export.

Command-line helper:

python scripts/run_xrd_example.py Kaolinite_GII_0.0487.pdb --two-theta 5 70 --save-output

Minimal API example:

from atomipy import import_auto, xrd

atoms, box = import_auto("Kaolinite_GII_0.0487.pdb")
two_theta, intensity, fig = xrd(
    atoms=atoms,
    Box=box,
    wavelength=1.54187,      # Cu K-alpha
    two_theta_range=(5, 70), # degrees
    angle_step=0.02,
    save_output=False,
)
fig.show()

CIF to PDB helper script

Use scripts/run_cif2pdb.py to prepare CIF/mmCIF structures for downstream MD workflows.

Typical usage:

# Keep CIF atom names, auto-expand symmetry, fuse overlaps, write PDB
python scripts/run_cif2pdb.py input.cif

# Add BVS protonation diagnostics
python scripts/run_cif2pdb.py input.cif --check-protonation

# Assign MINFF atom types/charges before writing PDB
python scripts/run_cif2pdb.py input.cif --assign-minff

What it supports:

• Symmetry expansion to full unit cell (default, opt-out with --no-expand-symmetry)
• Layered structure z-unfolding in fractional coordinates:
  • automatic detection (default, opt-out with --no-auto-layer-z-unfold)
  • manual override with --layer-z-unfold --z-split <value>
• Overlap fusion with --fuse-rmax and --fuse-criteria
• Optional BVS-based protonation-need check (--check-protonation)
• Optional MINFF typing (--assign-minff) while keeping original CIF atom names when disabled

Common Issues for Beginners

• ModuleNotFoundError: Make sure you've installed all required packages.
• File not found errors: Check that your file paths are correct and that the files exist.
• No module named 'atomipy': Make sure you've installed the package correctly.

Function Documentation

File I/O

• import_pdb(file_path): Import a PDB file, returning (atoms, Cell).
• import_gro(file_path): Import a Gromacs GRO file, returning (atoms, Box_dim). Coordinates and box are converted from nm to Å.
• import_xyz(file_path): Import an XYZ file, returning (atoms, Cell).
• import_cif(file_path, expand_symmetry=True): Import a CIF/mmCIF file using GEMMI, returning (atoms, Cell). Standardizes symmetry expansion for minerals.
• import_cjson(file_path, resname=None): Import a Chemical JSON (.cjson/.json) file (Avogadro2 / Open Chemistry format), returning (atoms, Cell). Carries 3D coordinates, bonds (with order), formal charges, and partial charges.
• import_pqr(file_path): Import a PQR file (charge and radius in place of occupancy/temp).
• import_poscar(file_path): Import VASP POSCAR/CONTCAR files (fractional or Cartesian).
• import_traj(file_path, top=None, stride=1, start=0, stop=None): Multi-frame importer for .pdb/.gro, GROMACS .xtc/.trr (via the optional libxdrfile, or a gmx trjconv fallback — pass top= a companion .gro/.pdb/.cif for atom names/types), and .dcd (built-in zero-dependency reader) / .nc/.h5/.lammpstrj (via optional mdtraj). stride/start/stop subsample frames.
• import_auto(file_path): Auto-detect format and import.
• write_pdb(atoms, Box, file_path, ...): Write atoms to a PDB file with formal charge and scale matrices.
• write_gro(atoms, Box, file_path): Write atoms to a Gromacs GRO file (coordinates in nm).
• write_xyz(atoms, Box, file_path): Write atoms to an XYZ file.
• write_cif(atoms, Box, file_path): Write atoms to a CIF file.
• write_pqr(atoms, Box, file_path): Write atoms to a PQR file.
• write_poscar(atoms, Box, file_path): Write atoms to a VASP POSCAR file.
• write_sdf(atoms, Box, file_path): Write atoms to an SDF file.
• write_cjson(atoms, Box, file_path, name=None): Write atoms to a Chemical JSON (.cjson) file (coordinates, atomic numbers, bonds, formal + partial charges, optional unit cell).
• write_traj(frames, file_path): Write multiple frames to a trajectory file.

Bundled Molecule Library

A curated library of ~428 small organic molecules (Chemical JSON), vendored from the Avogadro2 molecules library (BSD-3-Clause, © 2016 Geoffrey Hutchison, University of Pittsburgh) and filtered to molecules composed only of GAFF / OpenFF-parameterizable elements. Categories include amino acids, nucleobases, carbohydrates (furanose/pyranose), alcohols, fatty acids, steroids, aromatics, heteroaromatics, and most common functional groups. Stored under atomipy/structures/molecules/ with an index.json manifest.

• molecule_categories(): Sorted list of available category names.
• list_molecules(category=None): List molecules ({name, file, formula, natoms, category}), optionally filtered to one category.
• load_molecule(name_or_file, resname=None): Load a molecule by name ('L-alanine') or path ('amino_acids/L-alanine.cjson'), returning (atoms, Cell).

import atomipy as ap
ap.molecule_categories()             # ['alcohols', 'amino_acids', 'aromatics', ...]
ap.list_molecules('nucleobases')     # adenine, cytosine, guanine, thymine, uracil
atoms, cell = ap.load_molecule('L-tryptophan')   # 27 atoms, resname 'LTRY'

Force Field

• minff(atoms, Box, ffname='minff', rmaxlong=2.45, rmaxH=1.2, log=False, log_file=None): Assign MINFF forcefield specific atom types to each atom. Set log=True to generate structure statistics; optionally specify log_file for the output path.
• clayff(atoms, Box, ffname='clayff', rmaxlong=2.45, rmaxH=1.2, log=False, log_file=None): Assign CLAYFF forcefield specific atom types to each atom. Set log=True to generate structure statistics.
• load_forcefield(json_path, blocks=None, units='lammps'): Load non-bonded parameters from a JSON file in ffparams. Returns a dictionary suitable for write_lmp.
• write_n2t(atoms, Box=None, n2t_file=None, verbose=True): Generate a GROMACS n2t (atom name to type) file based on structural analysis, honouring periodic boundary conditions when a 1×3 Box, 1×6 Cell, or 1×9 Box_dim array is supplied and merging nearly identical environments

Molecular Topology

• write_itp(atoms, Box, file_path): Write a Gromacs topology file
• write_psf(atoms, Box, file_path): Write a NAMD topology file
• write_lmp(atoms, Box, file_path): Write a LAMMPS topology file

Atom Properties

• element(atoms): Assign chemical elements to atoms based on their atom types and residue names
• mass(): Returns a comprehensive dictionary of atomic masses for different elements (includes all elements from the periodic table)
• radius(): Returns a dictionary of van der Waals radii for different elements
• set_atomic_masses(atoms): Set mass attributes for each atom in the list based on element type
• com(atoms, add_to_atoms=True): Calculate the center of mass of a molecule or slab without PBC wrapping, optionally adding COM coordinates to each atom dictionary

Structure Analysis

• bond_angle(atoms, Box, ...): Compute bonds and angles. Automatically scales between Direct and Sparse methods based on system size.
• dist_matrix(atoms, Box): Calculate a full distance matrix (N x N) for small systems.
• cell_list_dist_matrix(atoms, Box): Full cell-list based distance analysis returning sparse arrays.
• neighbor_list_fast(atoms, Box, cutoff): Extreme-performance sparse neighbor list ($O(N)$), capable of handling > 1,000,000 atoms.
• calculate_rdf(atoms, Box, rmax, dr, typeA=None, typeB=None, return_cn=False): radial distribution function; with return_cn=True also returns the running coordination number n(r).
• coordination_number(atoms, Box, cutoff): Efficient per-atom coordination analysis via sparse indexing.
• closest_atom(atoms, reference, Box): Find the neighbor closest to a target point or atom.
• find_H2O(atoms, Box_dim=None, rmin=1.25): Automated water discovery and molecule ID assignment.
• get_structure_stats(atoms, Box, ...): Comprehensive geometry and topology diagnostics.

Trajectory analysis (operate on a list of (atoms, Box) frames, e.g. from import_traj):

• density_profile(atoms, Box, axis='z', nbins=100, atom_types=None, mode='number'): 1-D density profile (number atoms/ų, mass g/cm³, or charge e/ų). density_frames(frames, ...) ensemble-averages it.
• rdf_frames(frames, typeA=None, typeB=None, rmax=15.0, dr=0.1, return_cn=False): ensemble-averaged g(r) (and optional n(r)) over a trajectory.
• msd(frames, atom_types=None, dims='xyz', dt=1.0, origin_stride=1): mean-square displacement and self-diffusion D (Einstein); PBC-unwrapped, multiple time origins; dims selects 3D/2D('xy')/1D('z'). displacement_distribution(frames, ...) gives the van Hove self-part.
• vacf(frames, atom_types=None, dt=1.0, ...): velocity autocorrelation (finite-difference velocities), Green–Kubo diffusion D, and the power spectrum / vibrational DOS (cm⁻¹).
• find_hbonds(atoms, Box, donor_types=None, acceptor_types=None, donor_resnames=None, acceptor_resnames=None, r_cut=3.5, angle_cut=30.0) / hbonds_frames(frames, ...): geometric hydrogen-bond detection (GROMACS gmx-hbond convention) with per-molecule counts.

Coordinate Transformations

• orthogonal_to_triclinic(...): Convert coordinates between basis sets.
• triclinic_to_orthogonal(...): Restore orthogonal coordinates from triclinic boxes.
• cartesian_to_fractional(...): Project atoms into unit-cell fractional space.
• fractional_to_cartesian(...): Restore absolute coordinates in Angstroms.
• wrap(atoms, Box): Wrap all atoms into the primary simulation cell.
• unwrap_coordinates(atoms, Box): Reconstruct split molecules across periodic boundaries.
• get_cell_vectors(Box): Derives absolute lattice vectors (H matrix).
• scale(atoms, Box, scale_factors): Rescale coordinates and box together. A scalar (or length-1 vector) scales isotropically; a length-3 [sx, sy, sz] scales each axis (and cell edge) independently. Works for orthogonal and triclinic cells (angles preserved).

Diffraction

• xrd(atoms, Box, wavelength=1.54187, angle_step=0.02, two_theta_range=(2, 90), ...): Calculate and optionally plot/save an XRD powder pattern using the high-performance diffraction module (see scripts/run_xrd_example.py for CLI usage)
• get_orthogonal_box(Box): Get orthogonal Box dimensions from triclinic parameters
• replicate_system(atoms, Box, replicate=[1, 1, 1]): Create supercells by replicating in a, b, c directions

Structure Building

• substitute(atoms, Box, num_oct_subst, o1_type, o2_type, min_o2o2_dist, ...): Perform isomorphous substitution by replacing atom types (e.g., Al→Mg in octahedral sites). Returns (atoms, Box, None). Advanced features include:
  • Automatic detection of centrosymmetric structures for proper centering during substitution
  • Fallback to element-based matching when atom type names aren't found (using the element module)
  • Even distribution of substitutions across the structure
  • Support for both octahedral and tetrahedral substitutions with minimum distance constraints
  • Optional spatial limits for targeting substitutions to specific regions
• molecule(atoms, molid=1, resname=None): Assign molecule ID and optionally residue name to all atoms in a list. Useful for grouping atoms as a single molecular unit before topology generation.
• fuse_atoms(atoms, Box, rmax=0.5, criteria='average'): Removes overlapping atoms within a cluster. Designed for cleaning up disordered CIFs by keeping the atom with highest occupancy, preserving order, or computing an average geometric center (matching the MATLAB library behavior).
• merge(atoms1, atoms2, Box, type_mode='molid', atom_label=None, min_distance=None): Merge two atom lists by removing atoms from atoms2 that are too close to atoms1. Supports molecule-aware removal via type_mode='molid'.
• solvate(limits, density=1000.0, min_distance=2.0, max_solvent='max', solvent_type='spce', ...): Fill a region with water/solvent. The pre-equilibrated template is isotropically scaled to the target density, solvent–solvent overlaps (tile-seam/PBC clashes) are removed, and molecules get contiguous molids. max_solvent may be 'max', a 'shellN' shell, or an exact integer count (raises if the box can't hold it).
• ionize(ion_type, resname, limits, num_ions, ...): Add ions to a system within specified region limits.
• insert(molecule_atoms, limits, ...): Insert molecules into a system within specified region limits.
• slice(atoms, limits, remove_partial_molecules=True): Extract atoms within a region defined by limits.
• assign_resname(atoms, default_resname='MIN'): Assign residue names to atoms based on their types. Assigns 'SOL' to water atoms, 'ION' to ions, and the specified default name to other atoms.
• add_H_atom(atoms, Box, target_type, h_type='H', bond_length=0.96, coordination=1, max_h_per_atom=1): Add hydrogen atoms to under-coordinated atoms (e.g. protonating edge sites).
• adjust_H_atom(atoms, Box, h_type='H', neighbor_type='O', distance=0.96): Adjust bond lengths of hydrogen atoms, useful for fixing distorted bonds.
• adjust_Hw_atom(atoms, Box, water_resname='SOL', water_model='OPC3'): Repair water molecules by adding missing hydrogens and fixing geometry (OH distance and HOH angle) for models like OPC3/SPC/TIP3P.
• is_centrosymmetric_along_z(atoms, tolerance=0.1): Check if a structure is approximately centrosymmetric along the z-axis by analyzing the distribution of z-coordinates

Local GROMACS engine (atomipy.gromacs)

Run GROMACS directly from an atomipy system, reusing the same MINFF/CLAYFF topology writers. This powers the web-module's GROMACS Simulate node (including on a Colab GPU), but is usable standalone wherever a gmx binary is available.

• detect_gmx(gmx='gmx'): locate/validate a GROMACS install (a gmx binary, a GMXRC, or an install dir) and report its version.
• mdp(stage, ...) / build_defines(...): generate .mdp text for em/nvt/npt/md and the MINFF/CLAYFF -D defines.
• stage_run_dir(...) / stage_minff(workdir, defines=...): stage the run directory / copy the (define-aware) min.ff force field into it.
• write_freeze_ndx(workdir, n_frozen, n_total): write an index file ([System] + [frozen]) for the frozen Dummy FF framework, passed to grompp -n so freezegrps/freezedim hold it rigid (mdp(..., freeze_group='frozen') also sets comm-mode=None).
• run_stage(...), run_pipeline(...), run_local_gmx(workdir, top, gro, stages, ..., ndx=..., freeze_group=...): run grompp + mdrun for one or more stages; accepts a verbatim mdp_text to override the generated .mdp, and an ndx/freeze_group for a frozen framework.
• trjconv_to_pdb(workdir, tpr=..., xtc=..., out=..., pbc='atom', ...): convert an .xtc/.trr trajectory to a multi-frame PDB (box per frame) for viewing/analysis.
• trjconv(workdir, tpr=..., src=..., out=..., pbc='mol', ...): general format conversion via gmx trjconv — the output format follows the out extension (.xtc/.trr/.gro/.pdb).
• energy_timeseries(workdir, edr, terms=..., ...): parse .edr thermodynamics (potential/temperature/pressure/volume/density vs time) via gmx energy, returning a {time, series} dict (one series per requested term).

Cell Utilities

• cell_utils.Box_dim2Cell(Box_dim): Convert Box dimensions to Cell parameters
• cell_utils.Cell2Box_dim(Cell): Convert Cell parameters to Box dimensions

Charges

• get_formal_charge(element_or_type): Get the full formal oxidation state charge for a chemical element or atom type
• get_half_formal_charge(element_or_type): Get half the formal charge for an element or atom type (used historically in some clay modelling applications)
• assign_formal_charges(atoms): Assign formal charges to ions and water molecules based on residue names (SOL, ION) and element types
• charge_minff(atoms, Box, atom_labels=None, charges=None, resname=None): Assign MINFF charges to atoms based on coordination environment
• charge_clayff(atoms, Box, atom_labels=None, charges=None, resname=None): Assign CLAYFF charges to atoms based on coordination environment
• balance_charges(atoms, resname=None): Balance the total charge of the system by adjusting oxygen charges

Oxidation states (rules + charge balance)

guess_oxidation_states(atoms, total_charge=0, group_by='all', method='auto', write=True, key='oxidation_state', verbose=False) guesses per-atom oxidation states without needing accurate geometry — a parallel alternative to the geometry-based compute_bvs (bond-valence sums) and the force-field assign_formal_charges. Two engines:

• ionic — apply reliable oxidation states (group 1 → +1, group 2 → +2, Al/B/Ga/Sc/Y/La → +3, Si/Ge/Ti/Zr/Hf → +4, F → −1, O → −2, Cl/Br/I → −1, H → +1) then solve the remaining variable elements (C, N, S, Fe, Mn, …) so each charge group sums to its net charge (0 by default — a neutral lattice). Single variable elements are solved exactly (incl. fractional mixed valence, e.g. Fe₃O₄ → +8/3); multiple variable elements are resolved by enumerating common oxidation states for a charge-balanced, most-probable combination. Best for crystals / unit cells.
• electronegativity — the rigorous definition: each bond's electrons go to the more electronegative atom, so ox = Σ(±bond order) + formal charge. Needs connectivity (atom['bonds']/atom['neigh']). Best for molecules.

method='auto' (default) picks the electronegativity engine when bond orders are present, else the ionic engine. Results are written to atom['oxidation_state'].

import atomipy as ap
atoms, cell = ap.import_auto('UC_conf/Pyrophyllite_GII_0.0.pdb')
ap.guess_oxidation_states(atoms)        # Si=+4, Al=+3, O=-2, H=+1 (neutral)
atoms, _ = ap.load_molecule('L-alanine')
ap.guess_oxidation_states(atoms)        # methyl C=-3, alpha C=0, carboxyl C=+3, N=-3

The Dummy FF — frozen model for inorganics not covered by the built-in force fields

For an inorganic whose framework the built-in force fields (CLAYFF / MINFF) can't type (e.g. MnO, NiO, Cr₂O₃, Ag/As/Mn oxides), assign_dummy_mineral_params(atoms, charge_mode='pauling', metal_site='Alo') builds a crude qualitative model. Charges follow one of two modes:

• pauling (default) — each cation gets a Pauling effective charge q_eff = oxidation × [1 − exp(−¼(χ_O − χ_M)²)] (the bracket is the bond's fractional ionic character): Si +1.79, Al +1.70, Mg +1.36, Ti +2.38, Fe²⁺ +0.95 / Fe³⁺ +1.43. Hydrogen is fixed at +0.4, and the anion (O, F) charges follow the coordination-resolved formula q_O = oxidation + Σ_j (oxidation_j − partial_j)/CN_j over each coordinating cation/H, so each framework stays neutral (ap.pauling_effective_charge(ox, element) exposes the cation formula).
• half — legacy charge_scale × oxidation state for every atom.

Lennard-Jones parameters come from lj_mode (all van der Waals data is from UFF, Rappé 1992):

• shannon (default) — oxygen (any) gets the OPC3 water-oxygen LJ (σ=0.31743 nm, ε=0.68369 kJ/mol); hydrogen gets zero LJ; every other element M has its LJ minimum placed at the Shannon crystal M–O bond distance d_MO = r_M + r_O. Under Lorentz–Berthelot the M–O pair minimum is (r_min_M + r_min_O)/2, so r_min_M = 2·d_MO − r_min_O(OPC3) and σ_M = r_min_M / 2^(1/6). The well depth ε_M is the per-element UFF value, clamped to within one order of magnitude of the OPC3-oxygen ε. A very short bond (shorter than the OPC3 oxygen radius, e.g. tetrahedral Si⁴⁺) gives M no LJ (Coulomb only — a small, buried, shielded core). A pure-metal structure (no anions) uses per-element UFF radii instead. Needs Box for the coordination numbers used in the Shannon lookup.
• element — the Dummy FF's own per-element UFF LJ (σ = x_i/2^(1/6), ε = D_i·4.184) for every element including O/F/H. Each element gets its own size (Mn ≠ Al), no borrowing. ap.uff_lj(element) exposes the conversion.
• minff — borrow from MINFF: oxygen → OPC3-O, fluorine → F⁻, H → none, metals → a small buried site (default metal_site='Alo'). Stronger O–water attraction.

The framework is flagged frozen. Freezing means no bonded parameters are needed, so only nonbonded terms remain — the material interacts with water/solutes via electrostatics + LJ. In OpenMM the frozen particles get mass 0; in the local GROMACS engine they are held rigid via freezegrps/freezedim in the .mdp plus an index group written by ap.gromacs.write_freeze_ndx.

• assign_dummy_mineral_params(atoms, Box=None, charge_mode='pauling', lj_mode='shannon', metal_site='Alo', rmaxlong=2.45, rmaxH=1.2): assign per-atom charge/LJ + the frozen flag (in place). Pass Box (needed by shannon for coordination) and for the coordination-resolved oxygen charges. Returns (atoms, report).
• write_dummy_mineral_itp(atoms, 'dummy.itp', mol_name='DUM'): a self-contained .itp (own [ atomtypes ] + bond-free [ moleculetype ]) that #includes like an organic itp.
• write_dummy_system_top(atoms, box, out_top, out_gro, water_model='spce', organic_itps=None): a complete .top + .gro for a frozen framework plus organics, water and ions — multiple dummy minerals, several different organics (each #included), and Na/Cl etc. are all supported; returns (ordered_atoms, n_frozen) so the caller freezes the leading framework particles. Atoms are ordered framework → organics → ions → water (SOL last).

Run EM / NVT only (a frozen rigid body is incompatible with an NPT barostat). Intended for qualitative questions (wetting, ion adsorption), not quantitative energetics. In the web builder, select Forcefield → Dummy FF; the structure node's Preview & Validate scan flags elements with no built-in force-field type and steers you here.

Data Structure

All atomic information is stored in a list of dictionaries called atoms. Each atom dictionary contains the following fields:

• molid: Molecule ID
• index: Atom index
• resname: Residue name
• x, y, z: Coordinates
• vx, vy, vz: Velocities (if present)
• neigh: List of neighbor indices
• bonds: List of pairs (j, distance) where j is the index of a bonded atom
• angles: List of pairs ((j, k), angle) where j,k are indices of atoms forming an angle with the central atom
• element: Chemical element
• type: Atom type
• fftype: Force field specific atom type
• cn: Coordination number

The simulation Cell can be represented in three ways with the unified Box parameter:

• Orthogonal box: A 1x3 array [lx, ly, lz] for simple rectangular boxes
• Cell parameters: A 1x6 array [a, b, c, alpha, beta, gamma] (used in PDB files)
• Triclinic box: A 1x9 array [lx, ly, lz, 0, 0, xy, 0, xz, yz] (used in Gromacs GRO files for triclinic cells)

Conversion utilities Box_dim2Cell() and Cell2Box_dim() can be used to convert between these formats.

Common Workflow Examples

Here are some common workflows that demonstrate how to use Atomipy for specific tasks:

Box Parameter Handling

Atomipy now uses a standardized approach for handling simulation Box parameters across all functions:

• All functions accept a generalized Box parameter and normalise internally

• The Box parameter supports three formats:

  1. Orthogonal box: A 1x3 array [lx, ly, lz] for simple rectangular boxes
  2. Cell parameters: A 1x6 array [a, b, c, alpha, beta, gamma] for crystallographic notation
  3. Triclinic box: A 1x9 array [lx, ly, lz, 0, 0, xy, 0, xz, yz] using GROMACS triclinic Box format

• Conversion utilities are automatically applied internally based on the format provided

Example usage with the different formats:

# Using orthogonal Box format
Box_ortho = [30.0, 30.0, 30.0]  # lx, ly, lz in Angstroms
atoms, bonds, angles = ap.bond_angle(atoms, Box=Box_ortho)

# Using Cell parameters format
Box_cell = [30.0, 30.0, 30.0, 90.0, 90.0, 90.0]  # a, b, c, alpha, beta, gamma
atoms, bonds, angles = ap.bond_angle(atoms, Box=Box_cell)

# Using triclinic Box format
Box_triclinic = [30.0, 30.0, 30.0, 0.0, 0.0, 5.0, 0.0, 5.0, 2.0]  # GROMACS format
atoms, bonds, angles = ap.bond_angle(atoms, Box=Box_triclinic)

High-Performance Structural Analysis

For large systems (>100,000 atoms), use the optimized analysis functions which leverage sparse neighbor lists.

import atomipy as ap
import matplotlib.pyplot as plt

# Load a large mineral-water system
atoms, Box = ap.import_gro("large_system.gro")

# 1. Calculate Radial Distribution Function g(r)
# Only compute for Si-O pairs up to 15 Å
r, g_r = ap.calculate_rdf(
    atoms, Box, rmax=15.0, dr=0.1, 
    pair_types=(['Si'], ['O'])
)

plt.plot(r, g_r)
plt.xlabel("r (Å)")
plt.ylabel("g(r)")
plt.show()

# 2. Calculate coordination numbers
# Count how many Oxygen atoms are within 2.5 Å of each Silicon atom
counts = ap.coordination_number(
    atoms, Box, cutoff=2.5, 
    atom_types=['Si'], neighbor_types=['O']
)
print(f"Average Si-O coordination: {sum(counts)/len(counts):.2f}")

# 3. Fix split molecules (Unwrapping)
# If molecules are broken across periodic boundaries, cluster them back together
unwrapped = ap.unwrap_coordinates(atoms, Box)
ap.write_pdb(unwrapped, Box, "unwrapped.pdb")

Basic Structure Processing

import atomipy as ap

# Load structure
atoms, Box_dim = ap.import_gro("my_structure.gro")

# Assign elements
atoms = ap.element(atoms)

# Calculate bonds and angles
atoms, bonds, angles = ap.bond_angle(atoms, Box=Box_dim)

# Save processed structure
ap.write_gro(atoms, Box=Box_dim, file_path="processed.gro")

Calculating Center of Mass

import atomipy as ap

# Load structure
atoms, Box_dim = ap.import_gro("my_molecule.gro")

# Ensure atoms have element information (needed for mass assignment)
atoms = ap.element(atoms)

# Calculate center of mass and add it to each atom
com_coords = ap.com(atoms, add_to_atoms=True)
print(f"Center of mass coordinates: x={com_coords[0]:.3f}, y={com_coords[1]:.3f}, z={com_coords[2]:.3f} Å")

# Access the COM stored in each atom
print(f"COM from first atom: {atoms[0]['com']}")

# Calculate COM without adding to atoms
com_coords_only = ap.com(atoms, add_to_atoms=False)

Creating a Topology File for Molecular Dynamics

import atomipy as ap

# Load structure
atoms, Box = ap.import_gro("my_mineral.gro")

# Process the structure (elements, bonds, etc.)
atoms = ap.element(atoms)

# Assign forcefield atom types (choose either MINFF or CLAYFF)
# For MINFF:
ap.minff(atoms, Box=Box)
# OR for CLAYFF:
# ap.clayff(atoms, Box=Box)

# Write a topology file for different simulation programs
# For GROMACS:
ap.write_itp(atoms, Box=Box, file_path="topology.itp")

# For LAMMPS:
ap.write_lmp(atoms, Box=Box, file_path="topology.data") 

# For NAMD:
ap.write_psf(atoms, Box=Box, file_path="topology.psf")

Creating a Supercell

import atomipy as ap

# Load structure
atoms, Box = ap.import_gro("unit_cell.gro")

# Process atoms
atoms = ap.element(atoms)

# Create a 2x2x2 supercell
replicated_atoms, new_Box, new_Cell = ap.replicate_system(atoms, Box, replicate=[2, 2, 2])

# Save the supercell
ap.write_gro(replicated_atoms, Box=new_Box, file_path="supercell.gro")

Isomorphous Substitution in Clay Minerals

import atomipy as ap

# Load a clay mineral structure (e.g., pyrophyllite or montmorillonite)
atoms, Cell = ap.import_pdb("Pyrophyllite_GII_0.071.pdb")

# The substitute function now automatically detects if a structure is centrosymmetric
# and centers it appropriately during substitution (important for pyrophyllite)

# Perform octahedral substitution: Replace 16 Al atoms with Mg atoms
# Using the mineral-specific atom types (will automatically look up by element if type not found)
atoms, Cell, _ = ap.substitute(
    atoms, 
    Cell, 
    num_oct_subst=16, 
    o1_type='Alo',   # Uses 'Alo' atom type for octahedral Al
    o2_type='Mgo',   # Will become Mg in octahedral position
    min_o2o2_dist=5.2  # Minimum distance between Mg atoms in Angstroms
)

# Perform both octahedral and tetrahedral substitutions
# Replace 4 Al with Mg (octahedral) and 8 Si with Al (tetrahedral)
# If atom types aren't found, it will try to match by element instead
atoms, Cell, _ = ap.substitute(
    atoms, 
    Cell, 
    num_oct_subst=4, 
    o1_type='Al',    # Will find atoms by element if 'Al' type not found
    o2_type='Mgo', 
    min_o2o2_dist=5.5,
    num_tet_subst=8,
    t1_type='Si',
    t2_type='Alt',   # Will become Al in tetrahedral position
    min_t2t2_dist=5.5
)

# Limit substitutions to specific region (e.g., only z > 10 Å)
# The function ensures even distribution between upper and lower parts
atoms, Box, _ = ap.substitute(
    atoms, 
    Box, 
    num_oct_subst=4, 
    o1_type='Al', 
    o2_type='Mgo', 
    min_o2o2_dist=5.5,
    lo_limit=10.0,
    hi_limit=50.0,
    dimension=3  # 1=x, 2=y, 3=z (MATLAB-style indexing)
)

# Check if a structure is centrosymmetric along z (can be used independently)
is_symmetric = ap.is_centrosymmetric_along_z(atoms)
print(f"Is the structure centrosymmetric along z? {is_symmetric}")

# Save the substituted structure
ap.write_pdb(atoms, Box=Box_dim, file_path="substituted_clay.pdb")

Example Scripts

Example Scripts

The package includes example scripts that demonstrate comprehensive workflows for processing mineral structures using atomipy:

scripts/run_minff_atomi.py and scripts/run_clayff_atomi.py

These scripts demonstrate workflows for using MINFF and CLAYFF forcefields respectively. Both serve as excellent starting points for users new to the package.

What the script does:

1. Imports a structure file (scripts/run_minff_atomi.py uses Kaolinite_GII_0.0487.pdb, scripts/run_clayff_atomi.py uses Kaolinite_GII_0.0487.gro)
2. Assigns chemical elements to each atom using chemical knowledge-based rules
3. Creates a supercell by replicating the unit cell to a target size (about 30 Å in each dimension)
4. Saves the replicated structure in both GRO and PDB formats
5. Calculates bonds and angles based on distance criteria with periodic boundary awareness
6. Assigns specialized MINFF/CLAYFF atom types based on their chemical environment and coordination
7. Generates a molecular topology file (ITP) for use in molecular dynamics simulations
8. Writes the final structure with all assigned properties to output files

Running the scripts:

Simply execute python scripts/run_minff_atomi.py or python scripts/run_clayff_atomi.py from the command line. The first expects Kaolinite_GII_0.0487.pdb; the second expects Kaolinite_GII_0.0487.gro in the working directory.

Output files:

• replicated_structure.gro - The enlarged supercell in GROMACS format
• replicated_structure.pdb - The enlarged supercell in PDB format
• molecular_topology.itp - GROMACS topology file with bond and angle definitions
• preem.gro - Final structure with MINFF typing and calculated properties

Bond valence analysis

Compute bond valence sums (BVS) and the Global Instability Index (GII) using the IUCr parameter table (bvparm2020.cif):

import atomipy as ap

atoms, Box_dim = ap.import_gro("structure.gro")
params = ap.load_bv_params()  # optional, uses bvparm2020.cif by default
results, gii = ap.compute_bvs(atoms, Box_dim, params=params)
print("GII =", gii)
for atom in results[:5]:
    print(atom["index"], atom["element"], atom["bvs"], atom["expected_ox"], atom["delta"])

results contains per-atom bond valence sums, the inferred oxidation state, and the deviation (delta). Use ap.summarize_bvs(results) to get quick per-element averages and the worst-offending sites.

Ionic/crystal radii and bond estimates (Shannon radii)

Use the Revised Shannon radii table to fetch ionic/crystal radii and estimate bond distances:

import atomipy as ap

# Radii lookups
r_al = ap.get_radius("Al", 3, 6)                     # ionic radius
r_si_cryst = ap.get_radius("Si", 4, 4, prefer="crystal")

# Bond distance estimates (radius sums)
d_al_o = ap.bond_distance("Al", 3, 6, "O", -2, 4)
d_si_o = ap.bond_distance("Si", 4, 4, "O", -2, 4, use_crystal=True)

Other Examples

# Import the entire package
import atomipy as ap
import numpy as np

# Import a PDB file
atoms, Cell = ap.import_pdb("structure.pdb")

# Guess elements for all atoms
atoms = ap.element(atoms)

# Calculate bonds and angles
atoms, bonds, angles = ap.bond_angle(atoms, Box=Cell)

# Assign charges using one of the available methods
# For formal charges (ions and water):
atoms = ap.assign_formal_charges(atoms)
# For MINFF charges:
# atoms = ap.charge_minff(atoms, Box=Cell)
# For CLAYFF charges:
# atoms = ap.charge_clayff(atoms, Box=Cell)

# Convert to fractional coordinates
frac_coords, atoms = ap.cartesian_to_fractional(atoms, Box=Cell)

# Create a 2x2x1 supercell
replicated_atoms, new_Box, new_Cell = ap.replicate_system(atoms, Cell, replicate=[2, 2, 1])

# Convert from triclinic to orthogonal coordinates if needed
ortho_atoms, _, _ = ap.triclinic_to_orthogonal(replicated_atoms, Box=new_Box)


# Export to GRO format
ap.write_gro(replicated_atoms, Box=new_Box, file_path="structure.gro")

Molecular Dynamics Simulations with OpenMM

atomipy now supports setting up and running Molecular Dynamics (MD) simulations directly from Python using OpenMM! It automates the force field setup, system parameterization, platform selection, and execution of NVT/NPT molecular dynamics.

Key Features

• Automated Force Field Setup: Directly converts your mineral slab and solvated systems to OpenMM topology and system objects using the MINFF or CLAYFF force fields.
• Automatic Platform Selector: Intelligent GPU/CPU acceleration selector. It automatically detects and uses the fastest available platform on your machine (prioritizing CUDA, OpenCL, Metal, then falling back to CPU or Reference).
• Empty Box Safety: The system has robust edge-case handling. If a box has no atoms, it instantiates a dummy atom at (0, 0, 0) with zero charge and mass to guarantee visualizer and MD engine compatibility without throwing errors.
• Seamless GROMACS Integration: Generates top/itp files, validates bond/angle structures, handles periodic boundary wrapping, and executes the simulation without external dependencies.

Basic MD Simulation Example

Here is how you can set up and run a 5,000-step NVT molecular dynamics simulation:

import atomipy as ap
from atomipy.openmm_interface import load_minff_into_openmm, run_openmm_simulation

# 1. Load your mineral-water system
atoms, Box = ap.import_gro("system.gro")

# 2. Assign MINFF/CLAYFF force field types
atoms = ap.minff(atoms, Box=Box)

# 3. Write out files for GROMACS/OpenMM
ap.write_itp(atoms, Box=Box, file_path="system_minff.top.itp")
# write_top generates the main top file including the forcefield definition
ap.write_top(atoms, Box=Box, file_path="system_minff.top")

# 4. Initialize OpenMM Topology and System
# This function automatically handles special parameters (e.g. flexible water models OPC3/SPC/TIP3P, 
# dynamic mineral bond/angle parameters, and ion corrections)
topology, system, positions = load_minff_into_openmm(
    top_path="system_minff.top",
    gro_path="system.gro"
)

# 5. Run NVT MD Simulation
# This automatically selects the best available GPU/CPU acceleration platform (OpenCL, CUDA, etc.)
# and outputs a standard PDB trajectory
run_openmm_simulation(
    topology=topology,
    system=system,
    positions=positions,
    dt_fs=1.0,               # Time step (fs)
    temperature_k=298.15,    # Temperature (K)
    num_steps=5000,          # Steps
    report_interval=1000,    # Output frequency
    traj_pdb_path="trajectory.pdb",
    result_pdb_path="final_nvt.pdb"
)

---

Differences from atom MATLAB library

This Python implementation is designed to provide similar functionality to the MATLAB atom library while following Python's conventions and making use of NumPy for efficient numerical operations. The data structure is dictionary-based rather than struct-based, and the function interfaces are designed for Python's style.

Acknowledgements & External Tools

atomipy relies on and interoperates with several external projects, each the property of its respective authors and used under its own license:

• Numerics — NumPy, Numba (optional), tqdm.
• Structure & trajectory I/O — GEMMI (CIF/mmCIF, optional); GROMACS libxdrfile (.xtc/.trr) and MDTraj (optional binary trajectories). The bundled inorganic crystal library derives from the Avogadro collection.
• Simulation — OpenMM and GROMACS.
• Organic force fields — OpenFF Toolkit / Interchange (Sage, Parsley), openmmforcefields + ParmEd and/or ACPYPE + AmberTools antechamber (GAFF), with RDKit / Open Babel.
• Force fields — MINFF; CLAYFF (Cygan, Liang & Kalinichev, J. Phys. Chem. B 2004, 108, 1255); GAFF; OpenFF (Sage/Parsley).
• Ported from the MATLAB atom toolbox.

Disclaimer

atomipy and its bundled force-field implementations (e.g. MINFF, CLAYFF, and the Dummy FF) are provided as-is as a beta/research tool, with no warranty of any kind. There is no guarantee that the generated structures, topologies, force-field parameters, or analysis/simulation results are correct. The author accepts no responsibility or liability for the accuracy of atomipy or the force-field implementations, or for any use of their output. Users are responsible for verifying that all generated files and results are reasonable and suitable for their purposes. Use at your own risk.

License

This project is released under the MIT License.