"""
Utilities for bandshape calculations.
Copyright Schrodinger, LLC. All rights reserved.
"""
import glob
import math
import os
from collections import OrderedDict
from collections import namedtuple
import numpy
from scipy import constants
from schrodinger import structure
from schrodinger.application.jaguar.input import JaguarInput
from schrodinger.application.jaguar.output import JaguarOutput
from schrodinger.application.matsci import anharmonic
from schrodinger.application.matsci import \
jaguar_multistage_workflow_utils as jmswfu
from schrodinger.application.matsci import jobutils
from schrodinger.application.matsci import msutils
from schrodinger.application.matsci import spectra
from schrodinger.infra import mm
from schrodinger.infra import table
from schrodinger.structutils import analyze
from schrodinger.structutils import transform
from schrodinger.utils import fileutils
FAILED_SUBJOBS = 'failed_subjobs'
BS_SPECTRUM_FILE_TAG = 's_matsci_Bandshape_Spectrum_File'
BS_SPECTRUM_FILE_KEY = BS_SPECTRUM_FILE_TAG + '({temp}K)'
EXTRA_TAG = '_bandshape'
# only atoms with atomic number >= SOC_HEAVY_ATOM_Z have ZORA
# basis sets available
SOC_HEAVY_ATOM_Z = 18
ABSORPTION = 'absorption'
EMISSION = 'emission'
TRANSITION_KEY = 's_matsci_transition_type'
GAUSSIAN = spectra.GAUSSIAN
LORENTZIAN = spectra.LORENTZIAN
MOLAR_ABSORPTIVITY = 'Molar absorptivity'
NUMBER_OF_EMITTED_PHOTONS = 'Number of emitted photons'
SPECTRUM_TYPE_STARTER = 'SPECTRA_TYPE'
HUANG_RHYS_STARTER = 'Huang-Rhys factor of each mode:'
HUANG_RHYS_ENDER = 'Huang-Rhys factor total:'
TRANSITION_STARTER = 'Band shape peak initial and final vibrational states:'
TRANSITION_ENDER = 'Writing energy (electron volt) vs band shape intensity to file'
TRANSITION_INITIAL_STARTER = 'Initial state:'
TRANSITION_FINAL_STARTER = 'final state:'
TRANSFORM_STARTER = 'Transform from mass-weighted cartesian to normal mode'
# frequency in cm-1, normal_mode in Ang.
NormalMode = namedtuple('NormalMode', ['frequency', 'normal_mode'])
PRINT_PEAK_STATES_LABEL = 'Print information for intense transitions'
OUT_EXT = '_out'
SIMPLIFIED_EXT = '_simplified.out'
GROUND_OPT_RIN = '_ground_opt.01.in'
INITIAL_OPT_RIN = '_initial_opt.01.in'
FINAL_OPT_RIN = '_final_opt.01.in'
FINAL_REOPT_RIN = '_final_reopt.01.in'
# simple namespace object
[docs]class VibState(object):
[docs] def __init__(self, index, n_quantum, huang_rhys_f):
self.index = index
self.n_quantum = n_quantum
self.huang_rhys_f = huang_rhys_f
# simple namespace object, energy in eV
[docs]class Transition(object):
[docs] def __init__(self, index, energy, intensity, initial_vib_states,
final_vib_states, max_n_imaginary):
self.index = index
self.energy = energy
self.intensity = intensity
self.initial_vib_states = initial_vib_states
self.final_vib_states = final_vib_states
self.max_n_imaginary = max_n_imaginary
[docs] @staticmethod
def getSimplified(vib_states):
"""
Return a simplified representation of the transition.
:type vib_states: list
:param vib_states: contains VibState
:rtype: str
:return: the simplified representation of the transition
"""
transition = ''
for vib_state in vib_states:
if vib_state.n_quantum:
transition += f'{vib_state.index}({vib_state.n_quantum})'
return transition or '0'
[docs]def is_bs_spectrum_file_key(key):
"""
Return True if the given structure property key is a
bandshape spectrum file key.
:type key: str
:param key: the structure property key to check
:rtype: bool
:return: True if the given key is a bandshape spectrum file key
"""
return key.startswith(BS_SPECTRUM_FILE_TAG)
[docs]def get_temp_from_bs_spectrum_file_key(key):
"""
Return the temperature from the given bandshape spectrum
file structure property key.
:type key: str
:param key: the structure property key containing the
temperature
:rtype: float
:return: the temperature in K
"""
empty = BS_SPECTRUM_FILE_KEY.format(temp='')
return float(key.strip(empty))
[docs]def get_bs_spectrum_file_temps_keys(struct):
"""
Return the bandshape spectrum file structure property
keys and associated temperatures for the given structure.
:type struct: schrodinger.structure.Structure
:param struct: the structure with the property keys
:rtype: list
:return: pair tuples of temperatures in K and bandshape
spectrum file structure property keys
"""
temps_keys = []
for key in list(struct.property):
if is_bs_spectrum_file_key(key):
temp = get_temp_from_bs_spectrum_file_key(key)
temps_keys.append((temp, key))
return sorted(temps_keys, key=lambda x: x[0])
[docs]def get_transition_type(atable):
"""
Return the transition type of the given table.
:type atable: table.Table
:param atable: the table to check
:rtype: str or None
:return: module constant ABSORPTION or
EMISSION or None if it is unknown
"""
try:
atype = atable[TRANSITION_KEY]
except mm.MmException:
atype = None
return atype
[docs]class BandshapeSpectrumFile(spectra.SpectrumFile):
"""
Manage a bandshape spectrum file.
"""
ENERGY_KEY = 'r_j_Excitation_Energy_(eV)'
ENERGY_TITLE = 'Excitation Energy (eV)'
INTENSITY_KEY = 'r_j_Intensity'
INTENSITY_TITLE = 'Intensity'
SYMMETRY_KEY = 's_j_Symmetry'
SYMMETRY_TITLE = 'Symmetry'
SPECTRUM_KEY = 's_j_spectrum_type'
SPECTRUM_TITLE = 'Electronic Transition w/ Bandshape'
X_KEY = 's_j_x_label'
X_ALIAS = ENERGY_KEY
Y_KEY = 's_j_y_label'
Y_ALIAS = INTENSITY_KEY
HEADERS = OrderedDict([(SPECTRUM_KEY, SPECTRUM_TITLE), (X_KEY, X_ALIAS),
(Y_KEY, Y_ALIAS)])
COLUMNS = OrderedDict([(ENERGY_KEY, ENERGY_TITLE),
(INTENSITY_KEY, INTENSITY_TITLE),
(SYMMETRY_KEY, SYMMETRY_TITLE)])
EXT = '.uvv.spm'
# the following default is 400. cm-1 or 0.05 eV
DEFAULT_BIN_WIDTH = 400.
[docs] def __init__(self, data_file_name, spm_file_name, transition_type=None):
"""
Create an instance.
:type data_file_name: str
:param data_file_name: the text file containing
whitespace separated energies and intensities,
one pair per line, energies in eV
:type spm_file_name: str
:param: spm_file_name: the name of the `*spm` file to create
:type transition_type: str or None
:param transition_type: module constant ABSORPTION or
EMISSION or None if it is unknown
"""
super().__init__(data_file_name, spm_file_name)
self.transition_type = transition_type
def _buildTable(self):
"""
Build the table.
"""
super()._buildTable()
if self.transition_type:
mm.mmtable_set_string(self.table_obj, TRANSITION_KEY,
self.transition_type)
xidx, yidx, zidx = list(
map(self.table_obj.getColumnIndex, list(self.COLUMNS)))
self.table_obj.columnSetVisible(zidx, False)
for idx in range(1, len(self.table_obj) + 1):
self.table_obj[idx][zidx] = ''
[docs] @staticmethod
def isBandshapeTable(atable):
"""
Return True if the given table is a bandshape table.
:type atable: table.Table
:param atable: the table to check
:rtype: bool
:return: True if the given table is a bandshape table
"""
try:
atype = atable[BandshapeSpectrumFile.SPECTRUM_KEY]
except TypeError:
atype = ''
return atype == BandshapeSpectrumFile.SPECTRUM_TITLE
[docs] @staticmethod
def getTrimmedTable(filepath, bin_width=None):
"""
Return a table of trimmed data for the given file path
and bin width.
:type filepath: str
:param filepath: the file path to the file with the data
:type bin_width: float
:param bin_width: the bin width in wavenumbers (cm-1) used
for binning
:raise ValueError: if there isn't any data
:rtype: table.Table
:return: table of trimmed data
"""
atable = table.Table(filepath)
if not atable or not len(atable):
msg = ('File path {apath} has no data.').format(apath=filepath)
raise ValueError(msg)
# only trim bandshape tables
if not BandshapeSpectrumFile.isBandshapeTable(atable):
return atable
# convert from wavenumber (cm-1) to energy (eV)
if bin_width is None:
bin_width = BandshapeSpectrumFile.DEFAULT_BIN_WIDTH
bin_width *= spectra.WAVENUMBER_TO_EV
# bin energies
x_label = atable.columnGetDataName(1)
x_values = numpy.array(
[atable[idx][x_label] for idx in range(1,
len(atable) + 1)]) # in eV
x_bin_start = x_values[0]
x_bin_stop = x_values[-1]
num_bins = int(math.floor((x_bin_stop - x_bin_start) / bin_width))
if num_bins == 0:
num_bins = 1
x_bin_starts = numpy.linspace(x_bin_start,
x_bin_stop,
num=num_bins,
endpoint=False)
x_bin_idxs = numpy.digitize(x_values, x_bin_starts, right=False)
# collect (energy, intensity) points into bins
y_label = atable.columnGetDataName(2)
binned_points = {}
for idx, x_bin_idx in enumerate(x_bin_idxs, 1):
x_value = x_values[idx - 1]
y_value = atable[idx][y_label]
binned_points.setdefault(x_bin_idx, []).append((x_value, y_value))
# build trimmed table by summing intensities and intensity
# averaging the energies, just use the same table instance
# to prevent having to handle various other metadata
for idx in reversed(range(len(binned_points) + 1, len(atable) + 1)):
atable.deleteRow(idx)
for idx, points in enumerate(binned_points.values(), 1):
y_value = sum(list(zip(*points))[1])
x_value = sum([x * y / y_value for x, y in points])
atable[idx][x_label] = x_value
atable[idx][y_label] = y_value
return atable
[docs]def get_translated_xyz_vec(st):
"""
Return a vector of cartesian coordinates of the given structure which
have been translated so that the center of mass is at the (0, 0, 0)
origin.
:type st: `schrodinger.structure.Structure`
:param st: the structure
:rtype: numpy.array
:return: the translated cartesian coordinates
"""
x, y, z = -1 * analyze.center_of_mass(st)
transform.translate_structure(st, x=x, y=y, z=z)
return st.getXYZ().flatten()
[docs]def get_huang_rhys_fs(spectrum_type,
initial_vib_jo,
final_vib_jo,
geoms_base_name=None):
"""
Return the Huang-Rhys factors for the initial and final states.
:type spectrum_type: str
:param spectrum_type: the spectrum type, either module constant ABSORPTION
or EMISSION
:type initial_vib_jo: JaguarOutput
:param initial_vib_jo: the Jaguar output class for the initial state frequency
calculation
:type final_vib_jo: JaguarOutput
:param final_vib_jo: the Jaguar output class for the final state frequency
calculation
:type geoms_base_name: str or None
:param geoms_base_name: if given use geometries from Jaguar restart
.01.in files rather than from the input Jaguar objects
:rtype: list, list
:return: Huang-Rhys factors for initial and final states
"""
if spectrum_type == ABSORPTION:
ground_vib_jo = initial_vib_jo
excited_vib_jo = final_vib_jo
else:
ground_vib_jo = final_vib_jo
excited_vib_jo = initial_vib_jo
# see MATSCI-11715 - for vertical FC calculations the inital and final
# frequency calculations have identical geometries, so instead get
# the difference in geometry from the optimization files which
# differ by a single step of optimization
if geoms_base_name:
ground_opt_rin_fn = f'{geoms_base_name}{GROUND_OPT_RIN}'
if spectrum_type == ABSORPTION:
excited_opt_rin_fn = f'{geoms_base_name}{FINAL_OPT_RIN}'
else:
excited_opt_rin_fn = f'{geoms_base_name}{INITIAL_OPT_RIN}'
# for emission jobs there is the possibility that the user
# chose to reoptimize the ground state starting from the
# optimized excited state geometry, if so the following
# file will exist and should be used instead of the
# original ground state file
final_reopt_rin_fn = f'{geoms_base_name}{FINAL_REOPT_RIN}'
if os.path.exists(final_reopt_rin_fn):
ground_opt_rin_fn = final_reopt_rin_fn
ground_opt_jagobj = JaguarInput(ground_opt_rin_fn)
excited_opt_jagobj = JaguarInput(excited_opt_rin_fn)
else:
ground_opt_jagobj = ground_vib_jo
excited_opt_jagobj = excited_vib_jo
# modeled after huang_rhys_factor_calc in mmlibs/bandshape/bandshape_futil.f
# and https://pubs.rsc.org/en/content/articlehtml/2017/ra/c7ra00417f
xyzs = []
for jagobj in [ground_opt_jagobj, excited_opt_jagobj]:
out_file = f'{jagobj.filename}'
if os.path.exists(out_file) and not geoms_base_name:
out_st = anharmonic.get_st_jaguar_output(out_file,
allow_new_dummies=False)
else:
out_st = jagobj.getStructure()
xyz = get_translated_xyz_vec(out_st)
xyzs.append(xyz)
gxyz, exyz = xyzs
delta = exyz - gxyz
gtrans = get_mwc_to_nm_transform(ground_vib_jo)
etrans = get_mwc_to_nm_transform(excited_vib_jo)
if gtrans is None or etrans is None:
return [], []
# determine the index for the first normal mode with a real frequency
get_n_imaginary = lambda x: len(
[i for i in x.normal_mode if i.frequency < 0])
idx = max(get_n_imaginary(ground_vib_jo), get_n_imaginary(excited_vib_jo))
huang_rhys_ge_fs = []
for jagout, trans in [(ground_vib_jo, gtrans), (excited_vib_jo, etrans)]:
huang_rhys_fs = []
for jdx, normal_mode in enumerate(jagout.normal_mode):
if jdx < idx:
continue
eigenvector = trans[:, jdx]
proj = numpy.dot(eigenvector, delta) * constants.angstrom # m
angular_freq = 2 * constants.pi * constants.c * normal_mode.frequency * 100 # s^-1
mass = normal_mode.reduced_mass * constants.value(
'atomic mass unit-kilogram relationship') # kg
huang_rhys_f = constants.pi * mass * angular_freq * proj**2 / constants.h
huang_rhys_fs.append(huang_rhys_f)
huang_rhys_ge_fs.append(huang_rhys_fs)
huang_rhys_g_fs, huang_rhys_e_fs = huang_rhys_ge_fs
if spectrum_type == ABSORPTION:
return huang_rhys_g_fs, huang_rhys_e_fs
else:
return huang_rhys_e_fs, huang_rhys_g_fs
[docs]def get_transitions(bs_out_file_name, initial_jo=None, final_jo=None):
"""
Return a list of transitions from the given band shape output file.
:type bs_out_file_name: str
:param bs_out_file_name: the band shape output file name
:type initial_jo: JaguarOutput or None
:param initial_jo: the Jaguar output class for the initial state
:type final_jo: JaguarOutput or None
:param final_jo: the Jaguar output class for the final state
:rtype: list
:return: contains Transition
"""
spectrum_type = None
huang_rhys_section = False
huang_rhys_fs = []
transition_section = False
transition_final_section = False
transition_initial_section = False
transition_final_quanta = []
transition_initial_quanta = []
transitions = []
with open(bs_out_file_name, 'r') as afile:
# use generator because these files can be huge
for line in afile:
line = line.strip()
if line.startswith(SPECTRUM_TYPE_STARTER):
spectrum_type = line.split()[1]
continue
if line.startswith(HUANG_RHYS_STARTER):
huang_rhys_section = True
continue
if line.startswith(HUANG_RHYS_ENDER):
huang_rhys_section = False
continue
if huang_rhys_section:
huang_rhys_fs.extend(map(float, line.split()))
continue
if line.startswith(TRANSITION_STARTER):
transition_section = True
continue
if line.startswith(TRANSITION_ENDER):
transition_section = False
continue
if transition_section:
if line.startswith(TRANSITION_FINAL_STARTER):
transition_final_section = True
transition_initial_section = False
continue
if line.startswith(TRANSITION_INITIAL_STARTER):
transition_final_section = False
transition_initial_section = True
continue
if transition_final_section and line:
transition_final_quanta.extend(map(int, line.split()))
continue
if transition_initial_section and line:
transition_initial_quanta.extend(map(int, line.split()))
continue
if line:
index, energy, intensity = line.split()
index = int(index)
energy = float(energy) # eV
intensity = float(intensity)
continue
# an empty line signifies the beginning of the next transition block
# so reset variables and collect data
if not line:
transition_final_section = False
transition_initial_section = False
if transition_final_quanta and transition_initial_quanta:
# final and initial quanta have the same length, the maximum number
# of imaginary frequencies corresponds to whichever state has more
assert len(transition_final_quanta) == len(
transition_initial_quanta)
max_n_imaginary = len(transition_final_quanta) - len(
huang_rhys_fs)
final_vib_states = []
initial_vib_states = []
pairs = zip(transition_final_quanta,
transition_initial_quanta)
for idx, (final_n_quantum,
initial_n_quantum) in enumerate(pairs, 1):
jdx = idx - 1 - max_n_imaginary
if final_n_quantum:
final_vib_state = VibState(
index=idx,
n_quantum=final_n_quantum,
huang_rhys_f=huang_rhys_fs[jdx])
final_vib_states.append(final_vib_state)
if initial_n_quantum:
initial_vib_state = VibState(
index=idx,
n_quantum=initial_n_quantum,
huang_rhys_f=None)
initial_vib_states.append(initial_vib_state)
transitions.append(
Transition(index=index,
energy=energy,
intensity=intensity,
initial_vib_states=initial_vib_states,
final_vib_states=final_vib_states,
max_n_imaginary=max_n_imaginary))
transition_final_quanta = []
transition_initial_quanta = []
continue
if initial_jo and final_jo:
# at this point for each transition the final vib states have Huang-Rhys
# factors defined but not the initial vib states so just define both now
if not spectrum_type:
return []
geoms_base_name = get_bs_base_name(bs_out_file_name)
huang_rhys_i_fs, huang_rhys_f_fs = get_huang_rhys_fs(
spectrum_type,
initial_jo,
final_jo,
geoms_base_name=geoms_base_name)
if not numpy.any(huang_rhys_i_fs) or not numpy.any(huang_rhys_f_fs):
return []
for transition in transitions:
pairs = [(transition.initial_vib_states, huang_rhys_i_fs),
(transition.final_vib_states, huang_rhys_f_fs)]
for states, huang_rhys_fs in pairs:
for state in states:
idx = state.index - 1 - transition.max_n_imaginary
state.huang_rhys_f = huang_rhys_fs[idx]
return transitions
[docs]def get_bs_base_name(name):
"""
Return the base name of the given bandshape name.
:type name: str
:param name: the bandshape name
:rtype: str
:return: the base name
"""
assert EXTRA_TAG in name
return name.rsplit(EXTRA_TAG, maxsplit=1)[0]
[docs]def get_bs_ext(name):
"""
Return the extension of the given bandshape name.
:type name: str
:param name: the bandshape name
:rtype: str
:return: the extension
"""
assert EXTRA_TAG in name
return EXTRA_TAG + name.rsplit(EXTRA_TAG, maxsplit=1)[1]
[docs]def get_out_file_paths(source_path, patterns):
"""
Return the output file paths.
:type source_path: str
:param source_path: the source path to the output files
:type patterns: list
:param patterns: the glob patterns of the output files
:raise NormalModeAnalysisException: if there is an issue
:rtype: list
:return: contains the output file paths in the same order
as patterns
"""
file_paths = []
for pattern in patterns:
all_file_paths = glob.glob(os.path.join(source_path, pattern))
if len(all_file_paths) == 1:
file_paths.append(all_file_paths[0])
else:
msg = (f'A single output file {pattern} can not be found.')
raise NormalModeAnalysisException(msg)
return file_paths
[docs]class NormalModeAnalysisException(Exception):
pass
[docs]class NormalModeAnalysis(object):
"""
Manage normal mode analysis.
"""
TAG = 'normal_mode_analysis'
N_ATOMS_KEY = 'i_j_atom_total'
FREQ_KEY = 'r_j_Frequency'
FREQ_TITLE = 'Frequency'
SYMM_KEY = 's_j_Symmetry'
SYMM_TITLE = 'Symmetry'
INTENSITY_KEY = 'r_j_Intensity'
INTENSITY_TITLE = 'Intensity'
RAMAN_ACT_KEY = 'r_j_Raman_Act'
RAMAN_ACT_TITLE = 'Raman Act'
RAMAN_INT_KEY = 'r_j_Raman_Int'
RAMAN_INT_TITLE = 'Raman Int'
X_KEY = 'r_j_x%s'
X_TITLE = 'x%s'
Y_KEY = 'r_j_y%s'
Y_TITLE = 'y%s'
Z_KEY = 'r_j_z%s'
Z_TITLE = 'z%s'
# fixed columns
COLUMNS = OrderedDict([
(FREQ_KEY, FREQ_TITLE),
(SYMM_KEY, SYMM_TITLE),
(INTENSITY_KEY, INTENSITY_TITLE),
(RAMAN_ACT_KEY, RAMAN_ACT_TITLE),
(RAMAN_INT_KEY, RAMAN_INT_TITLE)
]) # yapf: disable
[docs] def __init__(self, bs_out_file_name):
"""
Create an instance.
:type bs_out_file_name: str
:param bs_out_file_name: the band shape output file name
:raise NormalModeAnalysisException: if there is an issue
"""
if not os.path.exists(bs_out_file_name):
msg = (f'Output file {bs_out_file_name} can not be found.')
raise NormalModeAnalysisException(msg)
source_path, file_name = os.path.split(bs_out_file_name)
self.base_name = get_bs_base_name(file_name)
initial_out_file_name, final_out_file_name = get_out_file_paths(
source_path, ['*_initial_vib.out', '*_final_vib.out'])
# the normal modes in the acutal *vib files have more precision
# than those logged to the *out files but use the latter for now
self.initial_jo = JaguarOutput(initial_out_file_name)
if self.initial_jo.status != JaguarOutput.OK or self.initial_jo.fatal_error:
msg = (
f'Output file {initial_out_file_name} comes from a failed job.')
raise NormalModeAnalysisException(msg)
self.final_jo = JaguarOutput(final_out_file_name)
if self.final_jo.status != JaguarOutput.OK or self.final_jo.fatal_error:
msg = (
f'Output file {final_out_file_name} comes from a failed job.')
raise NormalModeAnalysisException(msg)
# for now store transitions in memory, there can be very many
# transitions and very many contributing vibrational states
# per transition so re-reading from disk may be necessary in
# the future
self.transitions = get_transitions(bs_out_file_name,
initial_jo=self.initial_jo,
final_jo=self.final_jo)
[docs] def getVibStates(self, energy_start, energy_stop):
"""
Get all vibrational states within the given energy window
for the initial and final states.
:type energy_start: float
:param energy_start: the lower bound on the energy in eV
:type energy_stop: float
:param energy_stop: the upper bound on the energy in eV
:rtype: dict, dict
:return: initial and final vibrational states, keys are indices,
values are sets of VibState
"""
assert energy_start < energy_stop
istates, fstates = {}, {}
for transition in self.transitions:
if energy_start <= transition.energy <= energy_stop:
for state in transition.initial_vib_states:
istates.setdefault(state.index, set()).add(state)
for state in transition.final_vib_states:
fstates.setdefault(state.index, set()).add(state)
return istates, fstates
[docs] def getVibStateWeight(self, state):
"""
Return the weight of the given vibrational state.
:type state: VibState
:param state: the vibrational state
:rtype: float
:return: the weight
"""
# FIXME - unclear if these should come directly from intensity data
# (see MATSCI-7053) or indirectly from initial and final state
# Huang-Rhys factors
# discussion in https://pubs.rsc.org/en/content/articlehtml/2017/ra/c7ra00417f
# suggests that Huang-Rhys factors could approximate the weights sought here
return state.huang_rhys_f
[docs] def getAllVibData(self, all_states, jagout, include_n_quantum):
"""
Return all frequencies and normal modes from the given vibrational
states.
:type all_states: dict
:param all_states: vibrational states, keys are indices,
values are sets of VibState
:type jagout: JaguarOutput
:param jagout: the Jaguar output class
:type include_n_quantum: bool
:param include_n_quantum: whether to include effects due to
vibrational quantum numbers
:rtype: list
:return: contains NormalMode
"""
data = []
for idx, states in all_states.items():
if include_n_quantum:
weight = numpy.mean([state.n_quantum for state in states])
else:
weight = 1
obj = jagout.normal_mode[idx - 1]
mode = weight * obj.displacement
data.append(NormalMode(frequency=obj.frequency, normal_mode=mode))
return sorted(data, key=lambda x: x.frequency)
[docs] def getAvgVibData(self, all_states, jagout, include_n_quantum):
"""
Return an averaged frequency and normal mode from the given vibrational
states.
:type all_states: dict
:param all_states: vibrational states, keys are indices,
values are sets of VibState
:type jagout: JaguarOutput
:param jagout: the Jaguar output class
:type include_n_quantum: bool
:param include_n_quantum: whether to include effects due to
vibrational quantum numbers
:rtype: NormalMode
:return: the averaged NormalMode
"""
weights = []
avg_freq = 0
avg_mode = numpy.zeros((len(jagout.atom), 3))
for idx, states in all_states.items():
if include_n_quantum:
weight = numpy.mean([state.n_quantum for state in states])
else:
weight = 1
# all states have the same weight
state = list(states)[0]
weight *= self.getVibStateWeight(state)
weights.append(weight)
obj = jagout.normal_mode[idx - 1]
avg_freq += weight * obj.frequency
avg_mode += weight * obj.displacement
avg_freq /= numpy.sum(weights)
avg_mode /= numpy.sum(weights)
return NormalMode(frequency=avg_freq, normal_mode=avg_mode)
[docs] def getMaxVibData(self, all_states, jagout, include_n_quantum):
"""
Return the frequency and normal mode with the largest weight
from the given vibrational states.
:type all_states: dict
:param all_states: vibrational states, keys are indices,
values are sets of VibState
:type jagout: JaguarOutput
:param jagout: the Jaguar output class
:type include_n_quantum: bool
:param include_n_quantum: whether to include effects due to
vibrational quantum numbers
:rtype: NormalMode
:return: the max (largest weight) NormalMode
"""
pairs = []
for idx, states in all_states.items():
if include_n_quantum:
weight = numpy.mean([state.n_quantum for state in states])
else:
weight = 1
# all states have the same weight
state = list(states)[0]
pairs.append((idx, weight * self.getVibStateWeight(state)))
max_idx = max(pairs, key=lambda x: x[1])[0]
obj = jagout.normal_mode[max_idx - 1]
mode = weight * obj.displacement
return NormalMode(frequency=obj.frequency, normal_mode=mode)
[docs] def writeVibFile(self, file_name, n_atoms, normal_modes):
"""
Write the `*.vib` file for the given normal modes.
:type file_name: str
:param file_name: the file name
:type n_atoms: int
:param n_atoms: the number of atoms
:type normal_modes: list
:param normal_modes: contains NormalMode
"""
# emulate the *vib files obtained by running Jaguar
mm.mmtable_initialize(mm.error_handler)
mm.mmerr_level(mm.error_handler, mm.MMERR_OFF)
handle = mm.mmtable_new()
table_obj = table.Table(table_handle=handle)
mm.mmtable_set_integer(table_obj, self.N_ATOMS_KEY, n_atoms)
columns = self.COLUMNS.copy()
for idx in range(1, n_atoms + 1):
columns[self.X_KEY % idx] = self.X_TITLE % idx
columns[self.Y_KEY % idx] = self.Y_TITLE % idx
columns[self.Z_KEY % idx] = self.Z_TITLE % idx
for acol, aname in columns.items():
table_obj.insertColumn(mm.MMTABLE_END, acol)
idx = table_obj.getColumnIndex(acol)
table_obj.columnSetName(idx, aname)
# nothing should be visible except frequency
idxs = list(map(table_obj.getColumnIndex, list(columns)))
for idx in idxs[1:]:
table_obj.columnSetVisible(idx, False)
# tables are 1-indexed
for idx, normal_mode in enumerate(normal_modes, 1):
table_obj.insertRow(mm.MMTABLE_END, 0)
table_obj[idx][idxs[0]] = normal_mode.frequency
table_obj[idx][idxs[1]] = '' # symmetry
table_obj[idx][idxs[2]] = 0 # intensity
table_obj[idx][idxs[3]] = 0 # raman activity
table_obj[idx][idxs[4]] = 0 # raman intensity
for jdx, value in zip(idxs[5:], normal_mode.normal_mode.flatten()):
table_obj[idx][jdx] = value
# avoid catenation of sucessive runs
if os.path.exists(file_name):
fileutils.force_remove(file_name)
mm.mmtable_m2io_write(file_name, table_obj, True)
mm.mmerr_level(mm.error_handler, mm.MMERR_WARNING)
mm.mmtable_terminate()
[docs] def writeFiles(self,
base_name,
initial_normal_modes,
final_normal_modes,
title=None):
"""
Write normal mode analysis files.
:type base_name: str
:param base_name: a base name for the files
:type initial_normal_modes: list or None
:param initial_normal_modes: contains NormalMode for the initial
state
:type final_normal_modes: list or None
:param final_normal_modes: contains NormalMode for the final
state
:type title: str
:param title: a title to use for structures in the normal mode analysis
output files
"""
mae_file = base_name + OUT_EXT + '.mae'
i_vib_file = base_name + '_initial.vib'
f_vib_file = base_name + '_final.vib'
if title:
title = ' ' + title
else:
title = ''
hierarchy = [self.TAG]
ist = self.initial_jo.getStructure()
ist.title = 'initial' + title
msutils.set_project_group_hierarchy(ist, hierarchy)
fst = self.final_jo.getStructure()
fst.title = 'final' + title
msutils.set_project_group_hierarchy(fst, hierarchy)
assert ist.atom_total == fst.atom_total
sts = []
smap_dict = {}
if initial_normal_modes:
sts.append(ist)
smap_dict[i_vib_file] = len(sts)
self.writeVibFile(i_vib_file, ist.atom_total, initial_normal_modes)
if final_normal_modes:
sts.append(fst)
smap_dict[f_vib_file] = len(sts)
self.writeVibFile(f_vib_file, fst.atom_total, final_normal_modes)
structure.write_cts(sts, mae_file)
jmswfu.create_smap(base_name + OUT_EXT, mae_file, smap_dict=smap_dict)
[docs] def run(self,
energy_start,
energy_stop,
average=False,
maximum=False,
include_n_quantum=False,
base_name=None,
title=None):
"""
Run it.
:type energy_start: float
:param energy_start: the lower bound on the energy in eV
:type energy_stop: float
:param energy_stop: the upper bound on the energy in eV
:type average: bool
:param average: whether to average the normal modes
:type maximum: bool
:param maximum: whether to report only the normal mode with the
largest weight
:type include_n_quantum: bool
:param include_n_quantum: whether to include effects due to
vibrational quantum numbers
:type base_name: str
:param base_name: a base name to use for the normal mode analysis
output files
:type title: str
:param title: a title to use for structures in the normal mode analysis
output files
:raise NormalModeAnalysisException: if there is an issue
:rtype: str
:return: the name of the Maestro file written
"""
assert not (average and maximum)
if title:
name = title
else:
name = self.base_name
if not self.transitions:
msg = (f'For {name} there is no transition data.')
raise NormalModeAnalysisException(msg)
initial_vib_states, final_vib_states = self.getVibStates(
energy_start, energy_stop)
if not initial_vib_states and not final_vib_states:
msg = (f'For {name} no initial or final vibrational modes '
'participate to the intensity in this energy '
'window. This means that it is a pure |0> --> |0> '
'transition in this region.')
raise NormalModeAnalysisException(msg)
if average:
getter = lambda *args: [self.getAvgVibData(*args)]
elif maximum:
getter = lambda *args: [self.getMaxVibData(*args)]
else:
getter = self.getAllVibData
initial_normal_modes = final_normal_modes = None
if initial_vib_states:
initial_normal_modes = getter(initial_vib_states, self.initial_jo,
include_n_quantum)
if final_vib_states:
final_normal_modes = getter(final_vib_states, self.final_jo,
include_n_quantum)
base_name = base_name or self.base_name
base_name += '_' + self.TAG
self.writeFiles(base_name,
initial_normal_modes,
final_normal_modes,
title=title)
return base_name + OUT_EXT + '.mae'
[docs]def write_simplified_transitions_file(bs_out_file_name,
s_bs_out_file_name=None):
"""
Write a simplified transitions file.
:type bs_out_file_name: str
:param bs_out_file_name: the band shape output file name
:type s_bs_out_file_name: str or None
:param s_bs_out_file_name: the simplified band shape output file name,
if None then it will be determined from bs_out_file_name
"""
obj = NormalModeAnalysis(bs_out_file_name)
if not obj.transitions:
return
if not s_bs_out_file_name:
s_bs_out_file_name = bs_out_file_name.replace('.out', SIMPLIFIED_EXT)
with open(s_bs_out_file_name, 'w') as afile:
for transition in obj.transitions:
afile.write(f'{transition.index} {transition.energy} '
f'{transition.intensity}\n')
afile.write('Final State:\n')
state = transition.getSimplified(transition.final_vib_states)
afile.write(f'{state}\n')
afile.write('Initial State:\n')
state = transition.getSimplified(transition.initial_vib_states)
afile.write(f'{state}\n')
afile.write('\n')
jobutils.add_outfile_to_backend(s_bs_out_file_name)