Compute the phonons and the infrared spectrum

MyBigDFT comes with some classes implementing particular workflows of calculations. These workflows define a queue of jobs, that can easily be run sequentially, without having to worry about the Job context manager. They also generally define a particular post-processing procedure, run after all the BigDFT calculations in order to extract some meaningful imformation.

The example provided here shows how to obtain the phonons and the infrared spectrum of the CO molecule using the Phonons and InfraredSpectrum class.

Initialization

You first need to import some useful classes to define a ground state calculation as well as the Phonons and InfraredSpectrum classes:

from mybigdft import Job, Posinp, Atom, InputParams
from mybigdft.workflows import Phonons, InfraredSpectrum

The Phonons class

Let us first obtain the the normal modes of the CO molecule.

To get the energies of the Raman spectrum, one needs to find the eigenvalues of the dynamical matrix, that is closely related to the Hessian matrix. To build these matrices, one must find the derivatives of the forces when the position of each atom is translated (along each coordinate) by a small amount around their equilibrium position. To get a better precision on the derivative, each coodinate is translated positively and negatively. The number of BigDFT calculations therefore amounts to \(2 * 3 n_{at} = 6 n_{at}\), where \(n_{at}\) is the number of atoms (3 for the coordinates (\(x\), \(y\) and \(z\)), 2 for the number of calculations per coordinates).

Be careful: you must compute a geometry optimization of your system before actually computing the phonons, as you may get unphysical results. This part is not reproduced here. Similarly to the polarizability tensor calculation, you only need to provide a ground state job while default values are used to define the amplitude of the moves along each direction (you might want to set them via the translation_amplitudes argument).

Note: you can select order=1 when initializing the Phonons instance, so as to run only one calculation per coordinates (meaning \(3 n_{at}\) calculations will be performed, instead of \(6 n_{at}\))

CO_ref = """2  angstroem
free
C -8.61468162115724335E-22 -5.85961390944064004E-22 -4.99273481690253648E-03
O  6.11054044310156593E-22  5.85961390944064004E-22  1.12999273481690232E+00"""
pos = Posinp.from_string(CO_ref)
inp = InputParams({"dft": {"ixc": 11, "gnrm_cv": 1.e-5, "hgrids": [0.35]*3, "rmult": [8, 10]}})
gs = Job(posinp=pos, inputparams=inp, name='CO', run_dir='CO/phonons/')
phonons = Phonons(gs)
# The line above actually amounts to:
# phonons = Phonons(gs, translation_amplitudes=[0.45/64]*3, order=2)

To run the calculation, simply use the run method:

phonons.run(nmpi=2, nomp=2)

/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/CO/phonons/atom0000/x+
Logfile log-CO.yaml already exists!

/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/CO/phonons/atom0000/x-
Logfile log-CO.yaml already exists!

/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/CO/phonons/atom0000/y+
Logfile log-CO.yaml already exists!

/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/CO/phonons/atom0000/y-
Logfile log-CO.yaml already exists!

/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/CO/phonons/atom0000/z+
Logfile log-CO.yaml already exists!

/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/CO/phonons/atom0000/z-
Logfile log-CO.yaml already exists!

/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/CO/phonons/atom0001/x+
Logfile log-CO.yaml already exists!

/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/CO/phonons/atom0001/x-
Logfile log-CO.yaml already exists!

/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/CO/phonons/atom0001/y+
Logfile log-CO.yaml already exists!

/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/CO/phonons/atom0001/y-
Logfile log-CO.yaml already exists!

/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/CO/phonons/atom0001/z+
Logfile log-CO.yaml already exists!

/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/CO/phonons/atom0001/z-
Logfile log-CO.yaml already exists!

In the end, you can get the phonon energies via the energies attribute. They are stored in a dictionary whose keys are the units.

phonons.energies  # in cm^-1

array([ 2.11205229e+03,  5.82170475e-06, -5.11175466e+00, -5.11174626e+00,
        3.52753775e-06, -2.43602537e-03])

The InfraredSpectrum class

To get the intensities of the infrared spectrum, one must compute the derivative of the dipole moment along the normal modes. To that end, no extra DFT calculation is required, since one already has access to the dipole moments for each diplaced geometry used to compute the phonons.

An InfraredSpectrum instance is initialized by a Phonons instance.

Here, the phonons instance defined above is re-used, hence the UserWarning below: the calculations for the phonons were already performed before, so we are said to set force_run to True if we wish to run them again (but this is not our case).

ir = InfraredSpectrum(phonons)
ir.run(nmpi=2, nomp=2)

/Users/maximemoriniere/Documents/Python/MyBigDFT/mybigdft/workflows/workflow.py:118: UserWarning: Calculations already performed; set the argument 'force_run' to True to re-run them.
  warnings.warn(warning_msg, UserWarning)

In the end, we can check the phonon energies:

ir.energies  # in cm^-1

array([ 2.11205229e+03,  5.82170475e-06, -5.11175466e+00, -5.11174626e+00,
        3.52753775e-06, -2.43602537e-03])

The intensities associated to each normal mode are obtained via the intensities attribute:

ir.intensities  # in (D/A)^2/amu

array([1.55251144e+00, 1.47776864e-22, 4.07193761e-03, 4.07018828e-03,
       1.68939224e-22, 2.06873213e-12])

The carbon monoxyde molecule being linear, it has \(3 n_{at} - 5 = 1\) normal mode. Its energy, found using the GGA exchange-correlation is around 2112 cm\(^{-1}\). Following the statement of B. G. Johnson *et al.*, *J. Chem. Phys.* **98**, 5612 (1993), our result should be compared with the harmonic experimental value, which is reported to be 2170 cm\(^{-1}\) in this same reference. There is a rather good agreement! Note that the actual experimental value (without neglecting anharmocity) is 2143 cm\(^{-1}\).

An intensity of 61.2 km.mol\(^{-1}\) was reported in D. Bishop and L. Cheung, *J. Phys. Chem. Ref. Data* **11**, 119 (1982). Our result, converted in the same units, is 65.6 km.mol\(^{-1}\)), again in rather good agreement.

Note: 1 \((\text{D} / \unicode[serif]{xc5})^2 \text{amu}^{-1}\) = 42.255 km.mol\(^{-1}\).

Keep in mind that the present calculations cannot be considered as converged (the grid extension should be increased and the grid spacing decreased).