vibrational_spectra module

This module contains methods and classes relevant to computing frequencies and intensities of vibrational spectroscopy (IR, Raman) using double-harmonic approximation. This information can be later used for constructing vibrational spectra.

Tip

Only normal mode finite displacements NMFD/ENMFD can be used for computing intensities

Note

Currently, on IR intensities are implemented. Raman intensities would be implemented in future releases.

vibrational_spectra.add_broadening(freq, osc_str, line_profile='Lorentzian', line_param=10, step=10)

This function adds broadening to the line spectra (i.e. normal mode frequencies and intensities) to construct a vibrational spectra.

Arguments:

freq = A numpy 1D aray containing normal-mode frequencies(float).

osc_str = A numpy 1D array containing oscillator strength (float) computed using double-harmonic approximation

line_profile = A string: ‘Lorentzian’ (Default) or ‘Gaussian’

line_param = A float defining width of the profile

step = A float defining the size of the bin of the returned spectra

Returns:

A numpy array with 2 rows:

X-axis values (frequency) of the spectra as a numpy 1D float array

Y-axis values (intensity) of the spectra as a numpy 1D float array

vibrational_spectra.remove_dipole_jumps(dipole_moment, cell)

Removes the dipole moment jumps due to change in polarization lattice of a periodic system. The reference dipole moment is that of the geometry optimized configurations.

Arguments:

dipole_moment = numpy float array of shape (M,3) where M is the number of configurations

cell = numpy array of length 6. Lengths a,b,c and angles (degree):alpha, beta, gamma

Returns:

numpy float array of shape (M,3) containing corrected dipole moments

vibrational_spectra.difference(x, y, forward=True)

Calculate gradient using Finite Difference between 2 points. Finte difference can be calculated using forward or backward differences.

Arguments:

x = Array of displacement values

y = Array of functional values off shifted w.r.t x = 0

forward = if True then forward difference, else backward difference

Returns:

A numpy array containing gradients

vibrational_spectra.oscillator_strength(grad_dipole_x, grad_dipole_y, grad_dipole_z)

Calculate oscillator strengths from cartesian components of dipole-moment gradients

Arguments:

grad_dipole_x = Gradient of x-component of dipole moment

grad_dipole_y = Gradient of y-component of dipole moment

grad_dipole_z = Gradient of z-component of dipole moment

Returns:

A numpy array of oscillator strengths

vibrational_spectra.calc_ir_freq_osc_str(energies, dipole_total, disp_steps, ngrid=1, freq_unit='cm-1')

This function calculates the IR frequencies and oscillator strengths.

Arguments:

energies = A numpy float array of length 2*ngrid*M+1 containing displaced energies where M is the number of modes to be treated. The first element: energies[0] should be the energy of the optimized geometry

dipole_total = A (2*ngrid*M+1, 3) numpy float array of containing cartesian components of dipole moment. The first element: dipole_total[0] should be an array of length 3 containg cartesian components of the total dipole moment of the optimized geometry.

disp_steps = A numpy float array of length M containing symmetric displacement step along each normal mode.

ngrid = The number of displacement grid points. Allowed values: 1(Default), 2, 3, 4 .

freq_unit = The desired unit of frequency

Returns:

A numpy array of normal-mode frequencies (units desired by user) computed from the normal-mode hessian.

A numpy array of oscillator strengths.

If ngrid = 1; oscillator strength from forward diff; otherwise returns NONE.

If ngrid = 1; oscillator strength from backward diff; otherwise returns NONE.

For items 2-4, each array has 4 rows. First three rows are Cartesian X, Y and Z components of the oscillator strengths while the last row is their sum, i.e. total oscillator strengths.

class vibrational_spectra.ir(file_path, mode=None, deltax=None, deltae=None, ngrid=None, disp_indices=None, nmode_only=None, asr=None, diff_cut=None, ratio_cut=None, save_path='ir_freq_osc_str.npz')

Bases: object

This class contains methods for calculating IR frequencies and intensities and broadening methods to obtain a spectra. These quantities are calculated using normal mode finite displacements (NMFD/ENMFD/NMS/ENMS) and using finite-difference. The class can be instantiated using either (i) a normal mode displacement log file or (ii) pyepfd (or i-PI) XML restart file. Bellow are the description of various arguments it accepts, the objects it creates and the methods it contains.

Arguments:

file_path = This is a mandatory argument accepting a string. The string defines:

EITHER

(1) Full path to the log file (or output) while normal mode displacements (using NMFD, ENMFD, NMS or ENMS) are performed.

OR

(2) Full path to the pyepfd restart XML file obtained after finite difference (FD) frequency calculation, which was used as input for normal-mode displacements.

Other arguments are not mandatory as they have set default values but should be checked.

(A) Arguments if log file is supplied

disp_indices = A python list defining the points (from the displacement grid of a normal-mode sampling) that is to be used to calculate frequency and oscillator strength (intensity). This is needed only if ngrid is larger than 4 when you performed normal mode sampling using coord_util.ionic_mover class.

Tip

Suppose you used a normal mode sampling with ngrid=8, i.e, 16 symmetrically displaced points. The implemented finite difference formula accepts upto ngrid = 4. Therefore, you can chose the 8 closest points to supply this class using disp_indices = [5,6,7,8,9,10,11,12].

Note

Remember point 8 in this example is the negatively displaced configuration, while point 9 is positively displaced configuration. Check the normal mode displacement output (logfile) for further clarifications.

If you want to increase the displacement, you can use disp_indices = [1, 3, 5, 7, 10, 12, 14, 16]. Defining this, will automatically define ngrid = 4. If you want to use a lower value, say, ngrid = 2, you can set disp_indices = [5, 7, 10, 12]

(B) Arguments if only XML restart file is supplied

mode = Finite displacement mode

Options are:
(1) ‘NMFD’: Normal Mode Finite Displacement,

OR
(2) ‘ENMFD’: Energy-scaled Normal Mode Finite Displacement

deltax = A displacement (float) in atomic unit (Bohr)

deltae = Energy scaled displacement (float) in atomic_unit (Hartree)

ngrid = The number of displacement grid points. Allowed values: 1,2,3,4.

Warning

If these values are not supplied, then default values, which are: mode = 'ENMFD', deltax = 0.005, deltae = 0.001, ngrid = 1, would be used for calculating displacements.

Danger

If these values are not same to the ones supplied to coord_util.ionic_mover class while performing the displacements, it would lead to wrong frequencies and intensities.

(C) Argument irrespective of logfile or XML file

nmode_only = A python-list of normal-mode indices along which displacements are performed. Default is None, i.e., all modes are included.

Warning

If XML restart file is supplied as input then the supplied list must match to the one supplied to coord_util.ionic_mover while performing the displacements.

If log file is supplied, then the supplied list should be a subset of the sampled modes in the logfile.

asr = acoustic sum rule. Options are: (1) ‘none’, (2) ‘crystal’, (3) ‘poly’, (4) ‘lin’

Tip

See the documentation of coord_util.ionic_mover for more information.

Warning

If it is not supplied, either it would be read from the XML restart file, or set to default ‘none’ when log file is supplied.

Note

The following two cutoff’s are required to ensure numerical stability in periodic calculations. For a periodic system, the dipole moment is not a single vector but a lattice with several branches, and therefore, one must ensure that the dipole momentsof the displaced configurations entering the finite difference are in the same branch. If there is a jump, there would be considerable difference between the forward and backward difference formula and central difference won’t be numerically accurate. Using these cut offs we ensure, to use the correct difference formula while calculating the dipole moment gradient.

diff_cut = Default is None. Otherwise, a floating point number, defining the acceptable difference between forward and backward difference for oscillator strength when ngrid = 1. If the difference is less than the cutoff we accept the central difference.

ratio_cut = Default is None. Otherwise, a floating point number > 1, defining the range of acceptable ratio between forward and backward differences. If the ratio is within the acceptable range: (1, ratio_cut), central differance would be accepted.

save_path = full path to an numpy .npz file where frequencies and oscillator strength would be written if save() method is used. This file could be used to obtain the spectra using the add_broadening method described above.

arrange_energies_dipoles(energies, dipole_total)

This class method accepts numpy arrays containing energies and dipoles and arrange them to make it consistent with the displacements obtained using the other keys to instantiate the ir object.

Arguments:

energies = A 1D numpy float array of length M+1 containing the energies of the optimized geometry and the M displaced configurations

dipole_total = A 2D numpy float array of shape (M+1,3) containing the cartesian total dipole moment components of the the optimized geometry and the M displaced configurations

Returns:

Creats ir.energies and ir.dipole_total objects within the class.

calc_freq_osc_str(freq_unit='cm-1')

This class method calculates the frequency and oscilator strength using the the global method calc_ir_freq_osc_str() described above.

Arguments:

freq_unit = A string defining the desired unit of IR frequency

Returns the following objects:

omega = 1D numpy array of frequencies

osc_str = oscillator strengths using central difference

osc_str_f = oscillator strength using forward difference

osc_str_b = oscillator strength using forward difference

For osc_str, osc_str_f, and osc_str_b, each array has 4 rows. First three rows are Cartesian X, Y and Z components of the oscillator strengths while the last row is their sum, i.e. total oscillator strengths.

print_freq_osc_str(): This class method prints the frequency and oscillator strengths in the standard output

save()

This class method saves the omega_osc_str object as an numpy .npz file to the path defined by save_path key for the ir class.

Keys to access the data from the .npz files are:

omega = Frequencies

osc_str = Total oscillator strength

osc_str_x = Cartesian X component of the oscillator strength

osc_str_y = Cartesian Y component of the oscillator strength

osc_str_z = Cartesian Z component of the oscillator strength

calc_spectra(line_profile='Lorentzian', line_param=10, step=10, save_path='ir_spectra.npz')

This class method caculates the IR spectra from the class object omega_osc_str using the global method add_broadening().

Arguments:

line_profile = A string: ‘Lorentzian’ (Default) or ‘Gaussian’

line_param = A float defining width of the profile, Default value is 10.0

step = A float defining the size of the bin of the returned spectra, default value is 10.0

save_path = Full path to an .npz file where the spectra data would be saved. Default value is ir_spectra.npz.

Returns:

All of the following objects are numpy 2D arrays with 2 rows:

The first row contains frequencies (X-axis of the IR spectra).

The second row contains intensities (Y-axis of the spectra)

The created objects names (as well as keys to access the data from the *.npz. files) are:

spectra = Total IR spectra

spectra_x = IR spectra contributed by X-component of the oscillator strength

spectra_y = IR spectra contributed by Y-component of the oscillator strength

spectra_z = IR spectra contributed by Y-component of the oscillator strength

class vibrational_spectra.qbox_ir(file_path, qb_output_path, qb_prefix, n_qb_outfile=10000, n_qb_skip=0, corr_dipole_mom=True, mode=None, deltax=None, deltae=None, ngrid=None, disp_indices=None, nmode_only=None, asr=None, diff_cut=0.001, ratio_cut=3, save_path='ir_freq_osc_str.npz')

Bases: vibrational_spectra.ir

This class parses qbox outputs from NMFD/ENMFD/NMS/ENMS calculations and can calculate normal mode frequencies, IR oscillator strengths, IR spectra etc. Base class qbox_ir class is ir class described above. Therefore, qbox_ir contains all the methods that ir contains and their usage, keys arguments or internal objects are similar. Therefore, they wont be described again for qbox_ir. Below we will discuss the arguments that are specific to the qbox_ir class.

Arguments:

qb_output_path = This is a mandatory argument accepting a string defining the full directory path* to the qbox outputs containing NMFD/ENMFD/NMS/ENMS calculations

qb_prefix = This is another mandatory argument accepting a string defining the prefix of the qbox calculations.

Tip

If your qboux outputs are: enmfd-1.r, enmfd-2.r, enmfd-3.r,…, then your input would be qb_prefix = 'enmfd'.

n_qb_skip = The number of first qbox scf calculations (in each output file) to discard. The default value is 0.

Tip

Here is a scenario when we may need this argument by making it non-zero. For low band gap system, often the first SCF cannot be converged without set nempty <nonzero>. However, this should not be used while calculating dipole moment. In this case a work around is performing first SCF calculation with non-zero nempty value. After converging this calculations; we recompute dipole moment by setting the nempty to 0 for the first frame and subsequent frames. In this case, within one output, we should discard the first SCF calculation for the IR spectra calculation. We can do that by setting n_qb_skip = 1.

corr_dipole_mom = If true, dipole moment shifts (jumps) of the undisplaced structure with respect to the geometry optimized structure is corrected. Default is True.