{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Compute the phonons and the Raman spectrum\n", "\n", "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.\n", "\n", "The example provided here shows how to obtain the phonons and the Raman spectrum of the N$_2$ molecule using the [Phonons](https://mmoriniere.gitlab.io/MyBigDFT/phonons.html#mybigdft.workflows.phonons.Phonons) and [RamanSpectrum](https://mmoriniere.gitlab.io/MyBigDFT/raman.html#mybigdft.workflows.ramanspectrum.RamanSpectrum) class.\n", "\n", "The goal of these calculations is to find the normal modes of vibration (or phonons) of the system under consideration as well as their respective Raman intensities. A good reference for the underlying theory of molecular vibrations is the book *Molecular Vibrations* by Wilson *et al.* (1955) or the *Advances in Chemical Physics, vol. 67* by K.P. Lawley or *Spectra of Atoms and Molecules* by Peter F. Bernath (where the conversion between atomic units, SI units and the units used in the litterature is well explained in pages 316-317).\n", "\n", "For more details on how all these quantities are computed, you may read [A. Stirling, *J. Chem. Phys.* **104**, 1254 (1996)](http://dx.doi.org/10.1063/1.470783) and [D. Porezag and M. R.Pederson, *PRB* **54**, 7830-7836 (1996)](https://link.aps.org/doi/10.1103/PhysRevB.54.7830). To check for the polarizability tensor and, more importantly, its derivative, we used [K. M. Gough *et al.*, *Can. J. Chem.* **74**, 1139-1144 (1996)](https://doi.org/10.1139/v96-128) as a comparison." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initialization\n", "\n", "You first need to import some useful classes to define a ground state calculation as well as the `Phonons` and `RamanSpectrum` classes:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from mybigdft import Job, Posinp, Atom, InputParams\n", "from mybigdft.workflows import Phonons, RamanSpectrum" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The [Phonons](https://mmoriniere.gitlab.io/MyBigDFT/phonons.html#mybigdft.workflows.phonons.Phonons) class" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us first obtain the the normal modes of the N2 molecule. \n", "\n", "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).\n", "\n", "**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).\n", "\n", "**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}$)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "N2_ref = \"\"\"2 angstroem\n", "free\n", "N 3.571946174 3.571946174 3.620526682\n", "N 3.571946174 3.571946174 4.71401439\"\"\"\n", "pos = Posinp.from_string(N2_ref)\n", "gs = Job(posinp=pos, name='N2', run_dir='../../../tests/phonons_N2')\n", "phonons = Phonons(gs)\n", "# The line above actually amounts to:\n", "# phonons = Phonons(gs, translation_amplitudes=[0.45/64]*3, order=2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To run the calculation, simply use the run method:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/x-\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/y-\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/z-\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/x-\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/y-\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/z-\n", "Logfile log-N2.yaml already exists!\n", "\n" ] } ], "source": [ "phonons.run(nmpi=2, nomp=2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the end, you can get the phonon energies via the `energies` attribute. They are stored in a dictionary whose keys are the units." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 2.38694633e+03, 6.64733640e+00, 6.64737304e+00, -1.86019702e-02,\n", " 1.18120290e-05, -4.75490601e-06])" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "phonons.energies # values in cm^-1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The [RamanSpectrum](https://mmoriniere.gitlab.io/MyBigDFT/raman.html#mybigdft.workflows.ramanspectrum.RamanSpectrum) class\n", "\n", "To get the **intensities** (or **activities**) of the Raman spectrum, one must compute the **derivative of the polarizability tensor along the normal modes**. To that end, one must compute the polarizability tensor at each of the positons used to get the vibrational energies, and this means applying an external electric field along each coordinate. This leads to 3 extra calculations per atom, meaning that $18 n_{at}$ additional BigDFT standard calculations are required to obtain the intensities of the Raman spectrum intensities, leading to 24 $n_{at}$ calculations in total.\n", "\n", "A RamanSpectrum instance is initialized by a Phonons instance (you may also provide electric field amplitudes (via the `ef_amplitudes` argument) to set the amplitude of the electric field applied when computing the polarizability tensor). \n", "\n", "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).\n", "\n", "**Note**: You can set ``order=2`` upon initializing the RamanSpectrum class to run two calculations per electric field coordinates per atom, leading to $36 n_{at}$ extra calculations, instead of $18 n_{at}$." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/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.\n", " warnings.warn(warning_msg, UserWarning)\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/x+/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/x+/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/x+/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/x-\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/x-/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/x-/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/x-/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/y+/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/y+/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/y+/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/y-\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/y-/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/y-/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/y-/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/z+/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/z+/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/z+/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/z-\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/z-/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/z-/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0000/z-/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/x+/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/x+/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/x+/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/x-\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/x-/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/x-/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/x-/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/y+/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/y+/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/y+/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/y-\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/y-/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/y-/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/y-/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/z+/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/z+/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/z+/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/z-\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/z-/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/z-/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/tests/phonons_N2/atom0001/z-/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n" ] } ], "source": [ "raman = RamanSpectrum(phonons)\n", "raman.run(nmpi=2, nomp=2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the end, we can check the phonon energies:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 2.38694633e+03, 6.64733640e+00, 6.64737304e+00, -1.86019702e-02,\n", " 1.18120290e-05, -4.75490601e-06])" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "raman.energies # in cm^-1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The intensities associated to each normal mode are obtained via the `intensities` attribute:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([2.25644573e+01, 1.07011259e+00, 1.06300226e+00, 1.46246429e-09,\n", " 2.58350510e-21, 1.44359058e-22])" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "raman.intensities # in Ang^4.amu^-1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another quantity is the depolarization ratio:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0.09412174, 0.75 , 0.75 , 0.12054558, 0.1427555 ,\n", " 0.29955284])" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "raman.depolarization_ratios" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The nitrogen molecule being linear, it has $3 n_{at} - 5 = 1$ normal mode. Its energy, found using the LDA exchange-correlation is around 2387 cm$^{-1}$. Following the statement of [B. G. Johnson *et al.*, *J. Chem. Phys.* **98**, 5612 (1993)](http://aip.scitation.org/doi/10.1063/1.464906), our result should be compared with the harmonic experimental value, which is reported to be 2360 cm$^{-1}$ in this same reference. There is a rather good agreement! Note that the actual experimental value (without neglecting anharmocity) is 2331 cm$^{-1}$.\n", "\n", "A value of 0.11 for the depolarization ratio was also reported in [P. L. Polavarapu, *J. Chem. Phys.* **94**, 8106-8112 (1990)](http://pubs.acs.org/doi/abs/10.1021/j100384a024), which compares rather well to our 0.094. An energy of 2725 cm$^{-1}$ was also reported, along with an intensity of 26.0 $\\unicode[serif]{xC5}^4$ amu$^{-1}$ (close to our 22.6 $\\unicode[serif]{xC5}^4$ amu$^{-1}$).\n", "\n", "Keep in mind that default input parameters were used, so that the present calculations cannot be considered as converged (the grid extension should be increased and the grid spacing decreased)." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.2" } }, "nbformat": 4, "nbformat_minor": 2 }