From d8347091b1d19b8cedcb9437e2ca9ff783227e34 Mon Sep 17 00:00:00 2001 From: Matt Ernst Date: Wed, 14 Jun 2023 10:58:08 -0700 Subject: [PATCH 1/4] update from https://nwchemgit.github.io/Special_AWCforum/sp/id5655.html --- contrib/parsers/nw_spectrum.py | 1067 +++++++++++++++----------------- 1 file changed, 512 insertions(+), 555 deletions(-) diff --git a/contrib/parsers/nw_spectrum.py b/contrib/parsers/nw_spectrum.py index 9f0934ae03..34f0564f5c 100644 --- a/contrib/parsers/nw_spectrum.py +++ b/contrib/parsers/nw_spectrum.py @@ -1,586 +1,579 @@ +#!/usr/bin/env python ## ## nw_spectrum ## ## Kenneth Lopata -## Last modified: 2013-03-08 +## Last modified: 2015-02-21 ## ## Python script for parsing NWChem output for TDDFT/vspec excitation ## energies, and optionally Lorentzian broadenening the spectra. For ## online help run: ## -## nw_spectrum --help +## nw_spectrum.py --help ## ## +## Examples, using data from NWChem QA reference output: +## tddft: nw_spectrum.py -v -i QA/tests/tddft_h2o/tddft_h2o.out -o tddft.dat +## vspec: TODO, no test cases exist! +## dos: nw_spectrum.py -v -f dos -i QA/tests/tddft_h2o/tddft_h2o.out -o dos.dat +## +## To plot an example file, in gnuplot, do e.g. +## plot "tddft.dat" with lines import sys import textwrap from optparse import OptionParser -ver = "2.1" +ver = "2.2" pname = "nw_spectrum" -def check_version (): - """Check version and ensure new print and string formatting is - allowed. Raises and exception if not satisfied.""" +class SpectrumProcessor(object): + def __init__(self, opts, instream, outstream): + self.instream = instream + self.outstream = outstream + self.opts = opts + self.preprocess_check_opts() + self.check_version() - if sys.version_info < (2, 6): - raise Exception("This script requires python >= 2.6") + def write(self, line): + """Write a line of output to the output stream""" - try: - newstring = "oldstring {v}".format(v=3.14) - except: - raise Exception("This script requires string.format()") - - -def ev2au(e_ev): - return (1.0 / 27.2114) * e_ev - - -def au2ev (e_au): - return 27.2114 * e_au - - -def ev2nm(e_ev): - return 27.2114 * 2.0 * 2.99792 * 2.41888 * 3.14159265359 / e_ev - - -def determine_data_type (): - """Parses stdin to see what data type, then rewinds stdin. - Returns 'vspec', 'tddft', or raises and Exception if neither - found. It choses based on the first tag found in the input, so - files with multiple data will only find the 1st. To extract the - 2nd, manually specify the data format via the command line args.""" - - tag_tddft = "NWChem TDDFT Module" - tag_vspec = "DFT Virtual Spectrum" - - lines = sys.stdin.readlines() - - for line in lines: - if tag_tddft in line: - sys.stdin.seek(0) - return "tddft" - elif tag_vspec in line: - sys.stdin.seek(0) - return "vspec" - - raise Exception ("Failed to determine data format, please specify manually.") - - -def parse_input_vspec (opts): - """Parses input from vspec and returns excitation energies in the - form [energy, f], in eV and atomic units units, respectively.""" - - lines = sys.stdin.readlines () - - inside_data = False - roots = [] - for l in lines: + self.outstream.write(line + "\n") - if "" in l: - try: - ls = l.split() - tag = ls[0] - nexcite = int (ls[1]) - except: - raise Exception ("Failed to parse tag and number: {0}".format(l)) + def check_version(self): + """Check version and ensure new print and string formatting is + allowed. Raises and exception if not satisfied.""" - iexcite = 0 - inside_data = True - continue + if sys.version_info < (2, 6): + raise Exception("This script requires python >= 2.6") - if inside_data: - if "" in l: - inside_data = False + try: + newstring = "oldstring {v}".format(v=3.14) + except: + raise Exception("This script requires string.format()") + + + def ev2au(self, e_ev): + return (1.0 / 27.2114) * e_ev + + def au2ev(self, e_au): + return 27.2114 * e_au + + def ev2nm(self, e_ev): + return 27.2114 * 2.0 * 2.99792 * 2.41888 * 3.14159265359 / e_ev + + def determine_data_type(self): + """Parses input to see what data type, then rewinds input. + Returns 'vspec' or 'tddft', or raises an exception if neither + found. It chosses based on the first tag found in the input, so + files with multiple data will only find the 1st. To extract the + 2nd, manually specify the data format via the command line args.""" + + tag_tddft = "NWChem TDDFT Module" + tag_vspec = "DFT Virtual Spectrum" + + lines = self.instream.readlines() + + for line in lines: + if tag_tddft in line: + self.instream.seek(0) + return "tddft" + elif tag_vspec in line: + self.instream.seek(0) + return "vspec" + + raise Exception("Failed to determine data format, please specify manually.") + + + def parse_input_vspec(self): + """Parses input from vspec and returns excitation energies in the + form [energy, f], in eV and atomic units units, respectively.""" + + lines = self.instream.readlines() + + inside_data = False + roots = [] + for l in lines: + if "" in l: + try: + ls = l.split() + tag = ls[0] + nexcite = int (ls[1]) + except: + raise Exception("Failed to parse tag and number: {0}".format(l)) + + iexcite = 0 + inside_data = True continue -# break - - try: - line_split = l.strip().split() - n = int (line_split[0]) - occ = int (line_split[1]) #not used - virtual = int (line_split[2]) #not used - energy_ev = float (line_split[3]) - osc = float (line_split[7]) - except: - raise Exception ("Failed to parse data line: {0}".format(l)) - - iexcite = iexcite + 1 - - if n != iexcite: - raise Exception ("Expected excitation number {0}, found {1}".format(iexcite, n)) - - - if energy_ev < 0.0: - print ("{0} Warning: Ignored negative vpsec excitation: {1} eV, {2}".format(opts.cchar, energy_ev, osc)) - - if opts.verbose: - sys.stderr.write ("Warning: Ignored negative vpsec excitation: {0} eV, {1}\n".format(energy_ev, osc)) - else: - roots.append ([energy_ev, osc]) + if inside_data: + if "" in l: + inside_data = False + continue + # break -# if not inside_data: -# raise Exception ("Failed to find tag") - - if iexcite != nexcite: - print ("{0} Warning: Expected {1} excitations, found {2}".format(opts.cchar, nexcite, iexcite)) - - if opts.verbose: - sys.stderr.write ("Warning: Expected {0} excitations, found {1}\n".format(nexcite,iexcite)) - - - if opts.verbose: - sys.stderr.write ("{0}: Found {1} vspec excitations\n".format(pname, len(roots))) - - - return roots - - -def parse_input_evals (opts): - """Parses input for eigenvalues and return as a list""" - - start_tag = "DFT Final Molecular Orbital Analysis" - end_tag = "Task times cpu" - - inside = False - - lines = sys.stdin.readlines () - iline = -1 - evals = [] - - while True: - iline = iline + 1 - try: - line = lines[iline] - except: - break - - if start_tag in line: - inside = True - - if end_tag in line: - inside = False - - if inside and "Vector" in line and "Occ" in line: - line_strip = line.strip() - - try: - tagloc = line_strip.rfind("E=") - evalue_str = line_strip[tagloc+2:tagloc+15].replace("D", "E") # after E= ; replace D with E - evalue = float(evalue_str) - except: - raise Exception ("Failed to parse eigenvalue: {0}".format(line_strip)) - - eval_ev = au2ev (evalue) - evals.append(eval_ev) #store eigenvalue in eV - - - if opts.verbose: - sys.stderr.write ("{0}: Found {1} eigenvalues\n".format(pname, len(evals))) - - return evals - - -def bin_evals (opts, evals): - """Take eigenvalues and bins them, and return [energy, N], where N - is the number of eigenvalues in the energy bin centered around - energy""" - - ## they should be sorted but let's make sure - evals_sort = sorted(evals) - - emin = evals_sort[0] - emax = evals_sort[-1] - - de = (emax - emin) / opts.nbin - - #=== XXX HARDCODE === -# de = 0.01 -# opts.nbin = int((emax - emin)/de) + 1 - #=== - - dos_raw = [] - for ie in range(opts.nbin+1): - ecenter = emin + ie*de - eleft = ecenter - 0.5*de - eright = ecenter + 0.5*de - - count = 0 - for val in evals_sort: - if val >= eleft and val <= eright: - count = count + 1 - - dos_raw.append ([ecenter, count]) - - - ## check that total sum is number of eigenvalues - ntot = 0 - for d in dos_raw: - ntot = ntot + d[1] - - if ntot != len (evals): - raise Exception ("Inconsistent integrated DOS and number of eigenvalues: {0} vs {1}".format(ntot, len(evals))) - - return dos_raw - - -def parse_input_tddft (opts): - """Parses input for singlet TDDFT excitation energies. - Returns excitation energies in the form [energy, f], in eV and - atomic units units, respectively.""" - - start_tag = "Convergence criterion met" - end_tag = "Excited state energy" - max_osc_search = 10 #max number of lines after root energy to look for oscillator strength - - - inside = False #true when we are inside output block - - lines = sys.stdin.readlines() - - iline = -1 - roots = [] - - while True: - - iline = iline + 1 - try: - line = lines[iline] - except: - break - - if start_tag in line: - inside = True - - if end_tag in line: - inside = False - - if inside and "Root" in line and "eV" in line: - - line_strip = line.strip() - - # ## - # ## OLD WAY - # ## - # ## Note, I would ideally .split() the line directly, but you - # ## can have cases like Root356 (for three digit root numbers) - # ## so I have to do it this ugly way... - # ## - # try: - # line_start = line_strip[0:4] # should contain RootXXX - # line_data = line_strip[7:].split() - # line_n = line_strip[4:7] #contains integer root number - # line_e = line_data[-2] # should contain excitation energy - # line_ev_tag = line_data[-1] # should contain "eV" - # except: - # raise Exception ("Failed to parse data line for root: {0}".format(line_strip)) - - ## - ## NEW WAY after Niri changed formatting, e.g.: - ## Root 48 singlet a 57.454053695 a.u. 1563.4050 eV - ## - line_strip = line.strip() - line_split = line_strip.split() - try: - line_start = line_split[0] # contains "Root" - line_n = line_split[1] # contains root number (int) - line_ev_tag = line_split[-1] # contains "eV" - line_e = line_split[-2] # contains excitation energy in eV - except: - raise Exception ("Failed to parse data line for root: {0}".format(line_strip)) - - if line_start == "Root" and line_ev_tag == "eV": try: - n = int(line_n) - energy_ev = float(line_e) + line_split = l.strip().split() + n = int (line_split[0]) + occ = int (line_split[1]) #not used + virtual = int (line_split[2]) #not used + energy_ev = float (line_split[3]) + osc = float (line_split[7]) except: - raise Exception ("Failed to convert root values: {0}".format(line_strip)) - else: - raise Exception ("Unexpected format for root: {0}".format(line_strip)) + raise Exception("Failed to parse data line: {0}".format(l)) + + iexcite = iexcite + 1 + + if n != iexcite: + raise Exception("Expected excitation number {0}, found {1}".format(iexcite, n)) + + if energy_ev < 0.0: + self.write("{0} Warning: Ignored negative vpsec excitation: {1} eV, {2}".format(opts.cchar, energy_ev, osc)) + if self.opts.verbose: + sys.stderr.write ("Warning: Ignored negative vpsec excitation: {0} eV, {1}\n".format(energy_ev, osc)) + else: + roots.append([energy_ev, osc]) + + + # if not inside_data: + # raise Exception("Failed to find tag") + + if iexcite != nexcite: + self.write("{0} Warning: Expected {1} excitations, found {2}".format(self.opts.cchar, nexcite, iexcite)) + + if self.opts.verbose: + sys.stderr.write("Warning: Expected {0} excitations, found {1}\n".format(nexcite, iexcite)) + + + if self.opts.verbose: + sys.stderr.write("{0}: Found {1} vspec excitations\n".format(pname, len(roots))) + + return roots + + + def parse_input_evals(self): + """Parses input for eigenvalues and return as a list""" + + start_tag = "DFT Final Molecular Orbital Analysis" + end_tag = "Task times cpu" + + inside = False + + lines = self.instream.readlines() + iline = -1 + evals = [] + + while True: + iline = iline + 1 + try: + line = lines[iline] + except: + break + + if start_tag in line: + inside = True + + if end_tag in line: + inside = False + + if inside and "Vector" in line and "Occ" in line: + line_strip = line.strip() - - if line_start == "Root" and line_ev_tag == "eV": try: - n = int(line_n) - energy_ev = float(line_e) + tagloc = line_strip.rfind("E=") + evalue_str = line_strip[tagloc + 2 : tagloc + 15].replace("D", "E") # after E= ; replace D with E + evalue = float(evalue_str) except: - raise Exception ("Failed to convert root values: {0}".format(line_strip)) - else: - raise Exception ("Unexpected format for root: {0}".format(line_strip)) + raise Exception("Failed to parse eigenvalue: {0}".format(line_strip)) + + eval_ev = self.au2ev(evalue) + evals.append(eval_ev) #store eigenvalue in eV + + if self.opts.verbose: + sys.stderr.write("{0}: Found {1} eigenvalues\n".format(pname, len(evals))) + + return evals - ## - ## Now look for oscillator strength, which will be a few - ## lines down (though the exact position may vary it seems). - ## - ioscline = -1 - while True: - ioscline = ioscline + 1 - if ioscline >= max_osc_search: - raise Exception ("Failed to find oscillator strength after looking {0} lines.".format(ioscline)) - - oscline = lines[iline + ioscline].strip() - - if "Dipole Oscillator Strength" in oscline: + def bin_evals(self, evals): + """Take eigenvalues and bins them, and return [energy, N], where N + is the number of eigenvalues in the energy bin centered around + energy""" + + ## they should be sorted but let's make sure + evals_sort = sorted(evals) + + emin = evals_sort[0] + emax = evals_sort[-1] + + de = (emax - emin) / self.opts.nbin + + #=== XXX HARDCODE === + # de = 0.01 + # self.opts.nbin = int((emax - emin)/de) + 1 + #=== + + dos_raw = [] + for ie in range(self.opts.nbin + 1): + ecenter = emin + ie * de + eleft = ecenter - 0.5 * de + eright = ecenter + 0.5 * de + + count = 0 + for val in evals_sort: + if val >= eleft and val <= eright: + count += 1 + + dos_raw.append([ecenter, count]) + + + ## check that total sum is number of eigenvalues + ntot = 0 + for d in dos_raw: + ntot = ntot + d[1] + + if ntot != len(evals): + raise Exception("Inconsistent integrated DOS and number of eigenvalues: {0} vs {1}".format(ntot, len(evals))) + + return dos_raw + + + def parse_input_tddft(self): + """Parses input for singlet TDDFT excitation energies. + Returns excitation energies in the form [energy, f], in eV and + atomic units units, respectively.""" + + start_tag = "Convergence criterion met" + end_tag = "Excited state energy" + singlet_tag = "singlet excited" + triplet_tag = "triplet excited" + state = "singlet" + + max_osc_search = 10 #max number of lines after root energy to look for oscillator strength + + + inside = False #true when we are inside output block + + lines = self.instream.readlines() + + iline = -1 + roots = [] + + while True: + iline = iline + 1 + try: + line = lines[iline] + except: + break + + #Only extract singlet state data. By default TDDFT calculations + #produce both singlet and triplet data, and the (possibly + #meaningless) triplet state data can otherwise interfere with + #proper parsing of singlets. + if start_tag in line and state == "singlet": + inside = True + + if end_tag in line: + inside = False + + if singlet_tag in line: + state = "singlet" + + if triplet_tag in line: + state = "triplet" + + if inside and "Root" in line and "eV" in line: + line_strip = line.strip() + line_split = line_strip.split() + try: + line_start = line_split[0] # contains "Root" + line_n = line_split[1] # contains root number (int) + line_ev_tag = line_split[-1] # contains "eV" + line_e = line_split[-2] # contains excitation energy in eV + except: + raise Exception("Failed to parse data line for root: {0}".format(line_strip)) + + if line_start == "Root" and line_ev_tag == "eV": try: - osc_str = oscline.split() - osc = float (osc_str[3]) + n = int(line_n) + energy_ev = float(line_e) except: - raise Exception ("Failed to convert oscillator strength: {0}".format(oscline)) - break - - - ## do some final checks, then append to data - if energy_ev < 0.0: - raise Exception ("Invalid negative energy: {0}".format(energy_ev)) + raise Exception("Failed to convert root values: {0}".format(line_strip)) + else: + raise Exception("Unexpected format for root: {0}".format(line_strip)) - if osc < 0.0: - raise Exception ("Invalid negative oscillator strength: {0}".format(osc)) + if line_start == "Root" and line_ev_tag == "eV": + try: + n = int(line_n) + energy_ev = float(line_e) + except: + raise Exception("Failed to convert root values: {0}".format(line_strip)) + else: + raise Exception("Unexpected format for root: {0}".format(line_strip)) + ## + ## Now look for oscillator strength, which will be a few + ## lines down (though the exact position may vary it seems). + ## + ioscline = -1 + while True: + ioscline = ioscline + 1 + if ioscline >= max_osc_search: + raise Exception("Failed to find oscillator strength after looking {0} lines.".format(ioscline)) - roots.append([energy_ev, osc]) + oscline = lines[iline + ioscline].strip() + + if "Dipole Oscillator Strength" in oscline: + try: + osc_str = oscline.split() + osc = float(osc_str[3]) + except: + raise Exception("Failed to convert oscillator strength: {0}".format(oscline)) + break + + ## do some final checks, then append to data + if energy_ev < 0.0: + raise Exception("Invalid negative energy: {0}".format(energy_ev)) + + if osc < 0.0: + raise Exception("Invalid negative oscillator strength: {0}".format(osc)) + + roots.append([energy_ev, osc]) - nroots = len (roots) + nroots = len(roots) - if nroots < 1: - raise Exception ("Failed to find any TDDFT roots") - else: - if opts.header: - print ("{0} Successfully parsed {1} TDDFT singlet excitations".format(opts.cchar,nroots)) + if nroots < 1: + raise Exception("Failed to find any TDDFT roots") + else: + if self.opts.header: + self.write("{0} Successfully parsed {1} TDDFT singlet excitations".format(self.opts.cchar, nroots)) - if opts.verbose: - sys.stderr.write ("{0}: Found {1} TDDFT excitations\n".format(pname, nroots)) + if self.opts.verbose: + sys.stderr.write("{0}: Found {1} TDDFT excitations\n".format(pname, nroots)) - return roots + return roots + def make_energy_list(self, roots): + """Computes the list of spectrum energies, and potentially adjusts + peak widths""" -def make_energy_list (opts, roots): - """Computes the list of spectrum energies, and potentially adjusts - peak widths""" + epad = 20.0 * self.opts.width - epad = 20.0*opts.width - - emin = roots[0][0] - epad + emin = roots[0][0] - epad -# if emin < opts.width: -# emin = opts.width + #if emin < self.opts.width: + # emin = self.opts.width - emax = roots[-1][0] + epad - de = (emax - emin) / opts.npoints + emax = roots[-1][0] + epad + de = (emax - emin) / self.opts.npoints - ## Use width of at least two grid points - if opts.width < 2*de: - opts.width = 2*de - print ("{0} Warning: Forced broadening to be {1} eV".format(opts.cchar, opts.width)) + ## Use width of at least two grid points + if self.opts.width < 2 * de: + self.opts.width = 2 * de + if self.opts.verbose: + sys.stderr.write("Warning: Forced broadening to be {0} eV\n".format(self.opts.width)) - if opts.verbose: - sys.stderr.write ("Warning: Forced broadening to be {0} eV\n".format(opts.width)) - -# opts.width = max (opts.width, 2*de) + eout = [emin + ie * de for ie in range(self.opts.npoints)] + return eout - eout = [ emin + ie*de for ie in range(opts.npoints) ] - return eout + def gen_spectrum(self, energies, roots): + """Generator for making Lorentzian broadenend spectrum.""" + + ## cutoff radius + #cutoff = 15.0 * self.opts.width + cutoff = 20.0 * self.opts.width + #cutoff = 9999999.0 * self.opts.width -### OLD SLOWER WAY -# def lorentzian_broaden (opts, roots): -# """Broadens raw roots into spectrum and returns the result.""" + ## + ## L(w; w0, gamma) = gamma/pi * 1 / [(w-w0)^2 + gamma^2] + ## + ## + ## multiply by 0.5 as FWHM was supplied by gamma in + ## lorenzian is actually HWHM. + ## + gamma = 0.5 * self.opts.width + gamma_sqrd = gamma * gamma -# ## cutoff radius -# cutoff = 15.0*opts.width - -# ## -# ## multiply by 0.5 as FWHM was supplied by gamma in lorenzian is -# ## actually HWHM. -# ## -# gamma = 0.5 * opts.width -# gamma_sqrd = gamma*gamma + de = (energies[-1] - energies[0]) / (len(energies) - 1) + prefac = gamma / 3.14159265359 * de -# ## -# ## L(w; w0, gamma) = gamma/pi * 1 / [(w-w0)^2 + gamma^2] -# ## -# prefac = gamma/3.14159265359*de + ## + ## used for verbose output (ie see progress for very large runs) + ## + # checkpt_energies = [] + # ne = len(energies) + # for ie in range(ne): + # if ie % ne == 0: + # checkpt_energies.append(energies[ie]) -# # spectrum = [ [emin + ie*de, 0] for ie in range(opts.npoints)] - -# # for point in spectrum: -# # stot = 0.0 -# # for root in roots: -# # xx0 = point[0] - root[0] + for energy in energies: + # if self.opts.verbose: + # if energy in checkpt_energies: + # sys.stderr.write("XXX\n") + + + stot = 0.0 + for root in roots: + xx0 = energy - root[0] + + if abs(xx0) <= cutoff: + stot += root[1] / ( xx0 * xx0 + gamma_sqrd) #Lorentzian + + yield [energy, stot * prefac] + + def dump_header(self): + """Add header to output""" + + if self.opts.header: + self.write("{0} ================================".format(self.opts.cchar)) + self.write("{0} NWChem spectrum parser ver {1}".format(self.opts.cchar, ver)) + self.write("{0} ================================".format(self.opts.cchar)) + self.write("{0} ".format(self.opts.cchar)) + self.write("{0} Parser runtime options: {1}".format(self.opts.cchar, self.opts)) + self.write("{0}".format(self.opts.cchar)) + + + def dump_data(self, roots): + """Dumps output. This works for either lists of raw roots or broadened + spectra.""" + + if self.opts.verbose: + sys.stderr.write("{0}: Dumping data ... \n".format(pname)) + + if self.opts.units == "ev": + if self.opts.header: + self.write("{0}".format(self.opts.cchar)) + self.write("{c}{s1}{delim}{s2}".format(c=self.opts.cchar, s1=" Energy [eV] ", delim=self.opts.delim, s2=" Abs. [au] ")) + self.write("{0}----------------------------------".format(self.opts.cchar)) + + for root in roots: + self.write("{energy:.10e}{delim}{osc:.10e}".format(energy=root[0], delim=self.opts.delim, osc=root[1])) + + + elif self.opts.units == "au": + if self.opts.header: + self.write("{0}".format(self.opts.cchar)) + self.write("{c}{s1}{delim}{s2}".format(c=self.opts.cchar, s1=" Energy [au] ", delim=self.opts.delim, s2=" Abs. [au] ")) + self.write("{0}----------------------------------".format(self.opts.cchar)) + + for root in roots: + eout = self.ev2au(root[0]) + self.write("{energy:.10e}{delim}{osc:.10e}".format(energy=eout, delim=self.opts.delim, osc=root[1])) + + elif self.opts.units == "nm": + if self.opts.header: + self.write("{0}".format(self.opts.cchar)) + self.write("{c}{s1}{delim}{s2}".format(c=self.opts.cchar, s1=" Wavelen. [nm] ", delim=self.opts.delim, s2=" Abs. [au] ")) + self.write("{0}----------------------------------".format(self.opts.cchar)) + + # roots.reverse () #reorder so we output in increasing wavelength + #XXX SHOULD REVERSE + for root in roots: + eout = self.ev2nm(root[0]) + self.write("{energy:.10e}{delim}{osc:.10e}".format(energy=eout, delim=self.opts.delim, osc=root[1])) + + else: + raise Exception("Invalid unit: {0}".format(self.opts.units)) + + def preprocess_check_opts(self): + """Check options and replaces with sane values if needed, stores + all opts as purely lowercase.""" + + self.opts.units = self.opts.units.lower() + self.opts.datafmt = self.opts.datafmt.lower() + + if self.opts.units not in ("nm", "ev", "au"): + raise Exception("Invalid unit type: {0}".format(self.opts.units)) + + if self.opts.datafmt not in ("tddft", "vspec", "auto", "dos"): + raise Exception("Invalid data format: {0}".format(self.opts.datafmt)) + + if self.opts.npoints < 100: + raise Exception("Increase number of points to at least 100 (you asked for {0})".format(self.opts.npoints)) + + if self.opts.width < 0.0: + raise Exception("Peak width must be positive (you supplied {0})".format(self.opts.width)) + + def process_all(self): + if self.opts.header: + self.dump_header() + + if self.opts.datafmt == "auto": + self.opts.datafmt = self.determine_data_type() + if self.opts.datafmt == "tddft": + self.write("{0} The input appears to contain TDDFT data.".format(self.opts.cchar)) + elif self.opts.datafmt == "vspec": + self.write("{0} The input appears to contain vspec data.".format(self.opts.cchar)) + else: + raise Exception("Invalid data format: {0}".format(self.opts.datafmt)) + + ## parse raw data + if self.opts.datafmt == "tddft": + roots = self.parse_input_tddft() + elif self.opts.datafmt == "vspec": + roots = self.parse_input_vspec() + elif self.opts.datafmt == "dos": + evals = self.parse_input_evals() + roots = self.bin_evals(evals) + else: + raise Exception("Invalid data format supplied: {0}".format(self.opts.datafmt)) + + + ## make broadened spectrum if desired + if self.opts.makespec: + energies = self.make_energy_list(roots) + spectrum = self.gen_spectrum(energies, roots) + + if self.opts.verbose: + sys.stderr.write("{0}: Initialized spectrum [{1: .3f} : {2: .3f}] ({3} points)\n".format(pname, energies[0], energies[-1], len(energies))) + + if self.opts.header: + self.write("{0} Roots were Lorentzian broadened with width {1} eV ".format(self.opts.cchar, self.opts.width)) + self.write("{0} Spectrum generated with {1} points.".format(self.opts.cchar, self.opts.npoints)) + + else: + spectrum = roots + if self.opts.header: + self.write("{0} Roots not broadened".format(self.opts.cchar)) + self.write("{0} No spectrum generated: output is list of raw excitations".format(self.opts.cchar)) + + + ## finally, dump data to stdout + self.dump_data(spectrum) + + if self.opts.verbose: + sys.stderr.write ("{0}: Done\n".format(pname)) -# # if abs(xx0) <= cutoff: -# # stot += prefac * root[1] / ( xx0*xx0 + gamma_sqrd) #Lorentzian - -# # point[1] = stot - - -# # npoints_made = len(spectrum) -# # if npoints_made != opts.npoints: -# # raise Exception ("Spectrum should have {0} points, instead has {1}".format(opts.npoints, npoints_made)) - -# # if opts.header: -# # print ("{0} Lorentzian broadened roots with width {1} eV".format(opts.cchar, opts.width)) -# # print ("{0} Created spectrum with {1} points".format(opts.cchar, npoints_made)) - -# return spectrum - - - -## Faster way--use a generator -def gen_spectrum (opts, energies, roots): - """Generator for making Lorentzian broadenend spectrum.""" - - ## cutoff radius -# cutoff = 15.0*opts.width - cutoff = 20.0*opts.width -# cutoff = 9999999.0*opts.width - - - ## - ## L(w; w0, gamma) = gamma/pi * 1 / [(w-w0)^2 + gamma^2] - ## - ## - ## multiply by 0.5 as FWHM was supplied by gamma in - ## lorenzian is actually HWHM. - ## - gamma = 0.5 * opts.width - gamma_sqrd = gamma*gamma - - de = (energies[-1] - energies[0]) / (len(energies)-1) - prefac = gamma/3.14159265359*de - - - ## - ## used for verbose output (ie see progress for very large runs) - ## - # checkpt_energies = [] - # ne = len(energies) - # for ie in range(ne): - # if ie % ne == 0: - # checkpt_energies.append(energies[ie]) - - - for energy in energies: - - # if opts.verbose: - # if energy in checkpt_energies: - # sys.stderr.write ("XXX\n") - - - stot = 0.0 - for root in roots: - xx0 = energy - root[0] - - if abs(xx0) <= cutoff: - stot += root[1] / ( xx0*xx0 + gamma_sqrd) #Lorentzian - - yield [energy, stot*prefac] - - - - - -def dump_header (opts): - """Print header to stdout""" - - if opts.header: - print ("{0} ================================".format(opts.cchar)) - print ("{0} NWChem spectrum parser ver {1}".format(opts.cchar,ver)) - print ("{0} ================================".format(opts.cchar)) - print ("{0} ".format(opts.cchar)) - print ("{0} Parser runtime options: {1}".format(opts.cchar,opts)) - print ("{0}".format(opts.cchar)) - - -def dump_data (opts, roots): - """Dumps output to stdout. This works for either lists of raw roots - or broadened spectra.""" - - if opts.verbose: - sys.stderr.write ("{0}: Dumping data to stdout ... \n".format(pname)) - - if opts.units == "ev": - if opts.header: - print ("{0}".format(opts.cchar)) - print ("{c}{s1}{delim}{s2}".format(c=opts.cchar, s1=" Energy [eV] ", delim=opts.delim, s2=" Abs. [au] ")) - print ("{0}----------------------------------".format(opts.cchar)) - - for root in roots: - print ("{energy:.10e}{delim}{osc:.10e}".format(energy=root[0], delim=opts.delim, osc=root[1])) - - - elif opts.units == "au": - if opts.header: - print ("{0}".format(opts.cchar)) - print ("{c}{s1}{delim}{s2}".format(c=opts.cchar, s1=" Energy [au] ", delim=opts.delim, s2=" Abs. [au] ")) - print ("{0}----------------------------------".format(opts.cchar)) - - for root in roots: - eout = ev2au (root[0]) - print ("{energy:.10e}{delim}{osc:.10e}".format(energy=eout, delim=opts.delim, osc=root[1])) - - elif opts.units == "nm": - if opts.header: - print ("{0}".format(opts.cchar)) - print ("{c}{s1}{delim}{s2}".format(c=opts.cchar, s1=" Wavelen. [nm] ", delim=opts.delim, s2=" Abs. [au] ")) - print ("{0}----------------------------------".format(opts.cchar)) - -# roots.reverse () #reorder so we output in increasing wavelength -#XXX SHOULD REVERSE - for root in roots: - eout = ev2nm (root[0]) - print ("{energy:.10e}{delim}{osc:.10e}".format(energy=eout, delim=opts.delim, osc=root[1])) - - else: - raise Exception ("Invalid unit: {0}".format(opts.units)) - - -def preprocess_check_opts (opts): - """Check options and replaces with sane values if needed, stores - all opts as purely lowercase.""" - - opts.units = opts.units.lower() - opts.datafmt = opts.datafmt.lower() - - if opts.units != "nm" and opts.units != "ev" and opts.units != "au": - raise Exception ("Invalid unit type: {0}".format(opts.units)) - - if opts.datafmt != "tddft" and opts.datafmt != "vspec" and opts.datafmt != "auto" and opts.datafmt != "dos": - raise Exception ("Invalid data format: {0}".format(opts.datafmt)) - - if opts.npoints < 100: - raise Exception ("Increase number of points to at least 100 (you asked for {0})".format(opts.npoints)) - - if opts.width < 0.0: - raise Exception ("Peak width must be positive (you supplied {0})".format(opts.width)) def main(): - ## ## Parse command line options. Note we later make all lowercase, ## so from the user's point of view they are case-insensitive. ## usage = "%prog [options]\n\n" - desc = "Reads NWChem output from stdin, parses for the linear response TDDFT or DFT vspec excitations, and prints the absorption spectrum to stdout. It will optionally broaden peaks using a Lorentzian with FWHM of at least two energy/wavelength spacings. By default, it will automatically determine data format (tddft or vspec) and generate a broadened spectrum in eV." - desc_wrap = textwrap.wrap (desc,80) + desc = "Reads NWChem output from stdin or specified infile, parses for the linear response TDDFT or DFT vspec excitations, and prints the absorption spectrum to stdout or writes it to outfile. It will optionally broaden peaks using a Lorentzian with FWHM of at least two energy/wavelength spacings. By default, it will automatically determine data format (tddft or vspec) and generate a broadened spectrum in eV." + desc_wrap = textwrap.wrap(desc, 80) for s in desc_wrap: usage += s + "\n" example = 'Create absorption spectrum in nm named "spectrum.dat" from the NWChem output file "water.nwo" named spectrum.dat with peaks broadened by 0.3 eV and 5000 points in the spectrum.' - example_wrap = textwrap.wrap (example,80) - usage += "\nExample:\n\n\t"+"nw_spectrum -b0.3 -p5000 -wnm < water.nwo > spectrum.dat\n\n" + example_wrap = textwrap.wrap(example, 80) + usage += "\nExample:\n\n\t"+"nw_spectrum -b0.3 -p5000 -wnm < water.nwo > spectrum.dat\n\n\tor:\n\n\tnw_spectrum -b0.3 -p5000 -wnm -i water.nwo -o spectrum.dat\n\n" for s in example_wrap: usage += s + "\n" - parser = OptionParser(usage=usage) + parser.add_option("-i", "--infile", type="string", dest="infile", + help="input file (reads from stdin if unspecified)", metavar="FILE") + parser.add_option("-o", "--outfile", type="string", dest="outfile", + help="output file (writes to stdout if unspecified)", metavar="FILE") parser.add_option("-f", "--format", type="string", dest="datafmt", help="data file format: auto (default), tddft, vspec, dos", metavar="FMT") parser.add_option("-b", "--broad", type="float", dest="width", @@ -603,61 +596,25 @@ def main(): help="echo warnings and progress to stderr") - parser.set_defaults(width = 0.1, npoints = 2000, units="eV", datafmt="auto", delim=" ", makespec=True, header=True, cchar="#", nbin=20, verbose=False) + parser.set_defaults(width=0.1, npoints=2000, units="eV", datafmt="auto", + delim=" ", makespec=True, header=True, cchar="#", + nbin=20, verbose=False) (opts, args) = parser.parse_args() - preprocess_check_opts(opts) - check_version () - - if opts.header: - dump_header (opts) - - if opts.datafmt == "auto": - opts.datafmt = determine_data_type () - if opts.datafmt == "tddft": - print ("{0} The input appears to contain TDDFT data.".format(opts.cchar)) - elif opts.datafmt == "vspec": - print ("{0} The input appears to contain vspec data.".format(opts.cchar)) - else: - raise Exception ("Invalid data format: {0}".format(opts.datafmt)) - - - ## parse raw data - if opts.datafmt == "tddft": - roots = parse_input_tddft (opts) - elif opts.datafmt == "vspec": - roots = parse_input_vspec (opts) - elif opts.datafmt == "dos": - evals = parse_input_evals (opts) - roots = bin_evals (opts, evals) + if opts.infile: + instream = open(opts.infile, "r") else: - raise Exception ("Invalid data format supplied: {0}".format(opts.datafmt)) + instream = sys.stdin - - ## make broadened spectrum if desired - if opts.makespec: - energies = make_energy_list (opts, roots) - spectrum = gen_spectrum (opts, energies, roots) - - if opts.verbose: - sys.stderr.write ("{0}: Initialized spectrum [{1: .3f} : {2: .3f}] ({3} points)\n".format(pname, energies[0], energies[-1], len(energies))) - - if opts.header: - print ("{0} Roots were Lorentzian broadened with width {1} eV ".format(opts.cchar, opts.width)) - print ("{0} Spectrum generated with {1} points.".format(opts.cchar, opts.npoints)) - + if opts.outfile: + outstream = open(opts.outfile, "w") else: - spectrum = roots - if opts.header: - print ("{0} Roots not broadened".format(opts.cchar)) - print ("{0} No spectrum generated: output is list of raw excitations".format(opts.cchar)) - - - ## finally, dump data to stdout - dump_data (opts, spectrum) + outstream = sys.stdout - if opts.verbose: - sys.stderr.write ("{0}: Done\n".format(pname)) + S = SpectrumProcessor(opts, instream, outstream) + S.process_all() + + if __name__ == "__main__": From f90b3ac46377f252c7f70d9c83a1496c0836283a Mon Sep 17 00:00:00 2001 From: edoapra Date: Wed, 14 Jun 2023 11:09:13 -0700 Subject: [PATCH 2/4] fix error spotted by flake8 --- contrib/parsers/nw_spectrum.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/contrib/parsers/nw_spectrum.py b/contrib/parsers/nw_spectrum.py index 34f0564f5c..380b3b5b00 100644 --- a/contrib/parsers/nw_spectrum.py +++ b/contrib/parsers/nw_spectrum.py @@ -129,7 +129,7 @@ class SpectrumProcessor(object): raise Exception("Expected excitation number {0}, found {1}".format(iexcite, n)) if energy_ev < 0.0: - self.write("{0} Warning: Ignored negative vpsec excitation: {1} eV, {2}".format(opts.cchar, energy_ev, osc)) + self.write("{0} Warning: Ignored negative vpsec excitation: {1} eV, {2}".format(self.opts.cchar, energy_ev, osc)) if self.opts.verbose: sys.stderr.write ("Warning: Ignored negative vpsec excitation: {0} eV, {1}\n".format(energy_ev, osc)) else: From 479934e8b5e912d9e1e52e0b3d667ec0154418f6 Mon Sep 17 00:00:00 2001 From: Daniel Mejia-Rodriguez Date: Wed, 14 Jun 2023 12:30:10 -0700 Subject: [PATCH 3/4] Update nw_spectrum.py Close files at the end --- contrib/parsers/nw_spectrum.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/contrib/parsers/nw_spectrum.py b/contrib/parsers/nw_spectrum.py index 380b3b5b00..5b5ca641fb 100644 --- a/contrib/parsers/nw_spectrum.py +++ b/contrib/parsers/nw_spectrum.py @@ -276,6 +276,7 @@ class SpectrumProcessor(object): if end_tag in line: inside = False + break if singlet_tag in line: state = "singlet" @@ -614,6 +615,10 @@ def main(): S = SpectrumProcessor(opts, instream, outstream) S.process_all() + if opts.infile: + instream.close() + if opts.outfile: + outstream.close() From 1f941b9cd266aa68f886ae5848f20413e86d8902 Mon Sep 17 00:00:00 2001 From: edoapra Date: Wed, 14 Jun 2023 13:32:24 -0700 Subject: [PATCH 4/4] chmod [ci skip] --- contrib/parsers/nw_spectrum.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100644 => 100755 contrib/parsers/nw_spectrum.py diff --git a/contrib/parsers/nw_spectrum.py b/contrib/parsers/nw_spectrum.py old mode 100644 new mode 100755