From c8658b2a3515c1f5080a600ea53d87b78c055891 Mon Sep 17 00:00:00 2001 From: "Josef M. Gallmetzer" <64498081+galjos@users.noreply.github.com> Date: Sun, 23 Aug 2026 21:30:09 +0200 Subject: [PATCH] Build the inertia tensor from centred coordinates rotational_modes centred the coordinates into atom_coords_cm but passed the original, un-centred coordinates to inertia_tensor. Away from the origin that tensor picks up the parallel-axis contribution of the centre of mass, so its eigenvectors are no longer the principal axes of the molecule. For a linear molecule the molecular axis is then not an eigenvector and none of the three rotational displacement vectors comes out with zero norm, so the linear-molecule filter keeps three rotations instead of two. The consequence is that a diatomic away from the origin has six external modes for six degrees of freedom, the internal subspace is empty, and hessian_sign_factor short-circuits to +1.0 instead of detecting the negative Hessian convention. A CO molecule with a negative-convention Hessian at a generic position was reported at -1436.17 cm-1, an imaginary mode, while the same molecule centred on the origin gave the correct +1436.17 cm-1. Building the tensor from the centred coordinates restores the principal axes. For non-linear molecules the three rotational vectors still span the same subspace, so the water reference spectrum is unchanged. --- .../vibrational/vibrational_analysis.py | 4 ++- .../vibrational/test_vibrational_analysis.py | 31 +++++++++++++++++++ 2 files changed, 34 insertions(+), 1 deletion(-) diff --git a/PQAnalysis/analysis/vibrational/vibrational_analysis.py b/PQAnalysis/analysis/vibrational/vibrational_analysis.py index 5d365815..764a2630 100644 --- a/PQAnalysis/analysis/vibrational/vibrational_analysis.py +++ b/PQAnalysis/analysis/vibrational/vibrational_analysis.py @@ -270,7 +270,9 @@ def rotational_modes( Calculate normalized rotational modes. """ atom_coords_cm = center_to_com(atom_coords, atom_masses) - _, eigenvectors = np.linalg.eigh(inertia_tensor(atom_coords, atom_masses)) + _, eigenvectors = np.linalg.eigh( + inertia_tensor(atom_coords_cm, atom_masses) + ) x_frame = eigenvectors p_frame = atom_coords_cm @ eigenvectors diff --git a/tests/analysis/vibrational/test_vibrational_analysis.py b/tests/analysis/vibrational/test_vibrational_analysis.py index bcc724a4..bf3fa713 100644 --- a/tests/analysis/vibrational/test_vibrational_analysis.py +++ b/tests/analysis/vibrational/test_vibrational_analysis.py @@ -8,9 +8,11 @@ from PQAnalysis.analysis.vibrational.vibrational_analysis import ( calculate, hessian_sign_factor, + internal_subspace, mode_displacement, mass_weighted_hessian, read_hessian_file, + rotational_modes, select_mode_indices, symmetrize_addition, wavenumber, @@ -241,6 +243,35 @@ def test_hessian_sign_factor_options(self): exception.value ) == "hessian_sign must be auto, positive, negative, 1, or -1." + def test_linear_molecule_off_origin_keeps_two_rotations(self): + axis = np.array([1.0, 1.0, 1.0]) / np.sqrt(3.0) + masses = np.array([12.011, 15.999]) + coords = np.array([np.zeros(3), 1.128 * axis]) + coords = coords - np.sum( + coords * masses[:, None], axis=0 + ) / np.sum(masses) + shifted = coords + np.array([1.0, 2.0, 3.0]) + + block = 1200.0 * np.outer(axis, axis) + hessian = np.zeros((6, 6)) + hessian[0:3, 0:3] = block + hessian[3:6, 3:6] = block + hessian[0:3, 3:6] = -block + hessian[3:6, 0:3] = -block + hessian = -hessian + + assert rotational_modes(shifted, masses).shape[1] == 2 + assert internal_subspace(shifted, masses).shape == (6, 1) + assert hessian_sign_factor(shifted, masses, hessian, "auto") == -1.0 + + reference = calculate(masses, coords, hessian) + translated = calculate(masses, shifted, hessian) + + assert reference.wavenumbers[-1] == pytest.approx(1436.17, abs=1e-2) + assert translated.wavenumbers[-1] == pytest.approx( + reference.wavenumbers[-1], abs=1e-6 + ) + def test_wavenumber_units(self): eigenvalues = np.array([1.0, -1.0])