"""
Utilities for the anharmonic corrections workflow.
Copyright Schrodinger, LLC. All rights reserved.
"""
# the implementation used in this module is based on the following:
# Beste et al., J. Phys. Chem. A, 2007, 111, 12118
# Beste, Chem. Phys. Lett., 2010, 493, 200
import argparse
import copy
import os
import warnings
from collections import OrderedDict
from collections import namedtuple
from types import SimpleNamespace
import numpy
from matplotlib.backends.backend_agg import FigureCanvasAgg
from matplotlib.figure import Figure
from scipy import constants
from scipy import integrate
from scipy import optimize
from schrodinger import structure
from schrodinger.application.jaguar.input import JaguarInput
from schrodinger.application.jaguar.output import JaguarOutput
from schrodinger.application.matsci import \
jaguar_multistage_workflow_utils as jmswfu
from schrodinger.application.matsci import jaguar_restart
from schrodinger.application.matsci import jaguarworkflows
from schrodinger.application.matsci import jobutils
from schrodinger.application.matsci import parserutils
from schrodinger.application.matsci import msprops
from schrodinger.job import jobcontrol
from schrodinger.utils import fileutils
from schrodinger.utils import units
FLAG_MAE_FILE = '-mae_file'
FLAG_JAG_OUT_FILE = '-jag_out_file'
FLAG_JAG_RIN_FILE = '-jag_restart_in_file'
FLAG_MAX_FREQ = '-max_freq'
FLAG_F_DATA = '-factor_data'
FLAG_JAG_KWARGS = '-jaguar_kwargs'
FLAG_T_DATA = '-temperature_data'
FLAG_P_DATA = '-pressure_data'
FLAG_MAX_I_FREQ = '-max_i_freq'
FLAG_PLOT = '-plot'
FLAG_PROCESS_NO_ANHARMONICITIES = '-process_no_anharmonicities'
FLAG_TPP = parserutils.FLAG_TPP
# frequencies in wavenumbers (cm^-1)
DEFAULT_MAX_FREQ = 300
# unitless factors that multiply the normal mode displacement
DEFAULT_F_START = 0.5
DEFAULT_F_STEP = 1.0
DEFAULT_F_N_POINTS = 4
DEFAULT_JAGUAR_KWARGS = OrderedDict([('dftname', 'B3LYP'), ('basis', 'LACVP**'),
('igeopt', 1), ('molchg', 0),
('multip', 1)])
# K
DEFAULT_T_START = jmswfu.DEFAULT_TEMP_START
DEFAULT_T_STEP = jmswfu.DEFAULT_TEMP_STEP
DEFAULT_T_N_POINTS = jmswfu.DEFAULT_TEMP_N
# atm
DEFAULT_P_START = jmswfu.DEFAULT_PRESS_START
DEFAULT_P_STEP = jmswfu.DEFAULT_PRESS_STEP
DEFAULT_P_N_POINTS = jmswfu.DEFAULT_PRESS_N
# frequencies in wavenumbers (cm^-1), this handles small imaginary
# frequencies
DEFAULT_MAX_I_FREQ = 0
FREQUENCY_TAG = '_freq'
SINGLE_POINT_TAG = '_sp'
# hartree
ENERGY_ATTR = 'scf_energy'
# unitless
CORRECTION_THRESH = 1e-6
INFINITY_THRESH = 1e6
# kcal/mol/K
IDEAL_GAS_CONSTANT = (constants.R / constants.calorie) / constants.kilo
KG_PER_AMU = constants.physical_constants['atomic mass constant'][0]
KEY_STARTER = 'r_matsci_anharmonic_'
U_STARTER = KEY_STARTER + 'Total_Internal_Energy_(au)'
H_STARTER = KEY_STARTER + 'Total_Enthalpy_(au)'
G_STARTER = KEY_STARTER + 'Total_Free_Energy_(au)'
LNQVIB_KEY = KEY_STARTER + 'ln(Qvib)_{temp}K_{press}atm'
LNQ_KEY = KEY_STARTER + 'ln(Q)_{temp}K_{press}atm'
U_KEY = U_STARTER + jaguarworkflows.TEMP_PRESS_KEY_EXT
H_KEY = H_STARTER + jaguarworkflows.TEMP_PRESS_KEY_EXT
G_KEY = G_STARTER + jaguarworkflows.TEMP_PRESS_KEY_EXT
SeqData = namedtuple('SeqData', ['start', 'step', 'n_points'])
OUT_EXT = '_out.mae'
[docs]def get_seq_data(options, flag):
"""
Return a sequence data for the given flag.
:type options: argparse.Namespace
:param options: the options
:type flag: str
:param flag: the flag
:rtype: SeqData
:return: the sequence data
"""
# data was typed using parserutils.type_positive_float
start, step, n_points = jobutils.get_option(options, flag)
return SeqData(start=start, step=step, n_points=int(n_points))
[docs]def evaluate_f(x, deriv_idx, coeffs):
"""
Evaluate the nth derivative of a polynomial described by the given
coefficients.
:type x: float
:param x: the point at which to evaluate
:type deriv_idx: int
:param deriv_idx: indicates what derivative of the polynomial to
evaluate, 0 is the polynomial itself, 1 is the first derivative, etc.
:type coeffs: tuple
:param coeffs: the coefficents of the polynomial, for a mth order polynomial
must be of lenth m + 1
:rtype: float
:return: the evaluated value
"""
# f(x) = a_0 + a_1 * x + a_2 * x ** 2 + a_3 * x ** 3 + ...
assert deriv_idx < len(coeffs)
value = 0
for idx, coeff in enumerate(coeffs):
prod = numpy.prod([idx - jdx for jdx in range(deriv_idx)])
if not prod:
continue
value += coeff * prod * x**(idx - deriv_idx)
return value
[docs]def angular_freq_to_freq(angular_freq):
"""
Convert the given angular frequency to frequency.
:type angular_freq: float
:param angular_freq: the angular frequency in s**-1
:rtype: float
:return: the frequency in cm**-1
"""
return angular_freq / (2 * constants.pi * constants.c * 100)
[docs]def freq_to_angular_freq(freq):
"""
Convert the given frequency to angular frequency.
:type freq: float
:param freq: the frequency in cm**-1
:rtype: float
:return: the angular frequency in s**-1
"""
return freq * 2 * constants.pi * constants.c * 100
[docs]def plotter(x_min,
x_max,
x_e_min,
x_e_max,
x_step,
y_func,
file_name,
title,
y_label,
x_values=None):
"""
Plot the given function.
:type x_min: float
:param x_min: the minimum value on the x-axis
:type x_max: float
:param x_max: the maximum value on the x-axis
:type x_e_min: float
:param x_e_min: the minimum value on the extended x-axis
:type x_e_max: float
:param x_e_max: the maximum value on the extended x-axis
:type x_step: float
:param x_step: the step size to use on the x-axis
:type y_func: function
:param y_func: the function to use to obtain y-axis values
:type file_name: str
:param file_name: the file name used to write the plot image
:type title: str
:param title: the title for the plot image
:type y_label: str
:param y_label: the y-axis label for the plot image
:type x_values: list or None
:param x_values: if not None then contains x values for points
to show in the plot
"""
n_points = int((x_e_max - x_e_min) / x_step)
xs = numpy.linspace(x_e_min, x_e_max, num=n_points)
ys = numpy.array(list(map(y_func, xs)))
# prevent infinities which can happen for plotting
# since the range is extended beyond the fit range
if max(ys) > INFINITY_THRESH:
pairs = list(zip(xs, ys))
y_max = max(pairs, key=lambda x: x[1]
if x_min <= x[0] <= x_max else 0)[1]
pairs = [(x, y) for x, y in pairs if y <= 10 * y_max]
xs, ys = zip(*pairs)
title += ' Infinities removed.'
y_min, y_max = min(ys), max(ys)
# Agg rather than QTAgg is supposed to be default backend (MATSCI-9819) to prevent
# requiring graphics libraries in case running on a remote compute host,
# the following is taken from schrodinger.application.desmond.mplchart.get_xy_plot
figure = Figure(figsize=(8, 4), dpi=100)
canvas = FigureCanvasAgg(figure)
plot = figure.add_axes([.1, .15, .85, .75])
plot.set_xlabel('Factor')
plot.set_ylabel(y_label)
plot.set_xlim((x_e_min, x_e_max))
plot.set_ylim((y_min, y_max))
plot.set_title(title)
plot.plot(xs, ys, color='k')
if x_values:
y_values = list(map(y_func, x_values))
plot.plot(x_values, y_values, 'bo')
figure.savefig(file_name)
backend = jobcontrol.get_backend()
if backend:
backend.copyOutputFile(file_name)
[docs]def get_normal_modes(jagout, max_i_freq=-numpy.inf):
"""
Return the normal modes from the given JaguarOutput.
:type jagout: schrodinger.application.jaguar.output.JaguarOutput
:param jagout: the Jaguar output object
:type max_i_freq: float
:param max_i_freq: tolerate small imaginary frequencies less than this value in
wavenumbers (cm^-1)
:rtype: list
:return: contains pair tuples, (normal mode index (1-based),
schrodinger.application.jaguar.results.NormalMode)
"""
# for non-frequency jobs it was observed that jagout.normal_mode can
# be an instance of NormalMode rather than an empty list, None, etc.
if not isinstance(jagout.normal_mode, list):
return []
normal_modes = []
for idx, normal_mode in enumerate(jagout.normal_mode, 1):
if normal_mode.frequency < max_i_freq:
continue
else:
normal_modes.append((idx, normal_mode))
return normal_modes
[docs]def check_imaginary_frequencies(jag_out, jag_in, max_i_freq=DEFAULT_MAX_I_FREQ):
"""
Check imaginary frequencies.
:type jag_out: JaguarOutput
:param jag_out: the Jaguar output object
:type jag_in: JaguarInput
:param jag_in: the Jaguar input object
:type max_i_freq: float
:param max_i_freq: tolerate small imaginary frequencies less than this value in
wavenumbers (cm^-1)
:raise AnharmonicException: if there is an issue
"""
is_ts = jag_in.getValue('igeopt') == 2
n_tss = 0
for idx, normal_mode in get_normal_modes(jag_out):
if normal_mode.frequency < max_i_freq:
if not is_ts:
msg = (f'An imaginary frequency less than {max_i_freq} '
'has been found.')
raise AnharmonicException(msg)
else:
n_tss += 1
if n_tss > 1:
msg = ('Multiple imaginary frequencies less than '
f'{max_i_freq} have been found.')
raise AnharmonicException(msg)
if is_ts and not n_tss:
msg = (f'An imaginary frequency less than {max_i_freq} '
'has not been found.')
raise AnharmonicException(msg)
[docs]def get_st_jaguar_output(jagout_file_name, allow_new_dummies=False):
"""
Return a structure from the given Jaguar output file.
:type jagout_file_name: str
:param jagout_file_name: the name of a Jaguar output file
:type allow_new_dummies: bool
:param allow_new_dummies: whether to allow mmjag's lewis structure
build to possibly add new dummy atoms
:rtype: schrodinger.structure.Structure
:return: the structure
"""
jagout = JaguarOutput(jagout_file_name)
st = jagout.getStructure()
st_natoms = st.atom_total
file_natoms = jagout.atom_total
if st_natoms > file_natoms and not allow_new_dummies:
new_dummy_idxs = list(range(file_natoms + 1, st_natoms + 1))
st.deleteAtoms(new_dummy_idxs)
return st
[docs]class AnharmonicException(Exception):
pass
[docs]class AnharmonicPotentials(object):
[docs] def __init__(self,
st=None,
jagout_file_name=None,
jagrin_file_name=None,
max_freq=DEFAULT_MAX_FREQ,
factor_data=None,
jaguar_kwargs=DEFAULT_JAGUAR_KWARGS,
temperature_data=None,
pressure_data=None,
max_i_freq=DEFAULT_MAX_I_FREQ,
plot=False,
process_no_anharmonicities=False,
tpp=1,
logger=None,
robust=False):
"""
Create an instance.
:type st: `schrodinger.structure.Structure` or None
:param st: a structure for which to calculate anharmonic potentials or
None if using Jaguar frequency files directly
:type jagout_file_name: str or None
:param jagout_file_name: the name of a Jaguar frequency output file for
which to calculate anharmonic potentials or None if using an input
structure
:type jagrin_file_name: str or None
:param jagrin_file_name: the name of a Jaguar freqency restart input file
for which to calculate anharmonic potentials or None if using an input
structure
:type max_freq: float
:param max_freq: anharmonic potentials are created for normal modes with
harmonic frequencies less than this value in wavenumbers (cm^-1)
:type factor_data: SeqData or None
:param factor_data: unitless factor data for factors that multiply
a normal mode displacement, if None then the defaults are used,
the number of points is in the positive direction only, excluding
zero and the negative direction, for example using a value
of 4 in turn means 2 * 4 + 1 = 9 points total
:type jaguar_kwargs: dict
:param jaguar_kwargs: Jaguar &gen section keyword arguments, used only if the
anharmonic potentials are being calculated from an input structure rather
than directly from Jaguar frequency files
:type temperature_data: SeqData or None
:param temperature_data: temperature data in K, if None then the defaults
are used
:type pressure_data: SeqData or None
:param pressure_data: pressure data in atm, if None then the defaults
are used
:type max_i_freq: float
:param max_i_freq: tolerate small imaginary frequencies less than this value in
wavenumbers (cm^-1)
:type plot: bool
:param plot: if True then return plots of the potentials and particle densities
:type process_no_anharmonicities: bool
:param process_no_anharmonicities: if True then do not exit with an error if
the given max_freq results in zero normal modes to be treated anharmonically
:type tpp: int
:param tpp: the threads-per-process to use for running Jaguar calculations
:type logger: logging.Logger or None
:param logger: output logger or None if there isn't one
:param bool robust: If True, use the robust Jaguar driver to run Jaguar
jobs. If False, use Jaguar directly.
"""
jag_file_names = [jagout_file_name, jagrin_file_name]
assert (st and not any(jag_file_names)) or (not st and
all(jag_file_names))
self.st = st
if self.st:
jagout_file_name = self.st.property.get(msprops.JAGUAR_OUTPATH_PROP)
jagrin_file_name = self.st.property.get(
's_j_Jaguar_restart_input_file')
jag_file_names = [jagout_file_name, jagrin_file_name]
if all(jag_file_names):
self.jagout, self.jagin = self._getJaguarObjs(jagout_file_name,
jagrin_file_name,
check_freq=True)
if not self.st:
self.st = get_st_jaguar_output(jagout_file_name,
allow_new_dummies=False)
self.jagin.setStructure(self.st, set_molchg_multip=False)
base_name = fileutils.get_basename(jagout_file_name)
else:
self.jagout = None
self.jagin = None
base_name = self.st.title.strip()
if not base_name:
base_name = 'anharmonic'
self.base_name = jobutils.get_jobname(base_name)
self.max_freq = max_freq
self.factor_data = factor_data or SeqData(start=DEFAULT_F_START,
step=DEFAULT_F_STEP,
n_points=DEFAULT_F_N_POINTS)
self.jaguar_kwargs = jaguar_kwargs
self.temperature_data = temperature_data or SeqData(
start=DEFAULT_T_START,
step=DEFAULT_T_STEP,
n_points=DEFAULT_T_N_POINTS)
self.pressure_data = pressure_data or SeqData(
start=DEFAULT_P_START,
step=DEFAULT_P_STEP,
n_points=DEFAULT_P_N_POINTS)
self.max_i_freq = max_i_freq
self.plot = plot
self.process_no_anharmonicities = process_no_anharmonicities
self.tpp = tpp
self.logger = logger
self.robust = robust
options = argparse.Namespace(TPP=self.tpp)
self.job_dj = jobutils.create_queue(options=options,
host=jobutils.AUTOHOST)
self.potentials = {}
def _getJaguarObjs(self,
jagout_file_name,
jagrin_file_name,
check_freq=True):
"""
Return the Jaguar output and input objects for the given files.
:type jagout_file_name: str
:param jagout_file_name: the name of a Jaguar output file for
which to calculate anharmonic potentials
:type jagrin_file_name: str
:param jagrin_file_name: the name of a Jaguar restart input file
for which to calculate anharmonic potentials
:type check_freq: bool
:param check_freq: whether to ensure that the given Jaguar output file is
for a frequency job
:raise AnharmonicException: if there is an issue
:rtype: JaguarOutput, JaguarInput
:return: the Jaguar output and input objects
"""
jag_file_names = [jagout_file_name, jagrin_file_name]
if not all(map(os.path.exists, jag_file_names)):
msg = (
f'At least one of the Jaguar output files {jagout_file_name} '
f'and {jagrin_file_name} can not be found.')
raise AnharmonicException(msg)
jagout = JaguarOutput(jagout_file_name)
if jagout.status != JaguarOutput.OK or jagout.fatal_error:
msg = (f'The Jaguar output file {jagout_file_name} comes from '
'a failed job.')
raise AnharmonicException(msg)
# see MATSCI-12044 - allow atomic systems to fall through
st = jagout.getStructure()
if check_freq and st.atom_total > 1 and not get_normal_modes(jagout):
msg = (f'The Jaguar output file {jagout_file_name} must be from '
'a frequency job.')
raise AnharmonicException(msg)
return jagout, JaguarInput(input=jagrin_file_name)
[docs] def getJaguarJob(self, base_name, input_name):
"""
Get the job to run Jaguar with the given job base name and input
name. The job will run either Jaguar directly or the robust Jaguar
driver based on the current robust setting.
:param str base_name: The base name of this job
:param str input_name: The name of the input file
:rtype: `jobutils.RobustSubmissionJob`
:return: The job (suitable for JobDJ that will run Jaguar
"""
if self.robust:
cmd = [
jaguarworkflows.ROBUST_JAGUAR_DRIVER_PATH, '-JOBNAME',
f'{base_name}{jaguarworkflows.ROBUST_JOB_ENDING}'
]
else:
cmd = ['jaguar', 'run']
cmd += ['-TPP', str(self.tpp), input_name]
return jobutils.RobustSubmissionJob(cmd)
[docs] def runFrequencyJob(self):
"""
Run a Jaguar frequency job on the input structure.
:raise AnharmonicException: if there is an issue
"""
base_name = f'{self.base_name}{FREQUENCY_TAG}'
jagout_file_name = f'{base_name}.out'
jagrin_file_name = f'{base_name}.01.in'
jagrmae_file_name = f'{base_name}.01.mae'
# handle restarts
if jaguar_restart.needs_restart(base_name):
jagin_file_name = f'{base_name}.in'
# temperature and pressure data are not needed for potentials
# but will be needed for partition functions
jaguar_kwargs = self.jaguar_kwargs.copy()
jaguar_kwargs['ifreq'] = 1
jaguar_kwargs['tmpini'] = self.temperature_data.start
jaguar_kwargs['tmpstp'] = self.temperature_data.step
jaguar_kwargs['ntemp'] = self.temperature_data.n_points
jaguar_kwargs['press'] = self.pressure_data.start
jaguar_kwargs['press_step'] = self.pressure_data.step
jaguar_kwargs['npress'] = self.pressure_data.n_points
jagin = JaguarInput(structure=self.st, genkeys=jaguar_kwargs)
jagin.saveAs(jagin_file_name)
rsjob = self.getJaguarJob(base_name, jagin_file_name)
self.job_dj.addJob(rsjob)
self.job_dj.run()
if self.robust:
# The robust protocol may have run a number of sub-subjobs in
# multiple subdirectories, so make sure they all get added.
jobutils.add_subjob_files_to_backend(rsjob)
jaguarworkflows.add_jaguar_files_to_jc_backend(
base_name, others=['.01.smap', '.vib', '_vib.spm'])
try:
self.jagout, self.jagin = self._getJaguarObjs(jagout_file_name,
jagrin_file_name,
check_freq=True)
except AnharmonicException as err:
msg = (' The Jaguar frequency job died.')
raise AnharmonicException(str(err) + msg)
if os.path.exists(jagrmae_file_name):
self.st = structure.Structure.read(jagrmae_file_name)
else:
msg = (f'The Jaguar restart Maestro file {jagrmae_file_name} '
'can not be found.')
raise AnharmonicException(msg)
[docs] @staticmethod
def getFactors(factor_data):
"""
Return the factors.
:type factor_data: SeqData
:param factor_data: unitless factor data for factors that multiply
a normal mode displacement, the number of points is in the
positive direction only, excluding zero and the negative direction,
for example using a value of 4 in turn means 2 * 4 + 1 = 9 points
total
:rtype: tuple
:return: the factors
"""
# do not assume symmetry in the Jaguar energies and so do
# all factors, below are some values for a typical example system:
#
# factor SCF energy/H
# -3.5 -1152.50911446245
# -2.5 -1152.50987075683
# -1.5 -1152.51009115970
# -0.5 -1152.51011514460
# 0.0 -1152.51010128873
# 0.5 -1152.51010160525
# 1.5 -1152.51004642757
# 2.5 -1152.50977772989
# 3.5 -1152.50894987833
factors = numpy.array([
factor_data.start + i * factor_data.step
for i in range(factor_data.n_points)
])
return tuple(numpy.append(-numpy.flip(factors), factors))
[docs] def getExtendedFactors(self):
"""
Return the extended factors.
:rtype: tuple
:return: the extended factors
"""
# slightly extend the range of factors for plotting
start = self.factor_data.start
step = self.factor_data.step
n_points = 2 * self.factor_data.n_points
factor_data = SeqData(start=start, step=step, n_points=n_points)
return self.getFactors(factor_data)
[docs] def runSinglePointJobs(self):
"""
Run the Jaguar single point jobs from which to calculate the anharmonic
potentials.
:raise AnharmonicException: if there is an issue
"""
base_name = f'{self.base_name}{SINGLE_POINT_TAG}'
# Ang.
xyz = self.st.getXYZ()
st = self.st.copy()
jagin = copy.deepcopy(self.jagin)
jagin.setValue('igeopt', 0)
jagin.setValue('ifreq', 0)
# imaginary frequencies do not contribute to the vibrational
# partition function
all_rsjobs = []
for idx, normal_mode in get_normal_modes(self.jagout,
max_i_freq=self.max_i_freq):
if normal_mode.frequency > self.max_freq:
continue
base_names = []
for factor in self.getFactors(self.factor_data):
_factor = round(factor, 4)
base_name_p = f'{base_name}_{idx}_normal_mode_{_factor}_factor'
base_names.append(base_name_p)
# handle restarts
if not jaguar_restart.needs_restart(base_name_p):
continue
xyz_p = xyz + factor * normal_mode.displacement
st.setXYZ(xyz_p)
# use molchg and multip from the Jaguar output
jagin.setStructure(st, set_molchg_multip=False)
jagin_file_name = f'{base_name_p}.in'
jagin.saveAs(jagin_file_name)
rsjob = self.getJaguarJob(base_name_p, jagin_file_name)
all_rsjobs.append(rsjob)
self.job_dj.addJob(rsjob)
self.potentials[idx] = SimpleNamespace(base_names=tuple(base_names))
if not self.potentials:
msg = (f'There are no frequencies less than {self.max_freq}.')
if not self.process_no_anharmonicities:
raise AnharmonicException(msg)
elif self.logger:
self.logger.info(
f'{msg} Continuing with the harmonic approximation.')
# currently all single points use the zero point as the initial guess,
# however if negative curvature is observed the initial guesses could
# be followed, i.e. k from k - 1 rather than 0
if self.job_dj.total_added:
self.job_dj.run()
if self.robust:
for rsjob in all_rsjobs:
# The robust protocol may have run a number of sub-subjobs in
# multiple subdirectories, so make sure they all get added.
jobutils.add_subjob_files_to_backend(rsjob)
else:
for potential in self.potentials.values():
for base_name in potential.base_names:
jaguarworkflows.add_jaguar_files_to_jc_backend(base_name)
[docs] def collectEnergies(self):
"""
Update self.potentials with the Jaguar single point energies.
:raise AnharmonicException: if there is an issue
"""
for potential in self.potentials.values():
energies = []
for base_name in potential.base_names:
jagout_file_name = f'{base_name}.out'
jagrin_file_name = f'{base_name}.01.in'
try:
jagout, jagin = self._getJaguarObjs(jagout_file_name,
jagrin_file_name,
check_freq=False)
except AnharmonicException as err:
msg = (' The Jaguar single point job died.')
raise AnharmonicException(str(err) + msg)
energies.append(getattr(jagout, ENERGY_ATTR))
potential.energies = tuple(energies)
[docs] def getEnergies(self, idx):
"""
Return the energies (Hartree) used to build the potential for the given
normal mode.
:type idx: int
:param idx: the normal mode index, 1-based
:raise AnharmonicException: if there is an issue
:rtype: list
:return: the energies
"""
potential = self.potentials.get(idx)
if not potential:
msg = (
f'No anharmonic potential is available for normal mode {idx}.')
raise AnharmonicException(msg)
energy = getattr(self.jagout, ENERGY_ATTR)
energies = list(potential.energies)
energies.insert(self.factor_data.n_points, energy)
return energies
[docs] def collectFits(self):
"""
Update self.potentials with the anharmonic fit data.
:raise AnharmonicException: if there is an issue
"""
deriv_idx_0 = 0
deriv_idx_2 = 2
def real_evaluator(x, a0, a1, a2, a3, a4):
return evaluate_f(x, deriv_idx_0, (a0, a1, a2, a3, a4))
def imag_evaluator(x, a0, a1, a2, a3, a4):
f_0s = evaluate_f(x, deriv_idx_0, (a0, a1, a2, a3, a4))
f_2s = evaluate_f(x, deriv_idx_2, (a0, a1, a2, a3, a4))
# penalize negative curvature
penalties = []
for f_2 in f_2s:
if f_2 < 0:
with warnings.catch_warnings():
warnings.simplefilter('error', RuntimeWarning)
try:
expf = numpy.exp(-f_2)
except RuntimeWarning:
# MATSCI-9342 - use large number instead of inf
expf = 1e50
penalties.append(expf - 1)
else:
penalties.append(0)
return f_0s + numpy.array(penalties)
factors = list(self.getFactors(self.factor_data))
factors.insert(self.factor_data.n_points, 0)
for idx, potential in self.potentials.items():
energies = self.getEnergies(idx)
if self.jagout.normal_mode[idx - 1].frequency < 0:
evaluator = imag_evaluator
else:
evaluator = real_evaluator
with warnings.catch_warnings():
warnings.simplefilter('ignore', optimize.OptimizeWarning)
try:
popt, pcov = optimize.curve_fit(evaluator, factors,
energies)
except (TypeError, RuntimeError) as err:
msg = (
'The quartic polynomial fit of single point energies '
f'for normal mode {idx} has failed: {str(err)}')
raise AnharmonicException(msg)
potential.fit_coeffs = tuple(popt)
[docs] def evaluate_f(self, idx, factor, deriv_idx, convert_to_si=False):
"""
Evaluate the nth derivative of the anharmonic potential for the given normal
mode index.
:type idx: int
:param idx: the normal mode index, 1-based
:type factor: float
:param factor: the point at which to evaluate
:type deriv_idx: int
:param deriv_idx: indicates what derivative of the polynomial to
evaluate, 0 is the polynomial itself, 1 is the first derivative, etc.
:type convert_to_si: bool
:param convert_to_si: if True convert the returned value from units of
H/Ang.**deriv_idx to J/m**deriv_idx
:raise AnharmonicException: if there is an issue
:rtype: float
:return: the evaluated value in units of H/Ang.**deriv_idx or if
convert_to_si is True in units of J/m**deriv_idx
"""
# in principle the factor is unitless because it multiplies the normal
# mode displacement which is in units of Ang. however despite this
# the coefficients of the potential are still considered to be in
# units of Ang.
potential = self.potentials.get(idx)
if not potential:
msg = (
f'No anharmonic potential is available for normal mode {idx}.')
raise AnharmonicException(msg)
if convert_to_si:
conversion = units.JOULES_PER_HARTREE / constants.angstrom**deriv_idx
else:
conversion = 1
return conversion * evaluate_f(factor, deriv_idx, potential.fit_coeffs)
[docs] def getReducedMass1(self, idx):
"""
Return the reduced mass of the given normal mode using
the Jaguar definition.
:type idx: int
:param idx: the normal mode index, 1-based
:rtype: float
:return: the reduced mass in kg
"""
normal_mode = self.jagout.normal_mode[idx - 1]
vec = normal_mode.displacement
masses = [atom.atomic_weight for atom in self.st.atom]
reduced_mass = 0
for i in range(self.st.atom_total):
reduced_mass += masses[i]**3 * (vec[i, 0]**2 + vec[i, 1]**2 +
vec[i, 2]**2)**2
return KG_PER_AMU * reduced_mass
[docs] def getReducedMass2(self, idx):
"""
Return the reduced mass of the given normal mode using
the definition in the publications followed in this
module.
:type idx: int
:param idx: the normal mode index, 1-based
:rtype: float
:return: the reduced mass in kg
"""
normal_mode = self.jagout.normal_mode[idx - 1]
vec = normal_mode.displacement
reduced_mass = 1 / sum(i**2 for i in vec.flatten())
return KG_PER_AMU * reduced_mass
[docs] def getReducedMass(self, idx):
"""
Return the reduced mass of the given normal mode.
:type idx: int
:param idx: the normal mode index, 1-based
:rtype: float
:return: the reduced mass in kg
"""
# the reduced mass is not an observable, the reduced_mass_1 here
# is equivalent to normal_mode.reduced_mass (see jaguar-src/main/freq.f
# line 628 redmas variable) and is similar to the Mopac definition
# (http://openmopac.net/manual/reduced_masses.html), yet the papers
# followed in this module use reduced_mass_2, both are in units of amu,
# below are some values for a typical example system:
#
# reduced_mass_1 reduced_mass_2
# 0.84 1.92
# 0.16 1.06
# 1.16 2.26
return self.getReducedMass2(idx)
[docs] def collectAnharmonicFrequencies(self):
"""
Update self.potentials with the anharmonic frequencies in wavenumbers
(cm^-1).
:raise AnharmonicException: if there is an issue
"""
deriv_idx_0 = 0
deriv_idx_1 = 1
deriv_idx_2 = 2
x_0 = 0
for idx, potential in self.potentials.items():
# define a guess y_min using the points rather than the fit so
# that potentials can still be plotted prior to job failure
potential.y_min = min(self.getEnergies(idx))
# use Jaguar's reduced mass to be consistent with Jaguar
# frequencies
reduced_mass = self.getReducedMass1(idx) # kg
evaluator = lambda x: (self.evaluate_f(idx, x, deriv_idx_0),
self.evaluate_f(idx, x, deriv_idx_1))
result = optimize.minimize(evaluator, x_0, method='BFGS', jac=True)
if not result.success:
msg = (
'The minimization of the anharmonic potential needed for the '
f'anharmonic frequency of normal mode {idx} has failed: '
f'{result.message}')
raise AnharmonicException(msg)
potential.x_min = result.x[0]
# define the real y_min from the fit, do not rely on result.fun
# attribute because scipy can give an array instead of a float
potential.y_min = evaluator(potential.x_min)[0]
force_constant = self.evaluate_f(
idx, potential.x_min, deriv_idx_2,
convert_to_si=True) # J/m**2 (kg/s**2)
frequency = numpy.sqrt(abs(force_constant) / reduced_mass) # s**-1
frequency = angular_freq_to_freq(frequency) # cm**-1
potential.frequency = numpy.sign(force_constant) * frequency
[docs] def plotPotentials(self):
"""
Plot the potentials.
"""
# get lots of points by reducing the step size
x_step = self.factor_data.step / 100
# get range and extended range
factors = self.getExtendedFactors()
x_e_min, x_e_max = min(factors), max(factors)
factors = list(self.getFactors(self.factor_data))
factors.insert(self.factor_data.n_points, 0)
x_min, x_max = min(factors), max(factors)
deriv_idx = 0
for idx, potential in self.potentials.items():
if not getattr(potential, 'y_min', None):
continue
y_func = lambda x: 1000 * (self.evaluate_f(idx, x, deriv_idx) -
potential.y_min)
file_name = f'{self.base_name}_{idx}_normal_mode_potential.png'
title = f'Potential for normal mode {idx}.'
y_label = 'Potential / mH'
plotter(x_min,
x_max,
x_e_min,
x_e_max,
x_step,
y_func,
file_name,
title,
y_label,
x_values=factors)
[docs] def logCoefficientsTable(self):
"""
Log coefficients table.
"""
sep = 4 * ' '
col_1 = 'Normal mode'
col_2 = f'{sep}V_0/mH'
col_3 = f'{sep}V_1/mH/Ang.'
col_4 = f'{sep}V_2/mH/Ang.^2'
col_5 = f'{sep}V_3/mH/Ang.^3'
col_6 = f'{sep}V_4/mH/Ang.^4'
header_1 = ''.join([col_1, col_2, col_3, col_4, col_5, col_6]) + '\n'
header_2 = '-' * (len(header_1) - 1) + '\n'
self.logger.info(header_1)
self.logger.info(header_2)
for idx in sorted(self.potentials):
potential = self.potentials[idx]
a1, a2, a3, a4 = [
round(1000 * x, 3) for x in potential.fit_coeffs[1:]
]
a0 = round(1000 * (potential.fit_coeffs[0] - potential.y_min), 3)
fields = [f'{idx}'.ljust(len(col_1))]
fields.append(f'{a0}'.rjust(len(col_2)))
fields.append(f'{a1}'.rjust(len(col_3)))
fields.append(f'{a2}'.rjust(len(col_4)))
fields.append(f'{a3}'.rjust(len(col_5)))
fields.append(f'{a4}'.rjust(len(col_6)))
fields = ''.join(fields) + '\n'
self.logger.info(fields)
self.logger.info(header_2)
self.logger.info('')
[docs] def logFrequencyTable(self):
"""
Log frequency table.
"""
sep = 4 * ' '
col_1 = 'Normal mode'
col_2 = f'{sep}v_anharm/cm^-1'
col_3 = f'{sep}v_harm/cm^-1'
col_4 = f'{sep}Delta/cm^-1'
header_1 = ''.join([col_1, col_2, col_3, col_4]) + '\n'
header_2 = '-' * (len(header_1) - 1) + '\n'
self.logger.info(header_1)
self.logger.info(header_2)
for idx in sorted(self.potentials):
va = round(self.potentials[idx].frequency, 3)
vh = round(self.jagout.normal_mode[idx - 1].frequency, 3)
delta = round(va - vh, 3)
fields = [f'{idx}'.ljust(len(col_1))]
fields.append(f'{va}'.rjust(len(col_2)))
fields.append(f'{vh}'.rjust(len(col_3)))
fields.append(f'{delta}'.rjust(len(col_4)))
fields = ''.join(fields) + '\n'
self.logger.info(fields)
self.logger.info(header_2)
self.logger.info('')
[docs] def run(self):
"""
Calculate the anharmonic potentials.
"""
if not all([self.jagout, self.jagin]):
self.runFrequencyJob()
check_imaginary_frequencies(self.jagout,
self.jagin,
max_i_freq=self.max_i_freq)
self.runSinglePointJobs()
self.collectEnergies()
self.collectFits()
try:
self.collectAnharmonicFrequencies()
except AnharmonicException:
raise
finally:
if self.plot:
self.plotPotentials()
# potentially there should be a check on the curvature within
# the fit range here however so far no case has shown this to
# be necessary, curvature outside the fix range is plotted
# but not used in the model
if self.logger and self.potentials:
self.logCoefficientsTable()
self.logFrequencyTable()
[docs]class AnharmonicPartitionFunction(AnharmonicPotentials):
[docs] @staticmethod
def getBeta(temperature):
"""
Return beta.
:type temperature: float
:param temperature: the temperature in K
:rtype: float
:return: beta in 1/J
"""
return 1 / (constants.k * temperature)
[docs] def getClassicalParticleDensity(self, idx, temperature, factor):
"""
For the given normal mode return the classical particle density
evaluated at the given factor.
:type idx: int
:param idx: the normal mode index, 1-based
:type temperature: float
:param temperature: the temperature in K
:type factor: float
:param factor: the point at which to evaluate
:rtype: float
:return: the classical particle density in 1/Ang.
"""
mass = self.getReducedMass(idx) # kg
beta = self.getBeta(temperature) # 1/J
prefactor = numpy.sqrt(2 * constants.pi * mass / beta)
prefactor /= constants.h
prefactor *= constants.angstrom # 1/Ang.
rel_energy = self.evaluate_f(idx, factor, 0, convert_to_si=True) # J
rel_energy -= self.potentials[idx].y_min * units.JOULES_PER_HARTREE # J
return prefactor * numpy.exp(-beta * rel_energy) # 1/Ang.
[docs] def getCorrectionParticleDensity(self, idx, temperature, factor):
"""
For the given normal mode return the particle density
multiplicative correction evaluated at the given factor.
:type idx: int
:param idx: the normal mode index, 1-based
:type temperature: float
:param temperature: the temperature in K
:type factor: float
:param factor: the point at which to evaluate
:rtype: float
:return: the particle density multiplicative correction (unitless)
"""
mass = self.getReducedMass(idx) # kg
beta = self.getBeta(temperature) # 1/J
_lambda = constants.h**2 * beta**2
_lambda /= 8 * mass * constants.pi**2 # m**2/J
evaluator = lambda deriv_idx: self.evaluate_f(
idx, factor, deriv_idx, convert_to_si=True)
t1 = beta * evaluator(1)**2 / 12
t1 += -evaluator(2) / 6 # J/m**2
t2 = beta**2 * evaluator(1)**4 / 288
t2 += -11 * beta * evaluator(1)**2 * evaluator(2) / 360
t2 += evaluator(2)**2 / 40
t2 += evaluator(1) * evaluator(3) / 30
t2 += -evaluator(4) / (60 * beta) # J**2/m**4
correction = 1 + _lambda * t1 + _lambda**2 * t2 # unitless
return correction
[docs] def getParticleDensity(self, idx, temperature, factor):
"""
For the given normal mode return the particle density
evaluated at the given factor.
:type idx: int
:param idx: the normal mode index, 1-based
:type temperature: float
:param temperature: the temperature in K
:type factor: float
:param factor: the point at which to evaluate
:rtype: float
:return: the particle density in 1/Ang.
"""
particle_density = self.getClassicalParticleDensity(
idx, temperature, factor)
particle_density *= self.getCorrectionParticleDensity(
idx, temperature, factor)
return particle_density
[docs] def plotParticleDensity(self, idx, temperature):
"""
For the given normal mode plot the particle density.
:type idx: int
:param idx: the normal mode index, 1-based
:type temperature: float
:param temperature: the temperature in K
"""
# get lots of points by reducing the step size
x_step = self.factor_data.step / 100
# get range and extended range
factors = self.getFactors(self.factor_data)
x_min, x_max = min(factors), max(factors)
factors = self.getExtendedFactors()
x_e_min, x_e_max = min(factors), max(factors)
y_func = lambda x: self.getParticleDensity(idx, temperature, x)
file_name = (f'{self.base_name}_{idx}_normal_mode_{temperature}K'
'_particle_density.png')
title = f'Particle density for normal mode {idx} at {temperature}K.'
y_label = 'Particle Density / (1/Ang.)'
plotter(x_min, x_max, x_e_min, x_e_max, x_step, y_func, file_name,
title, y_label)
[docs] def checkCorrectionParticleDensity(self, idx, temperature):
"""
For the given normal mode check the particle density
multiplicative correction.
:type idx: int
:param idx: the normal mode index, 1-based
:type temperature: float
:param temperature: the temperature in K
:raise AnharmonicException: if there is an issue
"""
# the potential should be a small perturbation of a harmonic so
# the classical particle density should be well behaved but the
# perturbative correction to it should be checked here
# use the newton method for root finding because a bracketing
# interval in which the function changes sign can not be
# guaranteed, also newton does not guarantee finding a root
# so handle that
evaluator = lambda x: self.getCorrectionParticleDensity(
idx, temperature, x)
# checking with a single guess can be insufficient, but start at
# the middle
factors = self.getFactors(self.factor_data)
for factor in factors[self.factor_data.n_points:]:
for sign in [1, -1]:
factor *= sign
try:
root = optimize.newton(evaluator, factor, maxiter=1000)
except RuntimeError:
# no root has been found
continue
if evaluator(root) > CORRECTION_THRESH:
# a value has been found but it is not a root
continue
# a root has been found
msg = ('Detected divergence of the anharmonic correction '
f'for normal mode {idx} at {temperature}K. This may '
'occur at too low of a temperature for higher frequency '
'normal modes.')
raise AnharmonicException(msg)
[docs] def getAnharmonicVibPartitionFunctions(self, temperature):
"""
Return the ln of the anharmonic vibrational partition functions.
:type temperature: float
:param temperature: the temperature in K
:rtype: dict
:return: keys are normal mode indices, 1-based, values are
ln of the anharmonic vibrational partition functions
"""
# by not extending the integration range beyond the fit range,
# contributions from the tails of the particle densities for low
# frequency normal modes will be missed, typically these are small,
# this convention appears to be used in the papers that this module
# is based on
factors = self.getFactors(self.factor_data)
xmin, xmax = min(factors), max(factors)
# for x ranges larger than the range used in the fit (as in extrapolating
# rather than interpolating) the quadrature throws an AccuracyWarning, if the
# x range is near infinite then the error returned from the quadrature is
# significant (probably due to the limited number of iterations), if the x
# range is like that used here then the error is insignificant, so prevent
# the AccuracyWarning by increasing the maxiter from the default of 50 to 1000
lnz_a_vibs = {}
for idx in self.potentials.keys():
if idx in self.divergencies.get(temperature, []):
continue
evaluator = lambda x: self.getParticleDensity(idx, temperature, x)
val, err = integrate.quadrature(evaluator, xmin, xmax, maxiter=1000)
if self.logger:
err = round(err, 3)
self.logger.info('Quadrature error in q_vib_anharm for normal '
f'mode {idx} ({temperature}K): {err}')
lnz_a_vibs[idx] = numpy.log(val)
if self.logger and self.potentials:
self.logger.info('')
return lnz_a_vibs
[docs] def getHarmonicVibPartitionFunctions(self, temperature):
"""
Return the ln of the harmonic vibrational partition functions.
:type temperature: float
:param temperature: the temperature in K
:rtype: dict
:return: keys are normal mode indices, 1-based, values are
ln of the harmonic vibrational partition functions
"""
beta = self.getBeta(temperature) # 1/J
# imaginary frequencies do not contribute to the vibrational
# partition function
# by definition Jaguar does not include the ZPE prefactor in the vibrational
# partition function
# go ahead and get the harmonic partition functions of the normal modes that
# are being treated anharmonically so that a comparison can be logged
lnz_h_vibs = {}
for idx, normal_mode in get_normal_modes(self.jagout,
max_i_freq=self.max_i_freq):
angular_freq = freq_to_angular_freq(abs(
normal_mode.frequency)) # s**-1
val = 1 / (1 - numpy.exp(-beta * constants.hbar * angular_freq))
lnz_h_vibs[idx] = numpy.log(val)
return lnz_h_vibs
[docs] @staticmethod
def getVibPartitionFunction(lnz_a_vibs, lnz_h_vibs):
"""
Return the ln of the vibrational partition function.
:type lnz_a_vibs: dict
:type lnz_a_vibs: keys are normal mode indices, 1-based, values are
ln of the anharmonic vibrational partition functions
:type lnz_h_vibs: dict
:type lnz_h_vibs: keys are normal mode indices, 1-based, values are
ln of the harmonic vibrational partition functions
:rtype: float
:return: the ln of the vibrational partition function
"""
lnz_vib = sum(lnz_a_vibs.values())
lnz_vib += sum(lnz_h_vib for idx, lnz_h_vib in lnz_h_vibs.items()
if idx not in lnz_a_vibs)
return lnz_vib
[docs] def logLnQTable(self, temperature, lnz_a_vibs, lnz_h_vibs):
"""
Log lnQ table.
:type temperature: float
:param temperature: the temperature in K
:type lnz_a_vibs: dict
:type lnz_a_vibs: keys are normal mode indices, 1-based, values are
ln of the anharmonic vibrational partition functions
:type lnz_h_vibs: dict
:type lnz_h_vibs: keys are normal mode indices, 1-based, values are
ln of the harmonic vibrational partition functions
"""
sep = 4 * ' '
col_1 = 'Normal mode'
col_2 = f'{sep}ln(q_vib_anharm)'
col_3 = f'{sep}ln(q_vib_harm)'
col_4 = f'{sep}Delta'
col_5 = f'{sep}Temp./K'
header_1 = ''.join([col_1, col_2, col_3, col_4, col_5]) + '\n'
header_2 = '-' * (len(header_1) - 1) + '\n'
self.logger.info(header_1)
self.logger.info(header_2)
for idx in sorted(lnz_a_vibs):
lnz_a_vib = round(lnz_a_vibs[idx], 3)
lnz_h_vib = round(lnz_h_vibs[idx], 3)
delta = round(lnz_a_vib - lnz_h_vib, 3)
fields = [f'{idx}'.ljust(len(col_1))]
fields.append(f'{lnz_a_vib}'.rjust(len(col_2)))
fields.append(f'{lnz_h_vib}'.rjust(len(col_3)))
fields.append(f'{delta}'.rjust(len(col_4)))
fields.append(f'{temperature}'.rjust(len(col_5)))
fields = ''.join(fields) + '\n'
self.logger.info(fields)
lnz_a_vib = round(self.getVibPartitionFunction(lnz_a_vibs, lnz_h_vibs),
3)
lnz_h_vib = round(sum(lnz_h_vibs.values()), 3)
delta = round(lnz_a_vib - lnz_h_vib, 3)
self.logger.info(header_2)
self.logger.info(f'ln(q_vib_anharm): {lnz_a_vib}')
self.logger.info(f'ln(q_vib_harm): {lnz_h_vib}')
self.logger.info(f'Delta: {delta}')
self.logger.info('')
[docs] def setDivergencies(self):
"""
Set divergencies.
"""
self.divergencies = {}
# this means there are no normal modes with frequencies less than
# the anharmonic threshold so continue as entirely harmonic, this
# information has already been logged so just return
if not self.potentials:
return
# Jaguar thermo objects run over a matrix of unique temperatures and
# pressures, the vibrational partition function does not depend
# on pressure, so for each unique temperature plot the particle density
# for each normal mode that is being treated anharmonically, collect
# normal modes which at the given temperature feature a divergent particle
# density, rather than failing for temperatures that have divergencies
# log them and continue later with the harmonic approximation
temperatures = set()
for thermo in self.jagout.thermo:
if thermo.temp not in temperatures:
idxs = []
for idx in self.potentials.keys():
if self.plot:
self.plotParticleDensity(idx, thermo.temp)
try:
self.checkCorrectionParticleDensity(idx, thermo.temp)
except AnharmonicException:
idxs.append(idx)
if idxs:
self.divergencies[thermo.temp] = idxs
temperatures.add(thermo.temp)
if self.divergencies and self.logger:
self.logger.info(
'Detected divergence of the anharmonic correction '
'for the following normal modes at the given temperatures. This '
'may occur at too low of temperatures for higher frequency normal '
'modes.')
self.logger.info('')
for temperature in sorted(self.divergencies):
idxs = ','.join(map(str,
sorted(self.divergencies[temperature])))
self.logger.info(f'{temperature}K: {idxs}')
self.logger.info('')
bad_ts = [
k for k, v in self.divergencies.items()
if len(v) == len(self.potentials)
]
if bad_ts:
self.logger.info(
'For the following temperatures all normal modes '
'have divergent anharmonic corrections:')
self.logger.info('')
for bad_t in bad_ts:
self.logger.info(f'{bad_t}K')
self.logger.info('')
self.logger.info('Continuing with the harmonic approximation '
'for any divergent normal modes.')
self.logger.info('')
[docs] def run(self):
"""
Calculate the anharmonic partition function.
"""
super().run()
self.setDivergencies()
lnz_a_h_vibs = {}
for thermo in self.jagout.thermo:
if thermo.temp not in lnz_a_h_vibs:
# the vibrational partition function does not depend on
# pressure
lnz_a_vibs = self.getAnharmonicVibPartitionFunctions(
thermo.temp)
lnz_h_vibs = self.getHarmonicVibPartitionFunctions(thermo.temp)
if self.logger and self.potentials:
self.logLnQTable(thermo.temp, lnz_a_vibs, lnz_h_vibs)
# cache the partition functions
lnz_a_vib = self.getVibPartitionFunction(lnz_a_vibs, lnz_h_vibs)
lnz_h_vib = sum(lnz_h_vibs.values())
lnz_a_h_vibs[thermo.temp] = (lnz_a_vib, lnz_h_vib)
lnz_a_vib, lnz_h_vib = lnz_a_h_vibs[thermo.temp]
press = jaguarworkflows.format_pressure(thermo.press)
lnqvib_key = LNQVIB_KEY.format(temp=thermo.temp, press=press)
self.st.property[lnqvib_key] = lnz_a_vib
lnz_total = thermo.lnQ.total + lnz_a_vib - lnz_h_vib
lnq_key = LNQ_KEY.format(temp=thermo.temp, press=press)
self.st.property[lnq_key] = lnz_total
[docs]class AnharmonicThermochemicalProperties(AnharmonicPartitionFunction):
# the methods in this class follow from jaguar-src/main/futil.f
[docs] def getVibrationalTemperature(self, idx):
"""
Return the vibrational temperature of the given normal mode.
:type idx: int
:param idx: the normal mode index, 1-based
:rtype: float
:return: the vibrational temperature in K
"""
# Jaguar uses a default of freqcut = 10 cm^-1 which means that
# by default imaginary frequencies do not contribute to thermochemical
# properties and so their vibrational temperatures are not calculated
if idx in self.potentials:
freq = self.potentials[idx].frequency
else:
freq = self.jagout.normal_mode[idx - 1].frequency
angular_freq = freq_to_angular_freq(abs(freq)) # s**-1
vib_t = angular_freq * constants.h
vib_t /= 2 * constants.pi * constants.k # K
return vib_t
[docs] def getInternalEnergy(self, thermo):
"""
Return the internal energy.
:type thermo: schrodinger.application.jaguar.results.ThermoCollection
:param thermo: the thermo object
:rtype: float
:return: the internal energy in kcal/mol
"""
vib_u = 0
for idx, normal_mode in get_normal_modes(self.jagout,
max_i_freq=self.max_i_freq):
vib_t = self.getVibrationalTemperature(idx)
aexp = numpy.exp(-vib_t / thermo.temp)
vib_u += vib_t * aexp / (1 - aexp)
vib_u *= constants.k * units.KCAL_PER_MOL_PER_JOULE
tot_u = thermo.U.total - thermo.U.vibrational + vib_u
return tot_u
[docs] def getHeatCapacity(self, thermo):
"""
Return the heat capacity.
:type thermo: schrodinger.application.jaguar.results.ThermoCollection
:param thermo: the thermo object
:rtype: float
:return: the heat capacity in cal/(mol * K)
"""
vib_cv = 0
for idx, normal_mode in get_normal_modes(self.jagout,
max_i_freq=self.max_i_freq):
vib_t = self.getVibrationalTemperature(idx)
aexp = numpy.exp(-vib_t / thermo.temp)
vib_cv += (vib_t / thermo.temp)**2 * aexp / (1 - aexp)**2
vib_cv *= constants.k * units.KCAL_PER_MOL_PER_JOULE * 1000
tot_cv = thermo.Cv.total - thermo.Cv.vibrational + vib_cv
return tot_cv
[docs] def getEntropy(self, thermo):
"""
Return the entropy.
:type thermo: schrodinger.application.jaguar.results.ThermoCollection
:param thermo: the thermo object
:rtype: float
:return: the entropy in cal/(mol * K)
"""
vib_s = 0
for idx, normal_mode in get_normal_modes(self.jagout,
max_i_freq=self.max_i_freq):
vib_t = self.getVibrationalTemperature(idx)
aexp = numpy.exp(-vib_t / thermo.temp)
vib_s += vib_t * aexp / (thermo.temp * (1 - aexp))
vib_s -= numpy.log(1 - aexp)
vib_s *= constants.k * units.KCAL_PER_MOL_PER_JOULE * 1000
tot_s = thermo.S.total - thermo.S.vibrational + vib_s
return tot_s
[docs] def getEnthalpy(self, thermo):
"""
Return the enthalpy.
:type thermo: schrodinger.application.jaguar.results.ThermoCollection
:param thermo: the thermo object
:rtype: float
:return: the enthalpy in kcal/mol
"""
tot_u = self.getInternalEnergy(thermo)
vib_u = tot_u - (thermo.U.total - thermo.U.vibrational)
vib_h = vib_u
tot_h = thermo.H.total - thermo.H.vibrational + vib_h
return tot_h
[docs] def getGibbsFreeEnergy(self, thermo):
"""
Return the Gibbs free energy.
:type thermo: schrodinger.application.jaguar.results.ThermoCollection
:param thermo: the thermo object
:rtype: float
:return: the Gibbs free energy in kcal/mol
"""
tot_g = self.getEnthalpy(
thermo) - thermo.temp * self.getEntropy(thermo) / 1000
return tot_g
[docs] def logPropertyTable(self, thermo):
"""
Log property table.
:type thermo: schrodinger.application.jaguar.results.ThermoCollection
:param thermo: the thermo object
"""
sep = 4 * ' '
col_1 = f'Property{sep}{sep}'
col_2 = f'{sep}{sep}Anharm'
col_3 = f'{sep}{sep}Harm'
col_4 = f'{sep}Delta'
col_5 = f'{sep}Temp./K'
col_6 = f'{sep}Press./atm'
header_1 = ''.join([col_1, col_2, col_3, col_4, col_5, col_6]) + '\n'
header_2 = '-' * (len(header_1) - 1) + '\n'
self.logger.info(header_1)
self.logger.info(header_2)
tot_a_u = round(self.getInternalEnergy(thermo), 3)
tot_h_u = round(thermo.U.total, 3)
tot_a_cv = round(self.getHeatCapacity(thermo), 3)
tot_h_cv = round(thermo.Cv.total, 3)
tot_a_s = round(self.getEntropy(thermo), 3)
tot_h_s = round(thermo.S.total, 3)
tot_a_h = round(self.getEnthalpy(thermo), 3)
tot_h_h = round(thermo.H.total, 3)
tot_a_g = round(self.getGibbsFreeEnergy(thermo), 3)
tot_h_g = round(thermo.G.total, 3)
pairs = [(tot_a_u, tot_h_u), (tot_a_cv, tot_h_cv), (tot_a_s, tot_h_s),
(tot_a_h, tot_h_h), (tot_a_g, tot_h_g)]
props = [
'U/kcal/mol', 'Cv/cal/(mol*K)', 'S/cal/(mol*K)', 'H/kcal/mol',
'G/kcal/mol'
]
for prop, (tot_a, tot_h) in zip(props, pairs):
delta = round(tot_a - tot_h, 3)
fields = [f'{prop}'.ljust(len(col_1))]
fields.append(f'{tot_a}'.rjust(len(col_2)))
fields.append(f'{tot_h}'.rjust(len(col_3)))
fields.append(f'{delta}'.rjust(len(col_4)))
fields.append(f'{thermo.temp}'.rjust(len(col_5)))
fields.append(f'{thermo.press}'.rjust(len(col_6)))
fields = ''.join(fields) + '\n'
self.logger.info(fields)
self.logger.info(header_2)
self.logger.info('')
[docs] def run(self):
"""
Calculate the anharmonic thermochemical properties.
"""
super().run()
for thermo in self.jagout.thermo:
if self.logger and self.potentials:
self.logPropertyTable(thermo)
press = jaguarworkflows.format_pressure(thermo.press)
# internal energy, currently harmonic ZPE is used
tot_u = self.getInternalEnergy(thermo)
tot_u += self.jagout.zero_point_energy
tot_u *= units.HARTREE_PER_KCAL_PER_MOL
tot_u += getattr(self.jagout, ENERGY_ATTR)
u_key = U_KEY.format(temp=thermo.temp, press=press)
self.st.property[u_key] = tot_u
# enthalpy
tot_h = tot_u
tot_h += IDEAL_GAS_CONSTANT * units.HARTREE_PER_KCAL_PER_MOL * \
thermo.temp
h_key = H_KEY.format(temp=thermo.temp, press=press)
self.st.property[h_key] = tot_h
# Gibbs free energy
tot_g = tot_h
tot_g -= thermo.temp * self.getEntropy(
thermo) * units.HARTREE_PER_KCAL_PER_MOL / 1000
g_key = G_KEY.format(temp=thermo.temp, press=press)
self.st.property[g_key] = tot_g
mae_out_file = f'{self.base_name}{OUT_EXT}'
self.st.write(mae_out_file)
backend = None
jobutils.add_outfile_to_backend(mae_out_file, backend)