elph_classes module

This module contains all the methods / classes relavant to frozen phonon calculation

elph_classes.bose_einstein(omega, omega_unit='Ha', T=0.0)

This function calculates the Bose Einstein occupation

Arguments:

omega = A float number signifying mode-frequencies.

omega_unit = Unit of omega. Options: 'Ha' (default) OR 'eV'.

T = Temperature in K

Returns:

A float number. Bose Einstein occupation for the given mode at T.

elph_classes.coord_com(coord, mass, flatten=True)

This function retruns coordinates with respect to center of mass (com) and com coordinates.

Arguments:

coord = 3N dimensional cartesian coordinates for N atoms

mass = Numpy array of mass matrix.

Note

This function can accept mass matrix with length N or 3 N

Returns:

1. coord_com = coordinates with respect to c.o.m

AND

2. com = center of mass coordinate

elph_classes.write_dynmat(dynmat, out_file_name, header='#\n')

This function can write the dynamical matrix row-wise.

Arguments:

dynmat = dynamical matrix, a (3 N x 3 N ) numpy array; where N is the number of atoms.

out_file_name = path to the output file

header = A header string

elph_classes.reorder_dynmat(dynmat, original_ord, new_ord, mass, atoms)

This function reorders a given dynamical matrix by changing its atom sequence. This is sometimes necessary if the underlyfing first-principles code such as qbox reorders the atom sequences differently in two different calculations: (i) an MD and (ii) a normal mode calculation. Later, if we want to compare the normal-modes with the MD trajectory, the modes would be inconsistent due to atom-ordering. In that case dynamical matrix may be reordered to the atom-sequence we want. For example, in this case we want the dynamical matrix to be consistent with the sequence of the MD simulation.

Arguments:

dynmat = Dynamical matrix; a numpy array of 3 N x 3 N; where N is the number of atoms.

original_ord = A python list of tuples defining the first and last atom index in each block

new_order = A python list of tuples defining the first and last atom index (in terms of original indices) in new block

mass = Mass matrix in original atom order

atoms = List of atom-symbols in original order

Returns:

(1) New dynamical matrix with desired atom ordering.

(2) New mass matrix with desired atom ordering.

(3) New atom-symbol list with desired atom ordering.

class elph_classes.central_diff(x, y, order=2, ngrid=4, decreasing=True)

Bases: object

This class calculates the 1st and 2nd order central difference of any quantity (Orbital energies, total energies, forces) when displacements are symmetric.

Arguments:

x = A numpy column vector containing the width of displacement grids.

y = A numpy column vector containing values of a function (can be energy or forces) whose derivatives are to be calculated. this must be clustered in such a way that first ngrid contiguous values related to x(0), next ngrid contiguous values related to x(1), and so on and so forth.

order = Order of central difference. Allowed values: (1) 1: Gradients, OR (2) 2: Hessian.

ngrid = Symmetric displacement grid for central difference. Allowed values: 1, 2, 3 or 4.

Note

Actual number of points would be double as we do symmetric displacements.

decreasing = If True: the displacements are in decreasing order , i.e +2h,+h,-h,-2h. False: Otherwise.

class elph_classes.dm(dynmat, mass)

Bases: object

This class stores the essential informations of a dynamical matrix. It also diagonalizes the dynamical matrix and stores the normal mode vectors as objects. This is class is used as a base class for many classes in PyEPFD. It has methods to transform from normal mode cordinates to Cartesian coordinates and vice versa.

Arguments:

dynmat = dynamical matrix, a numpy array of 3 N x 3 N; where N is the number of atoms.

mass = (3 N x 3 N ) mass matrix

Objects:

dynmatrix = Here the input dynamical matrix object is saved. This is mass-weighted Hessian matrix.

w2 = 1-D numpy array of Eigenvalues in ascending order

V[:,i] = i-th eigenvector

omega = array of normal mode frequencies (atomic unit).

natoms = Number of atoms

massinv = A 1-D array of inverse mass i.e., (1/sqrt(mass))

hessian = Hessian matrix without mass-weighting

calc_hessian()

This method calculates Hessian matrix without mass weighting and stores them into an object name hessian.

Note

While dynamical matrix object dynmatrix is Mass-weighted Hessian, here hessian object is not mass weighted.

nm2cart_disp(nm_disp)

This method take a a 3N-dim column vector of normal mode displacements, nm_disp, and transforms it into a 3N-dim column vec containing cartesian displacements, cart_disp.

nm2cart_matrix(Mnm)

Argument:

Mnm = Any 3 N x 3 N matrix represented in normal modes

Returns:

Mcart = Mnm represented in cartesian

cart2nm_vec(cart_v, normed=False, mass_weight=True)

This function projects a 3 N dimensional cartesian vector (coordinate or force) into normal modes and return the coefficients along the normal modes as a numpy array of length 3 N.

Arguments:

cart_v = A numpy array of length 3 N containing a 3 N dimensional Cartesian Vector (coordinates or forces of atoms of a particular configuration).

normed = Should the normal-mode vector be normalized to 1? Allowed values: (1) True OR (2) False (Default).

mass_weight = Should the considered normal modes be mass-weighted? Allowed values: (1) True (Default) OR (2) False

apply_asr(opt_coord, asr='none')

This method can be used to apply acoustic sum rule to a dynamical matrix to project out the global rotation and/or translation from the dynamical matrix. After applying this method, two several new objects would be created within the dm class containing the refined dynamical matrix, eigenvectors, frequencies.

Arguments:

opt_coord = 3N-dimensional vector of cartesian coordinates of a optimized geometry/structure

asr = acoustic sum rule; Options are: (1) 'none' (Default): asr is not applied, OR (2) 'crystal' : For infinite systems / crystals, OR (3) 'poly' : For poly-atomic non-linear molecules, OR (4) 'lin': For any linear molecules

Returns:

New objects within dm class.

refdynmat = Refined dynamical matrix after applying ASR

refw2 = Refined eigen values in increasing order.

refV = Refined eigen vectors. A (3 N x 3 N ) numpy array.

refomega = Refined mode-frequencies. A numpy array of length 3 N.

write_dynmat(prefix)

Prefix of the output file name

calc_free_en(T=0.0, unit='Ha')

Returns Helmholtz Free Energy of vibration (A = U - TS) in chosen unit.

Note

If T = 0; it returns the zero-point vibrational energy. Electronic energy is not included here. Must be added if needed.

Arg:

T = Temperature in K

unit = Unit of output (Helmholz Free Energy).

class elph_classes.stoch_displacements(dynmat, mass, asr='none', temperature=0, algo='osap', ngrid=1, nmode_only=None)

Bases: elph_classes.dm

This class calculates the stochastic displacements along normal modes given an algorithm.

Arguments:

dynmat = Dynamical matrix

mass = mass matrix

asr = Acoustic sum rule; allowed values: 'none', 'crystal', 'lin', 'poly'

temperature = Temperature for which displacements are chosen. A float.

algo = algorithm

Options:

(1) 'os' : one-shot,

OR

(2) 'osap' : one-shot with anthetic pair

OR

(3) 'osr' : oneshot where displacement signs chosen at random

OR

(4) 'osrap' : Same as 'osr' but for each sampled point it’s antethetic pair would be included.

OR

(5) 'mc' : Monte Carlo where displacements are chosen at random from a Gaussian

OR

(6) 'mcap' : The same as 'mc' but for each sampled point it’s antethetic pair would be included

ngrid = number of integration grid (sample points);

Note

with antethetic pairs (‘ap’) actual number of single point calculation (sample points) will be 2 x ngrid.

nmode_only = A python list of normal modes over which only the displacements would be performed.

Returns:

A new object nmdisp would be created containing the displacements in normal modes.

class elph_classes.nm_sym_displacements(dynmat, mass, mode, deltax, deltae)

Bases: elph_classes.dm

This class computes the values of normal mode displacements (consistent with i-PI code).

Arguments:

dynmat = dynamical matrix, a numpy array of 3 N x 3 N; where N is the number of atoms.

mass = mass matrix

mode = Finite displacement mode.

Options are:

(1) 'nmfd' : Normal Mode Finite Displacement,

OR

(2) 'enmfd' : Energy-scaled Normal Mode Finite Displacement.

Warning

Unlike ionic_mover_class other modes: NMS/ENMS/SD etc are not allowed.

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

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

Returns:

Creates a new object displacements, a numply float array of length 3 N containing displacement grid-width for each normal mode.

class elph_classes.phonon_calculator(forces, mass, ngrid=1, mode='fd', deltax=0.005, dynmat=None, deltae=0.001)

Bases: object

This class accepts forces and displacements and computes dynamical matrix as Jacobian of force.

Arguments:

ngrid = Number of displacement grid points for central difference, allowed values: 1, 2, 3 or 4.

Note

Actual number of points would be double as we do symmetric displacements.

forces = A (2 x ngrid, 3 x N) numpy array. Each row represents 3N-force components for a specific displaced struc.

mass = mass matrix

mode = Options: (1) 'fd', (2) 'nmfd' OR (3) 'enmfd`.

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

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

Important

Use same value for the mode, deltax and deltae that was used for creating the displacements using ionic_mover_class.

dynmat = 3N x 3N dynamical matrix. Must be supplied for all modes except mode=FD.

Returns:

It creates two new objects:

1. dynmat = Dynamical Matrix

AND

2. hessian = Hessian Matrix

class elph_classes.epce_calculator(dynmat, mass, energy, eigval, logfile, asr=None, mode='enmfd', deltax=0.001, deltae=0.001, vib_freq_unit='cm-1', epce_unit='meV')

Bases: object

This class calculates the all important electron phonon related properties, e.g., elceton phonon coupling energies (epce) for each phonon mode, total zero point renormalization (zpr) of an orbital eigenvalue, or thermal average of eigen values.

Warning

Currently, stencils are not implemented for this class. To compute EPCEs and ZPR, one must choose ngrid=1 during finite displacements using ionic_mover_class and phonon calculations using phonon_claculator class.

Arguments:

dynmat = A numpy float array (3 N x 3 N ) of dynamical matrix, where N is the number of atoms

mass = A numpy float array of mass in atomic unit (length = 3 N )

energy = An float array of length (6 N + 1) of Born-Oppenheimer energies. The first array element (energy[0]) should be the energy of the equilibrium geometry. Later array elements correspond to the 2 displaced structures along each normal modes, starting from the 1st mode.

eigval = A (6 N + 1 x n) numpy float array where n denotes the number of orbitals.

logfile = logfile path; a string.

asr = acoustaic sum rule. A string.

Options are:

(1) ‘none’ (Default): asr is not applied,

OR

(2) ‘crystal’: For infinite systems / crystals,

OR

(3) ‘poly’: For poly-atomic non-linear molecules,

OR

(4) ‘lin’: For any linear molecules

mode = Finite displacement mode used for displacements. Options: 'nmfd' or 'enmfd'

deltax = Displacement (in Bohr) used while displacing the ions. A float number.

deltae = Energy scaled displacement (float) in atomic_unit (Hartree) used while displacing the ions.

vib_freq_unit = Output unit for vibrational frequencies. Options are: (1) 'Ha' , (2) 'cm-1', (3) 'eV', (4) 'meV', (5) 'kcal/mol', (6) 'kJ/mol', (7) 'K', (8) 'THz'. Default: 'cm-1'

epce_unit = Output unit for EPCE. Options are same as vib_freq_unit Default: 'meV'

Returns:

Creates the following objects within the class:

vib_freq = A numpy float array of length 3 N containing vibrational frequencies computed from Normal Mode Hessian.

epce = A (3 N -m , n ) numpy float array containing EPCEs of n orbitals m = 0/3/5/6 depending on asr = 'none', asr = 'crystal', asr = 'lin' and asr = 'poly', respectively.

zpr = A numpy float array of length n containing zero-point renormalizations of n orbitals.

renorm_eigval_at_temp(T)

This method calculates renormalized band energy, for example Kohn-Sham eigenvalues by performing a thermal average using the Bose Einstein factor. It returns renormalized eigenvalues in Ha.

Argument: T = Temperatute (a float) in K

Returns: A numpy float array of length n containing renormalized band energies for n bands/orbitals at temperature T.