Now vibrational analysis includes reduced masses, force constants and normalized modes (#413)

* Add force constants, reduced masses and normalized normal modes to vibrational analysis

* Formatting to make flake8 happy

examples/vibration_analysis.py

@@ -69,9 +69,11 @@ print(hessian.shape)
 # We are now ready to compute vibrational frequencies. The output has unit
 # cm^-1. Since there are in total 9 degree of freedom, there are in total 9
 # frequencies. Only the frequencies of the 3 vibrational modes are interesting.
-freq, modes = torchani.utils.vibrational_analysis(masses, hessian)
-freq = freq[-3:]
-modes = modes[-3:]
+# We output the modes as MDU (mass deweighted unnormalized), to compare with ASE.
+freq, modes, fconstants, rmasses = torchani.utils.vibrational_analysis(masses, hessian, mode_type='MDU')
+torch.set_printoptions(precision=3, sci_mode=False)
 
-print('Frequencies (cm^-1):', freq)
-print('Modes:', modes)
+print('Frequencies (cm^-1):', freq[6:])
+print('Force Constants (mDyne/A):', fconstants[6:])
+print('Reduced masses (AMU):', rmasses[6:])
+print('Modes:', modes[6:])

tests/test_vibrational.py

@@ -40,7 +40,7 @@ class TestVibrational(unittest.TestCase):
         coordinates = torch.from_numpy(molecule.get_positions()).unsqueeze(0).requires_grad_(True)
         _, energies = model((species, coordinates))
         hessian = torchani.utils.hessian(coordinates, energies=energies)
-        freq2, modes2 = torchani.utils.vibrational_analysis(masses[species], hessian)
+        freq2, modes2, _, _ = torchani.utils.vibrational_analysis(masses[species], hessian)
         freq2 = freq2[6:].float()
         modes2 = modes2[6:]
         ratio = freq2 / freq

torchani/utils.py

@@ -284,8 +284,31 @@ class FreqsModes(NamedTuple):
     modes: Tensor
 
 
-def vibrational_analysis(masses, hessian, unit='cm^-1'):
-    """Computing the vibrational wavenumbers from hessian."""
+class VibAnalysis(NamedTuple):
+    freqs: Tensor
+    modes: Tensor
+    fconstants: Tensor
+    rmasses: Tensor
+
+
+def vibrational_analysis(masses, hessian, mode_type='MDU', unit='cm^-1'):
+    """Computing the vibrational wavenumbers from hessian.
+
+    Note that normal modes in many popular software packages such as
+    Gaussian and ORCA are output as mass deweighted normalized (MDN).
+    Normal modes in ASE are output as mass deweighted unnormalized (MDU).
+    Some packages such as Psi4 let ychoose different normalizations.
+    Force constants and reduced masses are calculated as in Gaussian.
+
+    mode_type should be one of:
+    - MWN (mass weighted normalized)
+    - MDU (mass deweighted unnormalized)
+    - MDN (mass deweighted normalized)
+
+    MDU modes are orthogonal but not normalized, MDN modes are normalized
+    but not orthogonal. MWN modes are orthonormal, but they correspond
+    to mass weighted cartesian coordinates (x' = sqrt(m)x).
+    """
     if unit != 'cm^-1':
         raise ValueError('Only cm^-1 are supported right now')
     assert hessian.shape[0] == 1, 'Currently only supporting computing one molecule a time'
@@ -306,8 +329,25 @@ def vibrational_analysis(masses, hessian, unit='cm^-1'):
     frequencies = angular_frequencies / (2 * math.pi)
     # converting from sqrt(hartree / (amu * angstrom^2)) to cm^-1
     wavenumbers = frequencies * 17092
-    modes = (eigenvectors.t() * inv_sqrt_mass).reshape(frequencies.numel(), -1, 3)
-    return FreqsModes(wavenumbers, modes)
+
+    # Note that the normal modes are the COLUMNS of the eigenvectors matrix
+    mw_normalized = eigenvectors.t()
+    md_unnormalized = mw_normalized * inv_sqrt_mass
+    norm_factors = 1 / torch.norm(md_unnormalized, dim=1)  # units are sqrt(AMU)
+    md_normalized = md_unnormalized * norm_factors.unsqueeze(1)
+
+    rmasses = norm_factors**2  # units are AMU
+    # The conversion factor for Ha/(AMU*A^2) to mDyne/(A*AMU) is 4.3597482
+    fconstants = eigenvalues * rmasses * 4.3597482  # units are mDyne/A
+
+    if mode_type == 'MDN':
+        modes = (md_normalized).reshape(frequencies.numel(), -1, 3)
+    elif mode_type == 'MDU':
+        modes = (md_unnormalized).reshape(frequencies.numel(), -1, 3)
+    elif mode_type == 'MWN':
+        modes = (mw_normalized).reshape(frequencies.numel(), -1, 3)
+
+    return VibAnalysis(wavenumbers, modes, fconstants, rmasses)
 
 
 __all__ = ['pad', 'pad_atomic_properties', 'present_species', 'hessian',