Skip to content

kups.potential.classical.lennard_jones

Lennard-Jones potential implementations with tail corrections.

This module provides the Lennard-Jones 12-6 potential for van der Waals interactions:

\[ U(r) = 4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6\right] \]

Includes variants with smooth tail corrections and analytical long-range corrections for periodic systems. Supports Lorentz-Berthelot mixing rules for multi-component systems.

GlobalTailCorrectedLennardJonesParameters

Bases: LennardJonesParameters

Lennard-Jones parameters with analytical long-range correction.

Attributes:

Name Type Description
tail_corrected Array

Enable correction per species pair, shape (n_species, n_species).
Source code in src/kups/potential/classical/lennard_jones.py 
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
@dataclass
class GlobalTailCorrectedLennardJonesParameters(LennardJonesParameters):
    """Lennard-Jones parameters with analytical long-range correction.

    Attributes:
        tail_corrected: Enable correction per species pair, shape ``(n_species, n_species)``.
    """

    tail_corrected: Array  # (n_species, n_species) bool

    @classmethod
    @override
    def from_dict(
        cls,
        cutoff: float | Array,
        parameters: dict[str, tuple[float | None, float | None]],
        mixing_rule: MixingRule,
        tail_correction: bool = True,
    ) -> GlobalTailCorrectedLennardJonesParameters:
        """Create tail-corrected parameters from a dict of per-species values.

        Args:
            cutoff: Cutoff radius [Angstrom].
            parameters: Map from species label to ``(sigma, epsilon)`` pair.
                ``None`` values default to ``sigma=1.0``, ``epsilon=0.0``.
            mixing_rule: Combining rule for cross-species interactions.
            tail_correction: Whether to enable tail corrections for all
                non-zero epsilon pairs.
        """
        base = LennardJonesParameters.from_dict(cutoff, parameters, mixing_rule)
        mask = (
            base.epsilon > 0
            if tail_correction
            else jnp.zeros_like(base.epsilon, dtype=bool)
        )
        return cls(
            labels=base.labels,
            sigma=base.sigma,
            epsilon=base.epsilon,
            cutoff=base.cutoff,
            tail_corrected=mask,
        )

from_dict(cutoff, parameters, mixing_rule, tail_correction=True) classmethod

Create tail-corrected parameters from a dict of per-species values.

Parameters:

Name Type Description Default
cutoff float | Array

Cutoff radius [Angstrom].

 required
parameters dict[str, tuple[float | None, float | None]]

Map from species label to (sigma, epsilon) pair. None values default to sigma=1.0, epsilon=0.0.

 required
mixing_rule MixingRule

Combining rule for cross-species interactions.

 required
tail_correction bool

Whether to enable tail corrections for all non-zero epsilon pairs.

 True
Source code in src/kups/potential/classical/lennard_jones.py 
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
@classmethod
@override
def from_dict(
    cls,
    cutoff: float | Array,
    parameters: dict[str, tuple[float | None, float | None]],
    mixing_rule: MixingRule,
    tail_correction: bool = True,
) -> GlobalTailCorrectedLennardJonesParameters:
    """Create tail-corrected parameters from a dict of per-species values.

    Args:
        cutoff: Cutoff radius [Angstrom].
        parameters: Map from species label to ``(sigma, epsilon)`` pair.
            ``None`` values default to ``sigma=1.0``, ``epsilon=0.0``.
        mixing_rule: Combining rule for cross-species interactions.
        tail_correction: Whether to enable tail corrections for all
            non-zero epsilon pairs.
    """
    base = LennardJonesParameters.from_dict(cutoff, parameters, mixing_rule)
    mask = (
        base.epsilon > 0
        if tail_correction
        else jnp.zeros_like(base.epsilon, dtype=bool)
    )
    return cls(
        labels=base.labels,
        sigma=base.sigma,
        epsilon=base.epsilon,
        cutoff=base.cutoff,
        tail_corrected=mask,
    )

LennardJonesParameters

Lennard-Jones potential parameters.

Attributes:

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

Species labels as Index.
sigma Array

Length scale parameters [Ã…], shape (n_species, n_species).
epsilon Array

Energy well depths [energy units], shape (n_species, n_species).
cutoff Table[SystemId, Array]

Cutoff radius [Ã…], shape (n_systems,).
Source code in src/kups/potential/classical/lennard_jones.py 
 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
@dataclass
class LennardJonesParameters:
    """Lennard-Jones potential parameters.

    Attributes:
        labels: Species labels as ``Index``.
        sigma: Length scale parameters [Ã…], shape ``(n_species, n_species)``.
        epsilon: Energy well depths [energy units], shape ``(n_species, n_species)``.
        cutoff: Cutoff radius [Ã…], shape ``(n_systems,)``.
    """

    labels: tuple[Label, ...] = field(static=True)  # (n_species,)
    sigma: Array  # (n_species, n_species) float
    epsilon: Array  # (n_species, n_species) float
    cutoff: Table[SystemId, Array]  # (n_graphs,) float

    @classmethod
    def from_dict(
        cls,
        cutoff: float | Array,
        parameters: dict[str, tuple[float | None, float | None]],
        mixing_rule: MixingRule,
    ) -> LennardJonesParameters:
        """Create parameters from a dict of per-species values.

        Args:
            cutoff: Cutoff radius [Angstrom].
            parameters: Map from species label to ``(sigma, epsilon)`` pair.
                ``None`` values default to ``sigma=1.0``, ``epsilon=0.0``.
            mixing_rule: Combining rule for cross-species interactions.
        """
        labels = tuple(parameters.keys())
        raw = [(s or 1.0, e or 0.0) for s, e in parameters.values()]
        sigma, epsilon = jnp.asarray(raw).T
        cutoff_indexed = Table((SystemId(0),), jnp.array([cutoff]))
        match mixing_rule:
            case "lorentz_berthelot":
                return cls.from_lorentz_berthelot_mixing(
                    labels, sigma, epsilon, cutoff_indexed
                )
            case _ as unreachable:
                assert_never(unreachable)

    @classmethod
    def from_lorentz_berthelot_mixing(
        cls,
        labels: tuple[str, ...],
        sigma: Array,
        epsilon: Array,
        cutoff: Table[SystemId, Array],
    ) -> LennardJonesParameters:
        """Create parameters using Lorentz-Berthelot mixing rules.

        - σᵢⱼ = (σᵢ + σⱼ) / 2 (arithmetic mean)
        - εᵢⱼ = √(εᵢ × εⱼ) (geometric mean)
        """
        assert sigma.ndim == epsilon.ndim == 1
        sigma_matrix = (sigma[:, None] + sigma) / 2
        epsilon_matrix = jnp.sqrt(epsilon[:, None] * epsilon)
        return cls(tuple(map(Label, labels)), sigma_matrix, epsilon_matrix, cutoff)

from_dict(cutoff, parameters, mixing_rule) classmethod

Create parameters from a dict of per-species values.

Parameters:

Name Type Description Default
cutoff float | Array

Cutoff radius [Angstrom].

 required
parameters dict[str, tuple[float | None, float | None]]

Map from species label to (sigma, epsilon) pair. None values default to sigma=1.0, epsilon=0.0.

 required
mixing_rule MixingRule

Combining rule for cross-species interactions.

 required
Source code in src/kups/potential/classical/lennard_jones.py 
 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
@classmethod
def from_dict(
    cls,
    cutoff: float | Array,
    parameters: dict[str, tuple[float | None, float | None]],
    mixing_rule: MixingRule,
) -> LennardJonesParameters:
    """Create parameters from a dict of per-species values.

    Args:
        cutoff: Cutoff radius [Angstrom].
        parameters: Map from species label to ``(sigma, epsilon)`` pair.
            ``None`` values default to ``sigma=1.0``, ``epsilon=0.0``.
        mixing_rule: Combining rule for cross-species interactions.
    """
    labels = tuple(parameters.keys())
    raw = [(s or 1.0, e or 0.0) for s, e in parameters.values()]
    sigma, epsilon = jnp.asarray(raw).T
    cutoff_indexed = Table((SystemId(0),), jnp.array([cutoff]))
    match mixing_rule:
        case "lorentz_berthelot":
            return cls.from_lorentz_berthelot_mixing(
                labels, sigma, epsilon, cutoff_indexed
            )
        case _ as unreachable:
            assert_never(unreachable)

from_lorentz_berthelot_mixing(labels, sigma, epsilon, cutoff) classmethod

Create parameters using Lorentz-Berthelot mixing rules.

  • σᵢⱼ = (σᵢ + σⱼ) / 2 (arithmetic mean)
  • εᵢⱼ = √(εᵢ × εⱼ) (geometric mean)
Source code in src/kups/potential/classical/lennard_jones.py 
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
@classmethod
def from_lorentz_berthelot_mixing(
    cls,
    labels: tuple[str, ...],
    sigma: Array,
    epsilon: Array,
    cutoff: Table[SystemId, Array],
) -> LennardJonesParameters:
    """Create parameters using Lorentz-Berthelot mixing rules.

    - σᵢⱼ = (σᵢ + σⱼ) / 2 (arithmetic mean)
    - εᵢⱼ = √(εᵢ × εⱼ) (geometric mean)
    """
    assert sigma.ndim == epsilon.ndim == 1
    sigma_matrix = (sigma[:, None] + sigma) / 2
    epsilon_matrix = jnp.sqrt(epsilon[:, None] * epsilon)
    return cls(tuple(map(Label, labels)), sigma_matrix, epsilon_matrix, cutoff)

PairTailCorrectedLennardJonesParameters

Bases: LennardJonesParameters

Lennard-Jones parameters with smooth pairwise tail correction.

Attributes:

Name Type Description
truncation_radius Table[SystemId, Array]

Radius where smoothing begins [Ã…], shape (n_systems,).
Source code in src/kups/potential/classical/lennard_jones.py 
176
177
178
179
180
181
182
183
184
@dataclass
class PairTailCorrectedLennardJonesParameters(LennardJonesParameters):
    """Lennard-Jones parameters with smooth pairwise tail correction.

    Attributes:
        truncation_radius: Radius where smoothing begins [Ã…], shape ``(n_systems,)``.
    """

    truncation_radius: Table[SystemId, Array]  # (n_graphs,)

global_lennard_jones_tail_correction_energy(inp)

Compute analytical long-range tail correction energy.

Source code in src/kups/potential/classical/lennard_jones.py 
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
@jit
def global_lennard_jones_tail_correction_energy(
    inp: GlobalTailCorrectedLennardJonesInput,
) -> WithPatch[Table[SystemId, Energy], IdPatch[Any]]:
    """Compute analytical long-range tail correction energy."""
    density, _volume, term1, term2, tail_mask, n_graphs = (
        _global_tail_correction_common(inp)
    )
    sigma = inp.parameters.sigma
    epsilon = inp.parameters.epsilon
    result = (8 / 3) * jnp.pi * density * epsilon * sigma**3 * (term2 / 3 - term1)
    result *= tail_mask
    total_energies = Table.arange(result.sum(axis=(1, 2)), label=SystemId)
    assert len(total_energies) == n_graphs
    return WithPatch(total_energies, IdPatch[Any]())

global_lennard_jones_tail_correction_pressure(inp)

Compute analytical long-range tail correction for pressure.

Source code in src/kups/potential/classical/lennard_jones.py 
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
@jit
def global_lennard_jones_tail_correction_pressure(
    inp: GlobalTailCorrectedLennardJonesInput,
) -> WithPatch[Table[SystemId, Energy], IdPatch[Any]]:
    """Compute analytical long-range tail correction for pressure."""
    density, volume, term1, term2, tail_mask, n_graphs = _global_tail_correction_common(
        inp
    )
    sigma = inp.parameters.sigma
    epsilon = inp.parameters.epsilon
    result = (
        (16 / 3)
        * jnp.pi
        * density
        / volume
        * epsilon
        * sigma**3
        * (term2 / 3 * 2 - term1)
    )
    result *= tail_mask
    total_pressure = Table.arange(result.sum(axis=(1, 2)), label=SystemId)
    assert len(total_pressure) == n_graphs
    return WithPatch(total_pressure, IdPatch[Any]())

lennard_jones_edge_energy(inp)

Compute Lennard-Jones energy per edge.

Source code in src/kups/potential/classical/lennard_jones.py 
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
@jit
def lennard_jones_edge_energy(inp: LennardJonesInput) -> Array:
    """Compute Lennard-Jones energy per edge."""
    graph = inp.graph
    assert graph.edges.indices.shape[1] == 2
    sigma = inp.parameters.sigma
    epsilon = inp.parameters.epsilon
    assert sigma.ndim == 2 and sigma.shape[0] == sigma.shape[1]
    assert epsilon.ndim == 2 and epsilon.shape[0] == epsilon.shape[1]
    edg_species = graph.particles[graph.edges.indices].labels.indices_in(
        inp.parameters.labels
    )
    epsilon = epsilon[edg_species[:, 0], edg_species[:, 1]]
    sigma = sigma[edg_species[:, 0], edg_species[:, 1]]
    r2 = jnp.sum(graph.edge_shifts[:, 0] ** 2, axis=-1)
    c6 = (sigma**2 / r2) ** 3
    edge_energy = 4 * epsilon * (c6**2 - c6)
    batch = graph.edge_batch_mask.indices
    mask = r2 < jnp.pow(inp.parameters.cutoff.data, 2)[batch]
    return edge_energy * mask

lennard_jones_energy(inp)

Compute total Lennard-Jones energy per system.

Source code in src/kups/potential/classical/lennard_jones.py 
166
167
168
169
170
171
172
173
def lennard_jones_energy(
    inp: LennardJonesInput,
) -> WithPatch[Table[SystemId, Energy], IdPatch[Any]]:
    """Compute total Lennard-Jones energy per system."""
    graph = inp.graph
    edge_energy = lennard_jones_edge_energy(inp)
    total_energies = graph.edge_batch_mask.sum_over(edge_energy) / 2
    return WithPatch(total_energies, IdPatch[Any]())

make_global_lennard_jones_tail_correction_potential(particles_view, systems_view, parameter_view, gradient_lens, hessian_lens, hessian_idx_view, patch_idx_view=None, out_cache_lens=None)

Create analytical long-range tail correction for Lennard-Jones potential.

Source code in src/kups/potential/classical/lennard_jones.py 
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
def make_global_lennard_jones_tail_correction_potential[State, Gradients, Hessians](
    particles_view: View[State, Table[ParticleId, IsLJGraphParticles]],
    systems_view: View[State, Table[SystemId, HasCell[AnyPeriodicity]]],
    parameter_view: View[State, GlobalTailCorrectedLennardJonesParameters],
    gradient_lens: Lens[GCLJInp, 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, Patch[State]]:
    """Create analytical long-range tail correction for Lennard-Jones potential."""
    return PotentialFromEnergy(
        energy_fn=global_lennard_jones_tail_correction_energy,
        composer=FullGraphSumComposer(
            GraphConstructor(
                particles=particles_view,
                systems=systems_view,
                neighborlist=lambda _: EmptyNeighborList[Literal[0]](),
                probe=None,
            ),
            parameter_view=parameter_view,
        ),
        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,
    )

make_global_lennard_jones_tail_correction_pressure(particles_view, systems_view, parameter_view)

Create long-range pressure correction for Lennard-Jones systems.

Source code in src/kups/potential/classical/lennard_jones.py 
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
def make_global_lennard_jones_tail_correction_pressure[State](
    particles_view: View[State, Table[ParticleId, IsLJGraphParticles]],
    systems_view: View[State, Table[SystemId, HasCell[AnyPeriodicity]]],
    parameter_view: View[State, GlobalTailCorrectedLennardJonesParameters],
) -> StateProperty[State, Table[SystemId, Array]]:
    """Create long-range pressure correction for Lennard-Jones systems."""
    graph_constructor = GraphConstructor(
        particles=particles_view,
        systems=systems_view,
        neighborlist=lambda _: EmptyNeighborList[Literal[0]](),
        probe=None,
    )

    def pressure(key: Array, state: State) -> Table[SystemId, Array]:
        del key
        params = parameter_view(state)
        graph = graph_constructor(state, None)
        return global_lennard_jones_tail_correction_pressure(
            GraphPotentialInput(params, graph)
        ).data

    return pressure

make_lennard_jones_potential(particles_view, systems_view, neighborlist_view, parameter_view, probe, gradient_lens, hessian_lens, hessian_idx_view, patch_idx_view=None, out_cache_lens=None)

Create a standard Lennard-Jones potential with sharp cutoff.

Source code in src/kups/potential/classical/lennard_jones.py 
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
def make_lennard_jones_potential[
    State,
    Ptch: Patch[Any],
    Gradients,
    Hessians,
](
    particles_view: View[State, Table[ParticleId, IsLJGraphParticles]],
    systems_view: View[State, Table[SystemId, HasCell[AnyPeriodicity]]],
    neighborlist_view: View[State, NeighborList[Literal[2]]],
    parameter_view: View[State, LennardJonesParameters],
    probe: Probe[State, Ptch, IsGraphProbe[IsLJGraphParticles, Literal[2]]] | None,
    gradient_lens: Lens[LJRadiusInp, 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 a standard Lennard-Jones potential with sharp cutoff."""
    graph_fn = GraphConstructor(
        particles=particles_view,
        systems=systems_view,
        neighborlist=neighborlist_view,
        probe=probe,
    )
    composer = LocalGraphSumComposer(
        graph_constructor=graph_fn,
        parameter_view=parameter_view,
    )
    return PotentialFromEnergy(
        composer=composer,
        energy_fn=lennard_jones_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,
    )

make_pair_tail_corrected_lennard_jones_potential(particles_view, systems_view, neighborlist_view, parameter_view, probe, gradient_lens, hessian_lens, hessian_idx_view, patch_idx_view=None, out_cache_lens=None)

Create a Lennard-Jones potential with smooth pairwise tail correction.

Source code in src/kups/potential/classical/lennard_jones.py 
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
def make_pair_tail_corrected_lennard_jones_potential[
    State,
    Ptch: Patch[Any],
    Gradients,
    Hessians,
](
    particles_view: View[State, Table[ParticleId, IsLJGraphParticles]],
    systems_view: View[State, Table[SystemId, HasCell[AnyPeriodicity]]],
    neighborlist_view: View[State, NeighborList[Literal[2]]],
    parameter_view: View[State, PairTailCorrectedLennardJonesParameters],
    probe: Probe[State, Ptch, IsGraphProbe[IsLJGraphParticles, Literal[2]]] | None,
    gradient_lens: Lens[PCLJInp, 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 a Lennard-Jones potential with smooth pairwise tail correction."""
    radius_graph_fn = GraphConstructor(
        particles=particles_view,
        systems=systems_view,
        neighborlist=neighborlist_view,
        probe=probe,
    )
    composer = LocalGraphSumComposer(
        graph_constructor=radius_graph_fn,
        parameter_view=parameter_view,
    )
    return PotentialFromEnergy(
        composer=composer,
        energy_fn=pair_tail_corrected_lennard_jones_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,
    )

pair_tail_corrected_lennard_jones_energy(inp)

Compute Lennard-Jones energy with smooth pairwise tail correction.

Source code in src/kups/potential/classical/lennard_jones.py 
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
@jit
def pair_tail_corrected_lennard_jones_energy(
    inp: PairTailCorrectedLennardJonesInput,
) -> WithPatch[Table[SystemId, Energy], IdPatch[Any]]:
    """Compute Lennard-Jones energy with smooth pairwise tail correction."""
    graph = inp.graph
    r: Array = jnp.linalg.norm(graph.edge_shifts, axis=(-2, -1))
    edge_energy = lennard_jones_edge_energy(cast(LennardJonesInput, inp))

    batch = graph.edge_batch_mask
    r_tr = inp.parameters.truncation_radius[batch]
    r_cut = inp.parameters.cutoff[batch]
    mask = r > r_tr
    remove = r >= r_cut
    factor1 = ((r_cut**2) - r**2) ** 2
    factor2 = 2 * r**2 + (r_cut**2 - 3 * r_tr**2)
    div = (r_cut**2 - r_tr**2) ** 3
    corrected_edge_energy = jnp.where(
        mask, edge_energy * factor1 * factor2 / div, edge_energy
    )
    corrected_edge_energy = jnp.where(remove, 0.0, corrected_edge_energy)
    total_energies = batch.sum_over(corrected_edge_energy) / 2
    return WithPatch(total_energies, IdPatch[Any]())