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
• Project Structure
• Common Variables
  • Structure Containers
  • Atom dictionary fields
• Key Features
• Requirements
• Installation
  • Step 1: Install Python
  • Step 2: Install atomipy
  • Step 3: Verify Installation
• Configuration & Performance
  • Global Thresholds
  • Customizing Settings
• Getting Started for Python Beginners
  • Gromacs .n2t File Generation Example
  • Dihedrals and 1–4 pair detection
  • Import existing .itp files
  • MINFF atom typing & topology generation
  • Force Field Parameters (JSON)
  • XRD pattern simulation
  • CIF to PDB helper script
  • Common Issues for Beginners
• Function Documentation
  • File I/O
  • Bundled Molecule Library
  • Force Field
  • Molecular Topology
  • Atom Properties
  • Structure Analysis
  • Coordinate Transformations
  • Diffraction
  • Structure Building
  • Local GROMACS engine (atomipy.gromacs)
  • Cell Utilities
  • Charges
  • Oxidation states (rules + charge balance)
  • The Dummy FF — frozen model for inorganics not covered by the built-in force fields
  • Additional Functions
• Data Structure
• Common Workflow Examples
• Box Parameter Handling
  • High-Performance Structural Analysis
  • Basic Structure Processing
  • Calculating Center of Mass
  • Creating a Topology File for Molecular Dynamics
  • Creating a Supercell
  • Isomorphous Substitution in Clay Minerals
• Example Scripts
  • Bond valence analysis
  • Ionic/crystal radii and bond estimates (Shannon radii)
  • Other Examples
• Molecular Dynamics Simulations with OpenMM
  • What load_minff_into_openmm does
  • Example
• Differences from atom MATLAB library
• Citing
• Acknowledgements & External Tools
• Disclaimer
• License

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 (Holmboe, M.; Lavauzelle, N.; Swai, R. E. MINFF: A Family of Classical Force Fields for Minerals with Tunable Angle Terms. arXiv preprint 2026, doi:10.48550/arXiv.2610.18220). 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, min.ff)
│   ├── gromacs/             # Local GROMACS engine (detect gmx, .mdp, grompp/mdrun, trjconv, energy)
│   ├── structures/          # Bundled libraries of crystals and organic molecules
│   ├── topology/            # Topology interchange layer (README inside)
│   ├── __init__.py          # Package init, public API exports
│   ├── analysis.py          # High-performance analysis (RDF, CN, MSD, VACF, H-bonds, 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
│   ├── charge_formal.py     # Formal-charge assignment (ions, water)
│   ├── charge_minff.py      # MINFF charge assignment by coordination environment
│   ├── classify.py          # System-component classifier (mineral / water / ion / organic)
│   ├── composition.py       # System composition / atom classification (water/ion/organic)
│   ├── config.py            # Global settings (e.g. SPARSE_THRESHOLD)
│   ├── crystal_library.py   # Bundled inorganic crystal library (CIF)
│   ├── dcd.py               # Zero-dependency pure-Python DCD trajectory reader
│   ├── diffraction.py       # XRD pattern calculation
│   ├── distances.py         # O(N) Sparse & O(N²) Direct distance dispatchers
│   ├── distortion.py        # Phyllosilicate structural-distortion analysis
│   ├── dummy_mineral.py     # Dummy FF (frozen framework) for non-MINFF/CLAYFF inorganics
│   ├── element.py           # Chemical element assignment
│   ├── ffparams.py          # Force field parameter (JSON) loading and unit conversion
│   ├── forcefield.py        # MINFF/CLAYFF force field logic
│   ├── import_conf.py       # File importers (PDB, GRO, XYZ, CIF, etc.)
│   ├── import_top.py        # GROMACS .itp/.top readers
│   ├── import_topology.py   # Topology readers returning a `Topology` object
│   ├── lattice.py           # Parametric lattice / unit-cell builder
│   ├── mass.py              # Atomic masses
│   ├── merge_top.py         # Merge mineral + organic/water/ion topologies into one .top
│   ├── miller.py            # Miller-plane geometry
│   ├── minff_defines.py     # GROMACS defines that select a MINFF parameter set
│   ├── molecule_library.py  # Bundled organic molecule library (Chemical JSON)
│   ├── move.py              # Translate, rotate, place, center, bend
│   ├── node_cache.py        # Persist/restore workflow-node outputs (used by the web-module)
│   ├── openmm_interface.py  # OpenMM topology/system loader
│   ├── oxidation.py         # Oxidation-state guessing (Shannon/ionic)
│   ├── radius.py            # Atomic radii
│   ├── replicate.py         # Supercell replication
│   ├── resname.py           # Residue-name assignment
│   ├── shapes.py            # Region cuts (nanoparticles, nanowires, slabs)
│   ├── solvent.py           # Solvation & water molecule repairs
│   ├── surface.py           # Surface/slab builder for a Miller (hkl) face
│   ├── transform.py         # Coordinate transformations (frac/cart/wrap)
│   ├── write_conf.py        # File exporters (PDB, GRO, XYZ, etc.)
│   ├── write_top.py         # Topology exporters (ITP, PSF, LMP)
│   ├── write_topology.py    # Topology writers that consume a `Topology` object
│   └── xdrfile.py           # Optional ctypes wrapper for libxdrfile (.xtc/.trr)
├── 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 < 3,000 atoms: Use the Direct $O(N^2)$ method (faster for small systems due to NumPy vectorization).
• Systems ≥ 3,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 (the default is 3000 atoms)
ap.config.SPARSE_THRESHOLD = 10000

# Or force a specific method for a single call
atoms, cell = ap.import_pdb("Kaolinite_GII_0.0487.pdb")
atoms, bonds, angles = ap.bond_angle(atoms, Box=cell, dm_method='sparse')   # atoms, cell from e.g. ap.import_pdb(...)

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:

import atomipy as ap

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:

If you use MINFF, please cite the preprint: Holmboe, M., Lavauzelle, N., Swai, R. E. (2026). MINFF: A Family of Classical Force Fields for Minerals with Tunable Angle Terms. arXiv preprint. doi:10.48550/arXiv.2610.18220.

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):

import atomipy as ap

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=['MINFF_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*_all.json (Tailored MINFF for specific minerals, one file per angle force constant)

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, 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', 'aldehydes', 'alkanes', ...]
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=15.0, dr=0.1, atom_types=None, pair_types=None, typeA=None, typeB=None, return_cn=False): radial distribution function (pair_types=(['Si'], ['O']) or typeA/typeB select the pair; rmax is capped at half the shortest box edge); with return_cn=True also returns the running coordination number n(r).
• coordination_number(atoms, Box, cutoff=3.0, atom_types=None, neighbor_types=None): Efficient per-atom coordination analysis via sparse indexing.
• closest_atom(atoms, reference, Box=None): Find the neighbor closest to a target point or atom.
• find_H2O(atoms, Box_dim=None, rmin=1.3): 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, target_total_charge=None): Balance the total charge to an integer (or to target_total_charge) by distributing the excess or deficit evenly over the oxygen atoms

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('Pyrophyllite_GII_0.071.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(oxidation_state, 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 atomipy.gromacs.write_freeze_ndx (from atomipy import gromacs).

• 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, file_path, 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.

Additional Functions

Further public functions, grouped by task (see each function's docstring for all arguments):

Cutting and shaping

• cut_sphere(atoms, Box=None, radius=10.0, center=None, side='inside', ...) / cut_cylinder(atoms, Box=None, radius=10.0, axis='z', ...): carve a spherical nanoparticle or a nanowire (or the complementary pore).
• cut_miller(atoms, Box, h, k, l, level='auto', offset=0.0, side='below', ...) / cut_planes(atoms, Box, planes, ...): keep atoms on one side of a Miller plane, or inside the intersection of several half-spaces.
• make_slab(atoms, Box, hkl, layers=1, vacuum=0.0, gromacs_box=False, ...): build an oriented supercell/slab exposing the (hkl) face along z.
• miller_planes(h, k, l, Box=(1, 1, 1), ...), d_spacing(h, k, l, Box), hkil_to_hkl(h, k, i=None, l=0): Miller-plane geometry helpers.
• fit_box(atoms, padding=10.0, cubic=False, center=True) / condense(atoms, Box=None): size a box to the atoms' bounding box, or shrink a box and centre the atoms to remove vacuum.
• reduce_supercell(atoms, box, replicate=[1, 1, 1], merge=True, ...): reduce a replicated system back to one unit cell.

Editing atoms and atom lists

• delete_sites(atoms, atom_type=None, index=None, molid=None, x=None, y=None, z=None, logic='and', ...) (alias remove): delete atoms matching selection rules.
• translate(atoms, trans_vec, resname='all'): translate atoms by a vector.
• reorder(atoms, neworder, by=None), update(*atoms_list, molid=None, use_resname=True, force=False), join_and_reorder(*atom_lists): reorder atoms and renumber indices/molids when combining lists.
• create_grid(atom_type, density, limits, resname='ION', molid=None): fill a region with atoms on a grid at a target density.
• split_slabs_by_molid(atoms, base_resname=None): give each mineral molecule its own residue name so topologies can treat them separately.
• min_distances(atoms, Box, group_by='molid', n_pairs=10, cutoff=None): smallest distances between different molecules/residues.
• occupancy_atom(atoms, Box, rmax=1.0): occupancy of atomic sites within a radius.
• tetrahedral_rotation(atoms, Box, ...) / tetrahedral_rotation_files(files, ...): tetrahedral rotation angle of a phyllosilicate, for one structure or aggregated over many files.
• spc2tip4p(atoms, Box=None, om_dist=0.15) (alias tip3p2tip4p): convert three-site water to TIP4P by adding the M site.

Coordinates and boxes

• normalize_box(Box): normalise a 1×3, 1×6 or 1×9 box to both Box_dim and Cell.
• direct_cartesian_to_fractional(...) / direct_fractional_to_cartesian(...) / wrap_coordinates(...): direct coordinate conversions and wrapping.

Crystal library

• list_crystals(category=None), crystal_categories(), load_crystal(name_or_file, expand_symmetry=True): browse and load the bundled inorganic crystal library (CIF).
• lattice_types(): names of the preset lattices for the parametric lattice builder.
• import_mmcif(file_path, expand_symmetry=True): CIF/mmCIF import (same reader as import_cif).

Bond valence and diffraction helpers

• analyze_bvs(atoms, Box, ...) / conf2bvs(file_path, ...): one-call BVS/GII analysis of a structure or file, optionally writing a CSV.
• global_instability_index(results), add_hydrogens_bvs(atoms, Box, ...), load_shannon_radii(path=None): GII from BVS results, protonation of under-bonded oxygen, and the Shannon radii table.
• bragg_law(wavelength, value, mode='distance', n=1), atomic_scattering_factors(atom_label, wavelength, two_theta, b_iso=0.0), calculate_multiplicity(h_target, k_target, l_target, hkl_array): building blocks used by xrd.

Topology, force-field files and I/O

• write_auto(atoms, Box, file_path): choose the writer from the file extension.
• write_gmx_top(atoms, Box=None, file_path=None, ...): complete, self-contained GROMACS system .top (the top function of the write_top module; ap.write_top itself is the module).
• import_gro_coords(gro_path, atoms=None): read only coordinates from a .gro file into an existing atom list.
• import_gaff_top(top_path, itp_path=None): import an ACPYPE-generated GROMACS topology.
• merge_top_files(*top_gro_pairs, out_top, out_gro, defines=None, ...) / write_merged_top(...): merge several .top/.gro pairs into one system.
• get_mol_sequence(atoms) / get_mol_sequence_typed(atoms): the ordered GROMACS [ molecules ] sequence of an atom list.
• get_ffparams_dir(), list_ff_blocks(json_path): location of the bundled ffparams and the blocks available in a force-field JSON.
• have_xdrfile(), save_node_state(...), load_node_state(...): optional libxdrfile check, and the node cache used by the web-module.

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:

import atomipy as ap

atoms, _ = ap.import_pdb("Kaolinite_GII_0.0487.pdb")
atoms = ap.element(atoms)

# 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")
atoms = ap.element(atoms)   # the analysis functions select atoms by element/type

# 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) and build a supercell
atoms, Cell = ap.import_pdb("Pyrophyllite_GII_0.071.pdb")
atoms, Box, _ = ap.replicate_system(atoms, Cell, replicate=[6, 4, 1])

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

# Octahedral substitution: replace 6 Al atoms with Mg
# (atom types are matched first; if a type is not found it falls back to matching by element)
atoms, Box, _ = ap.substitute(
    atoms,
    Box,
    num_oct_subst=6,
    o1_type='Alo',   # octahedral Al
    o2_type='Mgo',   # becomes Mg in the octahedral position
    min_o2o2_dist=5.2  # minimum distance between Mg atoms in Angstroms
)

# Octahedral and tetrahedral substitutions in one call:
# replace 2 more Al with Mg (octahedral) and 4 Si with Al (tetrahedral)
atoms, Box, _ = ap.substitute(
    atoms,
    Box,
    num_oct_subst=2,
    o1_type='Alo',
    o2_type='Mgo',
    min_o2o2_dist=5.5,
    num_tet_subst=4,
    t1_type='Si',
    t2_type='Alt',   # becomes Al in the tetrahedral position
    min_t2t2_dist=5.5
)

# Limit substitutions to a region along one axis (here only the central octahedral sheet, z in [-2, 2] Å)
atoms, Box, _ = ap.substitute(
    atoms,
    Box,
    num_oct_subst=2,
    o1_type='Alo',
    o2_type='Mgo',
    min_o2o2_dist=5.2,
    lo_limit=-2.0,
    hi_limit=2.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, file_path="substituted_clay.pdb")

substitute prints whether it could place all requested substitutions (✓ ... successful! or ⚠ ... not optimal!); if the minimum-distance constraint cannot be met it places as many as it can, so check the printed summary.

Example Scripts

The package includes example scripts in scripts/ that demonstrate complete workflows for processing mineral structures with atomipy.

scripts/run_minff_atomi.py and scripts/run_clayff_atomi.py

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

What the scripts do:

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

Running the scripts:

Run python scripts/run_minff_atomi.py or python scripts/run_clayff_atomi.py from the repository root, where the input structures live (the first expects Kaolinite_GII_0.0487.pdb, the second Kaolinite_GII_0.0487.gro in the working directory).

Output files:

• replicated_structure.gro / replicated_structure.pdb - The enlarged supercell
• bonded_structure.gro - The supercell after the bond search
• minff_structure_stats.log (or clayff_structure_stats.log) - Structure statistics: dimensions, density, bonds, angles
• minff_topology.itp (or clayff_topology.itp) - GROMACS topology file with bond and angle definitions
• minff_topology.psf / minff_topology.data (or clayff_topology.psf / clayff_topology.data) - NAMD and LAMMPS topology files
• minff_atomtypes.n2t - GROMACS atom name to type file for gmx x2top (MINFF script only)
• preem.gro (and preem.pdb for the MINFF script) - Final structure with the assigned atom types and charges

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 writes self-contained GROMACS .top/.gro files for MINFF/CLAYFF systems and can load them straight into OpenMM with ap.load_minff_into_openmm, which returns the OpenMM Topology, System and positions. (The OpenMM package itself is optional and is only imported when this function is called.) You then run the simulation with ordinary OpenMM code.

What load_minff_into_openmm does

• Reads the atomipy-written .top + .gro pair with openmm.app.GromacsTopFile/GromacsGroFile.
• Activates the force-field variant you choose through defines (the same names as the GROMACS define = -D... line, e.g. MINFF_k500, OPC3).
• Uses MINFF defaults: flexible bonds/angles on the mineral body (constraints=None), PME electrostatics, a 1.0 nm cutoff, Lorentz–Berthelot mixing and the long-range LJ correction. Pass rigid_water=True for rigid water models.
• include_dir must point to the folder holding the MINFF .itp files (ap.get_ffparams_dir()).

Example

import atomipy as ap
from openmm import LangevinIntegrator, Platform
from openmm.app import Simulation, PDBReporter, StateDataReporter
from openmm.unit import kelvin, picosecond, femtoseconds

# 1. Load a mineral system and assign MINFF atom types/charges
atoms, Box = ap.import_gro("system.gro")
atoms = ap.minff(atoms, Box=Box)

# 2. Write a complete, self-contained GROMACS topology and the matching coordinates
ap.write_gmx_top(atoms, Box=Box, file_path="system_minff.top")
ap.write_gro(atoms, Box=Box, file_path="system_minff.gro")

# 3. Build the OpenMM Topology / System / positions
topology, system, positions = ap.load_minff_into_openmm(
    top_path="system_minff.top",
    gro_path="system_minff.gro",
    defines=["MINFF_k500", "OPC3_IOD_LM", "OPC3"],   # parameter set (angle k=500) and water model
    include_dir=ap.get_ffparams_dir(),
    nonbonded_cutoff_nm=0.9,    # must be less than half the shortest box edge
    constraints=None,           # MINFF needs flexible bonds and angles
    rigid_water=True,
)

# 4. Run the dynamics with plain OpenMM
integrator = LangevinIntegrator(298.15 * kelvin, 1 / picosecond, 1.0 * femtoseconds)
simulation = Simulation(topology, system, integrator)   # OpenMM picks the fastest available platform
simulation.context.setPositions(positions)
simulation.minimizeEnergy()
simulation.reporters.append(PDBReporter("trajectory.pdb", 1000))
simulation.reporters.append(StateDataReporter("log.csv", 1000, step=True, potentialEnergy=True, temperature=True))
simulation.step(5000)

scripts/run_openmm_minff.py shows the loading step (build, replicate, type, write, load into OpenMM and evaluate the initial energy). To run simulations without writing OpenMM code, use the visual builder at www.atomipy.io (OpenMM or a local GROMACS engine), or the atomipy.gromacs module described above for GROMACS.

---

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.

Citing

If you use atomipy, please cite it using CITATION.cff. If you use the MINFF force field (atom typing, charges, or the bundled parameters), please also cite:

Holmboe, M., Lavauzelle, N., Swai, R. E. (2026). MINFF: A Family of Classical Force Fields for Minerals with Tunable Angle Terms. arXiv preprint. doi:10.48550/arXiv.2610.18220

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, required); 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

The atomipy source code is released under the MIT License (© 2024–2026 Michael Holmboe).

Bundled third-party data keeps its own license and attribution — the MIT license does not re-license it. This covers molecular fragments from Avogadro2 under BSD 3-Clause, crystallographic structures from open databases such as COD/AMCSD, MINFF/CLAYFF force-field parameters from the published literature, the standard data and water model files taken unmodified from the GROMACS distribution under LGPL-2.1-or-later in atomipy/ffparams/min.ff/, and the tabulated reference data in atomipy/data/. See the NOTICE file for details.