← Back to branch summary

~ubuntu-branches/ubuntu/utopic/nwchem/utopic

« back to all changes in this revision

Viewing changes to contrib/parsers/nw_spectrum.py

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Daniel Leidert, Andreas Tille, Michael Banck
  • Date: 2013-07-04 12:14:55 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130704121455-5tvsx2qabor3nrui
Tags: 6.3-1
http://bugs.debian.org/696361
* New upstream release.
* Fixes anisotropic properties (Closes: #696361).
* New features include:
  + Multi-reference coupled cluster (MRCC) approaches
  + Hybrid DFT calculations with short-range HF 
  + New density-functionals including Minnesota (M08, M11) and HSE hybrid
    functionals
  + X-ray absorption spectroscopy (XAS) with TDDFT
  + Analytical gradients for the COSMO solvation model
  + Transition densities from TDDFT 
  + DFT+U and Electron-Transfer (ET) methods for plane wave calculations
  + Exploitation of space group symmetry in plane wave geometry optimizations
  + Local density of states (LDOS) collective variable added to Metadynamics
  + Various new XC functionals added for plane wave calculations, including
    hybrid and range-corrected ones
  + Electric field gradients with relativistic corrections 
  + Nudged Elastic Band optimization method
  + Updated basis sets and ECPs 

[ Daniel Leidert ]
* debian/watch: Fixed.

[ Andreas Tille ]
* debian/upstream: References

[ Michael Banck ]
* debian/upstream (Name): New field.
* debian/patches/02_makefile_flags.patch: Refreshed.
* debian/patches/06_statfs_kfreebsd.patch: Likewise.
* debian/patches/07_ga_target_force_linux.patch: Likewise.
* debian/patches/05_avoid_inline_assembler.patch: Removed, no longer needed.
* debian/patches/09_backported_6.1.1_fixes.patch: Likewise.
* debian/control (Build-Depends): Added gfortran-4.7 and gcc-4.7.
* debian/patches/10_force_gcc-4.7.patch: New patch, explicitly sets
  gfortran-4.7 and gcc-4.7, fixes test suite hang with gcc-4.8 (Closes:
  #701328, #713262).
* debian/testsuite: Added tests for COSMO analytical gradients and MRCC.
* debian/rules (MRCC_METHODS): New variable, required to enable MRCC methods.

Show diffs side-by-side

added added

removed removed

Lines of Context:
contrib/parsers/nw_spectrum.py
 
1
##
 
2
## nw_spectrum
 
3
##
 
4
## Kenneth Lopata
 
5
## Last modified: 2013-03-08
 
6
##
 
7
## Python script for parsing NWChem output for TDDFT/vspec excitation
 
8
## energies, and optionally Lorentzian broadenening the spectra.  For
 
9
## online help run:
 
10
##
 
11
## nw_spectrum --help
 
12
##
 
13
##
 
14
 
 
15
import sys
 
16
import textwrap
 
17
from optparse import OptionParser
 
18
 
 
19
ver = "2.1"
 
20
pname = "nw_spectrum"
 
21
 
 
22
 
 
23
def check_version ():
 
24
    """Check version and ensure new print and string formatting is
 
25
    allowed.  Raises and exception if not satisfied."""
 
26
 
 
27
    if sys.version_info < (2, 6):
 
28
        raise Exception("This script requires python >= 2.6")
 
29
 
 
30
    try:
 
31
        newstring = "oldstring {v}".format(v=3.14)
 
32
    except:
 
33
        raise Exception("This script requires string.format()")
 
34
 
 
35
 
 
36
def ev2au(e_ev):
 
37
    return (1.0 / 27.2114) * e_ev
 
38
 
 
39
 
 
40
def au2ev (e_au):
 
41
    return 27.2114 * e_au
 
42
 
 
43
 
 
44
def ev2nm(e_ev):
 
45
    return 27.2114 * 2.0 * 2.99792 * 2.41888 * 3.14159265359 / e_ev
 
46
 
 
47
 
 
48
def determine_data_type ():
 
49
    """Parses stdin to see what data type, then rewinds stdin.
 
50
    Returns 'vspec', 'tddft', or raises and Exception if neither
 
51
    found.  It choses based on the first tag found in the input, so
 
52
    files with multiple data will only find the 1st.  To extract the
 
53
    2nd, manually specify the data format via the command line args."""
 
54
    
 
55
    tag_tddft = "NWChem TDDFT Module"
 
56
    tag_vspec = "DFT Virtual Spectrum"
 
57
 
 
58
    lines = sys.stdin.readlines()
 
59
    
 
60
    for line in lines:
 
61
        if tag_tddft in line:
 
62
            sys.stdin.seek(0)
 
63
            return "tddft"
 
64
        elif tag_vspec in line:
 
65
            sys.stdin.seek(0)
 
66
            return "vspec"
 
67
 
 
68
    raise Exception ("Failed to determine data format, please specify manually.")
 
69
 
 
70
 
 
71
def parse_input_vspec (opts):
 
72
    """Parses input from vspec and returns excitation energies in the
 
73
    form [energy, f], in eV and atomic units units, respectively."""
 
74
 
 
75
    lines = sys.stdin.readlines ()
 
76
    
 
77
    inside_data = False
 
78
    roots = []
 
79
    for l in lines:
 
80
        
 
81
        if "<START>" in l:
 
82
            try:
 
83
                ls = l.split()
 
84
                tag = ls[0]
 
85
                nexcite = int (ls[1])
 
86
            except:
 
87
                raise Exception ("Failed to parse <START> tag and number: {0}".format(l))
 
88
 
 
89
            iexcite = 0
 
90
            inside_data = True
 
91
            continue
 
92
 
 
93
        if inside_data:
 
94
            if "<END>" in l:
 
95
                inside_data = False
 
96
                continue
 
97
#                break
 
98
            
 
99
            try:
 
100
                line_split = l.strip().split()
 
101
                n = int (line_split[0])
 
102
                occ = int (line_split[1])          #not used
 
103
                virtual = int (line_split[2])      #not used
 
104
                energy_ev = float (line_split[3])
 
105
                osc = float (line_split[7])
 
106
            except:
 
107
                raise Exception ("Failed to parse data line: {0}".format(l))
 
108
            
 
109
            iexcite = iexcite + 1
 
110
            
 
111
            if n != iexcite:
 
112
                raise Exception ("Expected excitation number {0}, found {1}".format(iexcite, n))
 
113
            
 
114
            
 
115
            if energy_ev < 0.0:
 
116
                print ("{0} Warning: Ignored negative vpsec excitation: {1} eV, {2}".format(opts.cchar, energy_ev, osc))
 
117
                
 
118
                if opts.verbose:
 
119
                    sys.stderr.write ("Warning: Ignored negative vpsec excitation: {0} eV, {1}\n".format(energy_ev, osc))
 
120
            else:
 
121
                roots.append ([energy_ev, osc])
 
122
 
 
123
 
 
124
#    if not inside_data:
 
125
#        raise Exception ("Failed to find <START> tag")
 
126
 
 
127
    if iexcite != nexcite:
 
128
        print ("{0} Warning: Expected {1} excitations, found {2}".format(opts.cchar, nexcite, iexcite))
 
129
        
 
130
        if opts.verbose:
 
131
            sys.stderr.write ("Warning: Expected {0} excitations, found {1}\n".format(nexcite,iexcite))
 
132
 
 
133
 
 
134
    if opts.verbose:
 
135
        sys.stderr.write ("{0}: Found {1} vspec excitations\n".format(pname, len(roots)))
 
136
            
 
137
            
 
138
    return roots
 
139
 
 
140
 
 
141
def parse_input_evals (opts):
 
142
    """Parses input for eigenvalues and return as a list"""
 
143
 
 
144
    start_tag = "DFT Final Molecular Orbital Analysis"
 
145
    end_tag = "Task  times  cpu"
 
146
 
 
147
    inside = False
 
148
    
 
149
    lines = sys.stdin.readlines ()
 
150
    iline = -1
 
151
    evals = []
 
152
 
 
153
    while True:
 
154
        iline = iline + 1
 
155
        try:
 
156
            line = lines[iline]
 
157
        except:
 
158
            break
 
159
 
 
160
        if start_tag in line:
 
161
            inside = True
 
162
 
 
163
        if end_tag in line:
 
164
            inside = False
 
165
 
 
166
        if inside and "Vector" in line and "Occ" in line:
 
167
            line_strip = line.strip()
 
168
 
 
169
            try:
 
170
                tagloc = line_strip.rfind("E=")
 
171
                evalue_str = line_strip[tagloc+2:tagloc+15].replace("D", "E")  # after E=  ; replace D with E
 
172
                evalue = float(evalue_str)
 
173
            except:
 
174
                raise Exception ("Failed to parse eigenvalue: {0}".format(line_strip))
 
175
 
 
176
            eval_ev = au2ev (evalue)
 
177
            evals.append(eval_ev)  #store eigenvalue in eV
 
178
 
 
179
 
 
180
    if opts.verbose:
 
181
        sys.stderr.write ("{0}: Found {1} eigenvalues\n".format(pname, len(evals)))
 
182
 
 
183
    return evals
 
184
 
 
185
 
 
186
def bin_evals (opts, evals):
 
187
    """Take eigenvalues and bins them, and return [energy, N], where N
 
188
    is the number of eigenvalues in the energy bin centered around
 
189
    energy"""
 
190
 
 
191
    ## they should be sorted but let's make sure
 
192
    evals_sort = sorted(evals)
 
193
 
 
194
    emin = evals_sort[0]
 
195
    emax = evals_sort[-1]
 
196
 
 
197
    de = (emax - emin) / opts.nbin
 
198
    
 
199
    #=== XXX HARDCODE ===
 
200
#    de = 0.01
 
201
#    opts.nbin = int((emax - emin)/de) + 1
 
202
    #=== 
 
203
 
 
204
    dos_raw = []
 
205
    for ie in range(opts.nbin+1):
 
206
        ecenter = emin + ie*de
 
207
        eleft = ecenter - 0.5*de
 
208
        eright = ecenter + 0.5*de
 
209
 
 
210
        count = 0
 
211
        for val in evals_sort:
 
212
            if val >= eleft and val <= eright:
 
213
                count = count + 1
 
214
 
 
215
        dos_raw.append ([ecenter, count])
 
216
 
 
217
 
 
218
    ## check that total sum is number of eigenvalues
 
219
    ntot = 0
 
220
    for d in dos_raw:
 
221
        ntot = ntot + d[1]
 
222
 
 
223
    if ntot != len (evals):
 
224
        raise Exception ("Inconsistent integrated DOS and number of eigenvalues: {0} vs {1}".format(ntot, len(evals)))
 
225
 
 
226
    return dos_raw
 
227
 
 
228
 
 
229
def parse_input_tddft (opts):
 
230
    """Parses input for singlet TDDFT excitation energies.
 
231
    Returns excitation energies in the form [energy, f], in eV and
 
232
    atomic units units, respectively."""
 
233
 
 
234
    start_tag = "Convergence criterion met"
 
235
    end_tag = "Excited state energy"
 
236
    max_osc_search = 10  #max number of lines after root energy to look for oscillator strength
 
237
 
 
238
 
 
239
    inside = False  #true when we are inside output block
 
240
    
 
241
    lines = sys.stdin.readlines()
 
242
 
 
243
    iline = -1
 
244
    roots = []
 
245
    
 
246
    while True:
 
247
        
 
248
        iline = iline + 1
 
249
        try:
 
250
            line = lines[iline]
 
251
        except:
 
252
            break
 
253
 
 
254
        if start_tag in line:
 
255
            inside = True
 
256
 
 
257
        if end_tag in line:
 
258
            inside = False
 
259
 
 
260
        if inside and "Root" in line and "eV" in line:
 
261
 
 
262
            line_strip = line.strip()
 
263
 
 
264
            ##
 
265
            ## Note, I would ideally .split() the line directly, but you
 
266
            ## can have cases like Root356 (for three digit root numbers)
 
267
            ## so I have to do it this ugly way...
 
268
            ##
 
269
            try:
 
270
                line_start = line_strip[0:4]          # should contain RootXXX
 
271
                line_data = line_strip[7:].split()   
 
272
                line_n = line_strip[4:7]              #contains integer root number
 
273
                line_e = line_data[-2]                # should contain excitation energy
 
274
                line_ev_tag = line_data[-1]           # should contain "eV"
 
275
            except:
 
276
                raise Exception ("Failed to parse data line for root: {0}".format(line_strip))
 
277
            
 
278
#            print (line_start, line_n, line_data)  #debugging
 
279
 
 
280
            if line_start == "Root" and line_ev_tag == "eV":
 
281
                try:
 
282
                    n = int(line_n)
 
283
                    energy_ev = float(line_e)
 
284
                except:
 
285
                    raise Exception ("Failed to convert root values: {0}".format(line_strip))
 
286
            else:
 
287
                raise Exception ("Unexpected format for root: {0}".format(line_strip))
 
288
 
 
289
 
 
290
            ##
 
291
            ## Now look for oscillator strength, which will be a few
 
292
            ## lines down (though the exact position may vary it seems).
 
293
            ##
 
294
            ioscline = -1
 
295
            while True:
 
296
                ioscline = ioscline + 1
 
297
                if ioscline >= max_osc_search:
 
298
                    raise Exception ("Failed to find oscillator strength after looking {0} lines.".format(ioscline))
 
299
                
 
300
                oscline = lines[iline + ioscline].strip()
 
301
                
 
302
                if "Dipole Oscillator Strength" in oscline:
 
303
                    try:
 
304
                        osc_str = oscline.split()
 
305
                        osc = float (osc_str[3])
 
306
                    except:
 
307
                        raise Exception ("Failed to convert oscillator strength: {0}".format(oscline))
 
308
                    break
 
309
            
 
310
            
 
311
            ## do some final checks, then append to data
 
312
            if energy_ev < 0.0:
 
313
                raise Exception ("Invalid negative energy: {0}".format(energy_ev))
 
314
 
 
315
            if osc < 0.0:
 
316
                raise Exception ("Invalid negative oscillator strength: {0}".format(osc))
 
317
 
 
318
            roots.append([energy_ev, osc])
 
319
 
 
320
 
 
321
    nroots = len (roots)
 
322
 
 
323
    if nroots < 1:
 
324
        raise Exception ("Failed to find any TDDFT roots")
 
325
    else:
 
326
        if opts.header:
 
327
            print ("{0} Successfully parsed {1} TDDFT singlet excitations".format(opts.cchar,nroots))
 
328
 
 
329
 
 
330
    if opts.verbose:
 
331
        sys.stderr.write ("{0}: Found {1} TDDFT excitations\n".format(pname, nroots))
 
332
 
 
333
    return roots
 
334
 
 
335
 
 
336
 
 
337
def make_energy_list (opts, roots):
 
338
    """Computes the list of spectrum energies, and potentially adjusts
 
339
    peak widths"""
 
340
 
 
341
    epad = 20.0*opts.width
 
342
    
 
343
    emin = roots[0][0] - epad
 
344
 
 
345
#    if emin < opts.width:
 
346
#        emin = opts.width
 
347
 
 
348
    emax = roots[-1][0] + epad
 
349
    de = (emax - emin) / opts.npoints
 
350
 
 
351
    ## Use width of at least two grid points
 
352
    if opts.width < 2*de:
 
353
        opts.width = 2*de
 
354
        print ("{0} Warning: Forced broadening to be {1} eV".format(opts.cchar, opts.width))
 
355
 
 
356
        if opts.verbose:
 
357
            sys.stderr.write ("Warning: Forced broadening to be {0} eV\n".format(opts.width))
 
358
    
 
359
#    opts.width = max (opts.width, 2*de)
 
360
 
 
361
    eout = [ emin + ie*de for ie in range(opts.npoints) ]
 
362
    return eout
 
363
 
 
364
 
 
365
### OLD SLOWER WAY
 
366
# def lorentzian_broaden (opts, roots):
 
367
#     """Broadens raw roots into spectrum and returns the result."""
 
368
 
 
369
#     ## cutoff radius
 
370
#     cutoff = 15.0*opts.width
 
371
 
 
372
#     ##
 
373
#     ## multiply by 0.5 as FWHM was supplied by gamma in lorenzian is
 
374
#     ## actually HWHM.
 
375
#     ##
 
376
#     gamma = 0.5 * opts.width
 
377
#     gamma_sqrd = gamma*gamma
 
378
 
 
379
 
 
380
#     ##
 
381
#     ## L(w; w0, gamma) = gamma/pi * 1 / [(w-w0)^2 + gamma^2]
 
382
#     ##
 
383
#     prefac = gamma/3.14159265359*de
 
384
 
 
385
 
 
386
#     # spectrum = [ [emin + ie*de, 0] for ie in range(opts.npoints)]
 
387
    
 
388
#     # for point in spectrum:
 
389
#     #     stot = 0.0
 
390
#     #     for root in roots:
 
391
#     #         xx0 = point[0] - root[0]
 
392
        
 
393
#     #         if abs(xx0) <= cutoff:
 
394
#     #             stot += prefac * root[1] / ( xx0*xx0 + gamma_sqrd)   #Lorentzian
 
395
            
 
396
#     #         point[1] = stot
 
397
 
 
398
 
 
399
#     # npoints_made = len(spectrum)
 
400
#     # if npoints_made != opts.npoints:
 
401
#     #     raise Exception ("Spectrum should have {0} points, instead has {1}".format(opts.npoints, npoints_made))
 
402
 
 
403
#     # if opts.header:
 
404
#     #     print ("{0} Lorentzian broadened roots with width {1} eV".format(opts.cchar, opts.width))
 
405
#     #     print ("{0} Created spectrum with {1} points".format(opts.cchar, npoints_made))
 
406
 
 
407
#     return spectrum
 
408
 
 
409
 
 
410
 
 
411
## Faster way--use a generator
 
412
def gen_spectrum (opts, energies, roots):
 
413
    """Generator for making Lorentzian broadenend spectrum."""
 
414
    
 
415
    ## cutoff radius
 
416
#    cutoff = 15.0*opts.width
 
417
    cutoff = 20.0*opts.width
 
418
#    cutoff = 9999999.0*opts.width
 
419
 
 
420
    
 
421
    ##
 
422
    ## L(w; w0, gamma) = gamma/pi * 1 / [(w-w0)^2 + gamma^2]
 
423
    ##
 
424
    ##
 
425
    ## multiply by 0.5 as FWHM was supplied by gamma in
 
426
    ## lorenzian is actually HWHM.
 
427
    ##
 
428
    gamma = 0.5 * opts.width
 
429
    gamma_sqrd = gamma*gamma
 
430
 
 
431
    de = (energies[-1] - energies[0]) / (len(energies)-1)
 
432
    prefac = gamma/3.14159265359*de
 
433
 
 
434
    
 
435
    ##
 
436
    ## used for verbose output (ie see progress for very large runs)
 
437
    ##
 
438
    # checkpt_energies = []
 
439
    # ne = len(energies)
 
440
    # for ie in range(ne):
 
441
    #     if ie % ne == 0:
 
442
    #         checkpt_energies.append(energies[ie])
 
443
 
 
444
 
 
445
    for energy in energies:
 
446
 
 
447
        # if opts.verbose:
 
448
        #     if energy in checkpt_energies:
 
449
        #         sys.stderr.write ("XXX\n")
 
450
            
 
451
        
 
452
        stot = 0.0
 
453
        for root in roots:
 
454
            xx0 = energy - root[0]
 
455
        
 
456
            if abs(xx0) <= cutoff:
 
457
                stot += root[1] / ( xx0*xx0 + gamma_sqrd)   #Lorentzian
 
458
            
 
459
        yield [energy, stot*prefac]
 
460
 
 
461
 
 
462
 
 
463
        
 
464
 
 
465
def dump_header (opts):
 
466
    """Print header to stdout"""
 
467
 
 
468
    if opts.header:
 
469
        print ("{0} ================================".format(opts.cchar))
 
470
        print ("{0}  NWChem spectrum parser ver {1}".format(opts.cchar,ver))
 
471
        print ("{0} ================================".format(opts.cchar))
 
472
        print ("{0} ".format(opts.cchar))
 
473
        print ("{0} Parser runtime options: {1}".format(opts.cchar,opts))
 
474
        print ("{0}".format(opts.cchar))
 
475
 
 
476
 
 
477
def dump_data (opts, roots):
 
478
    """Dumps output to stdout.  This works for either lists of raw roots
 
479
    or broadened spectra."""
 
480
 
 
481
    if opts.verbose:
 
482
        sys.stderr.write ("{0}: Dumping data to stdout ... \n".format(pname))
 
483
 
 
484
    if opts.units == "ev":
 
485
        if opts.header:
 
486
            print ("{0}".format(opts.cchar))
 
487
            print ("{c}{s1}{delim}{s2}".format(c=opts.cchar, s1="  Energy [eV]  ", delim=opts.delim, s2="   Abs. [au]   "))
 
488
            print ("{0}----------------------------------".format(opts.cchar))
 
489
            
 
490
        for root in roots:
 
491
            print ("{energy:.10e}{delim}{osc:.10e}".format(energy=root[0], delim=opts.delim, osc=root[1]))
 
492
 
 
493
 
 
494
    elif opts.units == "au":
 
495
        if opts.header:
 
496
            print ("{0}".format(opts.cchar))
 
497
            print ("{c}{s1}{delim}{s2}".format(c=opts.cchar, s1="  Energy [au]  ", delim=opts.delim, s2="   Abs. [au]   "))
 
498
            print ("{0}----------------------------------".format(opts.cchar))
 
499
            
 
500
        for root in roots:
 
501
            eout = ev2au (root[0])
 
502
            print ("{energy:.10e}{delim}{osc:.10e}".format(energy=eout, delim=opts.delim, osc=root[1]))
 
503
 
 
504
    elif opts.units == "nm":
 
505
        if opts.header:
 
506
            print ("{0}".format(opts.cchar))
 
507
            print ("{c}{s1}{delim}{s2}".format(c=opts.cchar, s1=" Wavelen. [nm] ", delim=opts.delim, s2="   Abs. [au]   "))
 
508
            print ("{0}----------------------------------".format(opts.cchar))
 
509
            
 
510
#        roots.reverse ()  #reorder so we output in increasing wavelength
 
511
#XXX SHOULD REVERSE            
 
512
        for root in roots:
 
513
            eout = ev2nm (root[0])
 
514
            print ("{energy:.10e}{delim}{osc:.10e}".format(energy=eout, delim=opts.delim, osc=root[1]))
 
515
 
 
516
    else:
 
517
        raise Exception ("Invalid unit: {0}".format(opts.units))        
 
518
        
 
519
 
 
520
def preprocess_check_opts (opts):
 
521
    """Check options and replaces with sane values if needed, stores
 
522
    all opts as purely lowercase."""
 
523
 
 
524
    opts.units = opts.units.lower()
 
525
    opts.datafmt = opts.datafmt.lower()
 
526
 
 
527
    if opts.units != "nm" and opts.units != "ev" and opts.units != "au":
 
528
        raise Exception ("Invalid unit type: {0}".format(opts.units))
 
529
 
 
530
    if opts.datafmt != "tddft" and opts.datafmt != "vspec" and opts.datafmt != "auto" and opts.datafmt != "dos":
 
531
        raise Exception ("Invalid data format: {0}".format(opts.datafmt))
 
532
        
 
533
    if opts.npoints < 100:
 
534
        raise Exception ("Increase number of points to at least 100 (you asked for {0})".format(opts.npoints))
 
535
 
 
536
    if opts.width < 0.0:
 
537
        raise Exception ("Peak width must be positive (you supplied {0})".format(opts.width))
 
538
 
 
539
 
 
540
def main():
 
541
    
 
542
    ##
 
543
    ## Parse command line options.  Note we later make all lowercase,
 
544
    ## so from the user's point of view they are case-insensitive.
 
545
    ##
 
546
    usage = "%prog [options]\n\n"
 
547
    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."
 
548
    desc_wrap = textwrap.wrap (desc,80)
 
549
    for s in desc_wrap:
 
550
        usage += s + "\n"
 
551
 
 
552
    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.'
 
553
    example_wrap = textwrap.wrap (example,80)
 
554
    usage += "\nExample:\n\n\t"+"nw_spectrum -b0.3 -p5000 -wnm < water.nwo > spectrum.dat\n\n"
 
555
    for s in example_wrap:
 
556
        usage += s + "\n"
 
557
 
 
558
 
 
559
    parser = OptionParser(usage=usage)
 
560
    parser.add_option("-f", "--format", type="string", dest="datafmt",
 
561
                      help="data file format: auto (default), tddft, vspec, dos", metavar="FMT")
 
562
    parser.add_option("-b", "--broad", type="float", dest="width", 
 
563
                      help="broaden peaks (FWHM) by WID eV (default 0.1 eV)", metavar="WID")
 
564
    parser.add_option("-n", "--nbin", type="int", dest="nbin",
 
565
                      help="number of eigenvalue bins for DOS calc (default 20)", metavar="NUM")
 
566
    parser.add_option("-p", "--points", type="int", dest="npoints",
 
567
                      help="create a spectrum with NUM points (default 2000)", metavar="NUM")
 
568
    parser.add_option("-w", "--units", type="string", dest="units",
 
569
                      help="units for frequency:  eV (default), au, nm", metavar="UNT")
 
570
    parser.add_option("-d", "--delim", type="string", dest="delim",
 
571
                      help="use STR as output separator (four spaces default)", metavar="STR")
 
572
    parser.add_option("-x", "--extract", action="store_false", dest="makespec",
 
573
                      help="extract unbroadened roots; do not make spectrum")
 
574
    parser.add_option("-C", "--clean", action="store_false", dest="header",
 
575
                      help="clean output; data only, no header or comments")
 
576
    parser.add_option("-c", "--comment", type="string", dest="cchar",
 
577
                      help="comment character for output ('#' default)", metavar="CHA")
 
578
    parser.add_option("-v", "--verbose", action="store_true", dest="verbose",
 
579
                      help="echo warnings and progress to stderr")
 
580
 
 
581
 
 
582
    parser.set_defaults(width = 0.1, npoints = 2000, units="eV", datafmt="auto", delim="    ", makespec=True, header=True, cchar="#", nbin=20, verbose=False)
 
583
    (opts, args) = parser.parse_args()
 
584
 
 
585
    preprocess_check_opts(opts)
 
586
    check_version ()
 
587
    
 
588
    if opts.header:
 
589
        dump_header (opts)
 
590
 
 
591
    if opts.datafmt == "auto":
 
592
        opts.datafmt = determine_data_type ()
 
593
        if opts.datafmt == "tddft":
 
594
            print ("{0} The input appears to contain TDDFT data.".format(opts.cchar))
 
595
        elif opts.datafmt == "vspec":
 
596
            print ("{0} The input appears to contain vspec data.".format(opts.cchar))
 
597
        else:
 
598
            raise Exception ("Invalid data format: {0}".format(opts.datafmt))
 
599
            
 
600
 
 
601
    ## parse raw data
 
602
    if opts.datafmt == "tddft":
 
603
        roots = parse_input_tddft (opts)
 
604
    elif opts.datafmt == "vspec":
 
605
        roots = parse_input_vspec (opts)
 
606
    elif opts.datafmt == "dos":
 
607
        evals = parse_input_evals (opts)
 
608
        roots = bin_evals (opts, evals)
 
609
    else:
 
610
        raise Exception ("Invalid data format supplied: {0}".format(opts.datafmt))
 
611
 
 
612
    
 
613
    ## make broadened spectrum if desired
 
614
    if opts.makespec:
 
615
        energies = make_energy_list (opts, roots)
 
616
        spectrum = gen_spectrum (opts, energies, roots)
 
617
 
 
618
        if opts.verbose:
 
619
            sys.stderr.write ("{0}: Initialized spectrum [{1: .3f} : {2: .3f}] ({3} points)\n".format(pname, energies[0], energies[-1], len(energies)))
 
620
 
 
621
        if opts.header:
 
622
            print ("{0} Roots were Lorentzian broadened with width {1} eV ".format(opts.cchar, opts.width))
 
623
            print ("{0} Spectrum generated with {1} points.".format(opts.cchar, opts.npoints))
 
624
            
 
625
    else:
 
626
        spectrum = roots
 
627
        if opts.header:
 
628
            print ("{0} Roots not broadened".format(opts.cchar))
 
629
            print ("{0} No spectrum generated: output is list of raw excitations".format(opts.cchar))
 
630
 
 
631
 
 
632
    ## finally, dump data to stdout
 
633
    dump_data (opts, spectrum)
 
634
        
 
635
    if opts.verbose:
 
636
        sys.stderr.write ("{0}: Done\n".format(pname))
 
637
 
 
638
 
 
639
if __name__ == "__main__":
 
640
    main()
 
641
 
 
642