API Reference

Thermochemistry

class panthera.thermochemistry.Thermochemistry(vibenergies, atoms, *args, **kwargs)[source]

Calculate thermochemistry in harmonic approximation

Parameters:
vibenergiesnumpy.array

Vibrational energies in Joules

atomsase.Atoms

Atoms obect

phasestr

Phase, should be either gas or solid

pointgroupstr
symmetrynumberstr

If pointgroup is specified symmetrynumber is obsolete, since it will be inferred from the pointgroup

get_enthalpy(T=273.15)[source]

Return the enthalpy H

Parameters:
Tfloat

Temperature in K

get_entropy(T=273.15)[source]

Return the entropy S

Parameters:
Tfloat

Temperature in K

get_heat_capacity(T=273.15)[source]

Heat capacity at constant pressure

get_internal_energy(T=273.15)[source]

Return the internal energy U

Parameters:
Tfloat

Temperature in K

get_qvibrational(T=273.15, uselog=True)[source]

Calculate the vibrational partition function at temperature T in kJ/mol

Parameters:
Tfloat

Temperature in K

uselogbool

When True return the natural logarithm of the partition function

Notes

\[q_{vib}(T) = \prod^{3N-6}_{i=1}\frac{1}{1 - \exp(-h\omega_{i}/k_{B}T)}\]
get_vibrational_energy(T=273.15)[source]

Calculate the vibational energy correction at temperature T in kJ/mol

Parameters:
Tfloat

Temperature in K

Notes

\[U_{vib}(T) = \frac{R}{k_{B}}\sum^{3N-6}_{i=1} \frac{h\omega_{i}}{\exp(h\omega_{i}/k_{B}T) - 1}\]
get_vibrational_entropy(T=273.15)[source]

Calculate the vibrational entropy at temperature T in kJ/mol

Parameters:
Tfloat

Temperature in K

Notes

\[S_{vib}(T) = R\sum^{3N-6}_{i=1}\left[ \frac{h\omega_{i}}{k_{B}T(\exp(h\omega_{i}/k_{B}T) - 1)} - \ln(1 - \exp(-h\omega_{i}/k_{B}T)) \right]\]
get_vibrational_heat_capacity(T=273.15)[source]

Return the heat capacity

Parameters:
Tfloat

Temperature in K

Notes

\[C_{p,vib}(T) = R\sum^{3N-6}_{i=1} \left(\frac{h\omega_{i}}{k_{B}T}\right)^{2} \frac{\exp(-h\omega_{i}/k_{B}T)}{\left[1 - \exp(-h\omega_{i}/k_{B}T)\right]^{2}}\]
get_zpve()[source]

Calculate the Zero Point Vibrational Energy (ZPVE) in kJ/mol

Notes

\[E_{\text{ZPV}} = \frac{1}{2}\sum^{3N-6}_{i=1} h\omega_{i}\]
summary(T=273.15, p=0.1)[source]

Print summary with the thermochemical data at temperature T in kJ/mol

Parameters:
Tfloat

Temperature in K

pfloat

Pressure in MPa

NormalModeBFGS

class panthera.nmrelaxation.NormalModeBFGS(atoms, phase, hessian=None, hessian_update='BFGS', logfile='-', trajectory=None, restart=None, proj_translations=True, proj_rotations=True, master=None, verbose=False)[source]

Normal mode optimizer with approximate hessian update

Parameters:
atomsase.Atoms

Atoms object with the structure to optimize

phasestr

Phase, should be either gas or solid

hessianarray_like (N, N)

Initial hessian matrix in eV/Angstrom^2

hessian_updatestr

Name of the approximate formula to udpate hessian, one of: BFGS, SR1, DFP

proj_translationsbool

If True translational degrees of freedom will be projected from the hessian

proj_rotationsbool

If True rotational degrees of freedom will be projected from the hessian

logfilestr

Name log the log file

trajectorystr

Name of the trajectory file

log(grad, grad_nm, step_nm)[source]

Print a line with convergence information

log_header()[source]

Header for the log with convergence information

run(fmax=0.05, steps=100000000)[source]

Run structure optimization algorithm.

This method will return when the forces on all individual atoms are less than fmax or when the number of steps exceeds steps.

step(grad)[source]

Calculate the step in cartesian coordinates based on the step in normal modes in the rational function approximation (RFO)

Args:
gradarray_like (N,)

Current gradient

update_hessian(coords, grad)[source]

Perform hessian update

Parameters:
coordsarray_like (N,)

Current coordinates as vector

gradarray_like (N,)

Current gradient

panthera.anharmonicity

Methods for solving the one dimentional vibrational eigenproblem

panthera.anharmonicity.anharmonic_frequencies(atoms, T, coeffs, modeinfo)[source]

Calculate the anharmonic frequencies

Parameters:
atomsase.Atoms

Atoms object

Tfloat

Temperature in K

coeffspandas.DataFrame
modeinfopandas.DataFrame
panthera.anharmonicity.factsqrt(m, n)[source]

Return a factorial like constant

Parameters:
mint

Argument of the series

nint

Length of the series

Notes

\[f(m, n) = \prod^{n - 1}_{i = 0} \sqrt{m - i}\]
panthera.anharmonicity.get_anh_state_functions(eigenvals, T)[source]

Calculate the internal energy U and entropy S for an anharmonic vibrational mode with eigenvalues eigvals at temperature T in kJ/mol

\[ \begin{align}\begin{aligned}U &= N_{A}\frac{\sum^{n}_{i=1} \epsilon_{i}\exp(\epsilon_{i}/k_{B}T) }{\sum^{n}_{i=1} \exp(\epsilon_{i}/k_{B}T)}\\S &= N_{A}k_{B}\log(\sum^{n}_{i=1} \exp(\epsilon_{i}/k_{B}T)) + \frac{N_{A}}{T}\frac{\sum^{n}_{i=1} \epsilon_{i}\exp(\epsilon_{i}/k_{B}T) }{\sum^{n}_{i=1} \exp(\epsilon_{i}/k_{B}T)}\end{aligned}\end{align} \]
Parameters:
eigenvalsnumpy.array

Eigenvalues of the anharmonic 1D Hamiltonian in Joules

Tfloat

Temperature in K

Returns:
(U, S)tuple of floats

Tuple with the internal energy and entropy in kJ/mol

panthera.anharmonicity.get_hamiltonian(rank, freq, mass, coeffs)[source]

Compose the Hamiltonian matrix for the anharmonic oscillator with the potential described by the sixth order polynomial.

Parameters:
rankint

Rank of the Hamiltonian matrix

freqfloat

Fundamental frequency in hartrees

massfloat

Reduced mass of the mode

coeffsarray

A one dimensional array with polynomial coeffients

Notes

\[H_{ij} = \left\langle \Psi_{i} \left| \hat{H} \right| \Psi_{j} \right\rangle\]

where

\[\hat{H} = -\frac{\hbar^2}{2}\frac{\partial^2}{\partial \boldsymbol{Q}^2} + \sum_{\mu=0}^{6}c_{\mu}\boldsymbol{Q}^{\mu}\]

and \(\Psi_{i}\) are the standard harmonic oscillator functions.

panthera.anharmonicity.harmonic_df(modeinfo, T)[source]

Calculate per mode contributions to the thermodynamic functions in the harmonic approximation

Parameters:
modeinfopandas.DataFrame
Tfloat

Temperature in K

Returns:
dfpandas.DataFrame
panthera.anharmonicity.merge_vibs(anh6, anh4, harmonic, verbose=False)[source]

Form a DataFrame with the per mode thermochemical contributions from three separate dataframes with sixth order polynomial fitted potentia, fourth order fitted potential and harmonic frequencies.

Parameters:
anh6pandas.DataFrame
anh4pandas.DataFrame
harmonicpandas.DataFrame
Returns:
dfpandas.DataFrame

panthera.displacements module

panthera.displacements.calculate_displacements(atoms, hessian, freqs, normal_modes, npoints=4, modes='all')[source]

Calculate displacements in internal coordinates

Parameters:
atomsase.Atoms

Atoms object with the equilibrium structure

hessianarray_like

Hessian matrix in atomic units

freqsarray_like

Frequencies (square roots of the hessian eigenvalues) in atomic units

normal_modesarray_like

Normal modes in atomic units

npointsint

Number of points to displace structure, the code will calculate 2*npoints displacements since + and - directions are taken

modesstr or list/tuple of ints, default ‘all’

Range of the modes for which the displacements will be calculated

Returns:
imagesdict of dicts

A nested (ordred) dictionary with the structures with mode, point as keys, where point is a number from -4, -3, -2, -1, 1, 2, 3, 4

mipandas.DataFrame

DataFrame with per mode characteristics, displacements, masses and vibrational population analysis

panthera.displacements.get_internals_and_bmatrix(atoms)[source]

internals is a numpy record array with ‘type’ and ‘value’ records bmatrix is a numpy array n_int x n_cart

Parameters:
atomsase.Atoms

Atoms object

panthera.displacements.get_modeinfo(hessian, freqs, ndof, Bmatrix_inv, Dmatrix, mwevecs, npoints, internals)[source]

Compose a DataFrame with information about the vibrations, each mode corresponds to a separate row

panthera.displacements.get_nvibdof(atoms, proj_rotations, proj_translations, phase, include_constr=False)[source]

Calculate the number of vibrational degrees of freedom

Parameters:
atomsase.Atoms
proj_translationsbool

If True translational degrees of freedom will be projected from the hessian

proj_rotationsbool

If True rotational degrees of freedom will be projected from the hessian

include_constrbool

If True the constraints will be included

Returns:
nvibdoffloat

Number of vibrational degrees of freedom

panthera.displacements.vib_population(hessian, freqs, Bmatrix_inv, Dmatrix, internals, mi)[source]

Calculate the vibrational population analysis

Parameters:
hessianarray_like

Hessian matrix

freqsarray_like

A vector of frequencies (square roots fof hessian eigenvalues)

Bmatrix_invarray_like

Inverse of the B matrix

Dmatrixarray_like

D matrix

internalsarray_like

Structured array with internal coordinates

mipandas.DataFrame

Modeinfo

Returns:
mipandas.DataFrame

Modeinfo DataFrame updated with columns with vibrational population analysis results

panthera.io module

Module providing functions for reading the input and other related files

panthera.io.get_symmetry_number(pointgroup)[source]

Return the symmetry number for a given point group

See also

C. J. Cramer, Essentials of Computational Chemistry, Theories and Models, 2nd Edition, p. 363

Parameters:
pointgroupstr

Symbol of the point group

panthera.io.parse_arguments()[source]

Parse the input/config file name from the command line, parse the config and return the parameters.

panthera.io.print_mode_thermo(df, info=False)[source]

After calculating all the anharmonic modes print the per mode themochemical functions

panthera.io.print_modeinfo(mi, output=None)[source]

Print the vibrational population data

Parameters:
mipandas.DataFrame
outputstr

Name of the file to store the printout, if None stdout will be used

panthera.io.read_bmatdat()[source]

Read the bmat.dat file with internal coordiantes and the B matrix produced by the original writeBmat code

Returns:
internals, Bmatrixtuple

Internal coordiantes and B matrix

panthera.io.read_em_freq(fname)[source]

Read the file fname with the frequencies, reduced masses and fitted fitted coefficients for the potential into a pandas DataFrame.

Parameters:
fnamestr

Name of the file with PES

panthera.io.read_pes(fname)[source]

Parse the file with the potential energy surface (PES) into a dict of numpy arrays with mode numbers as keys

Parameters:
fnamestr

Name of the file with PES

panthera.io.read_poscars(filename)[source]

Read POSCARs file with the displaced structures and return an OrderedDict with the Atoms objects

panthera.io.read_vasp_hessian(outcar='OUTCAR', symmetrize=True, convert_to_au=True, dof_labels=False)[source]

Parse the hessian from the VASP OUTCAR file into a numpy array

Parameters:
outcarstr

Name of the VASP output, default is OUTCAR

symmetrizebool

If True the hessian will be symmetrized

convert_to_aubool

If True convert the hessian to atomic units, in the other case hessian is returned in [eV/Angstrom**2]

dof_labelsbool, default is False

If True a list of labels corresponding to the degrees of freedom will also be returned

Returns:
hessiannumpy.array

Hessian matrix

Notes

Note

By default VASP prints negative hessian so the elements should be multiplied by -1 to restore the original hessian, this is done by default, hessian in the XML file is NOT symmetrized by default

panthera.io.read_vasp_hessian_xml(xml='vasprun.xml', convert_to_au=True, stripmass=True)[source]

Parse the hessian from the VASP vasprun.xml file into a numpy array

Parameters:
xmlstr

Name of the VASP output, default is vasprun.xml

convert_to_aubool

If True convert the hessian to atomic units, in the other case hessian is returned in [eV/Angstrom**2]

dof_labelsbool, default is False

If True a list of labels corresponding to the degrees of freedom will also be returned

stripmassbool

If True use VASP default masses to transform hessian to non-mass-weighted form

Returns:
hessiannumpy.array

Hessian matrix

Notes

Note

By default VASP prints negative hessian so the elements should be multiplied by -1 to restore the original hessian, this is done by default, hessian in the XML file is symmetrized by default

panthera.io.write_modes(filename='POSCARs')[source]

Convert a file with multiple geometries representing vibrational modes in POSCAR/CONTCAR format into trajectory files with modes.

panthera.nmrelaxation module

class panthera.nmrelaxation.NormalModeBFGS(atoms, phase, hessian=None, hessian_update='BFGS', logfile='-', trajectory=None, restart=None, proj_translations=True, proj_rotations=True, master=None, verbose=False)[source]

Bases: Optimizer, object

Normal mode optimizer with approximate hessian update

Parameters:
atomsase.Atoms

Atoms object with the structure to optimize

phasestr

Phase, should be either gas or solid

hessianarray_like (N, N)

Initial hessian matrix in eV/Angstrom^2

hessian_updatestr

Name of the approximate formula to udpate hessian, one of: BFGS, SR1, DFP

proj_translationsbool

If True translational degrees of freedom will be projected from the hessian

proj_rotationsbool

If True rotational degrees of freedom will be projected from the hessian

logfilestr

Name log the log file

trajectorystr

Name of the trajectory file

log(grad, grad_nm, step_nm)[source]

Print a line with convergence information

log_header()[source]

Header for the log with convergence information

read()[source]
run(fmax=0.05, steps=100000000)[source]

Run structure optimization algorithm.

This method will return when the forces on all individual atoms are less than fmax or when the number of steps exceeds steps.

step(grad)[source]

Calculate the step in cartesian coordinates based on the step in normal modes in the rational function approximation (RFO)

Args:
gradarray_like (N,)

Current gradient

update_hessian(coords, grad)[source]

Perform hessian update

Parameters:
coordsarray_like (N,)

Current coordinates as vector

gradarray_like (N,)

Current gradient

panthera.nmrelaxation.nmoptimize(atoms, hessian, calc, phase, proj_translations=True, proj_rotations=True, gtol=1e-05, verbose=False, hessian_update='BFGS', steps=100000)[source]

Relax the strcture using normal mode displacements

Parameters:
atomsase.Atoms

Atoms object with the structure to optimize

hessianarray_like

Hessian matrix in eV/Angstrom^2

calcase.Calculator

ASE Calcualtor instance to be used to calculate forces

phasestr

Phase, ‘solid’ or ‘gas’

gtolfloat, default=1.0e-5

Energy gradient threshold

hessian_updatestr

Approximate formula to update hessian, possible values are ‘BFGS’, ‘SR1’ and ‘DFP’

stepsint

Maximal number of iteration to be performed

verbosebool

If True additional debug information will be printed

Notes

Internally eV and Angstroms are used.

See also

Bour, P., & Keiderling, T. A. (2002). Partial optimization of molecular geometry in normal coordinates and use as a tool for simulation of vibrational spectra. The Journal of Chemical Physics, 117(9), 4126. doi:10.1063/1.1498468

panthera.nmrelaxation.update_hessian(grad, grad_old, dx, hessian, update='BFGS')[source]

Perform hessian update

Parameters:
gradarray_like (N,)

Current gradient

grad_oldarray_like (N,)

Previous gradient

dxarray_like (N,)

Step vector x_n - x_(n-1)

hessianarray_like (N, N)

Hessian matrix

updatestr

Name of the hessian update to perform, possible values are ‘BFGS’, ‘SR1’ and ‘DFP’

Returns:
hessianarray_like

Update hessian matrix

panthera.panthera module

panthera.pes module

panthera.pes.calculate_energies(images, calc, modes='all')[source]

Given a set of images as a nested OrderedDict of Atoms objects and a calculator, calculate the energy for each displaced structure

Parameters:
imagesOrderedDict

A nested OrderedDict of displaced Atoms objects

calccalculator instance

ASE calculator

modesstr or list

Mode for which the PES will be calculated

Returns:
energiespandas.DataFrame

DataFrame with the energies per displacement

panthera.pes.differentiate(displacements, energies, order=2)[source]

Calculate numerical detivatives using the central difference formula

Parameters:
displacementsarray_like
energiesDataFrame
orderint

Order of the derivative

Notes

Central difference coefficients taken from [1]

[1]

Fornberg, B. (1988). Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of Computation, 51(184), 699-699. doi:10.1090/S0025-5718-1988-0935077-0

panthera.pes.expandrange(modestr)[source]

Convert a comma separated string of indices and dash separated ranges into a list of integer indices

Parameters:
modestrstr
Returns:
indiceslist of ints

Examples

>>> from panthera.pes import expandrange
>>> s = "2,3,5-10,20,25-30"
>>> expandrange(s)
[2, 3, 5, 6, 7, 8, 9, 10, 20, 25, 26, 27, 28, 29, 30]
panthera.pes.fit_potentials(modeinfo, energies)[source]

Fit the potentials with 6th and 4th order polynomials

Parameters:
modeinfopandas.DataFrame

DataFrame with per mode characteristics, displacements, masses and a flag to mark it a mode is a stretching mode or not

energiespd.DataFrame

Energies per displacement

Returns:
out(coeffs6o, coeffs4o)

DataFrames with 6th and 4th polynomial coefficients fitted to the potential

panthera.pes.harmonic_potential(x, freq, mu)[source]

Calculate the harmonic potential

Parameters:
xfloat of numpy.array

Coordinate

mufloat

Reduced mass

freqfloat

Frequency in cm^-1

panthera.plotting module

Functions for plotting the each mode and PES fits

panthera.plotting.plotmode(mode, energies, mi, c6o, c4o, output=None)[source]

Plot a given mode

Parameters:
modeint

Mode number (indexed from 0)

energiespandas.DataFrame
mipandas.DataFrame

Modeinfo

c6opandas.DataFrame
c4opandas.DataFrame
outputstr

name o file to store the plot

panthera.plotting.plotmode_legacy(mode, pes, coeff6, coeff4, output=None)[source]

Plot a given mode using legacy files

panthera.vibrations module

panthera.vibrations.get_levicivita()[source]

Get the Levi_civita symemtric tensor

panthera.vibrations.harmonic_vibrational_analysis(hessian, atoms, proj_translations=True, proj_rotations=False, ascomplex=True, massau=True)[source]

Given a force constant matrix (hessian) perform the harmonic vibrational analysis, by calculating the eigevalues and eigenvectors of the mass weighted hessian. Additionally projection of the translational and rotational degrees of freedom can be performed by specifying proj_translations and proj_rotations argsuments.

Parameters:
hessianarray_like

Force constant (Hessian) matrix in atomic units, should be square and symmetric

atomsAtoms

ASE atoms object

proj_translationsbool

If True translational degrees of freedom will be projected from the hessian

proj_rotationsbool

If True rotational degrees of freedom will be projected from the hessian

massaubool

If True atomic units of mass will be used

ascomplexbool

If there are complex eigenvalues return the array as complex type otherwise make the complex values negative and return array of reals

Returns:
out(w, v)

Tuple of numpy arrays with hessian square roots of the eigevalues (frequencies) and eiegenvectors in atomic units, both sorted in descending order of eigenvalues

panthera.vibrations.project(atoms, hessian, ndof, proj_translations=True, proj_rotations=False, verbose=False)[source]

Project out the translational and/or rotational degrees of freedom from the hessian.

Parameters:
atomsase.Atoms

Atoms object

ndofint

Number of degrees of freedom

hessianarray_like

Hessian/force constant matrix

proj_translationsbool

If True translational degrees of freedom will be projected from the hessian

proj_rotationsbool

If True rotational degrees of freedom will be projected from the hessian

Returns:
proj_hessianarray_like

Hessian matrix with translational and/or rotational degrees of freedom projected out

panthera.vibrations.project_massweighted(args, atoms, ndof, hessian, verbose=False)[source]

Project translational and or rotatioanl dgrees of freedom from mass weighted hessian