TABLE OF CONTENTS
::pwtk::dmo::vibSpectrum
SYNOPSIS
proc ::pwtk::dmo::vibSpectrum {args} {
USAGE
vibSpectrum ?-b gauss|lorentz? ?-f FWHM? ?-t TERM? ?-xr XMIN:XMAX? ?-yr YMIN:YMAX? dynmatOutput
PURPOSE
Generate vibrational (IR and Raman) spectra from the dynmat.x output file. Raman spectrum is generated only if the corresponding data exists in the dynmat.x output file.
OPTIONS
- -b TYPE --- type of broadening (TYPE = gauss | lorentz)
- -f FWHM --- the value of FWHM
- -t TERM --- terminal to plot the spectrum to (any terminal returned by ::pwtk::gp::pwtk_terminals, e.g., eps, pdf, png...)
- -xr XMIN:XMAX --- the x-plotting range, specified as FROM:TO
- -yr YMIN:YMAX --- the y-plotting range, specified as FROM:TO
ARGUMENT
- dynmatOutput -- the dynmat.x output file
RETURN VALUE
The "head" of the created data ($head.dat), gnuplot ($head.gp), and image ($head.eps or $head.png ...) files. The "head" is obtained from the name of dynmatOutput.
SOURCE
# parse options set narg 1 set usage "?-b gauss|lorentz? ?-f fwhm? ?-t term? dynmatOutput" set options { {b.arg "gauss" "type of broadening (gauss|lorentz)"} {f.arg 5.0 "FWHM in cm^-1"} {t.arg eps "terminal to plot to"} {xr.arg {} "x range of the plot"} {yr.arg {} "y range of the plot"} } ::pwtk::parseOpt_ ::pwtk::checkOTypeStrict_ -f $opt(f) double number # -f set fwhm $opt(f) # -b switch -glob -- $opt(b) { gauss* - loren* { set broadening $opt(b) } default { ::pwtk::error "wrong broadening \"$opt(b)\", must be gauss or lorentz" 1 } } # -t set term [::pwtk::gp::gp2pwtk [::pwtk::gp::term_usable $opt(t)]] set dmo $args ::pwtk::fileMustExist $dmo "dynmat.x output" set head [regsub {^dynmat.} [file tail [file rootname $dmo]] {}] set nodata 1 set iread 0 set lraman 0 set nf 0 ::pwtk::print "Parsing dynmat.x output file: $dmo\n" # read freqs and IR & Raman intensities pwtk::lineread line $dmo { if { $iread } { ::pwtk::print $line set len [llength $line] if { $len >= 4 } { incr nf set freq($nf) [lindex $line 1] set ir($nf) [lindex $line 3] if { $lraman } { if { $len != 6 } { ::pwtk::error "error reading Raman activities in $dmo" 1 } set ra($nf) [lindex $line 4] } } else { set iread 0 } } elseif { [regexp {^# mode \[cm-1\]} $line] } { set iread 1 set nodata 0 if { [llength $line] > 5 } { set lraman 1 set head IR+Raman.$head } else { set head IR.$head } ::pwtk::print $line } } if { $nodata } { pwtk::warning "dynmat.x output file \"$dmo\" does not contain any IR/Raman data" return } set sigma [expr $fwhm / (2.0*sqrt(2.0*log(2.0)))] set gamma [expr $fwhm / 2.0] set xmax [expr $freq($nf) + min(100,10*$fwhm)] set x [expr max(0.0, $freq(1) - 10*$fwhm)] set dx [expr 0.1*$fwhm] # relative threshold (rth) used for cutting off the smearing, i.e., only # smeared frequencies that are close enough to a given wavenumber are taken into account set rth 1e-3 set fid [open $head.dat w] if { $broadening == "gauss" } { set str [format "(Gaussian), sigma = %.3f\n#" $sigma] set norm [expr 1.0/($sigma*sqrt($pwtk::tpi))] set xoff [expr $sigma * sqrt(-2*log($rth))] } else { set str [format "(Lorentzian), gamma = %.3f\n#" $gamma] set norm [expr 1.0/($pwtk::pi*$gamma)] set rth [expr 2*$rth]; # Lorentzian smearing is "broad at the bottom", reduce the threshold set xoff [expr $gamma*sqrt(1.0/$rth - 1)] } if { ! $lraman } { puts $fid "# smeared IR spectrum: FWHM = $fwhm $str" puts $fid "# freq(cm^-1) IR.activity" } else { puts $fid "# smeared IR & Raman spectra: FWHM = $fwhm $str" puts $fid "# freq(cm^-1) IR.activity Raman.activity" } # calculate the spectra ::pwtk::print "Processing ..." set imin 1 while 1 { set f_ir 0.0 set f_ra 0.0 for {set i $imin} {$i <= $nf} {incr i} { # take into account only frequencies that are close enough (+/-$xoff) to $x if { $freq($i) < $x - $xoff } { set imin $i } elseif { $freq($i) > $x + $xoff } { break } if { $broadening == "gauss" } { set f_ir [expr $f_ir + $ir($i)*exp(-0.5*(($x-$freq($i))/$sigma)**2)] if { $lraman } { set f_ra [expr $f_ra + $ra($i)*exp(-0.5*(($x-$freq($i))/$sigma)**2)] } } else { set f_ir [expr $f_ir + $ir($i)/(1.0 + (($x-$freq($i))/$gamma)**2)] if { $lraman } { set f_ra [expr $f_ra + $ra($i)/(1.0 + (($x-$freq($i))/$gamma)**2)] } } } if { ! $lraman } { puts $fid [format "%10.2f %12.4e" $x [expr $norm*$f_ir]] } else { puts $fid [format "%10.2f %12.4e %12.4e" $x [expr $norm*$f_ir] [expr $norm*$f_ra]] } set x [expr $x + $dx] if { $x > $xmax } break } close $fid # plot the spectra set p {} set p [expr { $lraman ? "-p" : "" }] set gp [::pwtk::gp::plot new {*}$p $head.$term] $gp options { xlabel {'Wavenumber / cm^{-1}'} ylabel {'IR intensity (arb.u.)' textcolor lt 2} key {top left} format.y '' format.y2 '' } if { $opt(xr) ne {} } { $gp options "xrange {\[$opt(xr)\]}" } if { $opt(yr) ne {} } { $gp options "yrange {\[$opt(yr)\]}\ny2range {\[$opt(yr)\]}" } else { $gp options "yrange {\[0:\]}\ny2range {\[0:\]}" } $gp unset grid ytics y2tics if { $lraman } { $gp options {y2label {'Raman activity (arb.u.)' textcolor lt 1 offset 1}} $gp plot "'$head.dat' u 1:2 axis x1y1 not w filledcurve y1=0.0 lt 2, \\ '$head.dat' u 1:3 axis x1y2 not w filledcurve y2=0.0 lt 1, \\ \\ '$head.dat' u 1:2 axis x1y1 title ' IR' w l lt 2 lw 2, \\ '$head.dat' u 1:3 axis x1y2 title ' Raman' w l lt 1 lw 2" $gp new_page $gp add "unset y2label" $gp plot "'$head.dat' u 1:2 axis x1y1 not w filledcurve y1=0.0 lt 2, \\ '$head.dat' u 1:2 axis x1y1 title ' IR' w l lt 2 lw 2" $gp new_page $gp options {ylabel {'Raman activity (arb.u.)' textcolor lt 1}} $gp plot "'$head.dat' u 1:3 axis x1y2 not w filledcurve y2=0.0 lt 1, \\ '$head.dat' u 1:3 axis x1y1 title ' Raman' w l lt 1 lw 2 " ::pwtk::print "Data for IR & Raman spectra written to : $head.dat" } else { $gp plot "'$head.dat' u 1:2 axis x1y1 not w filledcurve y1=0.0 lt 2, \\ '$head.dat' u 1:2 axis x1y1 title ' IR' w l lt 2 lw 2" ::pwtk::print "Data for IR spectrum written to : $head.dat" } $gp exec $gp destroy return $head }