Skip to content

Refactor angle force constant formula and add test suite#7

Open
Eashan-H wants to merge 2 commits intotheochem:mainfrom
Eashan-H:fix-angle-force-constant
Open

Refactor angle force constant formula and add test suite#7
Eashan-H wants to merge 2 commits intotheochem:mainfrom
Eashan-H:fix-angle-force-constant

Conversation

@Eashan-H
Copy link
Copy Markdown

Summary
This PR refactors the get_angle_force_constant method in ffprime/bond.py to strictly adhere to the standard Seminario (1996) method. The previous
implementation produced physically unrealistic values (~2500 kcal/mol/rad²) due to squared eigenvalue analysis and incorrect unit scaling. The corrected
code now yields standard molecular mechanics values (~10–25 kcal/mol/rad²) that match established literature results.

Physical & Mathematical Rationale
The derivation follows the Seminario method for calculating harmonic force constants directly from the QM Hessian. The core logic was updated to address
three critical issues:

  1. Eigenanalysis: Replaced the use of $H H^T$ (which effectively squared the force constants) with a direct eigenanalysis of the off-diagonal $3 \times
    3$ Hessian blocks ($H_{ij}$ and $H_{kj}$) connecting the central atom $j$ to terminal atoms $i$ and $k$.
  2. Directional Projections: Updated the projection of eigenvalues onto the bending directions ($u_{PA}$ and $u_{PC}$) using absolute weights rather than
    squared weights, as specified in the original paper.
  3. Harmonic Combination: Implemented the correct combination rule for the two "arms" of the angle. The angle force constant $k_{\theta}$ is derived via
    the harmonic mean of the projected Cartesian stiffness:
    $$\frac{1}{k_{\theta}} = \frac{r_{ji}^2}{k_{PA}} + \frac{r_{jk}^2}{k_{PC}}$$
    where $r$ is the equilibrium bond length and $k_{P}$ is the projected force constant.

Changes

  • ffprime/bond.py:
    • Refactored get_angle_force_constant with the corrected Seminario logic.
    • Updated get_force_constant and get_angle_force_constant to default the hess argument to self.hess_new, improving usability.
    • Improved docstrings with LaTeX equations and clear unit documentation ($kcal/mol/\text{\AA}^2$ to $kcal/mol/rad^2$).
    • Wrapped example execution code in if name == "main": to prevent side effects during library imports.
  • tests/test_bond.py:
    • Added a new pytest suite to verify the Bonded class.
    • Included accuracy tests that match the calculated results against the standard Seminario reference for the lig.log example.

Verification Results
The implementation was verified using the provided examples/lig.log and lig.fchk files:

  • Angle A1 (2-1-13):
    • Previous value: ~2498.5 kcal/mol/rad²
    • Corrected value: 14.70 kcal/mol/rad² (Matches reference)
  • Comprehensive Check: All 26 angles in the example system now fall within the physically expected range of 10.2 to 24.4 kcal/mol/rad².
  • Automated Tests: All 4 new test cases in tests/test_bond.py pass successfully.

Closes : Bug: Physically unrealistic angle force constants in get_angle_force_constant #6

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant