diff --git a/src/nomad_simulations/schema_packages/force_field.py b/src/nomad_simulations/schema_packages/force_field.py
new file mode 100644
index 00000000..d638cbcf
--- /dev/null
+++ b/src/nomad_simulations/schema_packages/force_field.py
@@ -0,0 +1,955 @@
+import numpy as np
+import pint
+import re
+import typing
+from scipy.interpolate import UnivariateSpline
+
+# from structlog.stdlib import BoundLogger
+from typing import Optional, List, Tuple
+from ase.dft.kpoints import monkhorst_pack, get_monkhorst_pack_size_and_offset
+
+from nomad.units import ureg
+from nomad.datamodel.data import ArchiveSection
+from nomad.datamodel.metainfo.annotations import ELNAnnotation
+from nomad.metainfo import (
+    URL,
+    Quantity,
+    SubSection,
+    MEnum,
+    Section,
+    Context,
+    JSON,
+)
+
+from nomad_simulations.schema_packages.model_system import ModelSystem
+from nomad_simulations.schema_packages.model_method import BaseModelMethod, ModelMethod
+from nomad_simulations.schema_packages.utils import is_not_representative
+
+MOL = 6.022140857e23
+
+
+class ParameterEntry(ArchiveSection):
+    """
+    Generic section defining a parameter name and value
+    """
+
+    name = Quantity(
+        type=str,
+        shape=[],
+        description="""
+        Name of the parameter.
+        """,
+    )
+
+    value = Quantity(
+        type=str,
+        shape=[],
+        description="""
+        Value of the parameter as a string.
+        """,
+    )
+
+    unit = Quantity(
+        type=str,
+        shape=[],
+        description="""
+        Unit of the parameter as a string.
+        """,
+    )
+
+
+#     # TODO add description quantity
+
+
+class Potential(BaseModelMethod):
+    """
+    Section containing information about an interaction potential.
+
+        name: str - potential name, can be as specific as needed
+        type: str - potential type, e.g., 'bond', 'angle', 'dihedral', 'improper dihedral', 'nonbonded'
+        functional_form: str - functional form of the potential, e.g., 'harmonic', 'Morse', 'Lennard-Jones'
+        external_reference: URL
+    """
+
+    parameters = SubSection(
+        sub_section=ParameterEntry.m_def,
+        repeats=True,
+        description="""
+        List of parameters for custom potentials.
+        """,
+    )
+
+    type = Quantity(
+        type=MEnum('bond', 'angle', 'dihedral', 'improper dihedral', 'nonbonded'),
+        shape=[],
+        description="""
+        Denotes the classification of the interaction.
+        """,
+    )
+
+    functional_form = Quantity(
+        type=str,
+        shape=[],
+        description="""
+        Specifies the functional form of the interaction potential, e.g., harmonic, Morse, Lennard-Jones, etc.
+        """,
+    )
+
+    n_interactions = Quantity(
+        type=np.dtype(np.int32),
+        shape=[],
+        description="""
+        Total number of interactions in the system for this potential.
+        """,
+    )
+
+    n_particles = Quantity(
+        type=np.int32,
+        shape=[],
+        description="""
+        Number of particles interacting via (each instance of) this potential.
+        """,
+    )
+
+    particle_labels = Quantity(
+        type=np.dtype(str),
+        shape=['n_interactions', 'n_particles'],
+        description="""
+        Labels of the particles for each instance of this potential, stored as a list of tuples.
+        """,
+    )
+
+    particle_indices = Quantity(
+        type=np.int32,
+        shape=['n_interactions', 'n_particles'],
+        description="""
+        Indices of the particles for each instance of this potential, stored as a list of tuples.
+        """,
+    )
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        # set the dimensions based on the particle indices, if stored
+        if not self.n_interactions:
+            self.n_interactions = (
+                len(self.particle_indices)
+                if self.particle_indices is not None
+                else None
+            )
+        if not self.n_particles:
+            self.n_interactions = (
+                len(self.particle_indices[0])
+                if self.particle_indices is not None
+                else None
+            )
+
+        # check the consistency of the dimensions of the particle indices and labels
+        if self.n_interactions and self.n_particles:
+            if self.particle_indices is not None:
+                assert len(self.particle_indices) == self.n_interactions
+                assert len(self.particle_indices[0]) == self.n_particles
+            if self.particle_labels is not None:
+                assert len(self.particle_labels) == self.n_interactions
+                assert len(self.particle_labels[0]) == self.n_particles
+
+
+class TabulatedPotential(Potential):
+    """
+    Abstract class for tabulated potentials. The value of the potential and/or force
+    is stored for a set of corresponding bin distances. The units for bins and forces
+    should be set in the individual subclasses.
+    """
+
+    bins = Quantity(
+        type=np.float64,
+        shape=[],
+        description="""
+        List of bin angles.
+        """,
+    )
+
+    energies = Quantity(
+        type=np.float64,
+        unit='J',
+        shape=[],
+        description="""
+        List of energy values associated with each bin.
+        """,
+    )
+
+    forces = Quantity(
+        type=np.float64,
+        shape=[],
+        description="""
+        List of force values associated with each bin.
+        """,
+    )
+
+    def compute_forces(self, bins, energies, smoothing_factor=None):
+        if smoothing_factor is None:
+            smoothing_factor = len(bins) - np.sqrt(2 * len(bins))
+
+        spline = UnivariateSpline(bins, energies, s=smoothing_factor)
+        forces = -1.0 * spline.derivative()(bins)
+
+        return forces
+
+    def compute_energies(self, bins, forces, smoothing_factor=None):
+        if smoothing_factor is None:
+            smoothing_factor = len(bins) - np.sqrt(2 * len(bins))
+
+        spline = UnivariateSpline(bins, forces, s=smoothing_factor)
+        energies = -1.0 * np.array([spline.integral(bins[0], x) for x in bins])
+        energies -= np.min(energies)
+
+        return energies
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        self.name = 'TabulatedPotential'
+        if not self.functional_form:
+            self.functional_form = 'tabulated'
+        elif self.functional_form != 'tabulated':
+            logger.warning(f'Incorrect functional form set for {self.name}.')
+
+        if self.bins is not None and self.energies is not None:
+            if len(self.bins) != len(self.energies):
+                logger.error(
+                    f'bins and energies values have different length in {self.name}'
+                )
+        if self.bins is not None and self.forces is not None:
+            if len(self.bins) != len(self.forces):
+                logger.error(f'bins and forces have different length in {self.name}')
+
+        if self.bins is not None:
+            smoothing_factor = len(self.bins) - np.sqrt(2 * len(self.bins))
+            tol = 1e-2
+            if self.forces is None and self.energies is not None:
+                try:
+                    # generate forces from energies numerically using spline
+                    self.forces = (
+                        self.compute_forces(
+                            self.bins.magnitude,
+                            self.energies.magnitude,
+                            smoothing_factor=smoothing_factor,
+                        )  # ? How do I deal with units here?
+                    )
+                    # re-derive energies to check consistency of the forces
+                    energies = (
+                        self.compute_energies(
+                            self.bins.magnitude,
+                            self.forces.magnitude,
+                            smoothing_factor=smoothing_factor,
+                        )
+                        * ureg.J
+                    )
+
+                    energies_diff = energies.to('kJ').magnitude * MOL - (
+                        self.energies.to('kJ').magnitude * MOL
+                        - np.min(self.energies.to('kJ').magnitude * MOL)
+                    )
+                    if np.all([x < tol for x in energies_diff]):
+                        logger.warning(
+                            f'Tabulated forces were generated from energies in {self.name},'
+                            f'with consistency errors less than tol={tol}. '
+                        )
+                    else:
+                        logger.warning(
+                            f'Unable to derive tabulated forces from energies in {self.name},'
+                            f'consistency errors were greater than tol={tol}.'
+                        )
+                        self.forces = None
+                except ValueError as e:
+                    logger.warning(
+                        f'Unable to derive tabulated forces from energies in {self.name},'
+                        f'Unkown error occurred in derivation: {e}'
+                    )
+
+            if self.forces is not None and self.energies is None:
+                try:
+                    # generated energies from forces numerically using spline
+                    self.energies = self.compute_energies(
+                        self.bins.magnitude,
+                        self.forces.magnitude,
+                        smoothing_factor=smoothing_factor,
+                    )
+                    # re-derive forces to check consistency of the energies
+                    forces = (
+                        self.compute_forces(
+                            self.bins.magnitude,
+                            self.energies.magnitude,
+                            smoothing_factor=smoothing_factor,
+                        )
+                        * ureg.J
+                        / self.bins.units
+                    )
+
+                    forces_diff = forces.to(f'kJ/{self.bins.units}').magnitude * MOL - (
+                        self.forces.to(f'kJ/{self.bins.units}').magnitude * MOL
+                    )
+                    if np.all([x < tol for x in forces_diff]):
+                        logger.warning(
+                            f'Tabulated energies were generated from forces in {self.name},'
+                            f'with consistency errors less than tol={tol}. '
+                        )
+                    else:
+                        logger.warning(
+                            f'Unable to derive tabulated energies from forces in {self.name},'
+                            f'consistency errors were greater than tol={tol}.'
+                        )
+                        self.energies = None
+                except ValueError as e:
+                    logger.warning(
+                        f'Unable to derive tabulated energies from forces in {self.name},'
+                        f'Unkown error occurred in derivation: {e}'
+                    )
+
+
+class PolynomialForceConstant(ParameterEntry):
+    """
+    Section defining a force constant for a potential of polynomial form.
+    """
+
+    name = Quantity(
+        type=str,
+        shape=[],
+        description="""
+        Name of the force constant.
+        """,
+    )
+
+    exponent = Quantity(
+        type=np.int32,
+        shape=[],
+        description="""
+        Exponent for this term in the polynomial.
+        """,
+    )
+
+    value = Quantity(
+        type=np.float64,
+        shape=[],
+        description="""
+        Value of the force constant.
+        """,
+    )
+
+
+class PolynomialPotential(Potential):
+    r"""
+    Abstract class for potentials with polynomial form:
+    $V(x) = [\left k_1 (x - x_0) + k_2 (x - x_0)^2 + x_3 (x - x_0)^3 + \dots + C$,
+    where $\{x_1, x_2, x_3 \dots}$ are the `force_constants` for each term in the polynomial
+    expression and $x_0$ is the `equilibrium_value` of $x$. $C$ is an arbitrary constant (not stored).
+    This class is intended to be used with another class specifying the potential type, e.g., BondPotential, AnglePotential, etc.
+    """
+
+    force_constants = SubSection(
+        sub_section=PolynomialForceConstant.m_def,
+        repeats=True,
+        description="""
+        List of force constants value and corresponding unit for polynomial potentials.
+        """,
+    )
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        self.name = 'PolynomialPotential'
+        if not self.functional_form:
+            self.functional_form = 'polynomial'
+        elif self.functional_form != 'polynomial':
+            logger.warning('Incorrect functional form set for PolynomialPotential.')
+
+
+class BondPotential(Potential):
+    """
+    Section containing information about bond potentials.
+
+    Suggested types are: harmonic, cubic, polynomial, Morse, fene, tabulated
+    """
+
+    equilibrium_value = Quantity(
+        type=np.float64,
+        unit='m',
+        shape=[],
+        description="""
+        Specifies the equilibrium bond distance.
+        """,
+    )
+
+    force_constant = Quantity(
+        type=np.float64,
+        shape=[],
+        unit='J / m **2',
+        description="""
+        Specifies the force constant of the bond potential.
+        """,
+    )
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        if not self.name:
+            self.name = 'BondPotential'
+        if not self.type:
+            self.type = 'bond'
+        elif self.type != 'bond':
+            logger.warning('Incorrect type set for BondPotential.')
+
+        if self.n_particles:
+            if self.n_particles != 2:
+                logger.warning('Incorrect number of particles set for BondPotential.')
+            else:
+                self.n_particles = 2
+
+
+class HarmonicBond(BondPotential):
+    r"""
+    Section containing information about a Harmonic bond potential:
+    $V(r) = \frac{1}{2} k_r (r - r_0)^2 + C$,
+    where $k_r$ is the `force_constant` and $r_0$ is the `equilibrium_value` of $r$.
+    $C$ is an arbitrary constant (not stored).
+    """
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        self.name = 'HarmonicBond'
+        if not self.functional_form:
+            self.functional_form = 'harmonic'
+        elif self.functional_form != 'harmonic':
+            logger.warning('Incorrect functional form set for HarmonicBond.')
+
+
+class CubicBond(BondPotential):
+    r"""
+    Section containing information about a Cubic bond potential:
+    $V(r) = \frac{1}{2} k_r (r - r_0)^2 + \frac{1}{3} k_c (r - r_0)^3 + C$,
+    where $k_r$ is the (harmonic) `force_constant`, $k_c$ is the `force_constant_cubic`,
+    and $r_0$ is the `equilibrium_value` of $r$.
+    C is an arbitrary constant (not stored).
+    """
+
+    force_constant_cubic = Quantity(
+        type=np.float64,
+        shape=[],
+        unit='J / m**3',
+        description="""
+        Specifies the cubic force constant of the bond potential.
+        """,
+    )
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        self.name = 'CubicBond'
+        if not self.functional_form:
+            self.functional_form = 'cubic'
+        elif self.functional_form != 'cubic':
+            logger.warning('Incorrect functional form set for CubicBond.')
+
+
+class PolynomialBond(PolynomialPotential, BondPotential):
+    """
+    Section containing information about a polynomial bond potential:
+    """
+
+    def __init__(self):
+        super().__init__()
+        docstring = PolynomialPotential.__doc__
+        pattern = r'\$V\(x\)(.*?)(\(not stored\)\.)'
+        match = re.search(pattern, docstring, re.DOTALL)
+        extracted_text = '<functional form missing>'
+        if match:
+            extracted_text = match.group().strip()  # .group(1).strip()
+        self.__doc__ = rf"""{self.__doc__} {extracted_text}.
+            Here the dependent variable of the potential, $x$, corresponds to the bond distance."""
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        self.name = 'PolynomialBond'
+
+
+class MorseBond(BondPotential):
+    r"""
+    Section containing information about a Morse potential:
+    $V(r) = D \left[ 1 - e^{-a(r - r_0)} \right]^2 + C$,
+    where $a = sqrt(k/2D)$ is the `well_steepness`, with `force constant` k.
+    D is the `well_depth`, and r_0 is the `equilibrium_value` of a.
+    C is an arbitrary constant (not stored).
+    """
+
+    well_depth = Quantity(
+        type=np.float64,
+        unit='J',
+        shape=[],
+        description="""
+        Specifies the depth of the potential well.
+        """,
+    )
+
+    well_steepness = Quantity(
+        type=np.float64,
+        unit='1/m',
+        shape=[],
+        description="""
+        Specifies the steepness of the potential well.
+        """,
+    )
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        self.name = 'MorseBond'
+        if not self.functional_form:
+            self.functional_form = 'morse'
+        elif self.functional_form != 'morse':
+            logger.warning('Incorrect functional form set for MorseBond.')
+
+        if self.well_depth is not None and self.well_steepness is not None:
+            self.force_constant = 2.0 * self.well_depth * self.well_steepness**2
+        elif self.well_depth is not None and self.force_constant is not None:
+            self.well_steepness = np.sqrt(self.force_constant / (2.0 * self.well_depth))
+
+
+class FeneBond(BondPotential):
+    r"""
+    Section containing information about a FENE potential:
+    $V(r) = -\frac{1}{2} k R_0^2 \ln \left[ 1 - \left( \frac{r - r_0}{R_0} \right)^2 \right] + C$,
+    $k$ is the `force_constant`, $r_0$ is the `equilibrium_value` of $r$, and $R_0$ is the
+    maximum allowable bond extension beyond $r_0$. $C$ is an arbitrary constant (not stored).
+    """
+
+    maximum_extension = Quantity(
+        type=np.float64,
+        unit='m',
+        shape=[],
+        description="""
+        Specifies the maximum extension beyond the equilibrium bond distance.
+        """,
+    )
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        self.name = 'FeneBond'
+        if not self.functional_form:
+            self.functional_form = 'fene'
+        elif self.functional_form != 'fene':
+            logger.warning('Incorrect functional form set for FeneBond.')
+
+
+class TabulatedBond(TabulatedPotential, BondPotential):
+    """
+    Section containing information about a tabulated bond potential.
+    The value of the potential and/or force is stored for a set of corresponding bin distances.
+    """
+
+    bins = Quantity(
+        type=np.float64,
+        unit='m',
+        shape=[],
+        description="""
+        List of bin distances.
+        """,
+    )
+
+    forces = Quantity(
+        type=np.float64,
+        unit='J/m',
+        shape=[],
+        description="""
+        List of force values associated with each bin.
+        """,
+    )
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        self.name = 'TabulatedBond'
+
+
+class AnglePotential(Potential):
+    """
+    Section containing information about angle potentials.
+
+    Suggested types are: harmonic, cosine, restricted_cosinse, fourier_series, urey_bradley, polynomial, tabulated
+    """
+
+    equilibrium_value = Quantity(
+        type=np.float64,
+        unit='degree',
+        shape=[],
+        description="""
+        Specifies the equilibrium angle.
+        """,
+    )
+
+    force_constant = Quantity(
+        type=np.float64,
+        shape=[],
+        unit='J / degree**2',
+        description="""
+        Specifies the force constant of the angle potential.
+        """,
+    )
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        self.name = 'AnglePotential'
+        if not self.type:
+            self.type = 'angle'
+        elif self.type != 'angle':
+            logger.warning('Incorrect type set for AnglePotential.')
+
+        if self.n_particles:
+            if self.n_particles != 3:
+                logger.warning('Incorrect number of particles set for AnglePotential.')
+            else:
+                self.n_particles = 3
+
+
+class HarmonicAngle(AnglePotential):
+    r"""
+    Section containing information about a Harmonic angle potential:
+    $V(\theta) = \frac{1}{2} k_\theta (\theta - \theta_0)^2 + C$,
+    where $k_\theta$ is the `force_constant` and $\theta_0$ is the `equilibrium_value` of
+    the angle $\theta$. $C$ is an arbitrary constant (not stored).
+    """
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        self.name = 'HarmonicAngle'
+        if not self.functional_form:
+            self.functional_form = 'harmonic'
+        elif self.functional_form != 'harmonic':
+            logger.warning('Incorrect functional form set for HarmonicAngle.')
+
+
+class CosineAngle(AnglePotential):
+    r"""
+    Section containing information about a Cosine angle potential:
+    $V(\theta) = \frac{1}{2} k_\theta \left[ \cos(\theta) - \cos(\theta_0) \right]^2 + C$,
+    where $k_\theta$ is the `force_constant` and $\theta_0$ is the `equilibrium_value` of
+    the angle $\theta$. $C$ is an arbitrary constant (not stored).
+
+    """
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        self.name = 'CosineAngle'
+        if not self.functional_form:
+            self.functional_form = 'cosine'
+        elif self.functional_form != 'cosine':
+            logger.warning('Incorrect functional form set for CosineAngle.')
+
+
+class RestrictedCosineAngle(AnglePotential):
+    r"""
+    Section containing information about a Restricted Cosine angle potential:
+    $V(\theta) = \frac{1}{2} k_\theta \frac{\left[ \cos(\theta) - \cos(\theta_0) \right]^2}{\sin^2(\theta)} + C$,
+    where $k_\theta$ is the `force_constant` and $\theta_0$ is the `equilibrium_value` of
+    the angle $\theta$. $C$ is an arbitrary constant (not stored).
+
+    """
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        self.name = 'RestrictedCosineAngle'
+        if not self.functional_form:
+            self.functional_form = 'restricted_cosine'
+        elif self.functional_form != 'restricted_cosine':
+            logger.warning('Incorrect functional form set for RestrictedCosineAngle.')
+
+
+# TODO I think this was intended for the dihedral, like RB function, could generalize or just sepecify as dihedral potential
+# TODO Similar to UB pot, just say AKA RB function
+class FourierSeriesAngle(AnglePotential):
+    r"""
+    Section containing information about a Fourier Series angle potential:
+    $V(\theta) = k_1 \left[ 1 + \cos(\theta - \theta_0) \right] +
+    k_2 \left[ 1 + \cos(2(\theta - \theta_0)) \right] + \dots + C$,
+    where $\{k_1, k_2, \dots}$ are the `fourier_force_constants` for each term in the series
+    and $\theta_0$ is the `equilibrium_value` of the angle $\theta$.
+    $C$ is an arbitrary constant (not stored).
+    """
+
+    fourier_force_constants = Quantity(
+        type=np.float64,
+        shape=['*'],
+        unit='J / degree**2',
+        description="""
+        Specifies the force constants for each term in the fourier series representation
+        of the angle potential. The force constants of any "missing" terms must be
+        explicitly set to zero.
+        """,
+    )
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        self.name = 'FourierSeriesAngle'
+        if not self.functional_form:
+            self.functional_form = 'fourier_series'
+        elif self.functional_form != 'fourier_series':
+            logger.warning('Incorrect functional form set for FourierSeriesAngle.')
+
+
+# TODO I think we should name these more generally and then say that AKA
+class UreyBradleyAngle(AnglePotential):
+    r"""
+    Section containing information about a Urey-Bradly angle potential:
+    $V(\theta) = \frac{1}{2} k_\theta (\theta - \theta_0)^2 +
+    \frac{1}{2} k_{13} (r_{13} - r_{13}^0)^2 + C$,
+    where where $k_\theta$ is the `force_constant` and $\theta_0$ is the `equilibrium_value` of
+    the angle $\theta$, as for a harmonic angle potential. $k_{13}$ is the `force_constant_UB`
+    for the 1-3 term, and $r_{13}^0$ is the `equilibrium_value_UB` of the 1-3 distance (i.e.,
+    the distance between the edge particles forming the 1-2-3 angle, $\theta$), $r_{13}$.
+    $C$ is an arbitrary constant (not stored).
+    """
+
+    equilibrium_value_UB = Quantity(
+        type=np.float64,
+        unit='m',
+        shape=[],
+        description="""
+        Specifies the equilibrium 1-3 distance.
+        """,
+    )
+
+    force_constant_UB = Quantity(
+        type=np.float64,
+        shape=[],
+        unit='J / m **2',
+        description="""
+        Specifies the force constant of the potential.
+        """,
+    )
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        self.name = 'UreyBradleyAngle'
+        if not self.functional_form:
+            self.functional_form = 'urey_bradley'
+        elif self.functional_form != 'urey_bradley':
+            logger.warning('Incorrect functional form set for UreyBradleyAngle.')
+
+
+class PolynomialAngle(PolynomialPotential, AnglePotential):
+    """
+    Section containing information about a polynomial angle potential:
+    """
+
+    def __init__(self):
+        super().__init__()
+        docstring = PolynomialPotential.__doc__
+        pattern = r'\$V\(x\)(.*?)(\(not stored\)\.)'
+        match = re.search(pattern, docstring, re.DOTALL)
+        extracted_text = '<functional form missing>'
+        if match:
+            extracted_text = match.group().strip()  # .group(1).strip()
+        self.__doc__ = rf"""{self.__doc__} {extracted_text}.
+            Here the dependent variable of the potential, $x$, corresponds to the angle between three particles."""
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        self.name = 'PolynomialAngle'
+
+
+class TabulatedAngle(AnglePotential, TabulatedPotential):
+    """
+    Section containing information about a tabulated bond potential. The value of the potential and/or force
+    is stored for a set of corresponding bin distances.
+    """
+
+    bins = Quantity(
+        type=np.float64,
+        unit='degree',
+        shape=[],
+        description="""
+        List of bin angles.
+        """,
+    )
+
+    forces = Quantity(
+        type=np.float64,
+        unit='J/degree',
+        shape=[],
+        description="""
+        List of force values associated with each bin.
+        """,
+    )
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        self.name = 'TabulatedAngle'
+
+
+class DihedralPotential(Potential):
+    """
+    Section containing information about dihedral potentials.
+
+    Suggested types are: fourier_series, tabulated
+
+    # ? Something about angle convention?
+    """
+
+    equilibrium_value = Quantity(
+        type=np.float64,
+        unit='degree',
+        shape=[],
+        description="""
+        Specifies the equilibrium dihedral angle.
+        """,
+    )
+
+    force_constant = Quantity(
+        type=np.float64,
+        shape=[],
+        unit='J / degree**2',
+        description="""
+        Specifies the force constant of the dihedral angle potential.
+        """,
+    )
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        self.name = 'DihedralPotential'
+        if not self.type:
+            self.type = 'dihedral'
+        elif self.type != 'dihedral':
+            logger.warning('Incorrect type set for DihedralPotential.')
+
+        if self.n_particles:
+            if self.n_particles != 4:
+                logger.warning(
+                    'Incorrect number of particles set for DihedralPotential.'
+                )
+            else:
+                self.n_particles = 4
+
+
+class TabulatedDihedral(DihedralPotential, TabulatedPotential):
+    """
+    Section containing information about a tabulated bond potential. The value of the potential and/or force
+    is stored for a set of corresponding bin distances.
+    """
+
+    bins = Quantity(
+        type=np.float64,
+        unit='degree',
+        shape=[],
+        description="""
+        List of bin dihedral angles.
+        """,
+    )
+
+    forces = Quantity(
+        type=np.float64,
+        unit='J/degree',
+        shape=[],
+        description="""
+        List of force values associated with each bin.
+        """,
+    )
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        self.name = 'TabulatedDihedral'
+
+
+class ForceField(ModelMethod):
+    """
+    Section containing the parameters of a (classical, particle-based) force field model.
+    Typical `numerical_settings` are ForceCalculations.
+    Lists of interactions by type and, if available, corresponding parameters can be given within `interactions`.
+    Additionally, a published model can be referenced with `reference`.
+    """
+
+    # name and external reference already defined in BaseModelMethod
+    # name = Quantity(
+    #     type=str,
+    #     shape=[],
+    #     description="""
+    #     Identifies the name of the model.
+    #     """,
+    # )
+
+    # reference = Quantity(
+    #     type=str,
+    #     shape=['0..*'],
+    #     description="""
+    #     List of references to the model e.g. DOI, URL.
+    #     """,
+    # )
+
+    kimid = Quantity(
+        type=URL,
+        description="""
+        Reference to a model stored on the OpenKim database.
+        """,
+        a_eln=ELNAnnotation(component='URLEditQuantity'),
+    )
+
+    #     interactions = SubSection(sub_section=Interactions.m_def, repeats=True)
+    contributions = SubSection(
+        sub_section=Potential.m_def,
+        repeats=True,
+        description="""
+        Contribution or sub-term of the total model Hamiltonian.
+        """,
+    )
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        self.name = 'ForceField'
+        logger.warning('in force field')
+
+
+## GROMACS BONDED INTERACTIONS
+
+# Fourth power potential -- Use PolynomialBond? (Could also use it for harmonic and cubic but I guess it's better to have simpler forms for these more common ones)
+
+# Linear Angle potential -- useful?
+
+# Bond-Angle cross term -- useful?
+
+# Quartic angle potential -- Use polynomial angle?
+
+# Improper dihedrals
+
+# Improper dihedrals: harmonic type
+
+# Improper dihedrals: periodic type
+
+# Proper dihedrals
+
+# Proper dihedrals: periodic type
+
+# Proper dihedral angles are defined according to the IUPAC/IUB convention, where
+
+# Proper dihedrals: Ryckaert-Bellemans function
+
+# Proper dihedrals: Combined bending-torsion potential
+
+# Bonded pair and 1-4 interactions
+# Most force fields do not use normal Lennard-Jones and Coulomb interactions for atoms separated by three bonds, the so-called 1-4 interactions. These interactions are still affected by the modified electronic distributions due to the chemical bonds and they are modified in the force field by the dihedral terms. For this reason the Lennard-Jones and Coulomb 1-4 interactions are often scaled down, by a fixed factor given by the force field. These factors can be supplied in the topology and the parameters can also be overriden per 1-4 interaction or atom type pair. The pair interactions can be used for any atom pair in a molecule, not only 1-4 pairs. The non-bonded interactions between such pairs should be excluded to avoid double interactions. Plain Lennard-Jones and Coulomb interactions are used which are not affected by the non-bonded interaction treatment and potential modifiers.
+
+# Tabulated bonded interaction functions
+
+
+# ! Need to survey Lammps and maybe openmm and check for other common potential types
diff --git a/src/nomad_simulations/schema_packages/model_method.py b/src/nomad_simulations/schema_packages/model_method.py
index c7b143ff..06b0878b 100644
--- a/src/nomad_simulations/schema_packages/model_method.py
+++ b/src/nomad_simulations/schema_packages/model_method.py
@@ -88,6 +88,15 @@ class ModelMethod(BaseModelMethod):
         """,
     )
 
+    # ? Not sure about this, as it is more linked to a workflow in a way, but one may apply different models to different systems stored within a single entry
+    model_system_ref = Quantity(
+        type=ModelSystem,
+        description="""
+        Reference to the `ModelSystem` section to which the model was applied.
+        """,
+        a_eln=ELNAnnotation(component='ReferenceEditQuantity'),
+    )
+
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
 
diff --git a/src/nomad_simulations/schema_packages/numerical_settings.py b/src/nomad_simulations/schema_packages/numerical_settings.py
index 515deea7..05c0a6fd 100644
--- a/src/nomad_simulations/schema_packages/numerical_settings.py
+++ b/src/nomad_simulations/schema_packages/numerical_settings.py
@@ -887,3 +887,86 @@ def __init__(self, m_def: 'Section' = None, m_context: 'Context' = None, **kwarg
 
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
+
+
+class ForceCalculations(NumericalSettings):
+    """
+    Section containing the parameters for force calculations according to a ForceField model.
+    """
+
+    vdw_cutoff = Quantity(
+        type=np.float64,
+        shape=[],
+        unit='m',
+        description="""
+        Cutoff for calculating VDW forces.
+        """,
+    )
+
+    coulomb_type = Quantity(
+        type=MEnum(
+            'cutoff',
+            'ewald',
+            'multilevel_summation',
+            'particle_mesh_ewald',
+            'particle_particle_particle_mesh',
+            'reaction_field',
+        ),
+        shape=[],
+        description="""
+        Method used for calculating long-ranged Coulomb forces.
+
+        Allowed values are:
+
+        | Method Name          | Description                               |
+
+        | ---------------------- | ----------------------------------------- |
+
+        | `""`                   | No thermostat               |
+
+        | `"Cutoff"`          | Simple cutoff scheme. |
+
+        | `"Ewald"` | Standard Ewald summation as described in any solid-state physics text. |
+
+        | `"Multi-Level Summation"` |  D. Hardy, J.E. Stone, and K. Schulten,
+        [Parallel. Comput. **35**, 164](https://doi.org/10.1016/j.parco.2008.12.005)|
+
+        | `"Particle-Mesh-Ewald"`        | T. Darden, D. York, and L. Pedersen,
+        [J. Chem. Phys. **98**, 10089 (1993)](https://doi.org/10.1063/1.464397) |
+
+        | `"Particle-Particle Particle-Mesh"` | See e.g. Hockney and Eastwood, Computer Simulation Using Particles,
+        Adam Hilger, NY (1989). |
+
+        | `"Reaction-Field"` | J.A. Barker and R.O. Watts,
+        [Mol. Phys. **26**, 789 (1973)](https://doi.org/10.1080/00268977300102101)|
+        """,
+    )
+
+    coulomb_cutoff = Quantity(
+        type=np.float64,
+        shape=[],
+        unit='m',
+        description="""
+        Cutoff for calculating short-ranged Coulomb forces.
+        """,
+    )
+
+    neighbor_update_frequency = Quantity(
+        type=int,
+        shape=[],
+        description="""
+        Number of timesteps between updating the neighbor list.
+        """,
+    )
+
+    neighbor_update_cutoff = Quantity(
+        type=np.float64,
+        shape=[],
+        unit='m',
+        description="""
+        The distance cutoff for determining the neighbor list.
+        """,
+    )
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
diff --git a/src/nomad_simulations/schema_packages/outputs.py b/src/nomad_simulations/schema_packages/outputs.py
index 1491fda3..95ca7c98 100644
--- a/src/nomad_simulations/schema_packages/outputs.py
+++ b/src/nomad_simulations/schema_packages/outputs.py
@@ -48,8 +48,6 @@ class Outputs(ArchiveSection):
     information if the output `is_derived` from another output section or directly parsed from the simulation output files.
     """
 
-    # TODO add time quantities
-
     normalizer_level = 2
 
     model_system_ref = Quantity(
diff --git a/src/nomad_simulations/schema_packages/properties/forces.py b/src/nomad_simulations/schema_packages/properties/forces.py
index cedf72f1..ef9f937f 100644
--- a/src/nomad_simulations/schema_packages/properties/forces.py
+++ b/src/nomad_simulations/schema_packages/properties/forces.py
@@ -26,12 +26,19 @@ class BaseForce(PhysicalProperty):
     """
 
     value = Quantity(
-        type=np.dtype(np.float64),
+        type=np.float64,
+        shape=['*'],
         unit='newton',
         description="""
         """,
     )
 
+    def __init__(
+        self, m_def: 'Section' = None, m_context: 'Context' = None, **kwargs
+    ) -> None:
+        super().__init__(m_def, m_context, **kwargs)
+        self.rank = [3]
+
     def normalize(self, archive: 'EntryArchive', logger: 'BoundLogger') -> None:
         super().normalize(archive, logger)
 
diff --git a/tests/test_force_field.py b/tests/test_force_field.py
new file mode 100644
index 00000000..951132aa
--- /dev/null
+++ b/tests/test_force_field.py
@@ -0,0 +1,889 @@
+import pytest
+import numpy as np
+from typing import List, Tuple, Dict, Any
+
+
+from nomad_simulations.schema_packages.general import Simulation
+
+# from nomad_simulations.schema_packages.method import ModelMethod
+from nomad_simulations.schema_packages.force_field import (
+    ForceField,
+    ParameterEntry,
+    Potential,
+    TabulatedPotential,
+    BondPotential,
+    HarmonicBond,
+    CubicBond,
+    PolynomialForceConstant,
+    PolynomialBond,
+    MorseBond,
+    FeneBond,
+    TabulatedBond,
+    AnglePotential,
+    HarmonicAngle,
+    CosineAngle,
+    RestrictedCosineAngle,
+    FourierSeriesAngle,
+    UreyBradleyAngle,
+    PolynomialAngle,
+    TabulatedAngle,
+)
+from nomad_simulations.schema_packages.numerical_settings import ForceCalculations
+
+from nomad.datamodel import EntryArchive
+
+# from structlog.stdlib import BoundLogger
+from . import logger
+from nomad.units import ureg
+
+
+MOL = 6.022140857e23
+
+
+def approx(value, abs=0, rel=1e-6):
+    return pytest.approx(value, abs=abs, rel=rel)
+
+
+def assert_dict_equal(d1, d2):
+    """
+    Recursively assert that two dictionaries are equal.
+
+    Args:
+        d1 (dict): First dictionary to compare.
+        d2 (dict): Second dictionary to compare.
+    """
+
+    assert isinstance(d1, dict), f'Expected dict, got {type(d1)}'
+    assert isinstance(d2, dict), f'Expected dict, got {type(d2)}'
+    assert d1.keys() == d2.keys(), f'Keys mismatch: {d1.keys()} != {d2.keys()}'
+
+    def compare_unknown(key, value1, value2):
+        assert value1 == None, f"Unknown types'{key}': {type(value1)} != {type(value2)}"
+
+    def compare_string(key, str1, str2):
+        assert str1 == str2, f"Value mismatch for key '{key}': {str1} != {str2}"
+
+    def compare_float(key, float1, float2):
+        if abs(float1) == float('inf'):
+            assert 'inf' == float2 if float1 > 0 else '-inf' == float2
+        else:
+            assert float1 == approx(
+                float2
+            ), f"Value mismatch for key '{key}': {float1} != {float2}"
+
+    def compare_arrays(key, arr1, arr2):
+        assert np.isclose(
+            arr1, arr2
+        ).all(), f"Value mismatch for key '{key}': {arr1} != {arr2}"
+
+    def compare_lists(key, l1, l2):
+        assert len(l1) == len(
+            l2
+        ), f"Length mismatch for key '{key}': {len(l1)} != {len(l2)}"
+
+        for i, l1_item in enumerate(l1):
+            if isinstance(l1_item, dict) and isinstance(l2[i], dict):
+                assert_dict_equal(l1_item, l2[i])
+            elif isinstance(l1_item, (str, bool)) and isinstance(l2[i], (str, bool)):
+                compare_string(f'{key}-{i}', l1_item, l2[i])
+            elif isinstance(l1_item, list) and isinstance(l2[i], list):
+                compare_lists(f'{key}-{i}', l1_item, l2[i])
+            elif isinstance(l1_item, np.ndarray) and isinstance(l2[i], np.ndarray):
+                compare_arrays(f'{key}-{i}', l1_item, l2[i])
+            elif isinstance(l1_item, (float, int)) and isinstance(l2[i], (float, int)):
+                compare_float(f'{key}-{i}', l1_item, l2[i])
+            else:
+                compare_unknown(f'{key}-{i}', l1_item, l2[i])
+
+    for key in d1:
+        print(f'key: {key}', d1[key], d2[key])
+        if isinstance(d1[key], dict) and isinstance(d2[key], dict):
+            assert_dict_equal(d1[key], d2[key])
+        elif isinstance(d1[key], (str, bool)) and isinstance(d2[key], (str, bool)):
+            compare_string(key, d1[key], d2[key])
+        elif isinstance(d1[key], list) and isinstance(d2[key], list):
+            compare_lists(key, d1[key], d2[key])
+        elif isinstance(d1[key], np.ndarray) and isinstance(d2[key], np.ndarray):
+            compare_arrays(key, d1[key], d2[key])
+        elif isinstance(d1[key], (float, int)) and isinstance(d2[key], (float, int)):
+            compare_float(key, d1[key], d2[key])
+        else:
+            compare_unknown(key, d1[key], d2[key])
+
+
+def get_simulation_template():
+    data = Simulation()
+    sec_FF = ForceField()
+    data.model_method.append(sec_FF)
+    sec_force_calculations = ForceCalculations()
+    data.model_method[0].numerical_settings.append(sec_force_calculations)
+
+    return data
+
+
+def populate_potential(
+    class_potential,
+    n_interactions=None,
+    n_particles=None,
+    particle_labels=None,
+    particle_indices=None,
+):
+    sec_potential = class_potential()
+    sec_potential.n_interactions = n_interactions
+    sec_potential.n_particles = n_particles
+    sec_potential.particle_indices = particle_indices
+    sec_potential.particle_labels = particle_labels
+
+    return sec_potential
+
+
+def populate_parameters(sec_potential, parameters):
+    for key, value in parameters.items():
+        if key == 'parameter_entries':
+            for entry in value:
+                sec_parameter = ParameterEntry()
+                sec_parameter.name = entry['name']
+                sec_parameter.value = entry['value']
+                sec_parameter.unit = entry['unit']
+                sec_potential.parameters.append(sec_parameter)
+        else:
+            setattr(sec_potential, key, value)
+
+
+# Test Data
+
+# BOND POTENTIALS
+
+# System: 3 x OH molecules
+#   particle number       particle label
+#   0                     O
+#   1                     H
+#   2                     O
+#   3                     H
+#   4                     O
+#   5                     H
+n_interactions = 3
+n_particles = 2
+particle_labels = [('O', 'H'), ('O', 'H'), ('O', 'H')]
+particle_indices = [(0, 1), (2, 3), (4, 5)]
+
+# harmonic
+results_harmonic_bond = {
+    'n_interactions': 3,
+    'n_particles': 2,
+    'particle_indices': [[0, 1], [2, 3], [4, 5]],
+    'particle_labels': [['O', 'H'], ['O', 'H'], ['O', 'H']],
+    'equilibrium_value': 9.6e-11,
+    'force_constant': 5.811886641495074,
+    'name': 'HarmonicBond',
+    'type': 'bond',
+    'functional_form': 'harmonic',
+}
+data_harmonic_bond = (
+    HarmonicBond,
+    n_interactions,
+    n_particles,
+    particle_labels,
+    particle_indices,
+    {
+        'equilibrium_value': 0.96 * ureg.angstrom,
+        'force_constant': 3500 * ureg.kJ / MOL / ureg.nanometer**2,
+    },
+    results_harmonic_bond,
+)
+
+
+# cubic
+results_cubic_bond = {
+    'n_interactions': 3,
+    'n_particles': 2,
+    'particle_indices': [[0, 1], [2, 3], [4, 5]],
+    'particle_labels': [['O', 'H'], ['O', 'H'], ['O', 'H']],
+    'equilibrium_value': 9.6e-11,
+    'force_constant': 8.302695202135819,
+    'force_constant_cubic': 332107808.0854328,
+    'name': 'CubicBond',
+    'type': 'bond',
+    'functional_form': 'cubic',
+}
+data_cubic_bond = (
+    CubicBond,
+    n_interactions,
+    n_particles,
+    particle_labels,
+    particle_indices,
+    {
+        'equilibrium_value': 0.96 * ureg.angstrom,
+        'force_constant': 5000 * ureg.kJ / MOL / ureg.nanometer**2,
+        'force_constant_cubic': 200 * ureg.kJ / MOL / ureg.nanometer**3,
+    },
+    results_cubic_bond,
+)
+
+# polynomial
+results_polynomial_bond = {
+    'n_interactions': 3,
+    'n_particles': 2,
+    'particle_indices': [[0, 1], [2, 3], [4, 5]],
+    'particle_labels': [['O', 'H'], ['O', 'H'], ['O', 'H']],
+    'equilibrium_value': 9.6e-11,
+    'force_constants': [
+        {
+            'name': 'k_2',
+            'value': 8.302695202135821e-21,
+            'unit': 'kilojoule / nanometer ** 2',
+            'exponent': 2,
+        },
+        {
+            'name': 'k_4',
+            'value': 3.4871319848970445e-21,
+            'unit': 'kilojoule / nanometer ** 4',
+            'exponent': 4,
+        },
+    ],
+    'name': 'PolynomialBond',
+    'type': 'bond',
+    'functional_form': 'polynomial',
+}
+data_polynomial_bond = (
+    PolynomialBond,
+    n_interactions,
+    n_particles,
+    particle_labels,
+    particle_indices,
+    {
+        'equilibrium_value': 0.96 * ureg.angstrom,
+        'force_constants': [
+            PolynomialForceConstant(
+                name='k_2',
+                exponent=2,
+                value=5000.0 / MOL,
+                unit=str(ureg.kJ / ureg.nanometer**2),
+            ),
+            PolynomialForceConstant(
+                name='k_4',
+                exponent=4,
+                value=2100.0 / MOL,
+                unit=str(ureg.kJ / ureg.nanometer**4),
+            ),
+        ],
+    },
+    results_polynomial_bond,
+)
+
+# morse
+results_morse_bond = {
+    'n_interactions': 3,
+    'n_particles': 2,
+    'particle_indices': [[0, 1], [2, 3], [4, 5]],
+    'particle_labels': [['O', 'H'], ['O', 'H'], ['O', 'H']],
+    'equilibrium_value': 9.6e-11,
+    'well_depth': 7.472425681922239e-18,
+    'well_steepness': 24999999999.999996,
+    'name': 'MorseBond',
+    'type': 'bond',
+    'functional_form': 'morse',
+    'force_constant': 9340.532102402796,
+}
+data_morse_bond = (
+    MorseBond,
+    n_interactions,
+    n_particles,
+    particle_labels,
+    particle_indices,
+    {
+        'equilibrium_value': 0.96 * ureg.angstrom,
+        'well_depth': 4500 * ureg.kJ / MOL,
+        'well_steepness': 25 * (1 / ureg.nanometer),
+    },
+    results_morse_bond,
+)
+
+# fene
+results_fene_bond = {
+    'n_interactions': 3,
+    'n_particles': 2,
+    'particle_indices': [[0, 1], [2, 3], [4, 5]],
+    'particle_labels': [['O', 'H'], ['O', 'H'], ['O', 'H']],
+    'equilibrium_value': 9.6e-11,
+    'maximum_extension': 5e-11,
+    'force_constant': 6.227021401601864,
+    'name': 'FeneBond',
+    'type': 'bond',
+    'functional_form': 'fene',
+}
+data_fene_bond = (
+    FeneBond,
+    n_interactions,
+    n_particles,
+    particle_labels,
+    particle_indices,
+    {
+        'equilibrium_value': 0.96 * ureg.angstrom,
+        'maximum_extension': 0.5 * ureg.angstrom,
+        'force_constant': 3750 * ureg.kJ / MOL / ureg.nanometer**2,
+    },
+    results_fene_bond,
+)
+
+# tabulated
+results_tabulated_bond = {
+    'n_interactions': 3,
+    'n_particles': 2,
+    'particle_indices': [[0, 1], [2, 3], [4, 5]],
+    'particle_labels': [['O', 'H'], ['O', 'H'], ['O', 'H']],
+    'bins': [
+        7.600e-11,
+        7.970e-11,
+        8.340e-11,
+        8.710e-11,
+        9.070e-11,
+        9.440e-11,
+        9.810e-11,
+        1.018e-10,
+        1.055e-10,
+        1.092e-10,
+        1.128e-10,
+        1.165e-10,
+        1.202e-10,
+        1.239e-10,
+        1.276e-10,
+        1.313e-10,
+        1.349e-10,
+        1.386e-10,
+        1.423e-10,
+        1.460e-10,
+    ],
+    'energies': [
+        1.32311751e-21,
+        8.81248069e-22,
+        5.28549577e-22,
+        2.65354139e-22,
+        9.18278089e-23,
+        8.30269520e-24,
+        1.47787975e-23,
+        1.11422170e-22,
+        2.98564919e-22,
+        5.76539155e-22,
+        9.45178822e-22,
+        1.40498208e-21,
+        1.95611499e-21,
+        2.59857754e-21,
+        3.33286791e-21,
+        4.15882003e-21,
+        5.07693206e-21,
+        6.08737007e-21,
+        7.19013405e-21,
+        8.38572215e-21,
+    ],
+    'forces': [
+        1.3216958010639183e-10,
+        1.0784837919852664e-10,
+        8.347183948070384e-11,
+        5.903996095292366e-11,
+        3.521528791034704e-11,
+        1.0674227406164319e-11,
+        -1.3922171907975959e-11,
+        -3.857391003207359e-11,
+        -6.32809869661286e-11,
+        -8.804340271014125e-11,
+        -1.1218967953067358e-10,
+        -1.3706127725108878e-10,
+        -1.6198821378146133e-10,
+        -1.8697048912179138e-10,
+        -2.1200810327207885e-10,
+        -2.3710105623232373e-10,
+        -2.6156893683081195e-10,
+        -2.867710717674596e-10,
+        -3.1202854551406456e-10,
+        -3.373413580706269e-10,
+    ],
+    'name': 'TabulatedBond',
+    'type': 'bond',
+    'functional_form': 'tabulated',
+}
+data_tabulated_bond = (
+    TabulatedBond,
+    n_interactions,
+    n_particles,
+    particle_labels,
+    particle_indices,
+    {
+        'bins': [
+            0.076,
+            0.0797,
+            0.0834,
+            0.0871,
+            0.0907,
+            0.0944,
+            0.0981,
+            0.1018,
+            0.1055,
+            0.1092,
+            0.1128,
+            0.1165,
+            0.1202,
+            0.1239,
+            0.1276,
+            0.1313,
+            0.1349,
+            0.1386,
+            0.1423,
+            0.146,
+        ]
+        * ureg.nanometer,
+        'energies': [  # ! pass in energies and test the auto-generation of forces
+            0.7968,
+            0.5307,
+            0.3183,
+            0.1598,
+            0.0553,
+            0.005,
+            0.0089,
+            0.0671,
+            0.1798,
+            0.3472,
+            0.5692,
+            0.8461,
+            1.178,
+            1.5649,
+            2.0071,
+            2.5045,
+            3.0574,
+            3.6659,
+            4.33,
+            5.05,
+        ]
+        * ureg.kJ
+        / MOL,
+    },
+    results_tabulated_bond,
+)
+
+# custom - LJ
+results_custom_bond = {
+    'n_interactions': 3,
+    'n_particles': 2,
+    'particle_indices': [[0, 1], [2, 3], [4, 5]],
+    'particle_labels': [['O', 'H'], ['O', 'H'], ['O', 'H']],
+    'name': 'BondPotential',
+    'parameters': [
+        {'name': 'epsilon', 'value': '2.5738355126621044e-25', 'unit': 'kilojoule'},
+        {'name': 'sigma', 'value': '0.96', 'unit': 'angstrom'},
+    ],
+}
+data_custom_bond = (
+    BondPotential,
+    n_interactions,
+    n_particles,
+    particle_labels,
+    particle_indices,
+    {
+        'parameter_entries': [
+            {'name': 'epsilon', 'value': 0.155 / MOL, 'unit': str(ureg.kJ)},
+            {'name': 'sigma', 'value': 0.96, 'unit': str(ureg.angstrom)},
+        ],
+    },
+    results_custom_bond,
+)
+
+# Angle POTENTIALS
+
+# System: 2 x H20 molecules
+#   particle number       particle label
+#   0                     O
+#   1                     H
+#   2                     H
+#   3                     O
+#   4                     H
+#   5                     H
+n_interactions = 2
+n_particles = 3
+particle_labels = [('O', 'H', 'H'), ('O', 'H', 'H')]
+particle_indices = [(0, 1, 2), (3, 4, 5)]
+
+# harmonic
+results_harmonic_angle = {
+    'n_interactions': 2,
+    'n_particles': 3,
+    'particle_indices': [[0, 1, 2], [3, 4, 5]],
+    'particle_labels': [['O', 'H', 'H'], ['O', 'H', 'H']],
+    'equilibrium_value': 104.45020605234907,
+    'force_constant': 3.7937183846251475e-23,
+    'name': 'HarmonicAngle',
+    'type': 'angle',
+    'functional_form': 'harmonic',
+}
+data_harmonic_angle = (
+    HarmonicAngle,
+    n_interactions,
+    n_particles,
+    particle_labels,
+    particle_indices,
+    {
+        'equilibrium_value': 1.823 * ureg.radian,
+        'force_constant': 75 * ureg.kJ / MOL / ureg.radian**2,
+    },
+    results_harmonic_angle,
+)
+
+# cosine
+results_cosine_angle = {
+    'n_interactions': 2,
+    'n_particles': 3,
+    'particle_indices': [[0, 1, 2], [3, 4, 5]],
+    'particle_labels': [['O', 'H', 'H'], ['O', 'H', 'H']],
+    'equilibrium_value': 104.45020605234907,
+    'force_constant': 3.7937183846251475e-23,
+    'name': 'CosineAngle',
+    'type': 'angle',
+    'functional_form': 'cosine',
+}
+data_cosine_angle = (
+    CosineAngle,
+    n_interactions,
+    n_particles,
+    particle_labels,
+    particle_indices,
+    {
+        'equilibrium_value': 1.823 * ureg.radian,
+        'force_constant': 75 * ureg.kJ / MOL / ureg.radian**2,
+    },
+    results_cosine_angle,
+)
+
+# restricted cosine
+results_restrictedcosine_angle = {
+    'n_interactions': 2,
+    'n_particles': 3,
+    'particle_indices': [[0, 1, 2], [3, 4, 5]],
+    'particle_labels': [['O', 'H', 'H'], ['O', 'H', 'H']],
+    'equilibrium_value': 104.45020605234907,
+    'force_constant': 8.70026082874034e-23,
+    'name': 'RestrictedCosineAngle',
+    'type': 'angle',
+    'functional_form': 'restricted_cosine',
+}
+data_restrictedcosine_angle = (
+    RestrictedCosineAngle,
+    n_interactions,
+    n_particles,
+    particle_labels,
+    particle_indices,
+    {
+        'equilibrium_value': 1.823 * ureg.radian,
+        'force_constant': 172 * ureg.kJ / MOL / ureg.radian**2,
+    },
+    results_restrictedcosine_angle,
+)
+
+# fourier_series
+results_fourier_angle = {
+    'n_interactions': 2,
+    'n_particles': 3,
+    'particle_indices': [[0, 1, 2], [3, 4, 5]],
+    'particle_labels': [['O', 'H', 'H'], ['O', 'H', 'H']],
+    'equilibrium_value': 104.45020605234907,
+    'fourier_force_constants': [5.0582911795001975e-23, 0.0, 1.2645727948750494e-23],
+    'name': 'FourierSeriesAngle',
+    'type': 'angle',
+    'functional_form': 'fourier_series',
+}
+data_fourier_angle = (
+    FourierSeriesAngle,
+    n_interactions,
+    n_particles,
+    particle_labels,
+    particle_indices,
+    {
+        'equilibrium_value': 1.823 * ureg.radian,
+        'fourier_force_constants': [100, 0, 25] * ureg.kJ / MOL / ureg.radian**2,
+    },
+    results_fourier_angle,
+)
+
+# urey_bradley
+results_ureybradley_angle = {
+    'n_interactions': 2,
+    'n_particles': 3,
+    'particle_indices': [[0, 1, 2], [3, 4, 5]],
+    'particle_labels': [['O', 'H', 'H'], ['O', 'H', 'H']],
+    'equilibrium_value': 104.45020605234907,
+    'force_constant': 3.7937183846251475e-23,
+    'equilibrium_value_UB': 1.5140000000000001e-10,
+    'force_constant_UB': 4.9816171212814925e-19,
+    'name': 'UreyBradleyAngle',
+    'type': 'angle',
+    'functional_form': 'urey_bradley',
+}
+data_ureybradley_angle = (
+    UreyBradleyAngle,
+    n_interactions,
+    n_particles,
+    particle_labels,
+    particle_indices,
+    {
+        'equilibrium_value': 1.823 * ureg.radian,
+        'force_constant': 75 * ureg.kJ / MOL / ureg.radian**2,
+        'equilibrium_value_UB': 1.514 * ureg.angstrom,
+        'force_constant_UB': 300 * ureg.kJ / MOL / ureg.m**2,
+    },
+    results_ureybradley_angle,
+)
+
+# polynomial
+results_polynomial_angle = {
+    'n_interactions': 2,
+    'n_particles': 3,
+    'particle_indices': [[0, 1, 2], [3, 4, 5]],
+    'particle_labels': [['O', 'H', 'H'], ['O', 'H', 'H']],
+    'equilibrium_value': 104.45020605234907,
+    'force_constants': [
+        {
+            'name': 'k_2',
+            'value': 1.245404280320373e-22,
+            'unit': 'kilojoule / radian ** 2',
+            'exponent': 2,
+        },
+        {
+            'name': 'k_3',
+            'value': -3.3210780808543285e-24,
+            'unit': 'kilojoule / radian ** 4',
+            'exponent': 3,
+        },
+        {
+            'name': 'k_4',
+            'value': 7.472425681922239e-24,
+            'unit': 'kilojoule / radian ** 4',
+            'exponent': 4,
+        },
+    ],
+    'name': 'PolynomialAngle',
+    'type': 'angle',
+    'functional_form': 'polynomial',
+}
+data_polynomial_angle = (
+    PolynomialAngle,
+    n_interactions,
+    n_particles,
+    particle_labels,
+    particle_indices,
+    {
+        'equilibrium_value': 1.823 * ureg.radian,
+        'force_constants': [
+            PolynomialForceConstant(
+                name='k_2',
+                exponent=2,
+                value=75 / MOL,
+                unit=str(ureg.kJ / ureg.radian**2),
+            ),
+            PolynomialForceConstant(
+                name='k_3',
+                exponent=3,
+                value=-2.0 / MOL,
+                unit=str(ureg.kJ / ureg.radian**4),
+            ),
+            PolynomialForceConstant(
+                name='k_4',
+                exponent=4,
+                value=4.5 / MOL,
+                unit=str(ureg.kJ / ureg.radian**4),
+            ),
+        ],
+    },
+    results_polynomial_angle,
+)
+
+# tabulated
+results_tabulated_angle = {
+    'n_interactions': 2,
+    'n_particles': 3,
+    'particle_indices': [[0, 1, 2], [3, 4, 5]],
+    'particle_labels': [['O', 'H', 'H'], ['O', 'H', 'H']],
+    'bins': [
+        92.99105014973262,
+        94.19727699638311,
+        95.40350384303362,
+        96.60973068968411,
+        97.81595810929241,
+        99.02218495594292,
+        100.22841180259343,
+        101.43463864924392,
+        102.64086549589443,
+        103.84709234254493,
+        105.05331976215322,
+        106.25954660880373,
+        107.46577345545423,
+        108.67200030210473,
+        109.87822714875524,
+        111.08445399540572,
+        112.29068141501403,
+        113.49690826166454,
+        114.70313510831502,
+        115.90936195496555,
+    ],
+    'energies': [
+        2.483908793147618e-21,
+        1.9871270690593166e-21,
+        1.5455433097527105e-21,
+        1.1591575152954348e-21,
+        8.279695415419479e-22,
+        5.51979703171419e-22,
+        3.311878298023737e-22,
+        1.6559392146862525e-22,
+        5.519797819553163e-23,
+        0.0,
+        0.0,
+        5.519797819553239e-23,
+        1.6559392146862525e-22,
+        3.311878298023748e-22,
+        5.519797031714216e-22,
+        8.2796954154194335e-22,
+        1.1591575152954338e-21,
+        1.545543309752715e-21,
+        1.98712706905931e-21,
+        2.483908793147623e-21,
+    ],
+    'forces': [
+        4.347281042004185e-22,
+        3.8896725108092945e-22,
+        3.432063979614403e-22,
+        2.974455448419512e-22,
+        2.516846920122808e-22,
+        2.0592383889279167e-22,
+        1.6016298577330257e-22,
+        1.1440213265381345e-22,
+        6.864127953432432e-23,
+        2.288042641483519e-23,
+        -2.288042641483519e-23,
+        -6.864127953432432e-23,
+        -1.1440213265381345e-22,
+        -1.6016298577330257e-22,
+        -2.0592383889279167e-22,
+        -2.516846920122808e-22,
+        -2.974455448419512e-22,
+        -3.432063979614403e-22,
+        -3.8896725108092945e-22,
+        -4.347281042004185e-22,
+    ],
+    'name': 'TabulatedAngle',
+    'type': 'angle',
+    'functional_form': 'tabulated',
+}
+data_tabulated_angle = (
+    TabulatedAngle,
+    n_interactions,
+    n_particles,
+    particle_labels,
+    particle_indices,
+    {
+        'bins': [
+            1.623,
+            1.64405263,
+            1.66510526,
+            1.68615789,
+            1.70721053,
+            1.72826316,
+            1.74931579,
+            1.77036842,
+            1.79142105,
+            1.81247368,
+            1.83352632,
+            1.85457895,
+            1.87563158,
+            1.89668421,
+            1.91773684,
+            1.93878947,
+            1.95984211,
+            1.98089474,
+            2.00194737,
+            2.023,
+        ]
+        * ureg.radian,
+        'forces': [  # ! This time input the forces and test the auto-generation of energies
+            15.0,
+            13.42105263,
+            11.84210526,
+            10.26315789,
+            8.68421053,
+            7.10526316,
+            5.52631579,
+            3.94736842,
+            2.36842105,
+            0.78947368,
+            -0.78947368,
+            -2.36842105,
+            -3.94736842,
+            -5.52631579,
+            -7.10526316,
+            -8.68421053,
+            -10.26315789,
+            -11.84210526,
+            -13.42105263,
+            -15.0,
+        ]
+        * ureg.kJ
+        / MOL
+        / ureg.radian,
+    },
+    results_tabulated_angle,
+)
+
+
+@pytest.mark.parametrize(
+    'potential_class, n_interactions, n_particles, particle_labels, particle_indices, parameters, results',
+    [
+        data_harmonic_bond,
+        data_cubic_bond,
+        data_polynomial_bond,
+        data_morse_bond,
+        data_fene_bond,
+        data_tabulated_bond,
+        data_custom_bond,
+        data_harmonic_angle,
+        data_cosine_angle,
+        data_restrictedcosine_angle,
+        data_fourier_angle,
+        data_ureybradley_angle,
+        data_polynomial_angle,
+        data_tabulated_angle,
+    ],
+)
+def test_potentials(
+    potential_class: Potential,
+    n_interactions: int,
+    n_particles: int,
+    particle_labels: List[Tuple[str, str]],
+    particle_indices: List[Tuple[int, int]],
+    parameters: Dict[str, Any],
+    results: Dict[str, Any],
+):
+    """_summary_
+
+    Args:
+        input (str): _description_
+        result (Dict[Any]): _description_
+    """
+
+    data = get_simulation_template()
+    sec_FF = data.model_method[0]
+    sec_potential = populate_potential(
+        potential_class,
+        n_interactions=n_interactions,
+        n_particles=n_particles,
+        particle_labels=particle_labels,
+        particle_indices=particle_indices,
+    )
+    populate_parameters(sec_potential, parameters)
+
+    sec_FF.contributions.append(sec_potential)
+    sec_FF.contributions[-1].normalize(EntryArchive, logger)  # BoundLogger)
+
+    potential_dict = sec_FF.contributions[-1].m_to_dict()
+    potential_dict = {  # ! The dev is required to add new results to the dictionary, this will not be caught by the test!
+        key: value for key, value in potential_dict.items() if key in results
+    }
+    print(potential_dict)
+    print(results)
+    # assert 1 == 2
+    assert_dict_equal(potential_dict, results)