Function calculate_constraint_force 

fn calculate_constraint_force(
    geometry: &Geometry,
    constraint: &Constraint,
    violation: f64,
) -> Vec<f64>

Calculates the constraint force (gradient) for a single geometric constraint.

This function computes the force that should be applied to atomic coordinates to reduce the constraint violation. The force points in the direction that will move the coordinate value toward its target.

Theoretical Background

The constraint force is the negative gradient of the constraint violation:

F_constraint = -∇[g(x) - target] = -∇g(x)

For constraint satisfaction, we want to minimize |g(x) - target|², so:

F = -2 × (g(x) - target) × ∇g(x)

Algorithm Implementation

Bond Constraints (Analytical)

For bond constraints, the gradient is computed analytically:

∇r₁₂ = (r₂ - r₁) / |r₂ - r₁|  (unit vector along bond)

Forces are applied equally and oppositely to the two atoms.

Angle/Dihedral Constraints (Finite Difference)

For complex constraints, finite difference approximation is used:

∂g/∂xᵢ ≈ [g(x + δeᵢ) - g(x)] / δ

where eᵢ is the unit vector in coordinate direction i.

Arguments

• geometry - Current molecular geometry
• constraint - The geometric constraint to compute forces for
• violation - Current constraint violation (from calculate_constraint_violation)

Returns

A vector of forces with length 3×N_atoms, where each group of 3 consecutive elements represents the [Fx, Fy, Fz] force components for one atom.

Force Scaling

The force magnitude is proportional to the violation:

|F| = k × |violation|

where k = 10.0 is an empirical force constant.

Implementation Details

Bond Force Calculation

1. Calculate bond vector and distance
2. Compute unit vector along bond
3. Apply force proportional to violation
4. Equal and opposite forces on the two atoms

Finite Difference Parameters

• Step size (δ): 0.001 Angstrom (compromise between accuracy and numerical stability)
• One-sided difference: Uses forward difference for simplicity
• Coordinate perturbation: Each Cartesian coordinate perturbed independently

Limitations

• Finite difference errors: Approximate gradients for angles/dihedrals
• Fixed step size: Not adaptive to constraint type or system size
• No second derivatives: Cannot account for constraint curvature
• Proportional control: Simple linear relationship between violation and force