Compute the vibrational polarizability tensor¶
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 vibrational polarizability tensor of the CO molecule by using the VibPolTensor class.
The vibrational polarizability tensor represents “the second largest contribution to the static second-order polarizability tensor” [1]. It is “usually due to field-induced atomic relaxations” [1]. “Vibrational polarization is due to a distortion of the vibrational motion of the molecule by the field and it exists for all molecules having infrared-active transitions.” [2]
Once the infrared spectrum of a sytem is known, the mean vibrational polarizability \(\overline{\alpha_{vib}}\) can be computed:
\(\overline{\alpha_{vib}} = 1.338378076~10^{3} \sum_{n=1}^{N_vib} \frac{I_{n}}{{E_n}^2}\) (in \(\unicode[serif]{xc5}^3\)),
where \(I_n\) is the infrared intensity of a normal mode (in km.mol\(^{-1}\)) and \(E_n\) is the corresponding energy (in cm\(^{-1}\)).
Initialization¶
You first need to import some useful classes to define a ground state calculation as well as the VibPolTensor class:
[1]:
from mybigdft import Job, Posinp, Atom, InputParams
from mybigdft.workflows import Phonons, InfraredSpectrum, VibPolTensor
A VibPolTensor instance is based on an Infrared instance, which in turns is based on a Phonons instance. Another important parameter is e_cut (in cm\(^{-1}\)), which corresponds to an energy cut-off: if phonons have a lower energy than e_cut, they are not considered in the computation of the vibrational polarizability tensor.
[2]:
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)
infrared = InfraredSpectrum(phonons)
vib_pt = VibPolTensor(infrared)
# The line above actually amounts to:
# vib_pt = VibPolTensor(infrared, e_cut=200)
You can then run the workflow without having to care about the context manager to run each calculation:
[3]:
vib_pt.run()
/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!
You can access the mean vibrational polarizability via the mean_polarizability attribute. The value is given in atomic units but can easily be converted to \(\unicode[serif]{xc5}^3\).
[4]:
from mybigdft.globals import B_TO_ANG
vib_pt.mean_polarizability*B_TO_ANG**3
[4]:
0.019682596858960125
This value is consistent with the experimental one reported in [2], that is 0.0178 \(\unicode[serif]{xc5}^3\).
References:
[2] Mark R. Pederson et al., *J. Chem. Theory Comput.* **1**, 4, (2005), pp 590–596