Bug: Physically unrealistic angle force constants in get_angle_force_constant

[Eashan-H]

Description

The current implementation of the Seminario method for calculating angle force constants in ffprime/bond.py produces physically unrealistic values. For a typical C–C–H angle in the provided lig.log example, the code currently yields ~2500 kcal/mol/rad², whereas standard molecular mechanics values for such interactions are expected to be in the range of 10–100 kcal/mol/rad².

Root Cause Analysis

Upon investigation, the implementation of get_angle_force_constant deviates from the standard Seminario (1996) formula in several critical ways:

1. Squared Eigenanalysis
   The code performs an eigenanalysis on H_ij @ H_ij.T. Since the Hessian H has units of force/length, this product effectively squares the force constants, leading to the observed magnitude error.

2. Incorrect Projections
   The logic uses squared projections (abs(dot)**2) of eigenvectors, while the Seminario paper specifies using the absolute value of the projection.

3. Unit Scaling and Combination
   The combination of the two bond arms does not follow the harmonic mean requirement:

   1 / k_theta = r_ji^2 / k_PA + r_jk^2 / k_PC
   

   Instead, it sums contributions and divides by the sum of r^2, which is dimensionally inconsistent for an angular force constant.

Proposed Fix

I propose refactoring the get_angle_force_constant method to strictly follow the Seminario (1996) derivation:

• Perform direct eigenanalysis on the off-diagonal Hessian blocks.
• Use absolute projections for directional stiffness.
• Implement the correct harmonic combination rule for the two arms of the angle.
• Add a dedicated test suite (tests/test_bond.py) to verify these values against reference calculations.

Verification

I verified a local fix using the files:

• examples/lig.log
• examples/lig.fchk

The corrected implementation yields 14.70 kcal/mol/rad² for the first angle (A1), matching the expected reference value.