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_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}}\]
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
- 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.
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 entropyS
for an anharmonic vibrational mode with eigenvalueseigvals
at temperatureT
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 originalwriteBmat
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.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
- 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.
- 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.plotting module¶
Functions for plotting the each mode and PES fits
panthera.vibrations module¶
- 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
andproj_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