Potential

The potential module defines potential energy functions and functional forms.

Base

Base classes for potential functions.

Potential

Base class for all potential functions in MolPy.

This class provides a template for defining potential functions that can be used in molecular simulations. Each potential class should implement calc_energy and calc_forces methods with specific parameters.

calc_energy

calc_energy(*args, **kwargs)

Calculate the potential energy.

Parameters

args: Arguments specific to the potential type *kwargs: Keyword arguments specific to the potential type

Returns

float The potential energy.

Source code in src/molpy/potential/base.py

def calc_energy(self, *args, **kwargs) -> float:
    """
    Calculate the potential energy.

    Parameters
    ----------
    *args: Arguments specific to the potential type
    **kwargs: Keyword arguments specific to the potential type

    Returns
    -------
    float
        The potential energy.
    """
    raise NotImplementedError("Subclasses must implement this method.")

calc_forces

calc_forces(*args, **kwargs)

Calculate the forces.

Parameters

args: Arguments specific to the potential type *kwargs: Keyword arguments specific to the potential type

Returns

np.ndarray An array of forces.

Source code in src/molpy/potential/base.py

def calc_forces(self, *args, **kwargs) -> np.ndarray:
    """
    Calculate the forces.

    Parameters
    ----------
    *args: Arguments specific to the potential type
    **kwargs: Keyword arguments specific to the potential type

    Returns
    -------
    np.ndarray
        An array of forces.
    """
    raise NotImplementedError("Subclasses must implement this method.")

Potentials

Bases: UserList[Potential]

Collection of potential functions.

This class provides a way to combine multiple potentials and calculate total energy and forces. However, since different potentials require different parameters, you need to call calc_energy and calc_forces for each potential separately and sum the results.

For a simpler interface, use the helper functions in potential.utils to extract data from Frame objects.

calc_energy

calc_energy(*args, **kwargs)

Calculate the total energy by summing energies from all potentials.

If a Frame object is passed as the first argument, automatically extracts the necessary data for each potential type.

Parameters

args: Arguments passed to each potential's calc_energy method. If first arg is a Frame, data is automatically extracted. *kwargs: Keyword arguments passed to each potential's calc_energy method

Returns

float The total energy.

Source code in src/molpy/potential/base.py

def calc_energy(self, *args, **kwargs) -> float:
    """
    Calculate the total energy by summing energies from all potentials.

    If a Frame object is passed as the first argument, automatically extracts
    the necessary data for each potential type.

    Parameters
    ----------
    *args: Arguments passed to each potential's calc_energy method.
           If first arg is a Frame, data is automatically extracted.
    **kwargs: Keyword arguments passed to each potential's calc_energy method

    Returns
    -------
    float
        The total energy.
    """
    # Check if first argument is a Frame
    if args and isinstance(args[0], Frame):
        from .utils import calc_energy_from_frame

        return sum(calc_energy_from_frame(pot, args[0]) for pot in self)

    # Otherwise, pass arguments directly
    return sum(pot.calc_energy(*args, **kwargs) for pot in self)

calc_forces

calc_forces(*args, **kwargs)

Calculate the total forces by summing forces from all potentials.

If a Frame object is passed as the first argument, automatically extracts the necessary data for each potential type.

Parameters

args: Arguments passed to each potential's calc_forces method. If first arg is a Frame, data is automatically extracted. *kwargs: Keyword arguments passed to each potential's calc_forces method

Returns

np.ndarray An array of total forces.

Source code in src/molpy/potential/base.py

def calc_forces(self, *args, **kwargs) -> np.ndarray:
    """
    Calculate the total forces by summing forces from all potentials.

    If a Frame object is passed as the first argument, automatically extracts
    the necessary data for each potential type.

    Parameters
    ----------
    *args: Arguments passed to each potential's calc_forces method.
           If first arg is a Frame, data is automatically extracted.
    **kwargs: Keyword arguments passed to each potential's calc_forces method

    Returns
    -------
    np.ndarray
        An array of total forces.
    """
    if not self:
        raise ValueError(
            "Cannot determine force shape: no potentials in collection"
        )

    # Check if first argument is a Frame
    if args and isinstance(args[0], Frame):
        from .utils import calc_forces_from_frame

        # Get forces from first potential to determine shape
        first_forces = calc_forces_from_frame(self[0], args[0])
        total_forces = np.zeros_like(first_forces)

        for pot in self:
            forces = calc_forces_from_frame(pot, args[0])
            total_forces += forces

        return total_forces

    # Otherwise, pass arguments directly
    # Get forces from first potential to determine shape
    first_forces = self[0].calc_forces(*args, **kwargs)
    total_forces = np.zeros_like(first_forces)

    for pot in self:
        forces = pot.calc_forces(*args, **kwargs)
        total_forces += forces

    return total_forces

Angle

Bond

Dihedral

Dihedral potentials.

DihedralOPLSStyle

DihedralOPLSStyle()

Bases: DihedralStyle

OPLS dihedral style with fixed name='opls'.

OPLS dihedral uses Ryckaert-Bellemans (RB) coefficients c0-c5. LAMMPS opls style expects k1-k4, which are computed from c0-c5 using analytical conversion according to GROMACS manual Eqs. 200-201.

Source code in src/molpy/potential/dihedral/opls.py

def __init__(self):
    super().__init__("opls")

def_type

def_type(itom, jtom, ktom, ltom, c0=0.0, c1=0.0, c2=0.0, c3=0.0, c4=0.0, c5=0.0, name='')

Define OPLS dihedral type.

Parameters:

Name	Type	Description	Default
itom	AtomType	First atom type	required
jtom	AtomType	Second atom type	required
ktom	AtomType	Third atom type	required
ltom	AtomType	Fourth atom type	required
c0-c5		OPLS Ryckaert-Bellemans coefficients	required
name	str	Optional name (defaults to itom-jtom-ktom-ltom)	''

Returns:

Type	Description
DihedralOPLSType	DihedralOPLSType instance

Source code in src/molpy/potential/dihedral/opls.py

def def_type(
    self,
    itom: AtomType,
    jtom: AtomType,
    ktom: AtomType,
    ltom: AtomType,
    c0: float = 0.0,
    c1: float = 0.0,
    c2: float = 0.0,
    c3: float = 0.0,
    c4: float = 0.0,
    c5: float = 0.0,
    name: str = "",
) -> DihedralOPLSType:
    """Define OPLS dihedral type.

    Args:
        itom: First atom type
        jtom: Second atom type
        ktom: Third atom type
        ltom: Fourth atom type
        c0-c5: OPLS Ryckaert-Bellemans coefficients
        name: Optional name (defaults to itom-jtom-ktom-ltom)

    Returns:
        DihedralOPLSType instance
    """
    if not name:
        name = f"{itom.name}-{jtom.name}-{ktom.name}-{ltom.name}"
    dt = DihedralOPLSType(name, itom, jtom, ktom, ltom, c0, c1, c2, c3, c4, c5)
    self.types.add(dt)
    return dt

to_lammps_params

to_lammps_params(dihedral_type)

Convert OPLS c0-c5 coefficients to LAMMPS k1-k4 format.

Note: In OPLS XML files, c1-c4 typically contain OPLS coefficients F1-F4 directly (not RB format). So we use c1-c4 directly as k1-k4.

For true RB format coefficients, use rb_to_opls() function instead.

Parameters:

Name	Type	Description	Default
dihedral_type	DihedralOPLSType	DihedralOPLSType with c0-c5 parameters	required

Returns:

Type	Description
list[float]	List of [k1, k2, k3, k4] for LAMMPS in kcal/mol

Source code in src/molpy/potential/dihedral/opls.py

def to_lammps_params(self, dihedral_type: DihedralOPLSType) -> list[float]:
    """Convert OPLS c0-c5 coefficients to LAMMPS k1-k4 format.

    Note: In OPLS XML files, c1-c4 typically contain OPLS coefficients F1-F4
    directly (not RB format). So we use c1-c4 directly as k1-k4.

    For true RB format coefficients, use rb_to_opls() function instead.

    Args:
        dihedral_type: DihedralOPLSType with c0-c5 parameters

    Returns:
        List of [k1, k2, k3, k4] for LAMMPS in kcal/mol
    """
    # OPLS XML stores F1-F4 in c1-c4 (already OPLS format)
    # Just return them directly (should already be in kcal/mol)
    return [
        dihedral_type.params.kwargs.get("c1", 0.0),
        dihedral_type.params.kwargs.get("c2", 0.0),
        dihedral_type.params.kwargs.get("c3", 0.0),
        dihedral_type.params.kwargs.get("c4", 0.0),
    ]

DihedralOPLSType

DihedralOPLSType(name, itom, jtom, ktom, ltom, c0=0.0, c1=0.0, c2=0.0, c3=0.0, c4=0.0, c5=0.0)

Bases: DihedralType

OPLS dihedral type with c0-c5 coefficients.

Parameters:

Name	Type	Description	Default
name	str	Type name	required
itom	AtomType	First atom type	required
jtom	AtomType	Second atom type	required
ktom	AtomType	Third atom type	required
ltom	AtomType	Fourth atom type	required
c0-c5		OPLS Ryckaert-Bellemans coefficients	required

Source code in src/molpy/potential/dihedral/opls.py

def __init__(
    self,
    name: str,
    itom: AtomType,
    jtom: AtomType,
    ktom: AtomType,
    ltom: AtomType,
    c0: float = 0.0,
    c1: float = 0.0,
    c2: float = 0.0,
    c3: float = 0.0,
    c4: float = 0.0,
    c5: float = 0.0,
):
    """
    Args:
        name: Type name
        itom: First atom type
        jtom: Second atom type
        ktom: Third atom type
        ltom: Fourth atom type
        c0-c5: OPLS Ryckaert-Bellemans coefficients
    """
    super().__init__(
        name,
        itom,
        jtom,
        ktom,
        ltom,
        c0=c0,
        c1=c1,
        c2=c2,
        c3=c3,
        c4=c4,
        c5=c5,
    )

Pair

CoulCut

Bases: PairPotential

Coulomb pair potential with cutoff.

The potential is defined as: V(r) = q_i * q_j / r

The force is: F(r) = q_i * q_j / r^3 * dr

calc_energy

calc_energy(r, pair_idx, charges)

Calculate Coulomb energy.

Parameters:

Name	Type	Description	Default
r	NDArray[floating]	Atom coordinates (shape: (n_atoms, 3))	required
pair_idx	NDArray[integer]	Pair indices (shape: (n_pairs, 2))	required
charges	NDArray[floating]	Atom charges (shape: (n_atoms,))	required

Returns:

Type	Description
float	Total Coulomb energy

Source code in src/molpy/potential/pair/coul.py

def calc_energy(
    self,
    r: NDArray[np.floating],
    pair_idx: NDArray[np.integer],
    charges: NDArray[np.floating],
) -> float:
    """
    Calculate Coulomb energy.

    Args:
        r: Atom coordinates (shape: (n_atoms, 3))
        pair_idx: Pair indices (shape: (n_pairs, 2))
        charges: Atom charges (shape: (n_atoms,))

    Returns:
        Total Coulomb energy
    """
    # Calculate distances
    dr = r[pair_idx[:, 1]] - r[pair_idx[:, 0]]
    dr_norm = np.linalg.norm(dr, axis=1)

    # Calculate energy
    energy = charges[pair_idx[:, 0]] * charges[pair_idx[:, 1]] / dr_norm

    return float(np.sum(energy))

calc_forces

calc_forces(r, pair_idx, charges)

Calculate Coulomb forces.

Parameters:

Name	Type	Description	Default
r	NDArray[floating]	Atom coordinates (shape: (n_atoms, 3))	required
pair_idx	NDArray[integer]	Pair indices (shape: (n_pairs, 2))	required
charges	NDArray[floating]	Atom charges (shape: (n_atoms,))	required

Returns:

Type	Description
NDArray[floating]	Array of forces on each atom (shape: (n_atoms, 3))

Source code in src/molpy/potential/pair/coul.py

def calc_forces(
    self,
    r: NDArray[np.floating],
    pair_idx: NDArray[np.integer],
    charges: NDArray[np.floating],
) -> NDArray[np.floating]:
    """
    Calculate Coulomb forces.

    Args:
        r: Atom coordinates (shape: (n_atoms, 3))
        pair_idx: Pair indices (shape: (n_pairs, 2))
        charges: Atom charges (shape: (n_atoms,))

    Returns:
        Array of forces on each atom (shape: (n_atoms, 3))
    """
    n_atoms = len(r)

    # Calculate distances
    dr = r[pair_idx[:, 1]] - r[pair_idx[:, 0]]
    dr_norm = np.linalg.norm(dr, axis=1, keepdims=True)

    # Calculate forces
    forces = charges[pair_idx[:, 0]] * charges[pair_idx[:, 1]] / dr_norm**3 * dr

    # Accumulate forces on atoms
    per_atom_forces = np.zeros((n_atoms, 3), dtype=np.float64)
    np.add.at(per_atom_forces, pair_idx[:, 0], -forces)
    np.add.at(per_atom_forces, pair_idx[:, 1], forces)

    return per_atom_forces

LJ126

LJ126(epsilon, sigma)

Bases: PairPotential

Lennard-Jones 12-6 pair potential with cutoff.

The potential is defined as: V(r) = 4 * ε * ((σ/r)^12 - (σ/r)^6)

The force is: F(r) = 24 * ε * (2 * (σ/r)^12 - (σ/r)^6) * dr / r^2

Attributes:

Name	Type	Description
epsilon		Depth of potential well for each atom type [energy]
sigma		Finite distance at which potential is zero [length]

Initialize LJ126 potential.

Parameters:

Name	Type	Description	Default
epsilon	float | NDArray[float64]	Depth of potential well, can be scalar or array for multiple types	required
sigma	float | NDArray[float64]	Finite distance at which potential is zero, can be scalar or array for multiple types	required

Source code in src/molpy/potential/pair/lj.py

def __init__(
    self,
    epsilon: float | NDArray[np.float64],
    sigma: float | NDArray[np.float64],
) -> None:
    """
    Initialize LJ126 potential.

    Args:
        epsilon: Depth of potential well, can be scalar or array for multiple types
        sigma: Finite distance at which potential is zero, can be scalar or array for multiple types
    """
    self.epsilon = np.array(epsilon, dtype=np.float64).reshape(-1, 1)
    self.sigma = np.array(sigma, dtype=np.float64).reshape(-1, 1)

    if self.epsilon.shape != self.sigma.shape:
        raise ValueError("epsilon and sigma must have the same shape")

calc_energy

calc_energy(dr, dr_norm, pair_types)

Calculate pair energy.

Parameters:

Name	Type	Description	Default
dr	NDArray[floating]	Pair displacement vectors (shape: (n_pairs, 3))	required
dr_norm	NDArray[floating]	Pair distances (shape: (n_pairs, 1) or (n_pairs,))	required
pair_types	NDArray[integer]	Pair types (shape: (n_pairs,))	required

Returns:

Type	Description
float	Total pair energy

Source code in src/molpy/potential/pair/lj.py

def calc_energy(
    self,
    dr: NDArray[np.floating],
    dr_norm: NDArray[np.floating],
    pair_types: NDArray[np.integer],
) -> float:
    """
    Calculate pair energy.

    Args:
        dr: Pair displacement vectors (shape: (n_pairs, 3))
        dr_norm: Pair distances (shape: (n_pairs, 1) or (n_pairs,))
        pair_types: Pair types (shape: (n_pairs,))

    Returns:
        Total pair energy
    """
    if len(pair_types) == 0:
        return 0.0

    if np.any(pair_types >= len(self.epsilon)) or np.any(pair_types < 0):
        raise ValueError(
            f"pair_types contains invalid indices. Must be in range [0, {len(self.epsilon)})"
        )

    # Ensure dr_norm has correct shape
    if dr_norm.ndim == 1:
        dr_norm = dr_norm[:, None]

    eps = self.epsilon[pair_types]
    sig = self.sigma[pair_types]

    # Calculate energy
    energy = 4 * eps * ((sig / dr_norm) ** 12 - (sig / dr_norm) ** 6)

    return float(np.sum(energy))

calc_forces

calc_forces(dr, dr_norm, pair_types, pair_idx, n_atoms)

Calculate pair forces.

Parameters:

Name	Type	Description	Default
dr	NDArray[floating]	Pair displacement vectors (shape: (n_pairs, 3))	required
dr_norm	NDArray[floating]	Pair distances (shape: (n_pairs, 1) or (n_pairs,))	required
pair_types	NDArray[integer]	Pair types (shape: (n_pairs,))	required
pair_idx	NDArray[integer]	Pair indices (shape: (n_pairs, 2))	required
n_atoms	int	Number of atoms	required

Returns:

Type	Description
NDArray[floating]	Array of forces on each atom (shape: (n_atoms, 3))

Source code in src/molpy/potential/pair/lj.py

def calc_forces(
    self,
    dr: NDArray[np.floating],
    dr_norm: NDArray[np.floating],
    pair_types: NDArray[np.integer],
    pair_idx: NDArray[np.integer],
    n_atoms: int,
) -> NDArray[np.floating]:
    """
    Calculate pair forces.

    Args:
        dr: Pair displacement vectors (shape: (n_pairs, 3))
        dr_norm: Pair distances (shape: (n_pairs, 1) or (n_pairs,))
        pair_types: Pair types (shape: (n_pairs,))
        pair_idx: Pair indices (shape: (n_pairs, 2))
        n_atoms: Number of atoms

    Returns:
        Array of forces on each atom (shape: (n_atoms, 3))
    """
    if len(pair_types) == 0:
        return np.zeros((n_atoms, 3), dtype=np.float64)

    if np.any(pair_types >= len(self.epsilon)) or np.any(pair_types < 0):
        raise ValueError(
            f"pair_types contains invalid indices. Must be in range [0, {len(self.epsilon)})"
        )

    # Ensure dr_norm has correct shape
    if dr_norm.ndim == 1:
        dr_norm = dr_norm[:, None]

    eps = self.epsilon[pair_types]
    sig = self.sigma[pair_types]

    # Calculate force magnitude
    force_mag = (
        24 * eps * (2 * (sig / dr_norm) ** 12 - (sig / dr_norm) ** 6) / (dr_norm**2)
    )

    # Calculate force vectors
    forces = force_mag * dr

    # Accumulate forces on atoms
    per_atom_forces = np.zeros((n_atoms, 3), dtype=np.float64)
    np.add.at(per_atom_forces, pair_idx[:, 0], -forces.squeeze())
    np.add.at(per_atom_forces, pair_idx[:, 1], forces.squeeze())

    return per_atom_forces

LJ126CoulLong

LJ126CoulLong(epsilon, sigma, charges, type_names)

Bases: PairPotential

Combined Lennard-Jones 12-6 and Coulomb long-range pair potential.

This is a composite potential that combines LJ and Coulomb interactions. Uses PairTypeIndexedArray internally for automatic combining rules.

V(r) = V_LJ(r) + V_Coul(r)

Initialize LJ126CoulLong potential.

Parameters:

Name	Type	Description	Default
epsilon	NDArray[float64]	Per-atom-type epsilon values (numpy array)	required
sigma	NDArray[float64]	Per-atom-type sigma values (numpy array)	required
charges	NDArray[float64]	Per-atom-type charges (numpy array)	required
type_names	list[str]	List of atom type names corresponding to array indices	required

Source code in src/molpy/potential/pair/lj.py

def __init__(
    self,
    epsilon: NDArray[np.float64],
    sigma: NDArray[np.float64],
    charges: NDArray[np.float64],
    type_names: list[str],
) -> None:
    """
    Initialize LJ126CoulLong potential.

    Args:
        epsilon: Per-atom-type epsilon values (numpy array)
        sigma: Per-atom-type sigma values (numpy array)
        charges: Per-atom-type charges (numpy array)
        type_names: List of atom type names corresponding to array indices
    """
    from molpy.potential.utils import TypeIndexedArray
    from molpy.potential.pair_params import PairTypeIndexedArray

    # Create dictionaries from arrays and type names
    epsilon_dict = {name: float(eps) for name, eps in zip(type_names, epsilon)}
    sigma_dict = {name: float(sig) for name, sig in zip(type_names, sigma)}
    charge_dict = {name: float(q) for name, q in zip(type_names, charges)}

    # Create TypeIndexedArrays with combining rules
    self.epsilon = PairTypeIndexedArray(epsilon_dict, combining_rule="geometric")
    self.sigma = PairTypeIndexedArray(sigma_dict, combining_rule="arithmetic")
    self.charges = TypeIndexedArray(charge_dict)

calc_energy

calc_energy(dr, dr_norm, pair_types_i, pair_types_j)

Calculate combined LJ + Coulomb energy.

Uses PairTypeIndexedArray to automatically apply combining rules.

Source code in src/molpy/potential/pair/lj.py

def calc_energy(
    self,
    dr: NDArray[np.floating],
    dr_norm: NDArray[np.floating],
    pair_types_i: NDArray,
    pair_types_j: NDArray,
) -> float:
    """
    Calculate combined LJ + Coulomb energy.

    Uses PairTypeIndexedArray to automatically apply combining rules.
    """
    if len(pair_types_i) == 0:
        return 0.0

    # Ensure dr_norm has correct shape
    if dr_norm.ndim == 1:
        dr_norm = dr_norm[:, None]

    # Use PairTypeIndexedArray pair indexing (automatic combining rules)
    pair_types = np.column_stack([pair_types_i, pair_types_j])
    eps = self.epsilon[pair_types]
    sig = self.sigma[pair_types]

    # Ensure correct shape for broadcasting
    if isinstance(eps, np.ndarray) and eps.ndim == 1:
        eps = eps[:, None]
    if isinstance(sig, np.ndarray) and sig.ndim == 1:
        sig = sig[:, None]

    # Calculate LJ energy
    lj_energy = 4 * eps * ((sig / dr_norm) ** 12 - (sig / dr_norm) ** 6)

    # Calculate Coulomb energy
    q_i = self.charges[pair_types_i]
    q_j = self.charges[pair_types_j]
    coul_energy = q_i * q_j / dr_norm.squeeze()

    if isinstance(coul_energy, np.ndarray) and coul_energy.ndim == 1:
        coul_energy = coul_energy[:, None]

    total_energy = lj_energy + coul_energy
    return float(np.sum(total_energy))

calc_forces

calc_forces(dr, dr_norm, pair_types_i, pair_types_j, pair_idx, n_atoms)

Calculate combined LJ + Coulomb forces.

Uses PairTypeIndexedArray to automatically apply combining rules.

Source code in src/molpy/potential/pair/lj.py

def calc_forces(
    self,
    dr: NDArray[np.floating],
    dr_norm: NDArray[np.floating],
    pair_types_i: NDArray,
    pair_types_j: NDArray,
    pair_idx: NDArray[np.integer],
    n_atoms: int,
) -> NDArray[np.floating]:
    """
    Calculate combined LJ + Coulomb forces.

    Uses PairTypeIndexedArray to automatically apply combining rules.
    """
    if len(pair_types_i) == 0:
        return np.zeros((n_atoms, 3), dtype=np.float64)

    # Ensure dr_norm has correct shape
    if dr_norm.ndim == 1:
        dr_norm = dr_norm[:, None]

    # Use PairTypeIndexedArray pair indexing
    pair_types = np.column_stack([pair_types_i, pair_types_j])
    eps = self.epsilon[pair_types]
    sig = self.sigma[pair_types]

    # Ensure correct shape for broadcasting
    if isinstance(eps, np.ndarray) and eps.ndim == 1:
        eps = eps[:, None]
    if isinstance(sig, np.ndarray) and sig.ndim == 1:
        sig = sig[:, None]

    # Calculate LJ force magnitude
    lj_force_mag = (
        24 * eps * (2 * (sig / dr_norm) ** 12 - (sig / dr_norm) ** 6) / (dr_norm**2)
    )

    # Calculate Coulomb force magnitude
    q_i = self.charges[pair_types_i]
    q_j = self.charges[pair_types_j]
    coul_force_mag = q_i * q_j / (dr_norm.squeeze() ** 3)

    if isinstance(coul_force_mag, np.ndarray) and coul_force_mag.ndim == 1:
        coul_force_mag = coul_force_mag[:, None]

    # Total force magnitude
    total_force_mag = lj_force_mag + coul_force_mag

    # Calculate force vectors
    forces = total_force_mag * dr

    # Accumulate forces on atoms
    per_atom_forces = np.zeros((n_atoms, 3), dtype=np.float64)
    np.add.at(per_atom_forces, pair_idx[:, 0], -forces.squeeze())
    np.add.at(per_atom_forces, pair_idx[:, 1], forces.squeeze())

    return per_atom_forces

PairCoulLongStyle

PairCoulLongStyle(cutoff=10.0)

Bases: PairStyle

Coulomb long-range pair style with fixed name='coul/long'.

Parameters:

Name	Type	Description	Default
cutoff	float	Cutoff distance in Angstroms (default: 10.0)	10.0

Source code in src/molpy/potential/pair/lj.py

def __init__(self, cutoff: float = 10.0):
    """
    Args:
        cutoff: Cutoff distance in Angstroms (default: 10.0)
    """
    super().__init__("coul/long", cutoff)

PairLJ126CoulCutStyle

PairLJ126CoulCutStyle(lj_cutoff=10.0, coul_cutoff=10.0, coulomb14scale=0.5, lj14scale=0.5)

Bases: PairStyle

Combined LJ 12-6 and Coulomb cut pair style.

This is a composite style that combines PairLJ126Style and Coulomb cut. The name is 'lj/cut/coul/cut' for LAMMPS compatibility.

Parameters:

Name	Type	Description	Default
lj_cutoff	float	LJ cutoff distance in Angstroms (default: 10.0)	10.0
coul_cutoff	float	Coulomb cutoff distance in Angstroms (default: 10.0)	10.0
coulomb14scale	float	1-4 Coulomb scaling factor (default: 0.5)	0.5
lj14scale	float	1-4 LJ scaling factor (default: 0.5)	0.5

Source code in src/molpy/potential/pair/lj.py

def __init__(
    self,
    lj_cutoff: float = 10.0,
    coul_cutoff: float = 10.0,
    coulomb14scale: float = 0.5,
    lj14scale: float = 0.5,
):
    """
    Args:
        lj_cutoff: LJ cutoff distance in Angstroms (default: 10.0)
        coul_cutoff: Coulomb cutoff distance in Angstroms (default: 10.0)
        coulomb14scale: 1-4 Coulomb scaling factor (default: 0.5)
        lj14scale: 1-4 LJ scaling factor (default: 0.5)
    """
    super().__init__("lj/cut/coul/cut", lj_cutoff, coul_cutoff)
    self.params.kwargs["coulomb14scale"] = coulomb14scale
    self.params.kwargs["lj14scale"] = lj14scale

def_type

def_type(itom, jtom=None, epsilon=0.0, sigma=0.0, charge=0.0, name='')

Define LJ 12-6 pair type (same as PairLJ126Style).

Parameters:

Name	Type	Description	Default
itom	AtomType	First atom type	required
jtom	AtomType | None	Second atom type (None for self-interaction)	None
epsilon	float	LJ epsilon parameter	0.0
sigma	float	LJ sigma parameter	0.0
charge	float	Atomic charge (optional)	0.0
name	str	Optional name	''

Returns:

Type	Description
PairLJ126Type	PairLJ126Type instance

Source code in src/molpy/potential/pair/lj.py

def def_type(
    self,
    itom: AtomType,
    jtom: AtomType | None = None,
    epsilon: float = 0.0,
    sigma: float = 0.0,
    charge: float = 0.0,
    name: str = "",
) -> PairLJ126Type:
    """Define LJ 12-6 pair type (same as PairLJ126Style).

    Args:
        itom: First atom type
        jtom: Second atom type (None for self-interaction)
        epsilon: LJ epsilon parameter
        sigma: LJ sigma parameter
        charge: Atomic charge (optional)
        name: Optional name

    Returns:
        PairLJ126Type instance
    """
    if jtom is None:
        jtom = itom
    if not name:
        name = itom.name if itom == jtom else f"{itom.name}-{jtom.name}"
    pt = PairLJ126Type(name, itom, jtom, epsilon, sigma, charge)
    self.types.add(pt)
    return pt

PairLJ126CoulLongStyle

PairLJ126CoulLongStyle(lj_cutoff=10.0, coul_cutoff=10.0, coulomb14scale=0.5, lj14scale=0.5)

Bases: PairStyle

Combined LJ 12-6 and Coulomb long pair style.

This is a composite style that combines PairLJ126Style and PairCoulLongStyle. The name is 'lj/cut/coul/long' for LAMMPS compatibility.

Parameters:

Name	Type	Description	Default
lj_cutoff	float	LJ cutoff distance in Angstroms (default: 10.0)	10.0
coul_cutoff	float	Coulomb cutoff distance in Angstroms (default: 10.0)	10.0
coulomb14scale	float	1-4 Coulomb scaling factor (default: 0.5)	0.5
lj14scale	float	1-4 LJ scaling factor (default: 0.5)	0.5

Source code in src/molpy/potential/pair/lj.py

def __init__(
    self,
    lj_cutoff: float = 10.0,
    coul_cutoff: float = 10.0,
    coulomb14scale: float = 0.5,
    lj14scale: float = 0.5,
):
    """
    Args:
        lj_cutoff: LJ cutoff distance in Angstroms (default: 10.0)
        coul_cutoff: Coulomb cutoff distance in Angstroms (default: 10.0)
        coulomb14scale: 1-4 Coulomb scaling factor (default: 0.5)
        lj14scale: 1-4 LJ scaling factor (default: 0.5)
    """
    super().__init__("lj/cut/coul/long", lj_cutoff, coul_cutoff)
    self.params.kwargs["coulomb14scale"] = coulomb14scale
    self.params.kwargs["lj14scale"] = lj14scale

def_type

def_type(itom, jtom=None, epsilon=0.0, sigma=0.0, charge=0.0, name='')

Define LJ 12-6 pair type (same as PairLJ126Style).

Parameters:

Name	Type	Description	Default
itom	AtomType	First atom type	required
jtom	AtomType | None	Second atom type (None for self-interaction)	None
epsilon	float	LJ epsilon parameter	0.0
sigma	float	LJ sigma parameter	0.0
charge	float	Atomic charge (optional)	0.0
name	str	Optional name	''

Returns:

Type	Description
PairLJ126Type	PairLJ126Type instance

Source code in src/molpy/potential/pair/lj.py

def def_type(
    self,
    itom: AtomType,
    jtom: AtomType | None = None,
    epsilon: float = 0.0,
    sigma: float = 0.0,
    charge: float = 0.0,
    name: str = "",
) -> PairLJ126Type:
    """Define LJ 12-6 pair type (same as PairLJ126Style).

    Args:
        itom: First atom type
        jtom: Second atom type (None for self-interaction)
        epsilon: LJ epsilon parameter
        sigma: LJ sigma parameter
        charge: Atomic charge (optional)
        name: Optional name

    Returns:
        PairLJ126Type instance
    """
    if jtom is None:
        jtom = itom
    if not name:
        name = itom.name if itom == jtom else f"{itom.name}-{jtom.name}"
    pt = PairLJ126Type(name, itom, jtom, epsilon, sigma, charge)
    self.types.add(pt)
    return pt

to_potential

to_potential()

Convert this style to a Potential object.

Source code in src/molpy/potential/pair/lj.py

def to_potential(self):
    """Convert this style to a Potential object."""
    from molpy.core.forcefield import PairType

    pair_types = list(self.types.bucket(PairType))
    if not pair_types:
        raise ValueError("No pair types defined in style")

    # Extract parameters as lists
    type_names = []
    epsilon_list = []
    sigma_list = []
    charge_list = []

    for pt in pair_types:
        epsilon = pt.params.kwargs.get("epsilon")
        sigma = pt.params.kwargs.get("sigma")
        charge = pt.params.kwargs.get("charge", 0.0)

        if epsilon is None or sigma is None:
            raise ValueError(
                f"PairType '{pt.name}' is missing required parameters: "
                f"epsilon={epsilon}, sigma={sigma}"
            )

        type_names.append(pt.itom.name)
        epsilon_list.append(epsilon)
        sigma_list.append(sigma)
        charge_list.append(charge)

    # Convert to numpy arrays
    epsilon_array = np.array(epsilon_list, dtype=np.float64)
    sigma_array = np.array(sigma_list, dtype=np.float64)
    charges_array = np.array(charge_list, dtype=np.float64)

    # Create Potential instance with numpy arrays
    return LJ126CoulLong(
        epsilon=epsilon_array,
        sigma=sigma_array,
        charges=charges_array,
        type_names=type_names,
    )

PairLJ126Style

PairLJ126Style(cutoff=10.0)

Bases: PairStyle

Lennard-Jones 12-6 pair style with fixed name='lj126'.

Parameters:

Name	Type	Description	Default
cutoff	float	Cutoff distance in Angstroms (default: 10.0)	10.0

Source code in src/molpy/potential/pair/lj.py

def __init__(self, cutoff: float = 10.0):
    """
    Args:
        cutoff: Cutoff distance in Angstroms (default: 10.0)
    """
    super().__init__("lj126", cutoff)

def_type

def_type(itom, jtom=None, epsilon=0.0, sigma=0.0, charge=0.0, name='')

Define LJ 12-6 pair type.

Parameters:

Name	Type	Description	Default
itom	AtomType	First atom type	required
jtom	AtomType | None	Second atom type (None for self-interaction)	None
epsilon	float	LJ epsilon parameter	0.0
sigma	float	LJ sigma parameter	0.0
charge	float	Atomic charge (optional)	0.0
name	str	Optional name	''

Returns:

Type	Description
PairLJ126Type	PairLJ126Type instance

Source code in src/molpy/potential/pair/lj.py

def def_type(
    self,
    itom: AtomType,
    jtom: AtomType | None = None,
    epsilon: float = 0.0,
    sigma: float = 0.0,
    charge: float = 0.0,
    name: str = "",
) -> PairLJ126Type:
    """Define LJ 12-6 pair type.

    Args:
        itom: First atom type
        jtom: Second atom type (None for self-interaction)
        epsilon: LJ epsilon parameter
        sigma: LJ sigma parameter
        charge: Atomic charge (optional)
        name: Optional name

    Returns:
        PairLJ126Type instance
    """
    if jtom is None:
        jtom = itom
    if not name:
        name = itom.name if itom == jtom else f"{itom.name}-{jtom.name}"
    pt = PairLJ126Type(name, itom, jtom, epsilon, sigma, charge)
    self.types.add(pt)
    return pt

PairLJ126Type

PairLJ126Type(name, itom, jtom=None, epsilon=0.0, sigma=0.0, charge=0.0)

Bases: PairType

Lennard-Jones 12-6 pair type with epsilon and sigma parameters.

Parameters:

Name	Type	Description	Default
name	str	Type name	required
itom	AtomType	First atom type	required
jtom	AtomType | None	Second atom type (None for self-interaction)	None
epsilon	float	LJ epsilon parameter	0.0
sigma	float	LJ sigma parameter	0.0
charge	float	Atomic charge (optional)	0.0

Source code in src/molpy/potential/pair/lj.py

def __init__(
    self,
    name: str,
    itom: AtomType,
    jtom: AtomType | None = None,
    epsilon: float = 0.0,
    sigma: float = 0.0,
    charge: float = 0.0,
):
    """
    Args:
        name: Type name
        itom: First atom type
        jtom: Second atom type (None for self-interaction)
        epsilon: LJ epsilon parameter
        sigma: LJ sigma parameter
        charge: Atomic charge (optional)
    """
    if jtom is None:
        jtom = itom
    super().__init__(name, itom, jtom, epsilon=epsilon, sigma=sigma, charge=charge)

PairPotential

Bases: Potential

Base class for pair potentials.

Pair Params

Specialized TypeIndexedArray for pair potentials with combining rules.

PairTypeIndexedArray

PairTypeIndexedArray(data, combining_rule='geometric')

Bases: TypeIndexedArray

Specialized TypeIndexedArray for pair potential parameters.

Handles combining rules (Lorentz-Berthelot, geometric, etc.) for computing pair parameters from individual atom type parameters.

For pair potentials, we store per-atom-type parameters (epsilon, sigma) and apply combining rules when indexing with two atom types.

Examples:

>>> # Create with per-atom-type parameters
>>> epsilon = PairTypeIndexedArray(
...     {'opls_135': 0.066, 'opls_140': 0.030},
...     combining_rule='geometric'
... )
>>>
>>> # Index with single type (self-interaction)
>>> epsilon['opls_135']  # 0.066
>>>
>>> # Index with two types (cross-interaction, uses combining rule)
>>> epsilon[('opls_135', 'opls_140')]  # sqrt(0.066 * 0.030) = 0.0445
>>>
>>> # Array indexing with pairs
>>> pairs = np.array([['opls_135', 'opls_140'], ['opls_140', 'opls_140']])
>>> epsilon[pairs]  # [0.0445, 0.030]

Initialize PairTypeIndexedArray.

Parameters:

Name	Type	Description	Default
data	dict[str, float] | NDArray[floating] | float	Dictionary mapping atom type names to parameters, or array	required
combining_rule	Literal['geometric', 'arithmetic', 'harmonic']	Rule for combining parameters: - 'geometric': sqrt(a * b) - for epsilon - 'arithmetic': (a + b) / 2 - for sigma - 'harmonic': 2 * a * b / (a + b) - alternative	'geometric'

Source code in src/molpy/potential/pair_params.py

def __init__(
    self,
    data: dict[str, float] | NDArray[np.floating] | float,
    combining_rule: Literal["geometric", "arithmetic", "harmonic"] = "geometric",
):
    """
    Initialize PairTypeIndexedArray.

    Args:
        data: Dictionary mapping atom type names to parameters, or array
        combining_rule: Rule for combining parameters:
            - 'geometric': sqrt(a * b) - for epsilon
            - 'arithmetic': (a + b) / 2 - for sigma
            - 'harmonic': 2 * a * b / (a + b) - alternative
    """
    super().__init__(data)
    self.combining_rule = combining_rule

get_pair_array

get_pair_array(type_pairs)

Get combined parameters for an array of type pairs.

Parameters:

Name	Type	Description	Default
type_pairs	NDArray	Array of shape (n_pairs, 2) with type labels or indices	required

Returns:

Type	Description
NDArray[floating]	Array of combined parameters for each pair

Source code in src/molpy/potential/pair_params.py

def get_pair_array(self, type_pairs: NDArray) -> NDArray[np.floating]:
    """
    Get combined parameters for an array of type pairs.

    Args:
        type_pairs: Array of shape (n_pairs, 2) with type labels or indices

    Returns:
        Array of combined parameters for each pair
    """
    return self[type_pairs]

get_pair_matrix

get_pair_matrix()

Get full pair parameter matrix for all type combinations.

Returns:

Type	Description
NDArray[floating]	Matrix of shape (n_types, n_types) with combined parameters
NDArray[floating]	for all pairs of atom types

Source code in src/molpy/potential/pair_params.py

def get_pair_matrix(self) -> NDArray[np.floating]:
    """
    Get full pair parameter matrix for all type combinations.

    Returns:
        Matrix of shape (n_types, n_types) with combined parameters
        for all pairs of atom types
    """
    if not self._use_labels:
        raise ValueError("Cannot create pair matrix without type labels")

    n_types = len(self._values)
    matrix = np.zeros((n_types, n_types))

    for i in range(n_types):
        for j in range(n_types):
            matrix[i, j] = self._combine(self._values[i], self._values[j])

    return matrix

create_lj_parameters

create_lj_parameters(epsilon_dict, sigma_dict)

Create LJ parameter arrays with standard Lorentz-Berthelot combining rules.

Parameters:

Name	Type	Description	Default
epsilon_dict	dict[str, float]	Per-atom-type epsilon values	required
sigma_dict	dict[str, float]	Per-atom-type sigma values	required

Returns:

Type	Description
tuple[PairTypeIndexedArray, PairTypeIndexedArray]	Tuple of (epsilon_array, sigma_array) with combining rules applied

Example

epsilon_dict = {'opls_135': 0.066, 'opls_140': 0.030} sigma_dict = {'opls_135': 3.5, 'opls_140': 2.5} epsilon, sigma = create_lj_parameters(epsilon_dict, sigma_dict)

Get cross-interaction parameters

eps_ij = epsilon[('opls_135', 'opls_140')] # sqrt(0.066 * 0.030) sig_ij = sigma[('opls_135', 'opls_140')] # (3.5 + 2.5) / 2

Source code in src/molpy/potential/pair_params.py

def create_lj_parameters(
    epsilon_dict: dict[str, float],
    sigma_dict: dict[str, float],
) -> tuple[PairTypeIndexedArray, PairTypeIndexedArray]:
    """
    Create LJ parameter arrays with standard Lorentz-Berthelot combining rules.

    Args:
        epsilon_dict: Per-atom-type epsilon values
        sigma_dict: Per-atom-type sigma values

    Returns:
        Tuple of (epsilon_array, sigma_array) with combining rules applied

    Example:
        >>> epsilon_dict = {'opls_135': 0.066, 'opls_140': 0.030}
        >>> sigma_dict = {'opls_135': 3.5, 'opls_140': 2.5}
        >>> epsilon, sigma = create_lj_parameters(epsilon_dict, sigma_dict)
        >>>
        >>> # Get cross-interaction parameters
        >>> eps_ij = epsilon[('opls_135', 'opls_140')]  # sqrt(0.066 * 0.030)
        >>> sig_ij = sigma[('opls_135', 'opls_140')]    # (3.5 + 2.5) / 2
    """
    # Epsilon uses geometric mean (Lorentz-Berthelot)
    epsilon = PairTypeIndexedArray(epsilon_dict, combining_rule="geometric")

    # Sigma uses arithmetic mean (Lorentz-Berthelot)
    sigma = PairTypeIndexedArray(sigma_dict, combining_rule="arithmetic")

    return epsilon, sigma

Utils

Utility functions and classes for potentials.

TypeIndexedArray

TypeIndexedArray(data)

Array-like container that supports both integer and string type name indexing.

This class allows potentials to accept either integer indices or string type labels for indexing parameters. It maintains an internal mapping between type names and indices.

Examples:

>>> # Create from dictionary (type name -> value)
>>> k = TypeIndexedArray({"CT-CT": 100.0, "CT-OH": 80.0})
>>> k[0]  # Access by integer index
100.0
>>> k["CT-CT"]  # Access by type name
100.0
>>> k[np.array([0, 1])]  # Array indexing with integers
array([100.,  80.])
>>> k[np.array(["CT-CT", "CT-OH"])]  # Array indexing with strings
array([100.,  80.])

>>> # Create from array (for backward compatibility)
>>> k = TypeIndexedArray(np.array([100.0, 80.0]))
>>> k[0]
100.0

Initialize TypeIndexedArray.

Parameters:

Name	Type	Description	Default
data	dict[str, float] | NDArray[floating] | float	Either a dictionary mapping type names to values, or an array/float for backward compatibility	required

Source code in src/molpy/potential/utils.py

def __init__(self, data: dict[str, float] | NDArray[np.floating] | float):
    """
    Initialize TypeIndexedArray.

    Args:
        data: Either a dictionary mapping type names to values, or an array/float
             for backward compatibility
    """
    if isinstance(data, dict):
        # Dictionary mode: maintain type name -> index mapping
        self._type_names = list(data.keys())
        self._type_to_idx = {name: idx for idx, name in enumerate(self._type_names)}
        self._values = np.array(
            [data[name] for name in self._type_names], dtype=np.float64
        )
        self._use_labels = True
    else:
        # Array mode: traditional integer indexing only
        self._values = np.array(data, dtype=np.float64)
        if self._values.ndim == 0:
            self._values = self._values.reshape(1)
        self._type_names = None
        self._type_to_idx = None
        self._use_labels = False

type_names property

type_names

Get the list of type names (if available).

values property

values

Get the underlying values array.

reshape

reshape(*args, **kwargs)

Reshape the values array (for backward compatibility).

Source code in src/molpy/potential/utils.py

def reshape(self, *args, **kwargs):
    """Reshape the values array (for backward compatibility)."""
    self._values = self._values.reshape(*args, **kwargs)
    return self

calc_energy_from_frame

calc_energy_from_frame(potential, frame)

Calculate energy from Frame for a potential.

This is a convenience function that extracts the necessary data from Frame and calls the potential's calc_energy method.

Parameters:

Name	Type	Description	Default
potential	Potential	Potential instance	required
frame	Frame	Frame containing the necessary data	required

Returns:

Type	Description
float	Potential energy

Raises:

Type	Description
TypeError	If potential type is not recognized
ValueError	If required data is missing from frame

Source code in src/molpy/potential/utils.py

def calc_energy_from_frame(potential: "Potential", frame: Frame) -> float:
    """
    Calculate energy from Frame for a potential.

    This is a convenience function that extracts the necessary data from Frame
    and calls the potential's calc_energy method.

    Args:
        potential: Potential instance
        frame: Frame containing the necessary data

    Returns:
        Potential energy

    Raises:
        TypeError: If potential type is not recognized
        ValueError: If required data is missing from frame
    """
    # Check if this is a wrapper (has a 'potential' attribute that wraps another potential)
    # Wrappers like BondPotentialWrapper/AnglePotentialWrapper handle data extraction internally
    if hasattr(potential, "potential") and callable(
        getattr(potential, "calc_energy", None)
    ):
        # Get the signature of calc_energy to check if it accepts only frame
        import inspect

        sig = inspect.signature(potential.calc_energy)
        params = list(sig.parameters.keys())
        # If calc_energy takes only 1 parameter (frame), call it directly
        if len(params) == 1:
            return potential.calc_energy(frame)

    # Check potential type by looking at its type attribute
    potential_type = getattr(potential, "type", None)

    if potential_type == "bond":
        r, bond_idx, bond_types = extract_bond_data(frame)
        return potential.calc_energy(r, bond_idx, bond_types)
    elif potential_type == "angle":
        r, angle_idx, angle_types = extract_angle_data(frame)
        return potential.calc_energy(r, angle_idx, angle_types)
    elif potential_type == "pair":
        # Check if it's LJ126CoulLong (combined) or separate LJ/Coul
        potential_name = getattr(potential, "name", "")
        if "lj" in potential_name.lower() and "coul" in potential_name.lower():
            # Combined LJ + Coulomb potential (like lj/cut/coul/long)
            dr, dr_norm, pair_types_i, pair_types_j, _, _ = extract_pair_data(frame)
            return potential.calc_energy(dr, dr_norm, pair_types_i, pair_types_j)
        elif "lj" in potential_name.lower():
            # Pure LJ potential
            dr, dr_norm, pair_types_i, pair_types_j, _, _ = extract_pair_data(frame)
            return potential.calc_energy(dr, dr_norm, pair_types_i, pair_types_j)
        elif "coul" in potential_name.lower():
            r, pair_idx, charges = extract_coul_data(frame)
            return potential.calc_energy(r, pair_idx, charges)
        else:
            raise TypeError(f"Unknown pair potential type: {potential_name}")
    else:
        raise TypeError(f"Unknown potential type: {potential_type}")

calc_energy_from_frame_multi

calc_energy_from_frame_multi(potentials, frame)

Calculate total energy from multiple potentials.

Source code in src/molpy/potential/utils.py

def calc_energy_from_frame_multi(potentials, frame: Frame) -> float:
    """Calculate total energy from multiple potentials."""
    total_energy = 0.0
    for potential in potentials:
        total_energy += calc_energy_from_frame(potential, frame)
    return total_energy

calc_forces_from_frame

calc_forces_from_frame(potential, frame)

Calculate forces from Frame for a potential.

This is a convenience function that extracts the necessary data from Frame and calls the potential's calc_forces method.

Parameters:

Name	Type	Description	Default
potential	Potential	Potential instance	required
frame	Frame	Frame containing the necessary data	required

Returns:

Type	Description
NDArray	Array of forces on each atom (shape: (n_atoms, 3))

Raises:

Type	Description
TypeError	If potential type is not recognized
ValueError	If required data is missing from frame

Source code in src/molpy/potential/utils.py

def calc_forces_from_frame(potential: "Potential", frame: Frame) -> NDArray:
    """
    Calculate forces from Frame for a potential.

    This is a convenience function that extracts the necessary data from Frame
    and calls the potential's calc_forces method.

    Args:
        potential: Potential instance
        frame: Frame containing the necessary data

    Returns:
        Array of forces on each atom (shape: (n_atoms, 3))

    Raises:
        TypeError: If potential type is not recognized
        ValueError: If required data is missing from frame
    """
    # Check if this is a wrapper (has a 'potential' attribute that wraps another potential)
    # Wrappers like BondPotentialWrapper/AnglePotentialWrapper handle data extraction internally
    if hasattr(potential, "potential") and callable(
        getattr(potential, "calc_forces", None)
    ):
        # Get the signature of calc_forces to check if it accepts only frame
        import inspect

        sig = inspect.signature(potential.calc_forces)
        params = list(sig.parameters.keys())
        # If calc_forces takes only 1 parameter (frame), call it directly
        if len(params) == 1:
            return potential.calc_forces(frame)

    # Check potential type by looking at its type attribute
    potential_type = getattr(potential, "type", None)

    if potential_type == "bond":
        r, bond_idx, bond_types = extract_bond_data(frame)
        return potential.calc_forces(r, bond_idx, bond_types)
    elif potential_type == "angle":
        r, angle_idx, angle_types = extract_angle_data(frame)
        return potential.calc_forces(r, angle_idx, angle_types)
    elif potential_type == "pair":
        # Check if it's LJ126CoulLong (combined) or separate LJ/Coul
        potential_name = getattr(potential, "name", "")
        if "lj" in potential_name.lower() and "coul" in potential_name.lower():
            # Combined LJ + Coulomb potential (like lj/cut/coul/long)
            dr, dr_norm, pair_types_i, pair_types_j, pair_idx, n_atoms = (
                extract_pair_data(frame)
            )
            return potential.calc_forces(
                dr, dr_norm, pair_types_i, pair_types_j, pair_idx, n_atoms
            )
        elif "lj" in potential_name.lower():
            # Pure LJ potential
            dr, dr_norm, pair_types_i, pair_types_j, pair_idx, n_atoms = (
                extract_pair_data(frame)
            )
            return potential.calc_forces(
                dr, dr_norm, pair_types_i, pair_types_j, pair_idx, n_atoms
            )
        elif "coul" in potential_name.lower():
            r, pair_idx, charges = extract_coul_data(frame)
            return potential.calc_forces(r, pair_idx, charges)
        else:
            raise TypeError(f"Unknown pair potential type: {potential_name}")
    else:
        raise TypeError(f"Unknown potential type: {potential_type}")

extract_angle_data

extract_angle_data(frame)

Extract angle data from Frame.

Returns:

Type	Description
NDArray	(r, angle_idx, angle_types) where:
NDArray	r: atom coordinates (n_atoms, 3)
NDArray	angle_idx: angle indices (n_angles, 3)
tuple[NDArray, NDArray, NDArray]	angle_types: angle types (n_angles,) - can be integers or strings

Source code in src/molpy/potential/utils.py

def extract_angle_data(frame: Frame) -> tuple[NDArray, NDArray, NDArray]:
    """Extract angle data from Frame.

    Returns:
        (r, angle_idx, angle_types) where:
        - r: atom coordinates (n_atoms, 3)
        - angle_idx: angle indices (n_angles, 3)
        - angle_types: angle types (n_angles,) - can be integers or strings
    """
    # Get coordinates from x, y, z fields (never use xyz)
    x = frame["atoms"]["x"]
    y = frame["atoms"]["y"]
    z = frame["atoms"]["z"]
    r = np.column_stack([x, y, z])
    angles = frame["angles"]
    angle_idx = angles[["atomi", "atomj", "atomk"]]
    angle_types = angles["type"]
    return r, angle_idx, angle_types

extract_bond_data

extract_bond_data(frame)

Extract bond data from Frame.

Returns:

Type	Description
NDArray	(r, bond_idx, bond_types) where:
NDArray	r: atom coordinates (n_atoms, 3)
NDArray	bond_idx: bond indices (n_bonds, 2)
tuple[NDArray, NDArray, NDArray]	bond_types: bond types (n_bonds,) - can be integers or strings

Source code in src/molpy/potential/utils.py

def extract_bond_data(frame: Frame) -> tuple[NDArray, NDArray, NDArray]:
    """Extract bond data from Frame.

    Returns:
        (r, bond_idx, bond_types) where:
        - r: atom coordinates (n_atoms, 3)
        - bond_idx: bond indices (n_bonds, 2)
        - bond_types: bond types (n_bonds,) - can be integers or strings
    """
    atoms = frame["atoms"]
    x = atoms["x"]
    y = atoms["y"]
    z = atoms["z"]
    r = np.column_stack([x, y, z])
    bonds = frame["bonds"]
    bond_idx = bonds[["atomi", "atomj"]]
    bond_types = bonds["type"]
    return r, bond_idx, bond_types

extract_coul_data

extract_coul_data(frame)

Extract Coulomb interaction data from Frame.

Source code in src/molpy/potential/utils.py

def extract_coul_data(frame: Frame):
    """Extract Coulomb interaction data from Frame."""
    # Implementation depends on Coulomb potential type
    raise NotImplementedError("Coulomb data extraction not yet implemented")

extract_pair_data

extract_pair_data(frame)

Extract pair interaction data from Frame.

Generates all pairwise interactions between atoms and calculates displacement vectors and distances.

Returns:

Type	Description
	(dr, dr_norm, pair_types_i, pair_types_j, pair_idx, n_atoms) where:
	dr: displacement vectors (n_pairs, 3)
	dr_norm: pair distances (n_pairs,)
	pair_types_i: atom types for first atom in each pair (n_pairs,)
	pair_types_j: atom types for second atom in each pair (n_pairs,)
	pair_idx: pair indices (n_pairs, 2)
	n_atoms: number of atoms

Source code in src/molpy/potential/utils.py

def extract_pair_data(frame: Frame):
    """Extract pair interaction data from Frame.

    Generates all pairwise interactions between atoms and calculates
    displacement vectors and distances.

    Returns:
        (dr, dr_norm, pair_types_i, pair_types_j, pair_idx, n_atoms) where:
        - dr: displacement vectors (n_pairs, 3)
        - dr_norm: pair distances (n_pairs,)
        - pair_types_i: atom types for first atom in each pair (n_pairs,)
        - pair_types_j: atom types for second atom in each pair (n_pairs,)
        - pair_idx: pair indices (n_pairs, 2)
        - n_atoms: number of atoms
    """
    # Get coordinates from x, y, z fields
    x = frame["atoms"]["x"]
    y = frame["atoms"]["y"]
    z = frame["atoms"]["z"]
    r = np.column_stack([x, y, z])
    n_atoms = len(r)

    # Get atom types
    atom_types = frame["atoms"]["type"]

    # Generate all pairs (i, j) where i < j
    # This avoids double counting and self-interactions
    pair_idx = []
    pair_types_i_list = []
    pair_types_j_list = []

    for i in range(n_atoms):
        for j in range(i + 1, n_atoms):
            pair_idx.append([i, j])
            # Store both atom types for proper parameter lookup
            pair_types_i_list.append(atom_types[i])
            pair_types_j_list.append(atom_types[j])

    pair_idx = np.array(pair_idx, dtype=np.int64)
    # Keep pair_types as-is (can be strings or integers)
    pair_types_i = np.array(pair_types_i_list)
    pair_types_j = np.array(pair_types_j_list)

    # Calculate displacement vectors
    if len(pair_idx) > 0:
        dr = r[pair_idx[:, 1]] - r[pair_idx[:, 0]]
        dr_norm = np.linalg.norm(dr, axis=1)
    else:
        dr = np.zeros((0, 3), dtype=np.float64)
        dr_norm = np.zeros(0, dtype=np.float64)

    return dr, dr_norm, pair_types_i, pair_types_j, pair_idx, n_atoms