API Reference

Thermochemistry

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

Calculate thermochemistry in harmonic approximation

Parameters:

vibenergies : numpy.array

Vibrational energies in Joules

atoms : ase.Atoms

Atoms obect

phase : str

Phase, should be either gas or solid

pointgroup : str

symmetrynumber : str

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:

T : float

Temperature in K

get_entropy(T=273.15)[source]

Return the entropy S

Parameters:

T : float

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:

T : float

Temperature in K

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

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

Parameters:

T : float

Temperature in K

uselog : bool

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:

T : float

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:

T : float

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:

T : float

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:

T : float

Temperature in K

p : float

Pressure in MPa

NormalModeBFGS

class panther.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:

atoms : ase.Atoms

Atoms object with the structure to optimize

phase : str

Phase, should be either gas or solid

hessian : array_like (N, N)

Initial hessian matrix in eV/Angstrom^2

hessian_update : str

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

proj_translations : bool

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

proj_rotations : bool

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

logfile : str

Name log the log file

trajectory : str

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:
grad : array_like (N,)
Current gradient
update_hessian(coords, grad)[source]

Perform hessian update

Parameters:

coords : array_like (N,)

Current coordinates as vector

grad : array_like (N,)

Current gradient

panther.anharmonicity

Methods for solving the one dimentional vibrational eigenproblem

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

Calculate the anharmonic frequencies

Parameters:

atoms : ase.Atoms

Atoms object

T : float

Temperature in K

coeffs : pandas.DataFrame

modeinfo : pandas.DataFrame

panther.anharmonicity.factsqrt(m, n)[source]

Return a factorial like constant

Parameters:

m : int

Argument of the series

n : int

Length of the series

Notes

\[f(m, n) = \prod^{n - 1}_{i = 0} \sqrt{m - i}\]
panther.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:

eigenvals : numpy.array

Eigenvalues of the anharmonic 1D Hamiltonian in Joules

T : float

Temperature in K

Returns:

(U, S) : tuple of floats

Tuple with the internal energy and entropy in kJ/mol

panther.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:

rank : int

Rank of the Hamiltonian matrix

freq : float

Fundamental frequency in hartrees

mass : float

Reduced mass of the mode

coeffs : array

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.

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

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

Parameters:

modeinfo : pandas.DataFrame

T : float

Temperature in K

Returns:

df : pandas.DataFrame

panther.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:

anh6 : pandas.DataFrame

anh4 : pandas.DataFrame

harmonic : pandas.DataFrame

Returns:

df : pandas.DataFrame

panther.displacements module

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

Calculate displacements in internal coordinates

Parameters:

atoms : ase.Atoms

Atoms object with the equilibrium structure

hessian : array_like

Hessian matrix in atomic units

freqs : array_like

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

normal_modes : array_like

Normal modes in atomic units

npoints : int

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

modes : str or list/tuple of ints, default ‘all’

Range of the modes for which the displacements will be calculated

Returns:

images : dict 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

mi : pandas.DataFrame

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

panther.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:

atoms : ase.Atoms

Atoms object

panther.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

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

Calculate the number of vibrational degrees of freedom

Parameters:

atoms : ase.Atoms

proj_translations : bool

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

proj_rotations : bool

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

include_constr : bool

If True the constraints will be included

Returns:

nvibdof : float

Number of vibrational degrees of freedom

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

Calculate the vibrational population analysis

Parameters:

hessian : array_like

Hessian matrix

freqs : array_like

A vector of frequencies (square roots fof hessian eigenvalues)

Bmatrix_inv : array_like

Inverse of the B matrix

Dmatrix : array_like

D matrix

internals : array_like

Structured array with internal coordinates

mi : pandas.DataFrame

Modeinfo

Returns:

mi : pandas.DataFrame

Modeinfo DataFrame updated with columns with vibrational population analysis results

panther.io module

Module providing functions for reading the input and other related files

panther.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:

pointgroup : str

Symbol of the point group

panther.io.parse_arguments()[source]

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

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

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

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

Print the vibrational population data

Parameters:

mi : pandas.DataFrame

output : str

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

panther.io.read_bmatdat()[source]

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

Returns:

internals, Bmatrix : tuple

Internal coordiantes and B matrix

panther.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:

fname : str

Name of the file with PES

panther.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:

fname : str

Name of the file with PES

panther.io.read_poscars(filename)[source]

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

panther.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:

outcar : str

Name of the VASP output, default is OUTCAR

symmetrize : bool

If True the hessian will be symmetrized

convert_to_au : bool

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

dof_labels : bool, default is False

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

Returns:

hessian : numpy.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

panther.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:

xml : str

Name of the VASP output, default is vasprun.xml

convert_to_au : bool

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

dof_labels : bool, default is False

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

stripmass : bool

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

Returns:

hessian : numpy.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

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

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

panther.nmrelaxation module

class panther.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: ase.optimize.optimize.Optimizer, object

Normal mode optimizer with approximate hessian update

Parameters:

atoms : ase.Atoms

Atoms object with the structure to optimize

phase : str

Phase, should be either gas or solid

hessian : array_like (N, N)

Initial hessian matrix in eV/Angstrom^2

hessian_update : str

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

proj_translations : bool

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

proj_rotations : bool

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

logfile : str

Name log the log file

trajectory : str

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:
grad : array_like (N,)
Current gradient
update_hessian(coords, grad)[source]

Perform hessian update

Parameters:

coords : array_like (N,)

Current coordinates as vector

grad : array_like (N,)

Current gradient

panther.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:

atoms : ase.Atoms

Atoms object with the structure to optimize

hessian : array_like

Hessian matrix in eV/Angstrom^2

calc : ase.Calculator

ASE Calcualtor instance to be used to calculate forces

phase : str

Phase, ‘solid’ or ‘gas’

gtol : float, default=1.0e-5

Energy gradient threshold

hessian_update : str

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

steps : int

Maximal number of iteration to be performed

verbose : bool

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

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

Perform hessian update

Parameters:

grad : array_like (N,)

Current gradient

grad_old : array_like (N,)

Previous gradient

dx : array_like (N,)

Step vector x_n - x_(n-1)

hessian : array_like (N, N)

Hessian matrix

update : str

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

Returns:

hessian : array_like

Update hessian matrix

panther.panther module

Python package for Anharmonic Thermochemistry

panther.panther.main()[source]

The main Thermo program

panther.panther.temperature_range(conditions)[source]

Calculate the temperature grid from the input values and return them as numpy array

Parameters:

conditions : dict

Variable for conditions read from the input/config

Returns:

temps : numpy.array

Array with the temperature grid

panther.pes module

panther.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:

images : OrderedDict

A nested OrderedDict of displaced Atoms objects

calc : calculator instance

ASE calculator

modes : str or list

Mode for which the PES will be calculated

Returns:

energies : pandas.DataFrame

DataFrame with the energies per displacement

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

Calculate numerical detivatives using the central difference formula

Parameters:

displacements : array_like

energies : DataFrame

order : int

Order of the derivative

Notes

Central difference coefficients taken from [R11]

[R11]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
panther.pes.expandrange(modestr)[source]

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

Parameters:modestr : str
Returns:indices : list of ints

Examples

>>> from panther.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]
panther.pes.fit_potentials(modeinfo, energies)[source]

Fit the potentials with 6th and 4th order polynomials

Parameters:

modeinfo : pandas.DataFrame

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

energies : pd.DataFrame

Energies per displacement

Returns:

out : (coeffs6o, coeffs4o)

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

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

Calculate the harmonic potential

Parameters:

x : float of numpy.array

Coordinate

mu : float

Reduced mass

freq : float

Frequency in cm^-1

panther.plotting module

Functions for plotting the each mode and PES fits

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

Plot a given mode

Parameters:

mode : int

Mode number (indexed from 0)

energies : pandas.DataFrame

mi : pandas.DataFrame

Modeinfo

c6o : pandas.DataFrame

c4o : pandas.DataFrame

output : str

name o file to store the plot

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

Plot a given mode using legacy files

panther.vibrations module

panther.vibrations.get_levicivita()[source]

Get the Levi_civita symemtric tensor

panther.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:

hessian : array_like

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

atoms : Atoms

ASE atoms object

proj_translations : bool

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

proj_rotations : bool

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

massau : bool

If True atomic units of mass will be used

ascomplex : bool

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

panther.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:

atoms : ase.Atoms

Atoms object

ndof : int

Number of degrees of freedom

hessian : array_like

Hessian/force constant matrix

proj_translations : bool

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

proj_rotations : bool

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

Returns:

proj_hessian : array_like

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

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

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