Fix low-mode vibrational analysis to use mass-weighted coordinates#5
Merged
Conversation
The GF-matrix method (Wilson, Decius & Cross) requires building the
rigid-body translation/rotation basis and projecting the Hessian in
mass-weighted Cartesian coordinates q = M^{1/2} x. The previous
implementation used unweighted Cartesian coordinates, causing heavy
atoms (Br, I, S) to artificially inflate mode eigenvalues and misrank
the true soft conformational modes.
Changes:
- _build_vibrational_basis: translation vectors now weighted by sqrt(m_i)
per atom; rotation vectors weighted by sqrt(m_i) * (r_i x e_k); the
near-zero threshold scales with sqrt(total_mass) accordingly.
- _select_low_modes: projects the mass-weighted Hessian F_mw = M^{-1/2} H M^{-1/2}
into the vibrational subspace, diagonalises, then unweights eigenvectors
back to Cartesian displacement vectors and renormalises.
- config.py docstring: corrects eigenvalue threshold units from kcal/mol/A2
to kcal/mol/A2/Da (mass-weighted) with calibration reference values.
Benchmark on 10 molecules spanning C/Cl/S/Br/I:
- C-only molecules: cosine similarity ~0.97-0.99 (fix is near no-op, as expected)
- Br/I molecules: cosine similarity 0.62-0.72 (modes substantially reranked)
- 2-iodothiophene: soft mode count corrected from 3 to 5
- Mode overhead is <1ms vs 4-12ms for Hessian computation (negligible)
Member
|
Hey @mooreneural, thanks for the PR. We originally did this intentionally w/ the non-MW Hessian, following the literature precedent, but we ran our own benchmarks on this PR and performance seems good for "normal" non–heavy element systems, so I think this is a real improvement. A couple notes:
I'll be happy to merge after this changes. |
Replace Unicode multiplication sign (x) with ASCII x in rotation comments to resolve Ruff RUF003. Update eigenvalue_threshold docstring in generate_low_mode_seeds to reflect mass-weighted units (kcal/mol/A2/Da) introduced by the GF-matrix fix.
Contributor
Author
|
Fixed both! Replaced the Unicode multiplication signs with ASCII x to clear the Ruff warnings, and updated the generate_low_mode_seeds docstring to reflect the mass-weighted units (kcal/mol/Ų/Da). Ready to merge when you are. |
Member
|
merged! ty for the contribution |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
The GF-matrix method (Wilson, Decius & Cross) requires building the rigid-body
translation/rotation basis and projecting the Hessian in mass-weighted
Cartesian coordinates q = M^{1/2} x. The previous implementation used
unweighted Cartesian coordinates, causing heavy atoms (Br, I, S) to artificially
inflate mode eigenvalues and misrank the true soft conformational modes.
Changes
_build_vibrational_basis: translation vectors now weighted by √m_i per atom;rotation vectors weighted by √m_i · (r_i × e_k); near-zero threshold scales
with √(total_mass) accordingly.
_select_low_modes: projects the mass-weighted Hessian F_mw = M^{-1/2} H M^{-1/2}into the vibrational subspace, diagonalises, then unweights eigenvectors back
to Cartesian displacement vectors and renormalises columns.
config.py: correctslow_mode_eigenvalue_thresholddocstring units fromkcal/mol/Ų to kcal/mol/Ų/Da (mass-weighted, proportional to ω²) with
calibration reference values for torsional vs. stretch modes.
Benchmark
Cosine similarity measures alignment between the old and new softest-mode vectors.
Values near 1.0 mean the fix is a near no-op; lower values indicate the mode
ordering was substantially corrected for that molecule.
For C-only molecules the fix is a near no-op (cosine ~0.99). For Br/I molecules
the soft-mode vectors are substantially corrected, iodothiophene's softest mode
previously assigned 59% of displacement to the I atom; the corrected mode assigns
~3%. Mode-select overhead is <1 ms vs. 4–12 ms for Hessian computation.