Force Fields¶

Prepare a molecule with 3D coordinates¶

Force-field evaluation and optimization require an existing 3D conformer. They do not generate coordinates or add hydrogens automatically. Supply coordinates in atom order, with explicit hydrogens when appropriate.

The following small, self-contained water example supplies coordinates directly; it is not a conformer-generation algorithm:

import cosmolkit as ck
import numpy as np

base = ck.mol_from_smiles("O").with_hydrogens()
builder = base.to_builder()
builder.add_3d_conformer([
    [0.000, 0.000, 0.000],  # O
    [0.960, 0.000, 0.000],  # H
    [-0.240, 0.930, 0.000], # H
])
mol = builder.build().with_assigned_valence()
conformer_id = mol.conformers_3d[0].id

An application can instead supply a molecule loaded from a 3D molecular file or produced by a conformer-generation API. Force-field calls below do not depend on how those coordinates were obtained.

MMFF parameters, energy, and gradient¶

if not mol.mmff_has_all_molecule_params():
    raise ValueError("MMFF parameters are unavailable for this molecule")

properties = mol.mmff_properties_with_params(
    ck.MmffPropertiesParams(
        mmff_variant=ck.MmffVariant.Mmff94,
        # mmff_variant="MMFF94",
    )
)
print(properties.atom_type(0))
print(properties.formal_charge(0))
print(properties.partial_charge(0))

evaluation = mol.mmff_energy_gradient_with_params(
    ck.MmffEvaluationParams(
        mmff_variant=ck.MmffVariant.Mmff94,
        # mmff_variant="MMFF94",
        conformer_id=conformer_id,
        non_bonded_threshold=100.0,
        ignore_interfragment_interactions=True,
    )
)
if evaluation is None:
    raise ValueError("MMFF typing did not produce a valid force field")

print(evaluation.energy)
gradient = np.asarray(evaluation.gradient, dtype=np.float64).reshape(-1, 3)
print(gradient.shape)

The short forms mol.mmff_properties() and mol.mmff_energy_gradient() use default parameters. The current gradient is a flat Python list in [dx0, dy0, dz0, dx1, dy1, dz1, ...] order; the reshape above is performed by the caller. Force is the negative energy gradient.

MMFF single- and multiple-conformer optimization¶

result = mol.with_mmff_optimized_with_params(
    ck.MmffOptimizationParams(
        mmff_variant=ck.MmffVariant.Mmff94,
        # mmff_variant="MMFF94",
        max_iterations=200,
        conformer_id=conformer_id,
        non_bonded_threshold=100.0,
    )
)
optimized = result.molecule
print(result.status_code)
print(result.needs_more)

# The single-conformer MMFF result has no energy() accessor.
final_evaluation = optimized.mmff_energy_gradient_with_params(
    ck.MmffEvaluationParams(conformer_id=conformer_id)
)
if final_evaluation is not None:
    print(final_evaluation.energy)

# Optimize every stored 3D conformer independently.
multiple = mol.with_mmff_optimized_conformers_with_params(
    ck.MmffConformerOptimizationParams(
        num_threads=1,
        max_iterations=200,
        mmff_variant=ck.MmffVariant.Mmff94,
        # mmff_variant="MMFF94",
        non_bonded_threshold=100.0,
    )
)
for row in multiple.conformer_results:
    print(row.status_code, row.needs_more, row.energy)
optimized_conformers = multiple.molecule

These value-style calls leave mol unchanged. Status 0 indicates convergence, 1 indicates that more iterations are needed, and MMFF status -1 indicates invalid parameterization. Therefore not result.needs_more alone is not a success check: inspect result.status_code as well.

Default single-conformer MMFF optimization uses 200 iterations and a non-bonded threshold of 100.0. Default multiple-conformer optimization uses 1000 iterations and a threshold of 10.0. The example explicitly chooses the same threshold for both workflows.

UFF evaluation and optimization¶

if not mol.uff_has_all_molecule_params():
    raise ValueError("UFF parameters are unavailable for this molecule")

evaluation = mol.uff_energy_gradient_with_params(
    ck.UffEvaluationParams(conformer_id=conformer_id, vdw_threshold=10.0)
)
print(evaluation.energy)
gradient = np.asarray(evaluation.gradient, dtype=np.float64).reshape(-1, 3)

result = mol.with_uff_optimized_with_params(
    ck.UffOptimizationParams(max_iterations=200, conformer_id=conformer_id)
)
optimized = result.molecule
print(result.status_code, result.needs_more, result.energy)

multiple = mol.with_uff_optimized_conformers_with_params(
    ck.UffConformerOptimizationParams(num_threads=1, max_iterations=200)
)
for row in multiple.conformer_results:
    print(row.conformer_id, row.status_code, row.energy)

The corresponding default calls are uff_energy_gradient(), with_uff_optimized(), and with_uff_optimized_conformers(). Default UFF optimization uses 1000 iterations; its van der Waals threshold is 10.0. Malformed inputs and evaluation failures propagate as Python exceptions; non-convergence is an optimization result, not an exception.

Interactive Python workflow¶

ff = mol.mmff_force_field(
    conformer_id=conformer_id,
    mmff_variant=ck.MmffVariant.Mmff94,
    # mmff_variant="MMFF94",
    non_bonded_threshold=100.0,
    ignore_interfragment_interactions=True,
)
# UFF exposes the same handle interface:
# ff = mol.uff_force_field(conformer_id=conformer_id, vdw_threshold=10.0)

x, y, z = ff.position(1)
ff.set_position_(1, (x + 0.2, y, z))
print(ff.energy())
gradient = np.asarray(ff.gradient(), dtype=np.float64)  # shape (N, 3)

ff.set_fixed_atoms_([0, 2])
full_gradient = ff.gradient_unconstrained()  # includes fixed-atom rows
outcome = ff.minimize_(max_iterations=20)
print(outcome.converged, outcome.iterations, outcome.energy)

# Another drag, followed by relaxation from the current coordinates.
ff.set_position_(1, (x + 0.4, y, z))
outcome = ff.minimize_(max_iterations=20)

positions = np.asarray(ff.positions, dtype=np.float64)  # independent snapshot
positions[:, 2] += 0.1
ff.set_positions_(positions)
ff.set_fixed_atoms_([])      # replace the fixed set with an empty set

evaluation = ff.energy_gradient()
print(evaluation.energy, np.asarray(evaluation.gradient).shape)

Factory and configuration surface¶

Python factories accept keyword-only configuration:

mmff = mol.mmff_force_field(
    conformer_id=conformer_id,
    mmff_variant=ck.MmffVariant.Mmff94,
    # mmff_variant="MMFF94",
    non_bonded_threshold=100.0, ignore_interfragment_interactions=True,
)

uff = mol.uff_force_field(
    conformer_id=conformer_id, vdw_threshold=10.0,
    ignore_interfragment_interactions=True,
)
print(mmff.energy(), uff.energy())

Reusable configuration remains available through parameter objects:

params = ck.MmffForceFieldParams(
    conformer_id=conformer_id,
    mmff_variant=ck.MmffVariant.Mmff94,
    # mmff_variant="MMFF94",
)
ff = mol.mmff_force_field(params)
outcome = ff.minimize_(
    ck.ForceFieldMinimizeParams(
        max_iterations=20, force_tolerance=1e-4, energy_tolerance=1e-6
    )
)

Python parameter-object and keyword calls use the same method, implementation, fields and defaults; do not combine the two call forms.

Distance, angle, torsion, and position constraints¶

All three parameterizations (UFF, MMFF94, MMFF94s) support the same four soft constraint contributions on a persistent handle. These operations require conformer support. Constraints add energy penalties; they do not clamp coordinates or replace the fixed-atom mask.

import cosmolkit as ck

mol = ck.mol_from_smiles("CCCC").with_hydrogens().with_3d_conformer()
field = mol.uff_force_field()
field.add_distance_constraint_(
    0, 1, relative=False,
    min_distance=1.3, max_distance=1.6, force_constant=100.0,
)
anchor = ck.ForceFieldPositionConstraintParams(
    max_displacement=0.2, force_constant=100.0,
)
anchor.force_constant = 200.0
field.add_position_constraint_(0, anchor)
outcome = field.minimize_(max_iterations=20)

The methods and their parameter types are:

• add_distance_constraint_(atom1, atom2, params): ForceFieldDistanceConstraintParams(relative, min_distance, max_distance, force_constant).

• add_angle_constraint_(atom1, atom2, atom3, params): ForceFieldAngleConstraintParams(relative, min_angle_deg, max_angle_deg, force_constant); atom2 is the angle’s vertex.

• add_torsion_constraint_(atom1, atom2, atom3, atom4, params): ForceFieldTorsionConstraintParams(relative, min_dihedral_deg, max_dihedral_deg, force_constant).

• add_position_constraint_(atom, params): ForceFieldPositionConstraintParams(max_displacement, force_constant).

Every listed field is required; no range, relative mode, or force constant is guessed. Constructors use keyword-only fields. A method accepts either one parameter instance or its keyword fields, not both.

Distances and displacements use angstroms; angle and torsion bounds use degrees. Force constants use the corresponding RDKit contribution’s native convention, without unit conversion. relative=True adds each bound to the distance, angle, or signed dihedral measured from the handle’s current coordinates when the constraint is added. A position constraint captures that atom’s current position as its anchor at addition time. Later setters or minimization do not recapture either baseline. Torsion bounds retain RDKit’s normalized wrapped interval semantics rather than being clipped to an ordinary ordered interval.

Successful calls append contributions in call order. Invalid indices, bounds, or angle ranges return typed ForceFieldError values without adding a term; the original molecule and independently created handles remain unchanged. InvalidAtomIndex reports the input slot in component and its value in atom_index. InvalidConstraintBounds and InvalidConstraintRange retain the concrete source cause instead of silently repairing the input.

The complete Python example covers all four constraints, fixed atoms, and explicit installation of optimized coordinates in a new molecule.

Handle and result API¶

See cosmolkit.MolecularForceField, cosmolkit.ForceFieldEnergyGradient and cosmolkit.ForceFieldMinimizeOutcome for generated signatures.

set_position_ accepts a three-component sequence. set_positions_ accepts nested numeric sequences or NumPy arrays of shape (N, 3); input values are copied into handle-owned storage. Python atom IDs are indices in the fixed atom order, not IDs transferable between unrelated molecules. Rust uses AtomId and borrowed coordinate slices; Python never receives a mutable view into the handle’s live storage.

Both keyword minimization and ForceFieldMinimizeParams expose energy_tolerance=1e-6, alongside max_iterations=200 and force_tolerance=1e-4. Rust constructs these parameters with ForceFieldMinimizeParams::new(max_iterations, force_tolerance, energy_tolerance). All three arguments reach the same source-backed optimizer. The pinned RDKit BFGSOpt.h explicitly ignores funcTol via RDUNUSED_PARAM(funcTol): changing energy tolerance alone therefore does not currently change convergence. CK preserves that behavior rather than inventing an energy stopping criterion.

Errors expose their typed discriminant through variant and preserve actual lower-level causes in __cause__. Payload fields are documented on the corresponding exception in the API reference. InvalidTolerance identifies the source optimizer’s rejected force tolerance (component == 0); energy tolerance is forwarded without a new CK-only constraint.

Ownership, updates, and numerical semantics¶

• The handle owns its coordinates and evaluation state. It can outlive the input molecule and never implicitly writes coordinates back to it.

• Factories select an existing 3D conformer. With conformer_id=None, the target follows the existing force-field default of selecting the first stored 3D conformer. No embedding or hydrogen addition is implicit.

• Topology, atom order, and parameterization are fixed at construction. Changing topology requires a new handle.

• Coordinate setters validate shape, atom indices, and finite values before committing. Invalid input leaves handle coordinates unchanged.

• Coordinate updates invalidate relevant distance caches automatically. Users do not call Initialize. A single-atom update avoids copying all positions or retyping atoms; this does not make energy evaluation O(1).

• positions and gradient() return independent nested lists. Use np.asarray(values, dtype=np.float64) for array operations. Editing a snapshot alone does not change the force field.

• gradient_unconstrained() returns full energy derivatives as independent rows of three floats, including fixed-atom rows. It omits only the fixed-atom zeroing performed by gradient(); it does not remove energy terms, change the fixed set or alter minimization. Physical force is the negative of this gradient.

• set_fixed_atoms_ replaces the complete fixed set. Fixed atoms remain in energy terms but have zero components in the constrained gradient and do not move during minimization. Explicit setters may move them to new anchors.

• Non-bonded contributions are selected at construction using the configured threshold and initial geometry. Updating coordinates refreshes evaluation caches, not the contribution list. Rebuild the handle if interactions excluded at construction must be reconsidered after a large displacement.

• Repeated minimization resumes from current coordinates, but each call starts a new optimizer history. Two 20-iteration calls are not promised to equal one 40-iteration call. Non-convergence is reported by converged=False.

Construction errors are MmffForceFieldError and UffForceFieldError; handle update, evaluation, and minimization errors use ForceFieldError with structured kinds and preserved causes. Invalid MMFF parameterization must fail construction rather than return a dummy handle.

The existing one-shot APIs remain useful for ordinary molecule evaluation and optimization. The persistent API adds an interactive evaluator; it does not replace the chemistry algorithms or introduce a second force-field kernel.