Skip to content

kups.potential.classical.dihedral

Dihedral/torsion potential as defined in UFF.

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

Functional form:

\[ U(\phi) = \frac{1}{2} V_\phi \left[1 - \cos(n \phi_0) \cos(n \phi)\right] \]

where:

  • \(\phi\) is the dihedral angle (angle between planes i-j-k and j-k-l)
  • \(V_\phi\) is the barrier height
  • \(n\) is the periodicity (typically 1, 2, 3, or 6)
  • \(\phi_0\) is the equilibrium dihedral angle

DihedralParameters

UFF dihedral/torsion potential parameters.

Attributes:

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

Species labels, shape (n_species,).
V Array

Barrier heights [energy], shape (n_species, n_species, n_species, n_species).
n Array

Periodicities, shape (n_species, n_species, n_species, n_species).
phi0 Array

Equilibrium dihedral angles [radians], shape (n_species, n_species, n_species, n_species).
Source code in src/kups/potential/classical/dihedral.py 
 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
129
130
131
132
133
134
135
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
163
164
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
@dataclass
class DihedralParameters:
    """UFF dihedral/torsion potential parameters.

    Attributes:
        labels: Species labels, shape `(n_species,)`.
        V: Barrier heights [energy], shape `(n_species, n_species, n_species, n_species)`.
        n: Periodicities, shape `(n_species, n_species, n_species, n_species)`.
        phi0: Equilibrium dihedral angles [radians],
            shape `(n_species, n_species, n_species, n_species)`.
    """

    labels: tuple[Label, ...] = field(static=True)  # (n_species,)
    V: Array  # (n_species, n_species, n_species, n_species)
    n: Array  # (n_species, n_species, n_species, n_species)
    phi0: Array  # (n_species, n_species, n_species, n_species)

    @classmethod
    def from_uff(
        cls,
        labels: tuple[str, ...],
        bond_angle: Array,
        torsion_sp3: Array,
        torsion_sp2: Array,
        group: Array,
        bond_order: Array | None = None,
        hybridization_tol: float = 5.0 * 3.141592653589793 / 180.0,
    ) -> "DihedralParameters":
        r"""Create dihedral parameters using UFF formulas.

        Computes torsion parameters based on central bond hybridization following
        the UFF paper (Rappe et al. 1992). For dihedral i-j-k-l, parameters depend
        on the central j-k bond and neighboring atoms.

        Args:
            labels: Species labels, shape `(n_species,)`.
            bond_angle: Natural valence angle [radians], shape `(n_species,)`.
            torsion_sp3: sp3 torsional barrier V [kcal/mol], shape `(n_species,)`.
            torsion_sp2: sp2 torsional parameter U [kcal/mol], shape `(n_species,)`.
            group: Periodic table group (1-18), shape `(n_species,)`.
            bond_order: Bond order for j-k bond, shape `(n_species, n_species)`.
            hybridization_tol: Tolerance for hybridization detection [radians].

        Returns:
            DihedralParameters with full 4D interaction matrices.
        """
        nt = bond_angle.shape[0]

        # Determine hybridization from bond angle
        is_sp2 = jnp.abs(bond_angle - jnp.radians(120.0)) < hybridization_tol
        is_sp3 = ~(jnp.abs(bond_angle - jnp.pi) < hybridization_tol) & ~is_sp2

        # Group 6 and main group detection
        is_group6 = group == 6
        is_main_group = group > 0

        # Broadcast to 4D for dihedral i-j-k-l
        is_sp2_i = is_sp2[:, None, None, None]
        is_sp3_j = is_sp3[None, :, None, None]
        is_sp3_k = is_sp3[None, None, :, None]
        is_sp2_j = is_sp2[None, :, None, None]
        is_sp2_k = is_sp2[None, None, :, None]
        is_sp2_l = is_sp2[None, None, None, :]

        is_group6_j = is_group6[None, :, None, None]
        is_group6_k = is_group6[None, None, :, None]
        is_main_j = is_main_group[None, :, None, None]
        is_main_k = is_main_group[None, None, :, None]

        # Torsional parameters broadcast to 4D
        V_j = torsion_sp3[None, :, None, None]
        V_k = torsion_sp3[None, None, :, None]
        U_j = torsion_sp2[None, :, None, None]
        U_k = torsion_sp2[None, None, :, None]

        # Bond order for central j-k bond
        if bond_order is None:
            BO_jk = jnp.ones((1, nt, nt, 1))
        else:
            BO_jk = bond_order[None, :, :, None]

        # V from Eq.17: V = 5 * sqrt(U_j * U_k) * (1 + 4.18 * ln(BO))
        V_eq17 = (
            5.0
            * jnp.sqrt(jnp.abs(U_j * U_k))
            * (1.0 + 4.18 * jnp.log(jnp.maximum(BO_jk, 1e-10)))
        )

        # Initialize arrays
        V_arr = jnp.zeros((nt, nt, nt, nt))
        n_arr = jnp.ones((nt, nt, nt, nt), dtype=jnp.int32)
        phi0_arr = jnp.zeros((nt, nt, nt, nt))

        # === Case 1: sp3-sp3 general ===
        sp3_sp3 = is_sp3_j & is_sp3_k
        V_arr = jnp.where(sp3_sp3, jnp.sqrt(jnp.abs(V_j * V_k)), V_arr)
        n_arr = jnp.where(sp3_sp3, 3, n_arr)
        phi0_arr = jnp.where(sp3_sp3, jnp.pi, phi0_arr)

        # === Case 2: sp3-sp3 Group 6 pair (override n and phi0) ===
        sp3_sp3_group6 = sp3_sp3 & is_group6_j & is_group6_k
        n_arr = jnp.where(sp3_sp3_group6, 2, n_arr)
        phi0_arr = jnp.where(sp3_sp3_group6, jnp.pi / 2, phi0_arr)

        # === Case 3: sp2-sp2 ===
        sp2_sp2 = is_sp2_j & is_sp2_k
        V_arr = jnp.where(sp2_sp2, V_eq17, V_arr)
        n_arr = jnp.where(sp2_sp2, 2, n_arr)
        phi0_arr = jnp.where(sp2_sp2, jnp.pi, phi0_arr)

        # === Case 4: sp3-sp2 general ===
        sp3_sp2 = (is_sp3_j & is_sp2_k) | (is_sp2_j & is_sp3_k)
        V_arr = jnp.where(sp3_sp2, 1.0, V_arr)
        n_arr = jnp.where(sp3_sp2, 6, n_arr)
        phi0_arr = jnp.where(sp3_sp2, 0.0, phi0_arr)

        # === Case 5: sp3-sp2 propene-like (sp2 bonded to another sp2) ===
        propene = (is_sp2_j & is_sp3_k & is_sp2_i) | (is_sp3_j & is_sp2_k & is_sp2_l)
        V_arr = jnp.where(propene, 2.0, V_arr)
        n_arr = jnp.where(propene, 3, n_arr)
        phi0_arr = jnp.where(propene, jnp.pi, phi0_arr)

        # === Case 6: sp3-sp2 with Group 6 sp3 ===
        group6_sp3_sp2 = (is_sp3_j & is_group6_j & is_sp2_k) | (
            is_sp2_j & is_sp3_k & is_group6_k
        )
        V_arr = jnp.where(group6_sp3_sp2, V_eq17, V_arr)
        n_arr = jnp.where(group6_sp3_sp2, 2, n_arr)
        phi0_arr = jnp.where(group6_sp3_sp2, jnp.pi / 2, phi0_arr)

        # === Case 7: Non-main-group (V=0) ===
        non_main = ~is_main_j | ~is_main_k
        V_arr = jnp.where(non_main, 0.0, V_arr)

        return cls(labels=tuple(map(Label, labels)), V=V_arr, n=n_arr, phi0=phi0_arr)

from_uff(labels, bond_angle, torsion_sp3, torsion_sp2, group, bond_order=None, hybridization_tol=5.0 * 3.141592653589793 / 180.0) classmethod

Create dihedral parameters using UFF formulas.

Computes torsion parameters based on central bond hybridization following the UFF paper (Rappe et al. 1992). For dihedral i-j-k-l, parameters depend on the central j-k bond and neighboring atoms.

Parameters:

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

Species labels, shape (n_species,).

 required
bond_angle Array

Natural valence angle [radians], shape (n_species,).

 required
torsion_sp3 Array

sp3 torsional barrier V [kcal/mol], shape (n_species,).

 required
torsion_sp2 Array

sp2 torsional parameter U [kcal/mol], shape (n_species,).

 required
group Array

Periodic table group (1-18), shape (n_species,).

 required
bond_order Array | None

Bond order for j-k bond, shape (n_species, n_species).

 None
hybridization_tol float

Tolerance for hybridization detection [radians].

 5.0 * 3.141592653589793 / 180.0

Returns:

Type Description
DihedralParameters

DihedralParameters with full 4D interaction matrices.
Source code in src/kups/potential/classical/dihedral.py 
 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
129
130
131
132
133
134
135
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
163
164
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
@classmethod
def from_uff(
    cls,
    labels: tuple[str, ...],
    bond_angle: Array,
    torsion_sp3: Array,
    torsion_sp2: Array,
    group: Array,
    bond_order: Array | None = None,
    hybridization_tol: float = 5.0 * 3.141592653589793 / 180.0,
) -> "DihedralParameters":
    r"""Create dihedral parameters using UFF formulas.

    Computes torsion parameters based on central bond hybridization following
    the UFF paper (Rappe et al. 1992). For dihedral i-j-k-l, parameters depend
    on the central j-k bond and neighboring atoms.

    Args:
        labels: Species labels, shape `(n_species,)`.
        bond_angle: Natural valence angle [radians], shape `(n_species,)`.
        torsion_sp3: sp3 torsional barrier V [kcal/mol], shape `(n_species,)`.
        torsion_sp2: sp2 torsional parameter U [kcal/mol], shape `(n_species,)`.
        group: Periodic table group (1-18), shape `(n_species,)`.
        bond_order: Bond order for j-k bond, shape `(n_species, n_species)`.
        hybridization_tol: Tolerance for hybridization detection [radians].

    Returns:
        DihedralParameters with full 4D interaction matrices.
    """
    nt = bond_angle.shape[0]

    # Determine hybridization from bond angle
    is_sp2 = jnp.abs(bond_angle - jnp.radians(120.0)) < hybridization_tol
    is_sp3 = ~(jnp.abs(bond_angle - jnp.pi) < hybridization_tol) & ~is_sp2

    # Group 6 and main group detection
    is_group6 = group == 6
    is_main_group = group > 0

    # Broadcast to 4D for dihedral i-j-k-l
    is_sp2_i = is_sp2[:, None, None, None]
    is_sp3_j = is_sp3[None, :, None, None]
    is_sp3_k = is_sp3[None, None, :, None]
    is_sp2_j = is_sp2[None, :, None, None]
    is_sp2_k = is_sp2[None, None, :, None]
    is_sp2_l = is_sp2[None, None, None, :]

    is_group6_j = is_group6[None, :, None, None]
    is_group6_k = is_group6[None, None, :, None]
    is_main_j = is_main_group[None, :, None, None]
    is_main_k = is_main_group[None, None, :, None]

    # Torsional parameters broadcast to 4D
    V_j = torsion_sp3[None, :, None, None]
    V_k = torsion_sp3[None, None, :, None]
    U_j = torsion_sp2[None, :, None, None]
    U_k = torsion_sp2[None, None, :, None]

    # Bond order for central j-k bond
    if bond_order is None:
        BO_jk = jnp.ones((1, nt, nt, 1))
    else:
        BO_jk = bond_order[None, :, :, None]

    # V from Eq.17: V = 5 * sqrt(U_j * U_k) * (1 + 4.18 * ln(BO))
    V_eq17 = (
        5.0
        * jnp.sqrt(jnp.abs(U_j * U_k))
        * (1.0 + 4.18 * jnp.log(jnp.maximum(BO_jk, 1e-10)))
    )

    # Initialize arrays
    V_arr = jnp.zeros((nt, nt, nt, nt))
    n_arr = jnp.ones((nt, nt, nt, nt), dtype=jnp.int32)
    phi0_arr = jnp.zeros((nt, nt, nt, nt))

    # === Case 1: sp3-sp3 general ===
    sp3_sp3 = is_sp3_j & is_sp3_k
    V_arr = jnp.where(sp3_sp3, jnp.sqrt(jnp.abs(V_j * V_k)), V_arr)
    n_arr = jnp.where(sp3_sp3, 3, n_arr)
    phi0_arr = jnp.where(sp3_sp3, jnp.pi, phi0_arr)

    # === Case 2: sp3-sp3 Group 6 pair (override n and phi0) ===
    sp3_sp3_group6 = sp3_sp3 & is_group6_j & is_group6_k
    n_arr = jnp.where(sp3_sp3_group6, 2, n_arr)
    phi0_arr = jnp.where(sp3_sp3_group6, jnp.pi / 2, phi0_arr)

    # === Case 3: sp2-sp2 ===
    sp2_sp2 = is_sp2_j & is_sp2_k
    V_arr = jnp.where(sp2_sp2, V_eq17, V_arr)
    n_arr = jnp.where(sp2_sp2, 2, n_arr)
    phi0_arr = jnp.where(sp2_sp2, jnp.pi, phi0_arr)

    # === Case 4: sp3-sp2 general ===
    sp3_sp2 = (is_sp3_j & is_sp2_k) | (is_sp2_j & is_sp3_k)
    V_arr = jnp.where(sp3_sp2, 1.0, V_arr)
    n_arr = jnp.where(sp3_sp2, 6, n_arr)
    phi0_arr = jnp.where(sp3_sp2, 0.0, phi0_arr)

    # === Case 5: sp3-sp2 propene-like (sp2 bonded to another sp2) ===
    propene = (is_sp2_j & is_sp3_k & is_sp2_i) | (is_sp3_j & is_sp2_k & is_sp2_l)
    V_arr = jnp.where(propene, 2.0, V_arr)
    n_arr = jnp.where(propene, 3, n_arr)
    phi0_arr = jnp.where(propene, jnp.pi, phi0_arr)

    # === Case 6: sp3-sp2 with Group 6 sp3 ===
    group6_sp3_sp2 = (is_sp3_j & is_group6_j & is_sp2_k) | (
        is_sp2_j & is_sp3_k & is_group6_k
    )
    V_arr = jnp.where(group6_sp3_sp2, V_eq17, V_arr)
    n_arr = jnp.where(group6_sp3_sp2, 2, n_arr)
    phi0_arr = jnp.where(group6_sp3_sp2, jnp.pi / 2, phi0_arr)

    # === Case 7: Non-main-group (V=0) ===
    non_main = ~is_main_j | ~is_main_k
    V_arr = jnp.where(non_main, 0.0, V_arr)

    return cls(labels=tuple(map(Label, labels)), V=V_arr, n=n_arr, phi0=phi0_arr)

IsBondedParticles

Bases: HasPositionsAndLabels, IsRadiusGraphPoints, Protocol

Particle data with positions, labels, and system index.

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

    ...

dihedral_energy(inp)

Compute UFF dihedral/torsion energy for all dihedrals.

Parameters:

Name Type Description Default
inp DihedralInput

Graph potential input with dihedral parameters

 required

Returns:

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

Total dihedral energy per system
Source code in src/kups/potential/classical/dihedral.py 
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
def dihedral_energy(
    inp: DihedralInput,
) -> WithPatch[Table[SystemId, Energy], IdPatch[Any]]:
    r"""Compute UFF dihedral/torsion energy for all dihedrals.

    Args:
        inp: Graph potential input with dihedral parameters

    Returns:
        Total dihedral energy per system
    """
    graph = inp.graph
    assert graph.edges.indices.indices.shape[1] == 4, (
        "Dihedral potential only supports quadruplet interactions (order=4)."
    )
    edg = graph.particles[graph.edges.indices]
    edg_species = edg.labels.indices_in(inp.parameters.labels)
    s0, s1, s2, s3 = (
        edg_species[:, 0],
        edg_species[:, 1],
        edg_species[:, 2],
        edg_species[:, 3],
    )

    V = inp.parameters.V[s0, s1, s2, s3]
    n = inp.parameters.n[s0, s1, s2, s3]
    phi0 = inp.parameters.phi0[s0, s1, s2, s3]

    # Edge shifts: all relative to first atom (i); convert to bond vectors
    r_ij = graph.edge_shifts[:, 0]  # j - i
    r_jk = graph.edge_shifts[:, 1] - r_ij  # k - j
    r_kl = graph.edge_shifts[:, 2] - graph.edge_shifts[:, 1]  # l - k

    # Compute normal vectors to planes
    n1 = jnp.cross(r_ij, r_jk)  # Normal to plane i-j-k
    n2 = jnp.cross(r_jk, r_kl)  # Normal to plane j-k-l

    # Normalize the normal vectors
    n1_norm = jnp.linalg.norm(n1, axis=-1, keepdims=True)
    n2_norm = jnp.linalg.norm(n2, axis=-1, keepdims=True)

    eps = 1e-10
    n1_normalized = n1 / (n1_norm + eps)
    n2_normalized = n2 / (n2_norm + eps)

    # Compute dihedral angle using atan2 for proper sign
    cos_phi = jnp.sum(n1_normalized * n2_normalized, axis=-1)
    r_jk_norm = jnp.linalg.norm(r_jk, axis=-1, keepdims=True)
    r_jk_normalized = r_jk / (r_jk_norm + eps)
    sin_phi = jnp.sum(
        jnp.cross(n1_normalized, n2_normalized) * r_jk_normalized, axis=-1
    )

    phi = jnp.arctan2(sin_phi, cos_phi)

    # UFF torsion energy: U(phi) = 1/2 V [1 - cos(n*phi0)*cos(n*phi)]
    edge_energy = 0.5 * V * (1.0 - jnp.cos(n * phi0) * jnp.cos(n * phi))

    total_energies = graph.edge_batch_mask.sum_over(edge_energy)
    return WithPatch(total_energies, IdPatch[Any]())

make_dihedral_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 UFF dihedral potential for explicitly defined dihedrals.

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 dihedral connectivity (quadruplets)

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

Extracts indexed system data (cell)

 required
parameter_view View[State, DihedralParameters]

Extracts DihedralParameters

 required
probe Probe[State, Ptch, IsGraphProbe[IsBondedParticles, Literal[4]]] | None

Graph probe for incremental particle and neighbor-list updates

 required
gradient_lens Lens[DihedralInput, 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, Ptch]

UFF dihedral Potential
Source code in src/kups/potential/classical/dihedral.py 
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
def make_dihedral_potential[
    State,
    Ptch: 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, DihedralParameters],
    probe: Probe[State, Ptch, IsGraphProbe[IsBondedParticles, Literal[4]]] | None,
    gradient_lens: Lens[DihedralInput, 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, Ptch]:
    """Create UFF dihedral potential for explicitly defined dihedrals.

    Args:
        particles_view: Extracts particle data (positions, species) with system index
        edge_indices_view: Extracts dihedral connectivity (quadruplets)
        systems_view: Extracts indexed system data (cell)
        parameter_view: Extracts [DihedralParameters][kups.potential.classical.dihedral.DihedralParameters]
        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:
        UFF dihedral [Potential][kups.core.potential.Potential]
    """
    graph_fn = GraphConstructor(
        particles=particles_view,
        systems=systems_view,
        neighborlist=lambda state: FixedEdgesNeighborList[Literal[4]](
            edge_indices_view(state)
        ),
        probe=probe,
    )
    composer = LocalGraphSumComposer(
        graph_constructor=graph_fn,
        parameter_view=parameter_view,
    )
    potential = PotentialFromEnergy(
        composer=composer,
        energy_fn=dihedral_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