Use Otter pair potentials in LAMMPS

tools/otter_lammps_md.py provides a reproducible path from an Otter effective ion–ion potential to a classical LAMMPS simulation. It supports one or many ionic species and performs the complete sequence

\[V_{ab}(r) \longrightarrow \text{LAMMPS table and input} \longrightarrow \text{NVT equilibration} \longrightarrow \text{NVE sampling} \longrightarrow g_{ab}(r),\ S_{ab}(k).\]

LAMMPS is an external optional program and is not installed by Otter. The default executable names are lmp and mpirun; both can be changed in the Python configuration. The driver has no command-line parameter parser: all physical and numerical settings remain visible in the calling script.

Which Otter potential is used?

The input is the QOZ effective ion–ion pair-potential matrix vij_r:

\[V_{ab}(r) = \texttt{ion["vij_r"][a,b,:]},\]

where \(a,b\) are ionic-species indices. This is not the electron–ion channel v_ie_r. Otter stores vij_r in Hartree on a radius grid in Bohr. The helper pair_potentials_from_otter accepts any of these sources:

  • a complete workflow result;

  • its result["ion"] dictionary;

  • a portable NPZ state produced by otter.save_plasma_state();

  • a prepared multicomponent QOZ object, before an HNC solve.

For \(N_s\) species, the helper verifies that vij_r has shape (N_s, N_s, N_r), is symmetric, and contains every unique pair. The pair order is (0,0), (0,1), ..., (N_s-1,N_s-1).

Complete mixture example from an Otter NPZ

The following example assumes that outputs/ch2_state.npz is a converged Otter state with species order C, H. Edit the settings at the beginning of the script and run the Python file from the repository root.

from pathlib import Path

import numpy as np

from tools.otter_lammps_md import (
    MDConfig,
    MDSpecies,
    pair_potentials_from_otter,
    run_otter_lammps_md,
)


# ------------------------- user settings -------------------------
STATE_PATH = Path("outputs/ch2_state.npz")
OUTPUT_DIR = Path("outputs/ch2_lammps_md")
TI_EV = 10.0

# Exact CH2 stoichiometry: 1024 C + 2048 H.
PARTICLE_COUNTS = (1024, 2048)
MPI_PROCESSES = 8

# Demonstration values, not universal defaults.  Validate them below.
TIMESTEP_PS = 1.0e-5
THERMOSTAT_DAMP_PS = 1.0e-3
NVT_STEPS = 10_000
NVE_STEPS = 100_000
# -----------------------------------------------------------------


with np.load(STATE_PATH, allow_pickle=False) as state:
    symbols = tuple(str(value) for value in state["species_symbols"])
    if symbols != ("C", "H"):
        raise ValueError(f"Expected C,H in the NPZ, found {symbols}.")

    # zbar is the charge used to construct the saved QOZ potential.
    zbar = np.asarray(state["zbar"], dtype=float)
    total_ion_density = float(np.sum(state["n_i_bohr3"]))
    potentials = pair_potentials_from_otter(state)


config = MDConfig(
    output_dir=OUTPUT_DIR,
    species=(
        MDSpecies("C", 12.011, PARTICLE_COUNTS[0], charge_e=float(zbar[0])),
        MDSpecies("H", 1.008, PARTICLE_COUNTS[1], charge_e=float(zbar[1])),
    ),
    ion_density_bohr3=total_ion_density,
    ion_temperature_ev=TI_EV,
    timestep_ps=TIMESTEP_PS,
    thermostat_damp_ps=THERMOSTAT_DAMP_PS,
    equilibration_steps=NVT_STEPS,
    production_steps=NVE_STEPS,
    rdf_bins=500,
    rdf_every=100,
    rdf_repeat=50,
    trajectory_every=5_000,
    table_points=8_192,
    r_min_bohr=0.1,
    k_bin_width_angstrom_inv=0.1,
    k_max_angstrom_inv=8.3,
    structure_factor_workers=8,
    mpi_processes=MPI_PROCESSES,
    lammps_executable="lmp",
    mpi_launcher="mpirun",
    reuse_completed=True,
)

result = run_otter_lammps_md(config, potentials)

print("pairs:", result["md_pair_labels"])
print("NVE relative energy drift:", result["md_nve_relative_energy_drift"])
print("saved:", OUTPUT_DIR / "md_results.npz")

Why charge_e=zbar is specified

With nonzero charge_e, the driver decomposes each Otter potential as

\[V_{ab}(r) = \frac{\bar Z_a\bar Z_b}{r} + \Delta V_{ab}(r).\]

LAMMPS evaluates the analytic Coulomb term, while its table contains only the finite remainder \(\Delta V_{ab}\). Their sum is the original Otter potential. This avoids interpolating a \(1/r\) singularity. The combined energy and force are shifted smoothly to zero at the finite MD cutoff.

The charge must match the zbar used by QOZ. Do not substitute zbar_aa_ws or zstar unless the pair potential was explicitly rebuilt with that alternative charge definition. Setting every charge_e to zero is supported, but then the complete potential, including its short-range core, is tabulated directly.

The cutoff is the smaller of 0.48 * box_length and the largest available Otter radius. The NPZ export window must therefore retain enough of vij_r for the intended box. Increase state_r_max_bohr when generating the Otter state if necessary.

In-memory and single-species inputs

No intermediate NPZ is required. A finished workflow result can be passed directly:

workflow = solve_plasma_workflow(otter_config)
potentials = pair_potentials_from_otter(workflow)

For one species, use one MDSpecies and its total ion density:

ion = workflow["ion"]
config = MDConfig(
    output_dir=Path("outputs/al_md"),
    species=(
        MDSpecies(
            "Al",
            26.9815385,
            2048,
            charge_e=float(ion["zbar"]),
        ),
    ),
    ion_density_bohr3=float(ion["n_i"]),
    ion_temperature_ev=15.0,
    timestep_ps=1.0e-5,
    thermostat_damp_ps=1.0e-3,
)
result = run_otter_lammps_md(config, potentials)

For a mixture study that should not spend time solving HNC first, extract the potential from the reusable QOZ preparation:

prepared = prepare_multicomponent_ion_structure_from_electronic_result(
    config,
    electronic_kind=electronic_kind,
    electronic_result=electronic_result,
)
potentials = pair_potentials_from_otter(prepared)

The masses, integer MD particle counts, ion temperature, and LAMMPS controls still belong in MDConfig. Particle counts must reproduce the intended number fractions; they do not alter the density supplied separately through ion_density_bohr3.

Generated LAMMPS files

Before execution, the driver writes:

LAMMPS input artifacts

File

Contents

atoms.data

Reproducible randomized species placement on an FCC starting lattice, exact particle counts, masses, charges, and the cubic periodic box.

pair_potentials.table

One table section for every unique pair, including energy and force.

in.otter_md

Complete metal-unit LAMMPS input: velocities, NVT equilibration, NVE production, block RDF accumulation, and trajectory output.

After a successful run, the directory also contains log.lammps, screen.log, rdf_blocks.dat, trajectory.lammpstrj, md_results.npz, and run_metadata.json. The metadata records the configuration, LAMMPS version, box and cutoff, wall time, and hashes of the input and output artifacts. With reuse_completed=True, results are reused only when the previous run completed normally and the three input hashes still match.

How \(g_{ab}(r)\) and \(S_{ab}(k)\) are obtained

LAMMPS accumulates each partial RDF only during the NVE production segment. md_gij_r is the mean of the saved RDF blocks and md_gij_block_sem is the standard error across those blocks. The columns are identified by md_pair_labels.

The partial structure factors are not obtained by Fourier transforming the finite-range RDF. For each saved NVE frame, the driver evaluates the periodic density modes

\[\rho_a(\mathbf{k}) = \sum_{j\in a} \exp(i\mathbf{k}\cdot\mathbf{r}_j), \qquad \mathbf{k}=\frac{2\pi}{L}\mathbf{n},\]

and the Ashcroft–Langreth partial estimator

\[S_{ab}(k) = \left\langle \frac{\operatorname{Re}[\rho_a(\mathbf{k})\rho_b^*(\mathbf{k})]} {\sqrt{N_aN_b}} \right\rangle_{\text{vectors in bin, frames}}.\]

The relevant output fields are:

labels = result["md_pair_labels"]

r = result["md_r_bohr"]
gab = result["md_gij_r"]
gab_sem = result["md_gij_block_sem"]

k = result["md_k_bohr_inv"]
Sab = result["md_sij_k"]
Sab_sem = result["md_sij_frame_sem"]
vectors_per_bin = result["md_vectors_per_k_bin"]

The reported \(S_{ab}\) SEM is the standard deviation of the per-frame shell averages divided by the square root of the saved-frame count. It does not correct for temporal autocorrelation, so it is a sampling diagnostic, not automatically a rigorous confidence interval. Small-\(k\) bins contain few periodic vectors and can have much larger finite-cell uncertainty; inspect md_vectors_per_k_bin before plotting them.

Numerical checks before accepting a run

At minimum, verify all of the following:

  • the NVT segment is discarded and observables use only NVE production;

  • md_nve_relative_energy_drift is acceptably small and decreases when the timestep is reduced;

  • the mean NVE temperature is compatible with the target state;

  • increasing NVT duration does not change the production averages;

  • increasing NVE duration reduces sampling noise without shifting the mean;

  • increasing particle number and the box length does not change the resolved peaks or the retained low-\(k\) region;

  • the MD cutoff lies in a region where the shifted potential is physically negligible;

  • every mixture pair appears exactly once in md_pair_labels.

These checks matter more than adopting one fixed timestep or step count. Hot light ions generally require a smaller timestep, while low-\(k\) structure generally requires a larger box and longer sampling. The driver automates data conversion and auditing, but it cannot choose a universally converged MD protocol for every material and thermodynamic state.

LAMMPS is described by Thompson et al. [2022].