.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/vibration_analysis.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_vibration_analysis.py: 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`. .. GENERATED FROM PYTHON SOURCE LINES 13-20 .. code-block:: default import ase import ase.optimize import torch import torchani import math .. GENERATED FROM PYTHON SOURCE LINES 21-22 Let's now manually specify the device we want TorchANI to run: .. GENERATED FROM PYTHON SOURCE LINES 22-25 .. code-block:: default device = torch.device('cpu') model = torchani.models.ANI1x(periodic_table_index=True).to(device).double() .. rst-class:: sphx-glr-script-out .. code-block:: none /opt/hostedtoolcache/Python/3.8.16/x64/lib/python3.8/site-packages/torchani-2.2.3-py3.8.egg/torchani/resources/ .. GENERATED FROM PYTHON SOURCE LINES 26-27 Let's first construct a water molecule and do structure optimization: .. GENERATED FROM PYTHON SOURCE LINES 27-37 .. code-block:: default 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) .. rst-class:: sphx-glr-script-out .. code-block:: none Step Time Energy fmax BFGS: 0 18:14:27 -2078.633392 0.6576 BFGS: 1 18:14:27 -2078.637244 0.1905 BFGS: 2 18:14:28 -2078.637801 0.0292 BFGS: 3 18:14:28 -2078.637827 0.0183 BFGS: 4 18:14:28 -2078.637849 0.0014 BFGS: 5 18:14:28 -2078.637849 0.0001 BFGS: 6 18:14:28 -2078.637849 0.0000 BFGS: 7 18:14:28 -2078.637849 0.0000 True .. GENERATED FROM PYTHON SOURCE LINES 38-40 Now let's extract coordinates and species from ASE to use it directly with TorchANI: .. GENERATED FROM PYTHON SOURCE LINES 40-43 .. code-block:: default 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) .. GENERATED FROM PYTHON SOURCE LINES 44-47 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: .. GENERATED FROM PYTHON SOURCE LINES 47-49 .. code-block:: default masses = torchani.utils.get_atomic_masses(species) .. GENERATED FROM PYTHON SOURCE LINES 50-53 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: .. GENERATED FROM PYTHON SOURCE LINES 53-55 .. code-block:: default energies = model((species, coordinates)).energies .. GENERATED FROM PYTHON SOURCE LINES 56-57 We can now use the energy graph to compute analytical Hessian matrix: .. GENERATED FROM PYTHON SOURCE LINES 57-59 .. code-block:: default hessian = torchani.utils.hessian(coordinates, energies=energies) .. GENERATED FROM PYTHON SOURCE LINES 60-62 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`. .. GENERATED FROM PYTHON SOURCE LINES 62-64 .. code-block:: default print(hessian.shape) .. rst-class:: sphx-glr-script-out .. code-block:: none torch.Size([1, 9, 9]) .. GENERATED FROM PYTHON SOURCE LINES 65-69 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. .. GENERATED FROM PYTHON SOURCE LINES 69-76 .. code-block:: default 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:]) .. rst-class:: sphx-glr-script-out .. code-block:: none 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) .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 1.777 seconds) .. _sphx_glr_download_examples_vibration_analysis.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: vibration_analysis.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: vibration_analysis.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_