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

ARGUMENT

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
}