Computing Vibrational Frequencies Using Analytical Hessian

TorchANI is able to use ASE interface to do structure optimization and vibration analysis, but the Hessian in ASE’s vibration analysis is computed numerically, which is slow and less accurate.

TorchANI therefore provide an interface to compute the Hessian matrix and do vibration analysis analytically, thanks to the super power of torch.autograd.

import ase
import ase.optimize
import torch
import torchani
import math

Let’s now manually specify the device we want TorchANI to run:

device = torch.device('cpu')
model = torchani.models.ANI1x(periodic_table_index=True).to(device).double()

/opt/hostedtoolcache/Python/3.8.17/x64/lib/python3.8/site-packages/torchani-2.2.3-py3.8.egg/torchani/resources/

Let’s first construct a water molecule and do structure optimization:

d = 0.9575
t = math.pi / 180 * 104.51
molecule = ase.Atoms('H2O', positions=[
    (d, 0, 0),
    (d * math.cos(t), d * math.sin(t), 0),
    (0, 0, 0),
], calculator=model.ase())
opt = ase.optimize.BFGS(molecule)
opt.run(fmax=1e-6)

      Step     Time          Energy         fmax
BFGS:    0 00:04:45    -2078.633392        0.6576
BFGS:    1 00:04:45    -2078.637244        0.1905
BFGS:    2 00:04:45    -2078.637801        0.0292
BFGS:    3 00:04:45    -2078.637827        0.0183
BFGS:    4 00:04:45    -2078.637849        0.0014
BFGS:    5 00:04:45    -2078.637849        0.0001
BFGS:    6 00:04:45    -2078.637849        0.0000
BFGS:    7 00:04:46    -2078.637849        0.0000

True

Now let’s extract coordinates and species from ASE to use it directly with TorchANI:

species = torch.tensor(molecule.get_atomic_numbers(), device=device, dtype=torch.long).unsqueeze(0)
coordinates = torch.from_numpy(molecule.get_positions()).unsqueeze(0).requires_grad_(True)

TorchANI needs the masses of elements in AMU to compute vibrations. The masses in AMU can be obtained from a tensor with atomic numbers by using this utility:

masses = torchani.utils.get_atomic_masses(species)

To do vibration analysis, we first need to generate a graph that computes energies from species and coordinates. The code to generate a graph of energy is the same as the code to compute energy:

energies = model((species, coordinates)).energies

We can now use the energy graph to compute analytical Hessian matrix:

hessian = torchani.utils.hessian(coordinates, energies=energies)

The Hessian matrix should have shape (1, 9, 9), where 1 means there is only one molecule to compute, 9 means 3 atoms * 3D space = 9 degree of freedom.

print(hessian.shape)

torch.Size([1, 9, 9])

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. 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[6:])
print('Force Constants (mDyne/A):', fconstants[6:])
print('Reduced masses (AMU):', rmasses[6:])
print('Modes:', modes[6:])

Frequencies (cm^-1): tensor([1913.200, 3882.404, 3911.185], dtype=torch.float64)
Force Constants (mDyne/A): tensor([2.339, 9.269, 9.743], dtype=torch.float64)
Reduced masses (AMU): tensor([1.085, 1.044, 1.081], dtype=torch.float64)
Modes: tensor([[[     0.015,      0.677,      0.000],
         [     0.652,      0.184,     -0.000],
         [    -0.042,     -0.054,      0.000]],

        [[     0.689,     -0.054,      0.000],
         [    -0.225,      0.654,      0.000],
         [    -0.029,     -0.038,      0.000]],

        [[    -0.678,     -0.006,      0.000],
         [    -0.164,      0.658,     -0.000],
         [     0.053,     -0.041,     -0.000]]], dtype=torch.float64)