Skip to content

kups.potential.classical.morse

Morse bond stretching potential.

Reference: Rappé et al. (1992) "UFF, a Full Periodic Table Force Field" J. Am. Chem. Soc. 114, 10024-10035. DOI: 10.1021/ja00051a040

Functional form:

\[ U(r) = D \left[1 - e^{-\alpha(r - r_0)}\right]^2 \]

More accurate than harmonic for large displacements with proper dissociation behavior. Near equilibrium, Morse approximates harmonic with force constant \(k = 2 D \alpha^2\).

IsBondedParticles

Bases: HasPositionsAndLabels, IsRadiusGraphPoints, Protocol

Particle data with positions, labels, and system index.

Source code in src/kups/potential/classical/morse.py 
57
58
59
60
61
@runtime_checkable
class IsBondedParticles(HasPositionsAndLabels, IsRadiusGraphPoints, Protocol):
    """Particle data with positions, labels, and system index."""

    ...

MorseBondParameters

Morse bond potential parameters.

Attributes:

Name Type Description
labels tuple[Label, ...]

Species labels, shape (n_species,).
r0 Array

Equilibrium bond lengths [Å], shape (n_species, n_species).
D Array

Bond dissociation energy (well depth), shape (n_species, n_species).
alpha Array

Width parameter [Å⁻¹], shape (n_species, n_species).
Source code in src/kups/potential/classical/morse.py 
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
@dataclass
class MorseBondParameters:
    r"""Morse bond potential parameters.

    Attributes:
        labels: Species labels, shape `(n_species,)`.
        r0: Equilibrium bond lengths [Å], shape `(n_species, n_species)`.
        D: Bond dissociation energy (well depth), shape `(n_species, n_species)`.
        alpha: Width parameter [Å⁻¹], shape `(n_species, n_species)`.
    """

    labels: tuple[Label, ...] = field(static=True)  # (n_species,)
    r0: Array  # (n_species, n_species)
    D: Array  # (n_species, n_species)
    alpha: Array  # (n_species, n_species)

    @classmethod
    def from_harmonic(
        cls, labels: tuple[str, ...], r0: Array, k: Array, D: Array
    ) -> "MorseBondParameters":
        r"""Create Morse parameters from harmonic force constant.

        Args:
            labels: Species labels, shape `(n_species,)`
            r0: Equilibrium bond lengths [Å], shape `(n_species, n_species)`
            k: Harmonic force constants [energy/Ų], shape `(n_species, n_species)`
            D: Bond dissociation energies, shape `(n_species, n_species)`

        Returns:
            MorseBondParameters with computed alpha values
        """
        alpha = jnp.sqrt(k / D)
        return cls(labels=tuple(map(Label, labels)), r0=r0, D=D, alpha=alpha)

    @classmethod
    def from_uff(
        cls,
        labels: tuple[str, ...],
        bond_radius: Array,
        electronegativity: Array,
        effective_charge: Array,
        dissociation_energy: Array,
    ) -> "MorseBondParameters":
        r"""Create Morse parameters using UFF bond length/force constant formulas.

        Args:
            labels: Species labels, shape `(n_species,)`
            bond_radius: Valence bond radii [Å], shape `(n_species,)`
            electronegativity: GMP electronegativity, shape `(n_species,)`
            effective_charge: Effective atomic charge Z*, shape `(n_species,)`
            dissociation_energy: Bond dissociation energy D, shape `(n_species, n_species)`

        Returns:
            MorseBondParameters with full interaction matrices
        """
        r0 = compute_uff_bond_length(bond_radius, electronegativity)
        Z_i, Z_j = effective_charge[:, None], effective_charge[None, :]

        # Force constant (Eq. 6): k = 664.12 * Z_i * Z_j / r_ij^3
        k = 664.12 * Z_i * Z_j / (r0**3)
        alpha = jnp.sqrt(k / (2.0 * dissociation_energy))

        return cls(
            labels=tuple(map(Label, labels)), r0=r0, D=dissociation_energy, alpha=alpha
        )

from_harmonic(labels, r0, k, D) classmethod

Create Morse parameters from harmonic force constant.

Parameters:

Name Type Description Default
labels tuple[str, ...]

Species labels, shape (n_species,)

 required
r0 Array

Equilibrium bond lengths [Å], shape (n_species, n_species)

 required
k Array

Harmonic force constants [energy/Ų], shape (n_species, n_species)

 required
D Array

Bond dissociation energies, shape (n_species, n_species)

 required

Returns:

Type Description
MorseBondParameters

MorseBondParameters with computed alpha values
Source code in src/kups/potential/classical/morse.py 
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
@classmethod
def from_harmonic(
    cls, labels: tuple[str, ...], r0: Array, k: Array, D: Array
) -> "MorseBondParameters":
    r"""Create Morse parameters from harmonic force constant.

    Args:
        labels: Species labels, shape `(n_species,)`
        r0: Equilibrium bond lengths [Å], shape `(n_species, n_species)`
        k: Harmonic force constants [energy/Ų], shape `(n_species, n_species)`
        D: Bond dissociation energies, shape `(n_species, n_species)`

    Returns:
        MorseBondParameters with computed alpha values
    """
    alpha = jnp.sqrt(k / D)
    return cls(labels=tuple(map(Label, labels)), r0=r0, D=D, alpha=alpha)

from_uff(labels, bond_radius, electronegativity, effective_charge, dissociation_energy) classmethod

Create Morse parameters using UFF bond length/force constant formulas.

Parameters:

Name Type Description Default
labels tuple[str, ...]

Species labels, shape (n_species,)

 required
bond_radius Array

Valence bond radii [Å], shape (n_species,)

 required
electronegativity Array

GMP electronegativity, shape (n_species,)

 required
effective_charge Array

Effective atomic charge Z*, shape (n_species,)

 required
dissociation_energy Array

Bond dissociation energy D, shape (n_species, n_species)

 required

Returns:

Type Description
MorseBondParameters

MorseBondParameters with full interaction matrices
Source code in src/kups/potential/classical/morse.py 
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
@classmethod
def from_uff(
    cls,
    labels: tuple[str, ...],
    bond_radius: Array,
    electronegativity: Array,
    effective_charge: Array,
    dissociation_energy: Array,
) -> "MorseBondParameters":
    r"""Create Morse parameters using UFF bond length/force constant formulas.

    Args:
        labels: Species labels, shape `(n_species,)`
        bond_radius: Valence bond radii [Å], shape `(n_species,)`
        electronegativity: GMP electronegativity, shape `(n_species,)`
        effective_charge: Effective atomic charge Z*, shape `(n_species,)`
        dissociation_energy: Bond dissociation energy D, shape `(n_species, n_species)`

    Returns:
        MorseBondParameters with full interaction matrices
    """
    r0 = compute_uff_bond_length(bond_radius, electronegativity)
    Z_i, Z_j = effective_charge[:, None], effective_charge[None, :]

    # Force constant (Eq. 6): k = 664.12 * Z_i * Z_j / r_ij^3
    k = 664.12 * Z_i * Z_j / (r0**3)
    alpha = jnp.sqrt(k / (2.0 * dissociation_energy))

    return cls(
        labels=tuple(map(Label, labels)), r0=r0, D=dissociation_energy, alpha=alpha
    )

make_morse_bond_potential(particles_view, edge_indices_view, systems_view, parameter_view, probe, gradient_lens, hessian_lens, hessian_idx_view, patch_idx_view=None, out_cache_lens=None)

Create Morse bond potential for explicitly defined bonds.

Parameters:

Name Type Description Default
particles_view View[State, Table[ParticleId, IsBondedParticles]]

Extracts particle data (positions, species) with system index

 required
edge_indices_view View[State, Index[ParticleId]]

Extracts bond connectivity

 required
systems_view View[State, Table[SystemId, HasCell[AnyPeriodicity]]]

Extracts indexed system data (cell)

 required
parameter_view View[State, MorseBondParameters]

Extracts MorseBondParameters

 required
probe Probe[State, P, IsGraphProbe[IsBondedParticles, Literal[2]]] | None

Graph probe for incremental particle and neighbor-list updates

 required
gradient_lens Lens[MorseBondInput, Gradients]

Specifies gradients to compute

 required
hessian_lens Lens[Gradients, Hessians]

Specifies Hessians to compute

 required
hessian_idx_view View[State, Hessians]

Hessian index structure

 required
patch_idx_view View[State, PotentialOut[Gradients, Hessians]] | None

Cached output index structure

 None
out_cache_lens Lens[State, PotentialOut[Gradients, Hessians]] | None

Cache location lens

 None

Returns:

Type Description
Potential[State, Gradients, Hessians, P]

Morse bond Potential
Source code in src/kups/potential/classical/morse.py 
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
def make_morse_bond_potential[
    State,
    P: Patch[Any],
    Gradients,
    Hessians,
](
    particles_view: View[State, Table[ParticleId, IsBondedParticles]],
    edge_indices_view: View[State, Index[ParticleId]],
    systems_view: View[State, Table[SystemId, HasCell[AnyPeriodicity]]],
    parameter_view: View[State, MorseBondParameters],
    probe: Probe[State, P, IsGraphProbe[IsBondedParticles, Literal[2]]] | None,
    gradient_lens: Lens[MorseBondInput, Gradients],
    hessian_lens: Lens[Gradients, Hessians],
    hessian_idx_view: View[State, Hessians],
    patch_idx_view: View[State, PotentialOut[Gradients, Hessians]] | None = None,
    out_cache_lens: Lens[State, PotentialOut[Gradients, Hessians]] | None = None,
) -> Potential[State, Gradients, Hessians, P]:
    """Create Morse bond potential for explicitly defined bonds.

    Args:
        particles_view: Extracts particle data (positions, species) with system index
        edge_indices_view: Extracts bond connectivity
        systems_view: Extracts indexed system data (cell)
        parameter_view: Extracts [MorseBondParameters][kups.potential.classical.morse.MorseBondParameters]
        probe: Graph probe for incremental particle and neighbor-list updates
        gradient_lens: Specifies gradients to compute
        hessian_lens: Specifies Hessians to compute
        hessian_idx_view: Hessian index structure
        patch_idx_view: Cached output index structure
        out_cache_lens: Cache location lens

    Returns:
        Morse bond [Potential][kups.core.potential.Potential]
    """
    graph_fn = GraphConstructor(
        particles=particles_view,
        systems=systems_view,
        neighborlist=lambda state: FixedEdgesNeighborList[Literal[2]](
            edge_indices_view(state)
        ),
        probe=probe,
    )
    composer = LocalGraphSumComposer(
        graph_constructor=graph_fn,
        parameter_view=parameter_view,
    )
    potential = PotentialFromEnergy(
        composer=composer,
        energy_fn=morse_bond_energy,
        gradient_lens=gradient_lens,
        hessian_lens=hessian_lens,
        hessian_idx_view=hessian_idx_view,
        cache_lens=out_cache_lens,
        patch_idx_view=patch_idx_view,
    )
    return potential

morse_bond_energy(inp)

Compute Morse bond energy for all bonds.

Calculates energy as \(D [1 - e^{-\alpha(r - r_0)}]^2\) for each bond.

Parameters:

Name Type Description Default
inp MorseBondInput

Graph potential input with Morse bond parameters

 required

Returns:

Type Description
WithPatch[Table[SystemId, Energy], IdPatch[Any]]

Total bond energy per system
Source code in src/kups/potential/classical/morse.py 
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
def morse_bond_energy(
    inp: MorseBondInput,
) -> WithPatch[Table[SystemId, Energy], IdPatch[Any]]:
    r"""Compute Morse bond energy for all bonds.

    Calculates energy as $D [1 - e^{-\alpha(r - r_0)}]^2$ for each bond.

    Args:
        inp: Graph potential input with Morse bond parameters

    Returns:
        Total bond energy per system
    """
    graph = inp.graph
    assert graph.edges.indices.indices.shape[1] == 2, (
        "Morse bond potential only supports pairwise interactions (order=2)."
    )
    edg_species = graph.particles[graph.edges.indices].labels.indices_in(
        inp.parameters.labels
    )
    r0 = inp.parameters.r0[edg_species[:, 0], edg_species[:, 1]]
    D = inp.parameters.D[edg_species[:, 0], edg_species[:, 1]]
    alpha = inp.parameters.alpha[edg_species[:, 0], edg_species[:, 1]]
    r = jnp.linalg.norm(graph.edge_shifts[:, 0], axis=-1)
    edge_energy = D * (1 - jnp.exp(-alpha * (r - r0))) ** 2
    total_energies = graph.edge_batch_mask.sum_over(edge_energy)
    return WithPatch(total_energies, IdPatch[Any]())