Metadata-Version: 2.4
Name: atomipy
Version: 0.98
Summary: A Python toolbox for molecular structure analysis and simulation with support for both orthogonal and triclinic periodic boundary conditions
Home-page: https://github.com/mholmboe/atomipy
Author: Michael Holmboe
Author-email: michael.holmboe@umu.se
License: MIT
Keywords: molecular dynamics,periodic boundary conditions,minerals,chemistry,physics,triclinic,distance calculations,bond detection
Classifier: Programming Language :: Python :: 3
Classifier: License :: OSI Approved :: MIT License
Classifier: Operating System :: OS Independent
Classifier: Topic :: Scientific/Engineering :: Chemistry
Classifier: Topic :: Scientific/Engineering :: Physics
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Classifier: Topic :: Scientific/Engineering :: Visualization
Requires-Python: >=3.7
Description-Content-Type: text/markdown
License-File: LICENSE
License-File: NOTICE
Requires-Dist: numpy>=1.18.0
Requires-Dist: tqdm>=4.45.0
Requires-Dist: numba>=0.50.0; python_version >= "3.6"
Requires-Dist: gemmi>=0.7.0
Provides-Extra: xrd
Requires-Dist: matplotlib>=3.3.0; extra == "xrd"
Requires-Dist: scipy>=1.5.0; extra == "xrd"
Provides-Extra: cif
Requires-Dist: gemmi>=0.7.0; extra == "cif"
Dynamic: author
Dynamic: author-email
Dynamic: classifier
Dynamic: description
Dynamic: description-content-type
Dynamic: home-page
Dynamic: keywords
Dynamic: license
Dynamic: license-file
Dynamic: provides-extra
Dynamic: requires-dist
Dynamic: requires-python
Dynamic: summary

# 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**](https://github.com/mholmboe/atom) library and can in particular be used to generate molecular topology files for the [**MINFF**](https://github.com/mholmboe/minff) forcefield with a streamlined Python interface. An online molecular system and topology builder running [**www.atomipy.io**](https://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**](https://github.com/mholmboe/atom/tree/master/ATOM_scripts_lecture/MINFF).

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](#overview)
- [Two ways to use atomipy](#two-ways-to-use-atomipy)
- [Project Structure](#project-structure)
- [Common Variables](#common-variables)
  - [Structure Containers](#structure-containers)
  - [Atom dictionary fields](#atom-dictionary-fields)
- [Key Features](#key-features)
- [Requirements](#requirements)
- [Installation](#installation)
  - [Step 1: Install Python](#step-1-install-python)
  - [Step 2: Install atomipy](#step-2-install-atomipy)
  - [Step 3: Verify Installation](#step-3-verify-installation)
- [Configuration & Performance](#configuration--performance)
  - [Global Thresholds](#global-thresholds)
  - [Customizing Settings](#customizing-settings)
- [Getting Started for Python Beginners](#getting-started-for-python-beginners)
  - [Gromacs .n2t File Generation Example](#gromacs-n2t-file-generation-example)
  - [Dihedrals and 1–4 pair detection](#dihedrals-and-14-pair-detection)
  - [Import existing .itp files](#import-existing-itp-files)
  - [MINFF atom typing & topology generation](#minff-atom-typing--topology-generation)
  - [Force Field Parameters (JSON)](#force-field-parameters-json)
  - [XRD pattern simulation](#xrd-pattern-simulation)
  - [CIF to PDB helper script](#cif-to-pdb-helper-script)
  - [Common Issues for Beginners](#common-issues-for-beginners)
- [Function Documentation](#function-documentation)
  - [File I/O](#file-io)
  - [Bundled Molecule Library](#bundled-molecule-library)
  - [Force Field](#force-field)
  - [Molecular Topology](#molecular-topology)
  - [Atom Properties](#atom-properties)
  - [Structure Analysis](#structure-analysis)
  - [Coordinate Transformations](#coordinate-transformations)
  - [Diffraction](#diffraction)
  - [Structure Building](#structure-building)
  - [Local GROMACS engine (`atomipy.gromacs`)](#local-gromacs-engine-atomipygromacs)
  - [Cell Utilities](#cell-utilities)
  - [Charges](#charges)
  - [Oxidation states (rules + charge balance)](#oxidation-states-rules--charge-balance)
  - [The Dummy FF — frozen model for inorganics not covered by the built-in force fields](#the-dummy-ff--frozen-model-for-inorganics-not-covered-by-the-built-in-force-fields)
  - [Additional Functions](#additional-functions)
- [Data Structure](#data-structure)
- [Common Workflow Examples](#common-workflow-examples)
- [Box Parameter Handling](#box-parameter-handling)
  - [High-Performance Structural Analysis](#high-performance-structural-analysis)
  - [Basic Structure Processing](#basic-structure-processing)
  - [Calculating Center of Mass](#calculating-center-of-mass)
  - [Creating a Topology File for Molecular Dynamics](#creating-a-topology-file-for-molecular-dynamics)
  - [Creating a Supercell](#creating-a-supercell)
  - [Isomorphous Substitution in Clay Minerals](#isomorphous-substitution-in-clay-minerals)
- [Example Scripts](#example-scripts)
  - [Bond valence analysis](#bond-valence-analysis)
  - [Ionic/crystal radii and bond estimates (Shannon radii)](#ioniccrystal-radii-and-bond-estimates-shannon-radii)
  - [Other Examples](#other-examples)
- [Molecular Dynamics Simulations with OpenMM](#molecular-dynamics-simulations-with-openmm)
  - [What `load_minff_into_openmm` does](#what-load_minff_into_openmm-does)
  - [Example](#example)
- [Differences from atom MATLAB library](#differences-from-atom-matlab-library)
- [Citing](#citing)
- [Acknowledgements & External Tools](#acknowledgements--external-tools)
- [Disclaimer](#disclaimer)
- [License](#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**](https://github.com/mholmboe/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](https://doi.org/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](https://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](https://www.atomipy.io)** (open source:
[atomipy-web-module](https://github.com/mholmboe/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

```text
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](https://www.python.org/downloads/) (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](https://git-scm.com/downloads)
3. Clone the repository and install:
   ```bash
   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:
   ```bash
   cd path/to/extracted/atomipy
   ```
5. Install the package and its dependencies:
   ```bash
   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:

```bash
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:

```python
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:

```python
# 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:
```bash
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:

```python
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:

```python
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:

```python
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:

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

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

```bash
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**](https://github.com/mholmboe/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](https://doi.org/10.48550/arXiv.2610.18220).

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

```python
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`:

```python
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:
```bash
python scripts/run_xrd_example.py Kaolinite_GII_0.0487.pdb --two-theta 5 70 --save-output
```

Minimal API example:
```python
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:
```bash
# 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](https://github.com/OpenChemistry/molecules)
(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)`.

```python
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']`.

```python
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 `#include`s 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 `#include`d), 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:

```python
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.

```python
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

```python
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

```python
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

```python
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

```python
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

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

```python
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:

```python
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

```python
# 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

```python
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](https://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`](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](https://doi.org/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](https://numpy.org/), [Numba](https://numba.pydata.org/) (optional), [tqdm](https://tqdm.github.io/).
- **Structure & trajectory I/O** — [GEMMI](https://gemmi.readthedocs.io/) (CIF/mmCIF, required); GROMACS `libxdrfile` (`.xtc`/`.trr`) and [MDTraj](https://www.mdtraj.org/) (optional binary trajectories). The bundled inorganic crystal library derives from the [Avogadro](https://avogadro.cc/) collection.
- **Simulation** — [OpenMM](https://openmm.org/) and [GROMACS](https://www.gromacs.org/).
- **Organic force fields** — [OpenFF Toolkit / Interchange](https://openforcefield.org/) (Sage, Parsley), [openmmforcefields](https://github.com/openmm/openmmforcefields) + [ParmEd](https://parmed.github.io/ParmEd/) and/or [ACPYPE](https://github.com/alanwilter/acpype) + AmberTools *antechamber* (GAFF), with [RDKit](https://www.rdkit.org/) / [Open Babel](https://openbabel.org/).
- **Force fields** — [MINFF](https://github.com/mholmboe/minff); CLAYFF (Cygan, Liang & Kalinichev, *J. Phys. Chem. B* **2004**, *108*, 1255); GAFF; OpenFF (Sage/Parsley).
- Ported from the MATLAB [atom](https://github.com/mholmboe/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](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](https://avogadro.cc/) 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](https://www.gromacs.org/) distribution under LGPL-2.1-or-later in `atomipy/ffparams/min.ff/`, and the tabulated reference data in `atomipy/data/`. See the [NOTICE](NOTICE) file for details.
