Force fields and potentials¶

A ForceField is a parameter database: typed tables for bonded and non-bonded interactions, bonds, angles, dihedrals, pairs. A Potential is the executable kernel that evaluates energies/forces from those parameters.

When to use You want to parameterize an Atomistic structure, assign types and parameters., You want to convert parameter tables into executable potentials for energy/force evaluation., and You want to load a standard FF, e.g., OPLS-AA or build a minimal FF from scratch..

How this connects to Atomistic / Topology Atomistic holds the structure, atoms, bonds and may also hold derived Angle / Dihedral entities., Topology, see topology.ipynb is the graph view used to derive angles/dihedrals; mol.get_topo, gen_angle=True, gen_dihe=True can populate them., and Force-field types, AtomType/BondType/AngleType/DihedralType/PairType are what you assign to those structural entities, typically via typifiers..

In this tutorial you will: Create an AtomisticForcefield, Define atom/bond/angle/dihedral/pair types, and Convert styles, or the full force field into potentials.

Two common workflows: Load & apply: load a standard FF and typify a structure, see ../api/typifier.md. and Define manually: build a small, custom FF from scratch, this notebook..

1. Creating a force field¶

Start with an empty force field. You typically want AtomisticForcefield for atomistic workflows because it has convenience constructors for common styles.

Two classes: ForceField: generic container for styles + typed parameter tables and AtomisticForcefield: atomistic-oriented ForceField with helpers like def_atomstyle, def_bondstyle, def_anglestyle, def_dihedralstyle, def_pairstyle.

In [1]:

import molpy as mp

# Method 1: Create base ForceField
ff = mp.ForceField(name="MyCustomFF", units="real")

# Method 2: Create AtomisticForcefield (recommended for atomistic simulations)
ff_atomistic = mp.AtomisticForcefield(name="MyAtomisticFF", units="real")

print("Base ForceField:")
print(f"  {ff}")
print(f"  Name: {ff.name}, Units: {ff.units}")

print("\nAtomisticForcefield:")
print(f"  {ff_atomistic}")
print(f"  Name: {ff_atomistic.name}, Units: {ff_atomistic.units}")

# Use AtomisticForcefield for the rest of the tutorial
ff = ff_atomistic

import molpy as mp # Method 1: Create base ForceField ff = mp.ForceField(name="MyCustomFF", units="real") # Method 2: Create AtomisticForcefield (recommended for atomistic simulations) ff_atomistic = mp.AtomisticForcefield(name="MyAtomisticFF", units="real") print("Base ForceField:") print(f" {ff}") print(f" Name: {ff.name}, Units: {ff.units}") print("\nAtomisticForcefield:") print(f" {ff_atomistic}") print(f" Name: {ff_atomistic.name}, Units: {ff_atomistic.units}") # Use AtomisticForcefield for the rest of the tutorial ff = ff_atomistic

Base ForceField:
  <ForceField: MyCustomFF>
  Name: MyCustomFF, Units: real

AtomisticForcefield:
  <AtomisticForcefield: MyAtomisticFF>
  Name: MyAtomisticFF, Units: real

2. Defining atom types¶

Atom types are the base layer of most force fields. A type typically stores: Name: unique identifier, e.g., CT, HC, OH, Mass: g/mol, Charge: elementary charge units, and Non-bonded parameters, optional: e.g., Lennard–Jones sigma / epsilon.

Atom types live inside an AtomStyle. You’ll define an AtomStyle, then add types with atom_style.def_type, ....

In [2]:

# Method 1: Using built-in style method (recommended for AtomisticForcefield)
atom_style = ff.def_atomstyle("full")  # 'full' includes bonds, angles, dihedrals, and charges

# Method 2: Using generic def_style() (works for any ForceField)
# atom_style = ff.def_style(mp.AtomStyle("full"))

# Use the style's def_type() method to create atom types
ct = atom_style.def_type("CT", mass=12.011, charge=-0.18)  # Aliphatic carbon
hc = atom_style.def_type("HC", mass=1.008, charge=0.06)    # Hydrogen on carbon

# Add Lennard-Jones parameters (optional)
# sigma in Angstroms, epsilon in kcal/mol
ct["sigma"] = 3.50
ct["epsilon"] = 0.066

hc["sigma"] = 2.50
hc["epsilon"] = 0.030

print(f"Defined {len(ff.get_types(mp.AtomType))} atom types:")
for atype in ff.get_types(mp.AtomType):
    mass = atype["mass"]
    charge = atype["charge"]
    print(f"  {atype.name}: mass={mass:.3f}, charge={charge:.3f}")

# Method 1: Using built-in style method (recommended for AtomisticForcefield) atom_style = ff.def_atomstyle("full") # 'full' includes bonds, angles, dihedrals, and charges # Method 2: Using generic def_style() (works for any ForceField) # atom_style = ff.def_style(mp.AtomStyle("full")) # Use the style's def_type() method to create atom types ct = atom_style.def_type("CT", mass=12.011, charge=-0.18) # Aliphatic carbon hc = atom_style.def_type("HC", mass=1.008, charge=0.06) # Hydrogen on carbon # Add Lennard-Jones parameters (optional) # sigma in Angstroms, epsilon in kcal/mol ct["sigma"] = 3.50 ct["epsilon"] = 0.066 hc["sigma"] = 2.50 hc["epsilon"] = 0.030 print(f"Defined {len(ff.get_types(mp.AtomType))} atom types:") for atype in ff.get_types(mp.AtomType): mass = atype["mass"] charge = atype["charge"] print(f" {atype.name}: mass={mass:.3f}, charge={charge:.3f}")

Defined 2 atom types:
  HC: mass=1.008, charge=0.060
  CT: mass=12.011, charge=-0.180

3. Built-in styles¶

MolPy provides a set of built-in interaction styles that describe how different parts of a molecular system interact energetically. These styles correspond to the standard interaction families used in classical force fields, such as bonded terms (bonds, angles, dihedrals) and non-bonded terms (pair interactions). Rather than hard-coding specific functional forms, MolPy represents each interaction family through a style object, which defines the mathematical form of the interaction and the parameters it requires.

Within an AtomisticForcefield, these styles are created and registered through convenience methods such as def_atomstyle, def_bondstyle, def_anglestyle, def_dihedralstyle, def_improperstyle, and def_pairstyle. Each call defines a reusable interaction template that can later be assigned to specific atoms or topological entities. This separation allows the same functional form to be reused across many bonds or angles while keeping the underlying topology independent of the force-field details.

In [3]:

# Method 1: Using built-in style method (recommended)
bond_style = ff.def_bondstyle("harmonic")

# Method 2: Using generic def_style() (alternative)
# bond_style = ff.def_style(mp.BondStyle("harmonic"))

# Use the style's def_type() method to create bond types
ct_hc_bond = bond_style.def_type(ct, hc, k=340.0, r0=1.09)   # C-H bond
ct_ct_bond = bond_style.def_type(ct, ct, k=268.0, r0=1.529)  # C-C bond

print(f"Defined {len(ff.get_types(mp.BondType))} bond types:")
for btype in ff.get_types(mp.BondType):
    k = btype["k"]
    r0 = btype["r0"]
    print(f"  {btype.name}: k={k:.1f} kcal/mol/Ų, r0={r0:.3f} Å")

# Method 1: Using built-in style method (recommended) bond_style = ff.def_bondstyle("harmonic") # Method 2: Using generic def_style() (alternative) # bond_style = ff.def_style(mp.BondStyle("harmonic")) # Use the style's def_type() method to create bond types ct_hc_bond = bond_style.def_type(ct, hc, k=340.0, r0=1.09) # C-H bond ct_ct_bond = bond_style.def_type(ct, ct, k=268.0, r0=1.529) # C-C bond print(f"Defined {len(ff.get_types(mp.BondType))} bond types:") for btype in ff.get_types(mp.BondType): k = btype["k"] r0 = btype["r0"] print(f" {btype.name}: k={k:.1f} kcal/mol/Ų, r0={r0:.3f} Å")

Defined 2 bond types:
  CT-HC: k=340.0 kcal/mol/Ų, r0=1.090 Å
  CT-CT: k=268.0 kcal/mol/Ų, r0=1.529 Å

4. Defining Arbitrary Styles¶

You can define styles with any name using the generic def_style() method. This is useful when: You want a custom style name, You're using a style not covered by built-in methods, and You're defining a custom Style class.

The style instance must be an instantiated Style object, not a class.

In [4]:

# Example 1: Using built-in method (recommended)
angle_style = ff.def_anglestyle("harmonic")

# Example 2: Using generic def_style() with any name
# angle_style = ff.def_style(mp.AngleStyle("harmonic"))
# angle_style = ff.def_style(mp.AngleStyle("my_custom_angle_style"))  # Custom name works too!

# Use the style's def_type() method to create angle types
hc_ct_hc = angle_style.def_type(hc, ct, hc, k=33.0, theta0=107.8)   # H-C-H angle
hc_ct_ct = angle_style.def_type(hc, ct, ct, k=37.5, theta0=110.7)   # H-C-C angle
ct_ct_ct = angle_style.def_type(ct, ct, ct, k=58.35, theta0=112.7)  # C-C-C angle

print(f"Defined {len(ff.get_types(mp.AngleType))} angle types:")
for atype in ff.get_types(mp.AngleType):
    k = atype["k"]
    theta0 = atype["theta0"]
    print(f"  {atype.name}: k={k:.2f} kcal/mol/rad², θ0={theta0:.1f}°")

# Example 3: Define pair style (non-bonded interactions)
pair_style = ff.def_pairstyle("lj126/cut")
pair_style.def_type(ct, epsilon=0.066, sigma=3.50)
pair_style.def_type(hc, epsilon=0.030, sigma=2.50)

print(f"\nDefined {len(ff.get_types(mp.PairType))} pair types:")
for ptype in ff.get_types(mp.PairType):
    epsilon = ptype["epsilon"]
    sigma = ptype["sigma"]
    print(f"  {ptype.name}: ε={epsilon:.3f} kcal/mol, σ={sigma:.2f} Å")

# Example 1: Using built-in method (recommended) angle_style = ff.def_anglestyle("harmonic") # Example 2: Using generic def_style() with any name # angle_style = ff.def_style(mp.AngleStyle("harmonic")) # angle_style = ff.def_style(mp.AngleStyle("my_custom_angle_style")) # Custom name works too! # Use the style's def_type() method to create angle types hc_ct_hc = angle_style.def_type(hc, ct, hc, k=33.0, theta0=107.8) # H-C-H angle hc_ct_ct = angle_style.def_type(hc, ct, ct, k=37.5, theta0=110.7) # H-C-C angle ct_ct_ct = angle_style.def_type(ct, ct, ct, k=58.35, theta0=112.7) # C-C-C angle print(f"Defined {len(ff.get_types(mp.AngleType))} angle types:") for atype in ff.get_types(mp.AngleType): k = atype["k"] theta0 = atype["theta0"] print(f" {atype.name}: k={k:.2f} kcal/mol/rad², θ0={theta0:.1f}°") # Example 3: Define pair style (non-bonded interactions) pair_style = ff.def_pairstyle("lj126/cut") pair_style.def_type(ct, epsilon=0.066, sigma=3.50) pair_style.def_type(hc, epsilon=0.030, sigma=2.50) print(f"\nDefined {len(ff.get_types(mp.PairType))} pair types:") for ptype in ff.get_types(mp.PairType): epsilon = ptype["epsilon"] sigma = ptype["sigma"] print(f" {ptype.name}: ε={epsilon:.3f} kcal/mol, σ={sigma:.2f} Å")

Defined 3 angle types:
  HC-CT-HC: k=33.00 kcal/mol/rad², θ0=107.8°
  CT-CT-CT: k=58.35 kcal/mol/rad², θ0=112.7°
  HC-CT-CT: k=37.50 kcal/mol/rad², θ0=110.7°

Defined 2 pair types:
  HC: ε=0.030 kcal/mol, σ=2.50 Å
  CT: ε=0.066 kcal/mol, σ=3.50 Å

5. Converting Style to Potential¶

Key concept: A Style contains parameter definitions, types, while a Potential is an executable energy function.

Requirements: The style must have a to_potential() method, The style name must match a registered Potential class, and All required parameters must be defined in the types.

Let's convert our bond style to a potential:

In [5]:

# Convert bond style to potential
bond_potential = bond_style.to_potential()

print(f"Bond style converted to potential:")
print(f"  Style: {bond_style}")
print(f"  Potential: {bond_potential}")
print(f"  Potential type: {type(bond_potential).__name__}")
print(f"  Potential has k: {hasattr(bond_potential, 'k')}")
print(f"  Potential has r0: {hasattr(bond_potential, 'r0')}")
print(f"  Number of bond types: {len(bond_potential.k)}")

# Convert angle style to potential
angle_potential = angle_style.to_potential()
print(f"\nAngle style converted to potential:")
print(f"  Potential: {angle_potential}")
print(f"  Potential type: {type(angle_potential).__name__}")

# Convert pair style to potential
pair_potential = pair_style.to_potential()
print(f"\nPair style converted to potential:")
print(f"  Potential: {pair_potential}")
print(f"  Potential type: {type(pair_potential).__name__}")

# Convert bond style to potential bond_potential = bond_style.to_potential() print(f"Bond style converted to potential:") print(f" Style: {bond_style}") print(f" Potential: {bond_potential}") print(f" Potential type: {type(bond_potential).__name__}") print(f" Potential has k: {hasattr(bond_potential, 'k')}") print(f" Potential has r0: {hasattr(bond_potential, 'r0')}") print(f" Number of bond types: {len(bond_potential.k)}") # Convert angle style to potential angle_potential = angle_style.to_potential() print(f"\nAngle style converted to potential:") print(f" Potential: {angle_potential}") print(f" Potential type: {type(angle_potential).__name__}") # Convert pair style to potential pair_potential = pair_style.to_potential() print(f"\nPair style converted to potential:") print(f" Potential: {pair_potential}") print(f" Potential type: {type(pair_potential).__name__}")

Bond style converted to potential:
  Style: <BondStyle: harmonic>
  Potential: <molpy.potential.bond.harmonic.BondHarmonic object at 0x7fe102a9aba0>
  Potential type: BondHarmonic
  Potential has k: True
  Potential has r0: True
  Number of bond types: 2

Angle style converted to potential:
  Potential: <molpy.potential.angle.harmonic.AngleHarmonic object at 0x7fe102a9ae40>
  Potential type: AngleHarmonic

Pair style converted to potential:
  Potential: <molpy.potential.pair.lj.LJ126 object at 0x7fe102a9af90>
  Potential type: LJ126

6. Converting ForceField to Potentials¶

Instead of converting each style individually, you can convert the entire ForceField to a Potentials collection:

Method:

potentials = ff.to_potentials()

This automatically: Finds all styles with to_potential() method, Converts each style to its corresponding Potential, and Returns a Potentials collection, list-like container.

Benefits: Single call converts all compatible styles, Handles errors gracefully, skips styles without registered potentials, and Returns a collection ready for energy calculations.

In [6]:

# Convert entire force field to potentials
potentials = ff.to_potentials()

print(f"Force Field: {ff.name}")
print(f"=" * 50)
print(f"Atom types:     {len(ff.get_types(mp.AtomType))}")
print(f"Bond types:     {len(ff.get_types(mp.BondType))}")
print(f"Angle types:    {len(ff.get_types(mp.AngleType))}")
print(f"Pair types:     {len(ff.get_types(mp.PairType))}")

# Show all styles
from molpy.core.forcefield import Style
all_styles = ff.styles.bucket(Style)
print(f"\nStyles defined: {[s.name for s in all_styles]}")

# Show converted potentials
print(f"\nPotentials created: {len(potentials)}")
for i, pot in enumerate(potentials):
    print(f"  {i+1}. {type(pot).__name__}: {pot}")

# Convert entire force field to potentials potentials = ff.to_potentials() print(f"Force Field: {ff.name}") print(f"=" * 50) print(f"Atom types: {len(ff.get_types(mp.AtomType))}") print(f"Bond types: {len(ff.get_types(mp.BondType))}") print(f"Angle types: {len(ff.get_types(mp.AngleType))}") print(f"Pair types: {len(ff.get_types(mp.PairType))}") # Show all styles from molpy.core.forcefield import Style all_styles = ff.styles.bucket(Style) print(f"\nStyles defined: {[s.name for s in all_styles]}") # Show converted potentials print(f"\nPotentials created: {len(potentials)}") for i, pot in enumerate(potentials): print(f" {i+1}. {type(pot).__name__}: {pot}")

Force Field: MyAtomisticFF
==================================================
Atom types:     2
Bond types:     2
Angle types:    3
Pair types:     2

Styles defined: ['full', 'harmonic', 'harmonic', 'lj126/cut']

Potentials created: 3
  1. BondHarmonic: <molpy.potential.bond.harmonic.BondHarmonic object at 0x7fe137bcfd90>
  2. AngleHarmonic: <molpy.potential.angle.harmonic.AngleHarmonic object at 0x7fe137bcfed0>
  3. LJ126: <molpy.potential.pair.lj.LJ126 object at 0x7fe1022d8050>

7. Defining Dihedral Parameters¶

Let's add a dihedral style to complete our force field:

OPLS style: $$E_{dihedral} = \frac{1}{2}[K_1, 1+\cos\phi + K_2, 1-\cos 2\phi + K_3, 1+\cos 3\phi + K_4, 1-\cos 4\phi]$$

Where: $\phi$ = dihedral angle and $K_i$ = force constants.

In [7]:

# Define dihedral style (opls for OPLS-style multi-term)
dihedral_style = ff.def_dihedralstyle("opls")

# Use the style's def_type() method to create dihedral types for C-C-C-C rotation
ct_ct_ct_ct = dihedral_style.def_type(
    ct, ct, ct, ct,
    K1=1.3, K2=-0.05, K3=0.2, K4=0.0
)

print(f"Defined {len(ff.get_types(mp.DihedralType))} dihedral types:")
for dtype in ff.get_types(mp.DihedralType):
    K1 = dtype["K1"]
    K2 = dtype["K2"]
    K3 = dtype["K3"]
    print(f"  {dtype.name}: K1={K1:.2f}, K2={K2:.2f}, K3={K3:.2f}")

# Define dihedral style (opls for OPLS-style multi-term) dihedral_style = ff.def_dihedralstyle("opls") # Use the style's def_type() method to create dihedral types for C-C-C-C rotation ct_ct_ct_ct = dihedral_style.def_type( ct, ct, ct, ct, K1=1.3, K2=-0.05, K3=0.2, K4=0.0 ) print(f"Defined {len(ff.get_types(mp.DihedralType))} dihedral types:") for dtype in ff.get_types(mp.DihedralType): K1 = dtype["K1"] K2 = dtype["K2"] K3 = dtype["K3"] print(f" {dtype.name}: K1={K1:.2f}, K2={K2:.2f}, K3={K3:.2f}")

Defined 1 dihedral types:
  CT-CT-CT-CT: K1=1.30, K2=-0.05, K3=0.20

8. Loading OPLS-AA Force Field¶

MolPy can load standard force fields from XML files. The built-in OPLS-AA force field is available:

Method:

from molpy.io.forcefield.xml import read_xml_forcefield, read_oplsaa_forcefield

# Load OPLS-AA
ff_opls = read_xml_forcefield, "oplsaa.xml"

# Or explicitly use OPLS-AA reader
ff_opls = read_oplsaa_forcefield, "oplsaa.xml"

What gets loaded: Atom types with mass, charge, and LJ parameters, Bond parameters, harmonic, Angle parameters, harmonic, Dihedral parameters, OPLS style, and Non-bonded parameters, LJ and Coulomb.

Unit conversions: OPLS-AA XML uses kJ/mol and nm and read_oplsaa_forcefield() converts to kcal/mol and Å for LAMMPS compatibility.

In [8]:

# Load OPLS-AA force field
from molpy.io.forcefield.xml import read_xml_forcefield

try:
    # Load built-in OPLS-AA (filename only - searches in molpy/data/forcefield/)
    ff_opls = read_xml_forcefield("oplsaa.xml")
    
    print(f"Loaded OPLS-AA force field: {ff_opls.name}")
    print(f"  Units: {ff_opls.units}")
    print(f"  Atom types: {len(ff_opls.get_types(mp.AtomType))}")
    print(f"  Bond types: {len(ff_opls.get_types(mp.BondType))}")
    print(f"  Angle types: {len(ff_opls.get_types(mp.AngleType))}")
    print(f"  Dihedral types: {len(ff_opls.get_types(mp.DihedralType))}")
    print(f"  Pair types: {len(ff_opls.get_types(mp.PairType))}")
    
    # Show some example atom types
    print(f"\nExample atom types:")
    for atype in list(ff_opls.get_types(mp.AtomType))[:5]:
        print(f"  {atype.name}: mass={atype.get('mass', 'N/A')}, charge={atype.get('charge', 'N/A')}")
    
    # Convert to potentials
    potentials_opls = ff_opls.to_potentials()
    print(f"\nConverted to {len(potentials_opls)} potentials")
    
except FileNotFoundError as e:
    print(f"OPLS-AA file not found: {e}")
    print("\nTo load a custom XML file:")
    print("  from pathlib import Path")
    print("  ff = read_xml_forcefield(Path('/path/to/forcefield.xml'))")

# Load OPLS-AA force field from molpy.io.forcefield.xml import read_xml_forcefield try: # Load built-in OPLS-AA (filename only - searches in molpy/data/forcefield/) ff_opls = read_xml_forcefield("oplsaa.xml") print(f"Loaded OPLS-AA force field: {ff_opls.name}") print(f" Units: {ff_opls.units}") print(f" Atom types: {len(ff_opls.get_types(mp.AtomType))}") print(f" Bond types: {len(ff_opls.get_types(mp.BondType))}") print(f" Angle types: {len(ff_opls.get_types(mp.AngleType))}") print(f" Dihedral types: {len(ff_opls.get_types(mp.DihedralType))}") print(f" Pair types: {len(ff_opls.get_types(mp.PairType))}") # Show some example atom types print(f"\nExample atom types:") for atype in list(ff_opls.get_types(mp.AtomType))[:5]: print(f" {atype.name}: mass={atype.get('mass', 'N/A')}, charge={atype.get('charge', 'N/A')}") # Convert to potentials potentials_opls = ff_opls.to_potentials() print(f"\nConverted to {len(potentials_opls)} potentials") except FileNotFoundError as e: print(f"OPLS-AA file not found: {e}") print("\nTo load a custom XML file:") print(" from pathlib import Path") print(" ff = read_xml_forcefield(Path('/path/to/forcefield.xml'))")

2025-12-30 16:08:16,332 - molpy.potential.dihedral.opls - WARNING - RB coefficients do not lie on the ideal 4-term OPLS manifold (C0+C1+C2+C3+C4 = 10.041600, expected ≈ 0). Conversion will preserve forces and relative energies exactly, but will introduce a constant energy offset of ΔE = 10.041600 kJ/mol. This does not affect MD simulations.

Loaded OPLS-AA force field: OPLS-AA
  Units: real
  Atom types: 892
  Bond types: 272
  Angle types: 838
  Dihedral types: 888
  Pair types: 825

Example atom types:
  opls_473: mass=15.9994, charge=N/A
  H4: mass=N/A, charge=N/A
  opls_648: mass=12.011, charge=N/A
  opls_066: mass=16.043, charge=N/A
  opls_325: mass=1.008, charge=N/A

Converted to 3 potentials

9. Defining Custom Potential, Style, and Type¶

MolPy's force field system is extensible - you can define custom Potential, Style, and Type classes for novel interactions.

When to extend: Implementing new potential forms, e.g., polarizable models, reactive FFs, Custom bonded interactions, e.g., cross-terms, CMAP, Special constraints, e.g., virtual sites, Drude oscillators, and Coarse-grained force fields with custom bead types.

Extension pattern: Define custom Type class, inherits from Type or specialized type, Define custom Style class, inherits from Style or specialized style, Define custom Potential class, inherits from Potential, uses KernelMeta for auto-registration, and Implement to_potential() method in Style to create Potential instances.

In [9]:

# Example: Define a custom angle-bond cross-term (Class II force field style)
import numpy as np
from numpy.typing import NDArray
from molpy.core.forcefield import Type, Style, AtomType
from molpy.potential.base import Potential, KernelMeta

# Step 1: Define custom Type
class AngleRadialType(Type):
    """Custom type for angle-bond cross-term (like in Class II FFs)"""
    def __init__(self, name: str, itom: AtomType, jtom: AtomType, ktom: AtomType, **kwargs):
        super().__init__(name, **kwargs)
        self.itom = itom
        self.jtom = jtom  # Central atom
        self.ktom = ktom
        
    def __repr__(self):
        return f"<AngleRadialType: {self.itom.name}-{self.jtom.name}-{self.ktom.name}>"

# Step 2: Define custom Style
class AngleRadialStyle(Style):
    """Style for angle-bond cross-terms"""
    def def_type(self, itom: AtomType, jtom: AtomType, ktom: AtomType, name: str = "", **kwargs):
        """Define angle-radial coupling type
        
        Args:
            itom: First atom type
            jtom: Central atom type
            ktom: Third atom type
            name: Optional name (defaults to itom-jtom-ktom)
            **kwargs: Parameters (e.g., k_theta_r for coupling constant)
        """
        if not name:
            name = f"{itom.name}-{jtom.name}-{ktom.name}"
        art = AngleRadialType(name, itom, jtom, ktom, **kwargs)
        self.types.add(art)
        return art
    
    def to_potential(self):
        """Convert style to potential (optional - for energy calculations)"""
        # This is optional - only needed if you want to use the potential for calculations
        # For now, we'll just return None to show the pattern
        return None

# Step 3: Define custom Potential (optional - for energy calculations)
class AngleRadialPotential(Potential, metaclass=KernelMeta):
    """Custom potential for angle-bond cross-terms"""
    name = "angle_radial"
    type = "angle_radial"  # Registers in ForceField._kernel_registry["angle_radial"]
    
    def __init__(self, k_theta_r: NDArray[np.floating] | float):
        """Initialize angle-radial potential.
        
        Args:
            k_theta_r: Coupling constant
        """
        self.k_theta_r = np.array(k_theta_r, dtype=np.float64)
    
    def calc_energy(self, *args, **kwargs) -> float:
        """Calculate energy (implementation depends on your potential form)"""
        # Placeholder - implement your energy calculation here
        return 0.0
    
    def calc_forces(self, *args, **kwargs) -> NDArray:
        """Calculate forces (implementation depends on your potential form)"""
        # Placeholder - implement your force calculation here
        return np.zeros((1, 3))

# Step 4: Use in ForceField
angle_radial_style = ff.def_style(AngleRadialStyle("angle_radial"))
custom_term = angle_radial_style.def_type(ct, ct, ct, k_theta_r=5.0)

print(f"Custom type added: {custom_term}")
print(f"  Parameter k_theta_r: {custom_term['k_theta_r']}")
print(f"\nForce field styles: {list(ff.styles.bucket(Style))}")

# Verify potential registration
print(f"\nRegistered potentials in registry:")
from molpy.core.forcefield import ForceField
if "angle_radial" in ForceField._kernel_registry:
    print(f"  angle_radial: {list(ForceField._kernel_registry['angle_radial'].keys())}")

print(f"\nThis demonstrates extensibility for:")
print("  - Class II force fields (COMPASS, PCFF)")
print("  - Coarse-grained models (Martini, SDK)")
print("  - Reactive force fields (ReaxFF-like terms)")
print("  - Polarizable models (Drude, AMOEBA)")
print("  - Machine learning potentials (custom descriptors)")

# Example: Define a custom angle-bond cross-term (Class II force field style) import numpy as np from numpy.typing import NDArray from molpy.core.forcefield import Type, Style, AtomType from molpy.potential.base import Potential, KernelMeta # Step 1: Define custom Type class AngleRadialType(Type): """Custom type for angle-bond cross-term (like in Class II FFs)""" def __init__(self, name: str, itom: AtomType, jtom: AtomType, ktom: AtomType, **kwargs): super().__init__(name, **kwargs) self.itom = itom self.jtom = jtom # Central atom self.ktom = ktom def __repr__(self): return f"" # Step 2: Define custom Style class AngleRadialStyle(Style): """Style for angle-bond cross-terms""" def def_type(self, itom: AtomType, jtom: AtomType, ktom: AtomType, name: str = "", **kwargs): """Define angle-radial coupling type Args: itom: First atom type jtom: Central atom type ktom: Third atom type name: Optional name (defaults to itom-jtom-ktom) **kwargs: Parameters (e.g., k_theta_r for coupling constant) """ if not name: name = f"{itom.name}-{jtom.name}-{ktom.name}" art = AngleRadialType(name, itom, jtom, ktom, **kwargs) self.types.add(art) return art def to_potential(self): """Convert style to potential (optional - for energy calculations)""" # This is optional - only needed if you want to use the potential for calculations # For now, we'll just return None to show the pattern return None # Step 3: Define custom Potential (optional - for energy calculations) class AngleRadialPotential(Potential, metaclass=KernelMeta): """Custom potential for angle-bond cross-terms""" name = "angle_radial" type = "angle_radial" # Registers in ForceField._kernel_registry["angle_radial"] def __init__(self, k_theta_r: NDArray[np.floating] | float): """Initialize angle-radial potential. Args: k_theta_r: Coupling constant """ self.k_theta_r = np.array(k_theta_r, dtype=np.float64) def calc_energy(self, *args, **kwargs) -> float: """Calculate energy (implementation depends on your potential form)""" # Placeholder - implement your energy calculation here return 0.0 def calc_forces(self, *args, **kwargs) -> NDArray: """Calculate forces (implementation depends on your potential form)""" # Placeholder - implement your force calculation here return np.zeros((1, 3)) # Step 4: Use in ForceField angle_radial_style = ff.def_style(AngleRadialStyle("angle_radial")) custom_term = angle_radial_style.def_type(ct, ct, ct, k_theta_r=5.0) print(f"Custom type added: {custom_term}") print(f" Parameter k_theta_r: {custom_term['k_theta_r']}") print(f"\nForce field styles: {list(ff.styles.bucket(Style))}") # Verify potential registration print(f"\nRegistered potentials in registry:") from molpy.core.forcefield import ForceField if "angle_radial" in ForceField._kernel_registry: print(f" angle_radial: {list(ForceField._kernel_registry['angle_radial'].keys())}") print(f"\nThis demonstrates extensibility for:") print(" - Class II force fields (COMPASS, PCFF)") print(" - Coarse-grained models (Martini, SDK)") print(" - Reactive force fields (ReaxFF-like terms)") print(" - Polarizable models (Drude, AMOEBA)") print(" - Machine learning potentials (custom descriptors)")

Custom type added: <AngleRadialType: CT-CT-CT>
  Parameter k_theta_r: 5.0

Force field styles: [<AtomStyle: full>, <BondStyle: harmonic>, <AngleStyle: harmonic>, <PairStyle: lj126/cut>, <DihedralStyle: opls>, <AngleRadialStyle: angle_radial>]

Registered potentials in registry:
  angle_radial: ['angle_radial']

This demonstrates extensibility for:
  - Class II force fields (COMPASS, PCFF)
  - Coarse-grained models (Martini, SDK)
  - Reactive force fields (ReaxFF-like terms)
  - Polarizable models (Drude, AMOEBA)
  - Machine learning potentials (custom descriptors)

10. Summary and Best Practices¶

Let's inspect what we've created and review key concepts:

In [10]:

print(f"Force Field Summary: {ff.name}")
print("=" * 60)
print(f"Atom types:     {len(ff.get_types(mp.AtomType))}")
print(f"Bond types:     {len(ff.get_types(mp.BondType))}")
print(f"Angle types:    {len(ff.get_types(mp.AngleType))}")
print(f"Dihedral types: {len(ff.get_types(mp.DihedralType))}")
print(f"Pair types:     {len(ff.get_types(mp.PairType))}")

# Show all styles
all_styles = ff.styles.bucket(Style)
print(f"\nStyles defined ({len(all_styles)}):")
for style in all_styles:
    n_types = len(style.types.bucket(Type))
    print(f"  {style.__class__.__name__}({style.name}): {n_types} types")

# Show potentials
potentials = ff.to_potentials()
print(f"\nPotentials created: {len(potentials)}")
for pot in potentials:
    print(f"  {type(pot).__name__}")

print(f"Force Field Summary: {ff.name}") print("=" * 60) print(f"Atom types: {len(ff.get_types(mp.AtomType))}") print(f"Bond types: {len(ff.get_types(mp.BondType))}") print(f"Angle types: {len(ff.get_types(mp.AngleType))}") print(f"Dihedral types: {len(ff.get_types(mp.DihedralType))}") print(f"Pair types: {len(ff.get_types(mp.PairType))}") # Show all styles all_styles = ff.styles.bucket(Style) print(f"\nStyles defined ({len(all_styles)}):") for style in all_styles: n_types = len(style.types.bucket(Type)) print(f" {style.__class__.__name__}({style.name}): {n_types} types") # Show potentials potentials = ff.to_potentials() print(f"\nPotentials created: {len(potentials)}") for pot in potentials: print(f" {type(pot).__name__}")

Force Field Summary: MyAtomisticFF
============================================================
Atom types:     2
Bond types:     2
Angle types:    3
Dihedral types: 1
Pair types:     2

Styles defined (6):
  AtomStyle(full): 2 types
  BondStyle(harmonic): 2 types
  AngleStyle(harmonic): 3 types
  PairStyle(lj126/cut): 2 types
  DihedralStyle(opls): 1 types
  AngleRadialStyle(angle_radial): 1 types

Potentials created: 3
  BondHarmonic
  AngleHarmonic
  LJ126

Key takeaways¶

In MolPy, a ForceField represents parameter data, not executable code. It is a structured container that organizes interaction styles together with their typed parameter tables. A force field defines what parameters exist and how they are grouped, but it does not perform numerical evaluation by itself.

A Style represents a single interaction family within the force field, such as bonds, angles, dihedrals, or pairwise interactions. Each style specifies a functional form—for example, a harmonic bond or a Lennard–Jones pair interaction—and acts as a namespace for the parameter types that belong to that form. Styles are defined explicitly so that one functional form can be reused consistently across many interactions.

A Type is one concrete parameter set inside a style. For example, within a harmonic bond style, different bond types correspond to different equilibrium lengths and force constants, such as a specific C–H bond or C–C bond. Types are purely declarative: they store numbers and identifiers, but they do not execute any computation.

A Potential is the executable numerical kernel derived from a style, or from the entire force field. Converting a style or force field into a potential compiles the stored parameter data into a form that can be evaluated efficiently during simulation. This separation ensures that force-field definition remains declarative and editable, while numerical execution is handled only at the final stage.

The minimal workflow reflects this separation clearly:

# create a force field as a parameter container
ff = mp.AtomisticForcefield(name="MyFF", units="real")

# define interaction styles and their parameter spaces
atom_style = ff.def_atomstyle("full")
bond_style = ff.def_bondstyle("harmonic")
angle_style = ff.def_anglestyle("harmonic")
pair_style = ff.def_pairstyle("lj126/cut")

# convert parameter definitions into executable potentials
bond_pot = bond_style.to_potential()
pots = ff.to_potentials()

In this design, force fields describe what interactions exist, styles define how those interactions are modeled, types store the actual numerical parameters, and potentials are the point where everything becomes executable.

Next steps¶

Typing workflows, including typing rules and their application, are documented in typing. Molecular systems are first built or edited at the Atomistic graph level, from which topology is derived automatically. The structure is then typified and associated with a ForceField, after which it is exported through the Frame/I/O layer into simulator-ready formats. Simulation and downstream analysis operate on these exported representations, completing the pipeline:

Build Atomistic → Derive Topology → Typify & ForceField → Export → Simulate → Analyze

This separation keeps structural editing, parameter definition, numerical execution, and analysis cleanly decoupled while allowing them to integrate into a coherent end-to-end workflow.