"""
Post-processing of ORCA EPR outputs: g-value extraction, tabulation,
spectral comparison, and bond dissociation enthalpy.
"""
from PlottingTool import StandardPlotter
import matplotlib.pyplot as plt
from glob import glob
from eproducer.miscellaneous import output_terminal
import numpy as np
from eproducer.best_fit import TheoMatchExpe
import pandas as pd
from scipy.stats import pearsonr
from scipy.optimize import minimize
from eproducer.g_valsetup import ext_xyz_from_npz
# add2executable
[docs]
def several_gvals(names: str, experiment: np.ndarray, path: str = './'):
"""
Extract several gvals and compare them with experimental values.
Parameters
==========
names: str
names of the orca output files separated by commas.
Eg: 'cc-pVDZ.out,epr_info.out'
experiment: numpy.array
expected values according to the experiment [gx, gy, gz] in
decreasing order.
path: str. Default='./'
path to the directory containing the output files from orca.
Return
======
(list) one array per orca file with [gx, gy, gz, max absolute error,
mean absolute error] as strings.
"""
methods = names.split(',')
res = []
for method in methods:
values = extract_gvals(path + '/' + method)
errors = np.abs(experiment - values)
maxerror = np.max(errors)
avg = np.mean(errors)
values = np.append(values.astype('<U100'),
["{:.7f}".format(round(maxerror, 7)),
"{:.7f}".format(round(avg, 7))])
res.append(values)
return res
# add2executable
[docs]
def create_table_per_method(data: np.ndarray,
experiment: list,
methods: list) -> str:
"""
For a given system, creates a table with the different methods used to
compute the g-values.
Parameters
==========
data: numpy.array
g-values and errors in the shape [meth][gx, gy, gz, max absolute
error, mean absolute error], as returned by several_gvals.
experiment: numpy.array
expected values obtained from the exteriment. [gx, gy, gz] in
decreassing order.
methods: list
list of names of the methods used. they has to be in the same order
than data.
Return
======
(str) Table in LaTeX format of g-values with the different methods and
the associated max error and average error. The table is also printed.
"""
assert len(data) == len(methods), "There is not the same amount of " + \
"methods than g-values."
# object dtype, otherwise np.insert truncates the labels to the width of
# the g-value strings
res = np.array(data, dtype=object).T
header = np.array(["", "\\text{experiments}"] +
[f"\\text{{{meth}}}" for meth in methods])
rows_labels = ["g_x", "g_y", "g_z", "\\text{max absolute error}",
"\\text{Mean absolute error}"]
table = np.insert(res, 0, np.append(experiment.astype(str), ["", ""]),
axis=1)
table = np.insert(table, 0, rows_labels, axis=1)
table = np.insert(table, 0, header, axis=0)
# Calculate column widths
column_widths = [max(len(str(item)) for item in col) for col in zip(*table)]
# Print the data rows
table_str = []
for row in table:
row_data = [f"{col:{width}}" for col, width in zip(row, column_widths)]
line = " & ".join(row_data) + " \\\\ \\hline"
print(line)
table_str.append(line)
return "\n".join(table_str)
# add2executable
[docs]
def create_table_per_system(experiment: list,
systems: list,
orca_file: str='epr_info.out') -> str:
"""
Creates a markdown table with the g-values of several systems computed
with the same method, and their errors with respect to the experiment.
Parameters
==========
experiment: list
expected values obtained from the exteriment. [gx, gy, gz] in
decreassing order.
systems: list
paths to the directories of the systems, ending with '/'. Each
directory contains orca_file, and the name of the system in the table
links to <system>vmd_image.md.
orca_file: str. Default='epr_info.out'
name of the orca file containing the g-values.
Return
======
(str) Table of g-values with the different systems and the associated
max error and average error.
"""
# extract data
data = []
for subsys in systems:
gvals = several_gvals(names=orca_file,
experiment=experiment,
path=subsys)
data.append(gvals[0])
data = np.array(data)
# names
subsys_names = systems.copy()
for i in range(len(systems)):
subsys_names[i] = f"[{subsys_names[i].split('/')[-2]}]({systems[i]}vmd_image.md)"
# header
header = np.array(["System", "$g_x$", "$g_y$", "$g_z$", "Max absolute error",
"Mean absolute error"])
# experiment
experiment = np.append(np.array(experiment).astype('<U100'), ["", ""])
experiment = np.insert(experiment, 0, "Experiment")
# arrange resuls
res = np.insert(np.array(data), 0, subsys_names, axis=1)
# set table
table = np.insert(res, 0, experiment, axis=0)
table = np.insert(table, 0, header, axis=0)
# Calculate column widths
column_widths = [max(len(str(item)) for item in col) for col in zip(*table)]
# Print the data rows
table_str = []
for row in table:
row_data = [f"{col:{width}}" for col, width in zip(row, column_widths)]
table_str.append("| " + " | ".join(row_data) + " |")
header_widths = table_str[0].split('|')[1:-1]
header_sep = "".join(["|" + "-" * len(width) for width in header_widths ]) + "|"
table_str.insert(1, header_sep)
return "\n".join(table_str)
# add2executable
# add2executable
[docs]
def add_gval2npz(path):
"""
Add g-values to the npz files in the entry 'g-values', based on the files
with the shape <name>_<index>_epr.out next to <name>.npz. The g-values of
each file are stored in the row <index> of 'g-values', which is created
with zeros (one row per radical in 'heavy_atom_missing_idxs') if it does
not exist.
Parameters
==========
path: str
path to the directory containing the npz files and the epr.out files,
the npz file itself or a .dat file with a list of npz files.
"""
if path[-4:] == '.npz':
npz_files = [path]
elif path[-4:] == '.dat':
npz_files = np.loadtxt(path, dtype=str)
else:
npz_files = glob(path + '/*.npz')
npz_files.sort()
for file in npz_files:
infofil = np.load(file)
info = {key: infofil[key] for key in infofil.files}
infofil.close()
if not 'g-values' in info.keys():
info['g-values'] = np.zeros((len(info['heavy_atom_missing_idxs']), 3))
eprs = glob(file[:-4] + '_*_epr.out')
if len(eprs) == 0:
print(f'{file} <--- {file} no epr.', flush=True)
continue
for epr in eprs:
print('enter', flush=True)
try:
idx = int(epr.split('_')[-2])
info['g-values'][idx] = extract_gvals(epr)
print(f'{epr} <--- {file} worked.', flush=True)
except:
print(f'{epr} <--- {file} failed.', flush=True)
continue
np.savez(file, **info)
[docs]
class GvalComparison:
"""
Table (pandas.DataFrame stored in GvalComparison.df) to compare the
g-values and spectra of several candidate molecules with each other and
with experiments. Each method adds one or more columns to the table.
Parameters
==========
molecules: list
paths to the directories of the candidate molecules. They are stored
in the column 'molecule'.
"""
def __init__(self, molecules):
"""
Creates the table with the column 'molecule'. See the class docstring
for the parameters.
"""
self.tme = None
self.df = pd.DataFrame({'molecule': molecules})
[docs]
def add_gvals(self, output_files):
"""
Extracts the g-values of each molecule from its orca output and stores
them in the column 'gval'.
Parameters
==========
output_files: str or list
orca output files, one per molecule. If a str is given, it is taken
as the name of the file inside the directory of each molecule, i.e.
<molecule>/<output_files>.
Return
======
(list) g-values [gx, gy, gz] of each molecule.
"""
if isinstance(output_files, str):
output_files = [i + '/' + output_files
for i in self.df['molecule']]
gvals = [extract_gvals(out_file) for out_file in output_files]
self.df['gval'] = gvals
return gvals
[docs]
def d_target_gval(self, gval_ref, column_name):
"""
Computes the euclidean distance between the g-values of each molecule
and a reference, and stores it in a new column. Requires the column
'gval' (see GvalComparison.add_gvals).
Parameters
==========
gval_ref: list
reference g-values [gx, gy, gz], e.g. from the experiment.
column_name: str
name of the new column.
Return
======
(pandas.DataFrame) updated table.
"""
delta_fuction = lambda x: np.linalg.norm(np.array(x) - gval_ref)
self.df[column_name] = self.df['gval'].apply(delta_fuction)
return self.df
[docs]
def add_spectra(self, spect_files, experiment=None):
"""
Loads the computed spectrum of each molecule, interpolated onto the field
of the experiment (see TheoMatchExpe), and stores it in the column
'spectrum'. It also adds the column 'peak_loc'.
Parameters
==========
spect_files: str or list
.dat files (field vs intensity) of the computed spectra, one per
molecule. If a str is given, it is taken as the name of the file
inside the directory of each molecule, i.e. <molecule>/<spect_files>.
experiment: str. Default=None
.dat file of the experimental spectrum whose field is used for the
interpolation. If None, the field of the first spectrum is used.
Return
======
(numpy.array) [n_molecules x n_points] interpolated spectra.
"""
if isinstance(spect_files, str):
spect_files = [i + '/' + spect_files
for i in self.df['molecule']]
self.tme = TheoMatchExpe(files=spect_files, experiment_file=experiment)
self.df['spectrum'] = [np.array(row) for row in self.tme.intensities]
self.add_peak_loc()
return self.tme.intensities
[docs]
def add_peak_loc(self):
"""
Stores the field where each spectrum has its maximum in the column
'peak_loc'. Requires the column 'spectrum' (see
GvalComparison.add_spectra).
Return
======
(list) field of the maximum of each spectrum.
"""
ds = []
for spec in self.df['spectrum']:
ipeak, _ = self.peak_loc(self.tme.fieldexp, spec)
ds.append(ipeak)
self.df['peak_loc'] = ds
return ds
[docs]
def fit_to(self, experiment_file, **kwargs):
"""
Fits the spectra of the molecules to an experimental spectrum, keeping
only the molecules that contribute the most (see
TheoMatchExpe.gradual_cleaning). Requires GvalComparison.add_spectra to
be called first.
Parameters
==========
experiment_file: str
.dat file of the experimental spectrum.
**kwargs:
passed to TheoMatchExpe.gradual_cleaning, e.g. threshold, steps or
shifting_width.
Return
======
(tuple) files of the selected spectra and their percentages in the fit.
"""
self.tme.fieldexp, self.tme.intensexp = np.loadtxt(experiment_file,
unpack=True)
self.tme.intensexp = self.tme.intensexp / max(self.tme.intensexp)
self.tme.gradual_cleaning(**kwargs)
return self.tme.files, self.tme.percentages
[docs]
def lc_max_corr(self, experiment, column_name):
"""
Finds the coefficients of the linear combination of the spectra that
maximizes the Pearson correlation with an experimental spectrum, and
stores them in a new column. Requires the column 'spectrum' (see
GvalComparison.add_spectra).
Parameters
==========
experiment: numpy.array
experimental intensities, defined on the same field values as the
spectra in the column 'spectrum'.
column_name: str
name of the new column with the coefficients.
Return
======
(tuple) coefficients and the maximum correlation reached.
Note
====
The optimization starts from random coefficients, so the result can
change between runs.
"""
def lc_corr(coefs, functions, experiment):
"""
Negative Pearson correlation between the linear combination of functions
(weighted by coefs) and the experiment. It is the function minimized.
"""
linear_combination = np.sum(coefs * functions, axis=0)
return -pearsonr(linear_combination, experiment)[0]
x0 = np.random.rand(len(self.df['spectrum']))
result = minimize(lc_corr, x0, args=(self.df['spectrum'].to_numpy(),
experiment))
self.df[column_name] = result.x
return result.x, -result.fun
[docs]
def individual_corr(self, experiment, column_name):
"""
Computes the Pearson correlation between each spectrum and an
experimental spectrum, and stores it in a new column. Requires the column
'spectrum' (see GvalComparison.add_spectra).
Parameters
==========
experiment: numpy.array
[2 x n_points] experimental field and intensities. The spectra are
interpolated onto this field before computing the correlation.
column_name: str
name of the new column.
Return
======
(list) correlation of each spectrum with the experiment.
"""
corr = []
for spec in self.df['spectrum']:
intensity = np.interp(experiment[0], self.tme.fieldexp, spec)
corr.append(pearsonr(intensity, experiment[1])[0])
self.df[column_name] = corr
return corr
[docs]
def peak_loc(self, x, y):
"""
Finds the maximum of a spectrum.
Parameters
==========
x: numpy.array
field values.
y: numpy.array
intensities.
Return
======
(tuple) field where the maximum is located and the maximum intensity.
"""
ymax = np.max(y)
i = np.where(y == ymax)[0][0]
return x[i], ymax
[docs]
def pp_dist(self, reference, column_name):
"""
Computes the distance between the peak of each spectrum and a reference
peak, and stores it in a new column. Requires the column 'peak_loc' (see
GvalComparison.add_spectra).
Parameters
==========
reference: float, numpy.array or str
field of the reference peak, or a [2 x n_points] spectrum (field and
intensities) whose maximum is taken as reference, or the path to a
.dat file with that spectrum.
column_name: str
name of the new column.
Return
======
(pandas.Series) peak_loc - reference for each molecule.
"""
if isinstance(reference, str):
reference = np.loadtxt(reference, unpack=True)
if len(np.array(reference).shape) != 0:
reference, _ = self.peak_loc(reference[0], reference[1])
self.df[column_name] = self.df['peak_loc'].apply(lambda x: x - reference)
return self.df[column_name]
def _extract_enthalpy(file):
"""
Extract the numeric "Total Enthalpy" value from an orca output file.
Parameters
----------
file : str
Path to an orca output file to search for a line containing
'Total Enthalpy'.
Returns
-------
(float) The enthalpy value parsed from the fourth whitespace-separated
field on the matching line.
"""
e = output_terminal(f"grep 'Total Enthalpy' {file}" +
" | awk '{print $4}' ", print_output=False)
return float(e)
# add2executable
[docs]
def compute_bde(reactants=None, products=None):
"""
Compute the bond dissociation enthalpy (BDE) for a reaction as the
difference between the total enthalpy of products (P) and reactants (R),
such that BDE = H_P - H_R.
Parameters
----------
reactants : iterable
An iterable (e.g., list or tuple) of reactant species. Each element
should be an orca output file path with computed enthalpy.
products : iterable
An iterable (e.g., list or tuple) of product species. Each element
should be an orca output file path with computed enthalpy.
Returns
-------
(float) The BDE computed. The units are the same as in the orca outputs.
"""
if reactants is None:
raise ValueError("Non reactant recognized. See the documentation of" +
"this function (compute_bde).")
if products is None:
raise ValueError("Non products recognized. See the documentation of" +
"this function (compute_bde).")
e_react = 0
for react in reactants:
e_react += _extract_enthalpy(react)
e_produ = 0
for produ in products:
e_produ += _extract_enthalpy(produ)
return e_produ - e_react