coord_util module

This module contains methods and classes dealing with atomic coordinates.

coord_util.abc2h(abc): Returns a lattice vector (in input units) matrix in upper triangular form i.e., each columns of the matrix represent the unit vectors given a description in terms of the lattice vector lengths and the angles(degree) in between.

coord_util.h2abc(h): Returns abc(a,b,c,alpha,beta,gamma) when cell vectors (h) in upper triangular matrix are provided. Return unit: angles in degree, lengths in supplied unit

coord_util.grep(file_path, pattern, cols): It greps a line from a file and then reads specic columns from that line and returns it as an array.

Arguments:

file_path = A string specifying a path of the file to be read

pattern = A pattern (string) that we look for

cols = A tuple of columns that we want to read from the grepped line

coord_util.quaternion_fit(ref_coord, coord, mass): This function fits a quaternion that rotates a set of 3N-coordinates into a a set of ref coordinates.

See also

Simon K. Kearsley, Acta Cryst. (1989). A45, 208-210

Arguments:

ref_coord = A vector of 3N cartesian coordinates for the reference structure

coord = A vector of 3N cartesian coordinates for the rotated rigid body

mass = A (3 N ) numpy array of mass

Returns:

A (3 x 3) rigid-body rotation matrix

coord_util.rotate_vector(rot_matrix, vector): Performs a rigid body rotation of a 3N dimensional vector given a rotation matrix, and returns resultant 3N dimensional vector.

Arguments:

rot_matrix = A (3 x 3) rigid-body rotation matrix

vector = A 3N dimensional vector

coord_util.remove_trans_rot(ref_coord, coord, mass, rotation=True, forces=None, wrt_com=False)

Rotates coordinates and forces(optional) so that the coordinates superimposed to a reference coordinates by rigid body rotation. This is to be used to remove rotation from a trajectory.

Arguments:

ref_coord = A 3N-dimensional vector of cartesian coordinates defining the reference configuration of the rigid body

coord = A 3N-dimensional vector of cartesian coordinates defining the rotated rigid body

mass = A (3N x 3N) matrix of mass of 3N particles.

rotation = Do you want to remove angular momentum/rotation? Allowed values are: (1) True (Default): Angular momentum would be removed from coord w.r.t ref_coord OR (2) False: Angular momentum would not be removed.

forces (optional) = None (Default) OR A 3N-dimensional vector of forces in rotated reference frame. If supplied, the forces would be rotated back to ref_coord frame

wrt_com = Do you want to obtain the coordinates with respect to center of mass (c.o.m)? Allowed values are: (1) True: Coordinates with respect to c.o.m would be returned OR (2) False (Default): Coordinates with respect to original refence point would be returned.

Returns:

A 3N-dimensional vector of Cartesian coordinates after back-rotation

A 3N-dimensional vector of Cartesian forces after back-rotation if forces is not None

class coord_util.xyz(file_path, io='r', atoms=None, xyz_unit='atomic_unit', cell_unit='atomic_unit', quantity='pos', reorder_seq=None, append=False)

Bases: object

This class contains methods related to xyz file. It can read a given xyz file or given atoms and coordinates it can create an xyz file.

Arguments:

file_path = path of the XYZ file

io = input/output status. Allowed values are: (1) ‘r’: In this mode the class would read an xyz file and store data; OR (2) ‘w’ = In this mode the class would write an xyz file.

atoms = A python list containing atom symbols. This is mandatory while io = 'w'.

xyz_unit = Position unit. Allowed values are: (1) ‘atomic_unit’ (default), OR (2) ‘angstrom’. This argument is optional for io = ‘w’ mode; while in io = 'r' mode this argument is not used.

cell_unit = Cell parameter unit. Allowed values are: (1) ‘atomic_unit’ (default), OR (2) ‘angstrom’. This argument is optional for io = 'w' mode; while in io = 'r' mode this argument is not used.

quantity = Allowed values are: (1) ‘pos’ (default): positions, OR (2) ‘mom’: momentums, OR (3) ‘force’: forces, OR (4) ‘vel’: velocities This is only valid for io = 'r' mode.

reorder_seq (optional) = A python list of atom symbols defining the order in which the information read from an xyz file would be saved. Only valid for io = 'r' mode.

append = If True, then it would append the new configurations if the xyz file exists. Valid only for io = 'w' mode

get_frame(xyz_lines=None, frame=1): This method returns coordinates and cells of a specified frame in atomic_unit

Arguments:

xyz_lines = A python list containing lines from an xyz file. Default value is None in which case all lines of the given xyz file is included.

frame = An integer defining the frame number from the given xyz trajectory which we would like to save.

write(cell, coord, append=False): This method writes the atoms & coordinates in an XYZ file with output unit angstrom cell parameters would be written in abcABC format.

Arguments:

cell = A numpy array of length 6 or 9 containing the cell-parameters. If the length of the array is 6, it will assume that the supplied format is abcABC (lengths and angles). Otherwise, it will assume that all three cell vectors are supplied.

coord = A 3N-dimensional vector of cartesian coordinates.

class coord_util.ionic_mover(atoms, opt_coord, mode, deltax=0.005, deltae=0.001, dynmat=None, mass=None, ngrid=1, temperature=0, asr='none', algo='osap', nmode_only=None)

Bases: object

This class moves the ions(nuclei) according to FD: Finite Displacement in Cartesian, OR NMFD: Normal Mode Finite Displacement, OR ENMFD: Energy-scaled Normal Mode Finite Displacement, OR SD: Stochastic Displacement.

Arguments:

atoms = An N-dimensional character array containing atom Symbols

opt_coord = 3N-dim column vector of cartesian coordinates of optimized geometry/struc.

mode = Finite displacement mode.

Options are:
(1) ‘FD’: Finite Displacement in Cartesian ,

OR
(2) ‘NMFD’: Normal Mode Finite Displacement,

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

OR
(4) ‘SD’: Stochastic Displacement,

OR
(5) ‘NMS’: Normal Mode Sampling,

OR
(6) ‘ENMS’: Energy-scaled Normal Mode Sampling

Tip

More Information
FD must be used for computing Cartesian normal-mode frequencies. The dynamical matrix obtained in this way, can be used for NMFD/ENMFD/SD/NMS/ENMS calculations. NMFD/ENMFD can be used to refine Cartesian frequencies or electron-phonon renormalizations. SD can be used to compute zero-point energy, vibrational energy including nuclear-quantum effects or electronic properties at a finite T. NMS/ENMS can be used to scan/sample a few specific normal modes and plot how energies/band gaps/ HOMO / LUMO changes along these few normal modes.

Note

For FD/NMFD/ENMFD/NMS/ENMS symmetric displacement is used.

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

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

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

mass = Mass-matrix

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

temperature = A float. Default value is 0. It is only used for mode = SD.

algo = Algorithms to use for mode = SD.

Options are:

(1) ‘OS’: One-Shot method proposed by Zacharius and Giustino,

See also

Zacharius, M.; Giustino, F. One-shot calculation of temperature-dependent optical spectra and phonon-induced band-gap renormalization. Phys. Rev. B. 2016, 94, 075125

OR
(2) ‘OSAP’: (Default): One-Shot method but with anti-thetic pair.

OR
(3) ‘OSR’: One-shot sampling of random signs as described in the following reference [1] or thermal line sampling in reference [2]

See also

[1] Karsai, F.; Engel, M.; Flage-Larsen, E; Kresse, G. Electron-phonon coupling in semiconductors within the GW approximation. New J. Phys. 2018, 20, 123008.

[2] Monserrat, B. Vibrational averages along thermal lines. Phys. Rev. B 2016, 93, 014302.

OR
(4) ‘OSRAP’: Same as OSR but including the antithetic-pairs.

OR
(5) ‘MC’: Monte-Carlo sampling along normal modes assuming harmonic probability density

OR
(6) ‘MCAP’: Same as MC but including the antithetic pairs.

ngrid = The number of displacement grid points.

Allowed values:
(A) 1-4 for mode = FD/NMFD/ENMFD. Default value is 1.

(B) Any positive integer is allowed for mode = NMS/ENMS.

(C) 1 for mode = SD and algo = OS/OSAP.

(D) Any positive integer for mode = SD and algo = OSR/OSRAP/MC/MCAP.

nmode_only (Optional) = A python-list of normal-mode indices along which displacements are performed.

Returns:
None. Stores the displaced coordinates within the object coord_util.ionic_mover.disp_coord

Note

In the set of displaced coordinates, undistorted structure is always included. Since FD/NMFD/NMS/ENMS employes symmetric displacement, for each normal mode additional 2*ngrid number of displaced coordinates would ve created.

For mode = SD, if algo = OS: only 1 additional displaced coordinate would be created. If algo = OSAP: only 2 additional displaced coordinates would be created. If algo = OSR/MC: ngrid number of additional displaced coordinates would be created. If algo = OSRAP/MCAP: 2*ngrid number of additional displaced coordinates would be created.

Important

MPI4Py Parallelization is available for mode = SD and mode = (E)NMS/(E)NMFD.

(1) When mode = SD, the code is parallelized over ngrid. Therefore the number of employed mpi processes should preferably be a divisor of ngrid and should not be larger than 1/2 of ngrid. For algo=OS/OSAP by definition ngrid=1. Therefore, mpiprocess should not give any scaling.

(2) When mode = (E)NMS/(E)NMFD, the code is parallelized over the number of modes to be sampled, defined by either nmode_only list (see above). If this is not set then all modes are taken into account. Therefor number of employed mpi processes should preferably be a divisor of no. of modes to sample and should not be larger than the half of it.

class coord_util.qbox(file_path, io='w', atoms=None, run_cmd=None, reorder=True, reorder_seq=None)

Bases: object

This class has methods related to qbox inputs and outputs

Arguments:

file_path = Path to the qbox output/input file.

io = input/output status.

Options are:
(1) ‘r’: In this mode the class would read a qbox output file and store the data;

OR
(2) ‘w’ = In this mode the class would write a qbox input file.

atoms = A list of atom symbols needed for io = 'w' mode

run_cmd = A list of qbox run commands

Caution

Qbox reorders the atom sequence and print the results with reordered atom sequence. This creates problem while post-processing the files, for example, obtaining normal modes. The next two arguments help to interpret and store qbox output data by nullyfying the reordering done by qbox and store them according to the atoms sequence user needs.

reorder = This key is valid while io = 'r'. Do you want to reorder the atom sequence?

options are:

True (Default): Reorders according to given sequence.

OR

False: Stores everything according to qbox output sequence.

reorder_seq = A python list of atom symbols specifying the atom sequence. Only valid if reorder = True.

Tip

More Information on how to use. For example reorder_seq = ['C','H','O','N'] would store all informations on Carbon, Hydrogen, Oxygen and Nitrogen atoms sequentially. For the default value None, atoms would be ordered according to the input sequence.

write(cell, coord)

This method writes the atoms & coordinates in a qbox input file. This method can be used only if qbox class is initiated with io = 'w'.

Arguments:

cell = cell vectors in atomic_unit. A numpy array of length 9.

coord = coordinates in atomic_unit. A numpy array of length 3*N*.

Returns
qbox input file(s) at the path specified by file_path while intitiating the qbox class.

getv(quantity): This method gets the vectors (forces , xyz positions, or dipole”) from a a qbox output and returns it as a 2D vector where each rows are 3N (forces, xyz) coordinates or 3 dipole moments of a snapshot/configuration. This function is only valid is the qbox class is initiated with io = 'r'.

Arguments:

quantity = '<force>' OR '<position>' OR '<dipole_ion>' OR '<dipole_el>' OR '<dipole_total>'

getenergy(): This method gets the energies from a a qbox output and returns it as a numpy array.

coord_util.write_nmode(atoms, cell_v, opt_coord, mode_v, mode_freq, file_path='dynmat', fmt='axsf'): It takes atoms, cell vectors, their positions (optimized) and normal mode (mass scaled) coordinates and normal mode frequencies as arguments and writes them into either an (i) animated xsf (.axsf) file that can be visualized with Xcrysden program OR (ii) An normal mode file (.nmd) that can be visualized with VMD program OR (iii) molden file (.molden) that can be visualized with molden program.

Arguments:

atoms = A list of atom-symbols. An character array with length N, where N is the number of atoms.

cell_v = Cell vectors in atomic_unit. A numpy float array with length 9.

opt_coord = Cartesian coordinates of optimized geometry in atomic units. A numpy array with length 3*N*.

mode_v = Normal mode vectors. A 3*N* x 3*N* dimensional numpy array.

mode_freq = A numpy array (len = N) of mode frequencies (in atomic unit).

file_path = Full directory path + file prefix of the file where normal mode info would be written

fmt = A format of the normal mode file. Options are: (1) 'axsf', OR (2) 'nmd', OR (3) 'molden'

Returns:

An .axsf/.nmd/.molden file with the specified directory path + file prefix.

coord_util.xyz2qe(xyzdata_path, pw_opt_path, frames, pw_path)

This function converts an xyz trajectory files into an Quantum Espresso (QE) input deck.

Arguments:

xyzdata_path = Full path to the xyzdata file

pw_opt_path = Full path to the file containing all pw options

frames = A tuple or an integer of frame indices. If tuple then indices of (start, end) or (start, end, inc) where inc is the increment

pw_path = Path to the directory where pw_inputs would be written and saved

Returns:

At the given pw_path, it creates several files: pw1.in, pw2.in,…, pw<n>.in for each frames in the trajectory

Warning

This function will be removed in the later releases. Therefore, for converting xyz files to QE input decks, it is recommended to use the QE class (see below) as shown in Tutorial 2.1.

class coord_util.qe(path, frames, io='r', pw_opt_path=None, pw_path='PWINPs', frame_prefix='frame', pw_prefix='pwscf', bands=None, spin=0)

Bases: object

This class has methods to convert xyz files to quantum espresso inputs and parse quantum espresso outputs

Arguments:

io = input/output status.

Options are:
(1) ‘r’: In this mode the class would read series of quantum espresso output files (from single-point calculations) and store the data;

OR
(2) ‘w’ = In this mode the class would write a series of quantum espresso input files each performing single-point calculations

path = (1) For io = 'w' mode, this is the path to the xyz file containing displaced

configurations (FD/ENMFD/SD etc)

OR
(2) For io = 'r' mode, this should be the root directory path where qe single- point calculations will be performed.

frames = A tuple or an integer of frame indices.
If tuple then indices of (start, end) or (start, end, inc) where inc is the increment. Same for both io = 'w'/'r'.

pw_opt_path = Full path to the file containing all pw input options
needed for io='w' mode.

pw_path = Path to the directory, required for io = 'w' mode.
Here pw input files would be written and saved.

frame_prefix = Active for io = 'r' mode.
A string defining the prefix of the directories where QE single point calculations are performed. The default value is 'frame', meaning within the root directory defined by path variable, there would be directories frame-1, frame-2,…, frame-n, where actual QE calculation and their outputs exist.

bands = An integrer or a tuple of integers specifying the frame
indices. If tuple then indices of (start, end) or (start, end, inc) where inc is the increment. It is an optional argument for io='r' mode. If used, the band energies will be extracted and written to a .dat file. For example, the name of the .dat file will be

is_0_orb-2-5.dat, if

spin = 0 (see below), and bands = (2,5,1) are chosen.

spin = An integer specifying the spin. It is an optional argument for
io='r' mode to decide band energies for which spin would be written to the .dat file. Options are: 0 = Spin-up (default), 1 = Spin-down (spin-polarized only), 2 = Both spins (spin-polarized only).

pw_prefix = A string that is used as quantum espresso prefix in a calculation.
This would create a file with name pw_prefix.xml, which will pe parsed by the class. The default value is 'pwscf'. This argument is active only for io='r' mode.