lobster.tcl

::pwtk::lobster::sumdos

::pwtk::lobster::sumdos

SYNOPSIS

proc ::pwtk::lobster::sumdos {atomList {doscar DOSCAR} } {

USAGE

   ::pwtk::lobster::sumdos  atomList  ?doscar?

ARGUMENTS

atomList ... list of atoms to sum, specified in the ::pwtk::parseRangeString syntax (example: 1-4,5,7)
doscar ... (optional, default = DOSCAR) lobster DOSCAR file, e.g., DOSCAR.lobster, DOSCAR.LSO.lobster, DOSCAR.LCFO.lobster; can be abbreviated, e.g., DOS stands for DOSCAR.lobster, DOS.LCFO stands for DOSCAR.LCFO.lobster, etc.

PURPOSE

Sum atomic/fragment PDOSes. If compatible PDOSes are summed, both atomic/fragment PDOS and obital-projected PDOSes will be summed, otherwise only atomic/fragment PDOSes are summed. In the first case, the output filename has the orbatm*.dat (or orgfrag*.dat) suffix, wheras in the second case

RETURN VALUE

The name of the output file containing summed DOS, which has:

.orbatom*.dat or .orbfrag*.dat suffix (* stands for the compressed atom list) when both atomic/fragment and obital-projected PDOSes were summed
.atom*.dat or .frag*.dat suffix when only atomic/fragment PDOSes were summed

SOURCE

    set doscar [::pwtk::lobster::lobsterfile_ $doscar]    
    set head   [file rootname $doscar]
    set lcfo   [expr { [regexp LCFO $head] ? 1 : 0 }]
    
    array set info [::pwtk::lobster::DOSCAR_info_ $doscar]
    if { $lcfo } {
        set type frag
        set types fragments
        # for LCFO, nat == number of fragments
        set info(nat) $info(nproj)
    } else {
        set type atom
        set types atoms
    }
        
    set atmL  [join [regsub -all , [regsub -all { *\- *} $atomList -] { }] ,]
    set atmL  [regsub {\-$} $atmL -$info(nat)]    
    set atomList [::pwtk::parseRangeString $atomList $info(nat)]

    if { [llength $atomList] == 0 } {
        # DOSCAR.atom.dat
        return $head.$type.dat
    } elseif { [llength $atomList] == 1 } {
        # DOSCAR.orbatomX.dat
        return $head.orb${type}$atomList.dat        
    }

    array set info [::pwtk::lobster::DOSCAR_dat_ $doscar]
    
    # check if obitals can be summed
    
    set orb 1
    set out $head.orb${type}$atmL.dat

    foreach ia $atomList {
        set ind([incr np]) $ia
    }
    for {set ip 1} {$ip < $np} {incr ip} {
        set ia  $ind($ip)
        set ia1 $ind([expr $ip + 1])
        if { $info(proj$ia) ne $info(proj$ia1) } {
            set orb 0
            set out $head.${type}$atmL.dat
            break
        }
    }

    # atom No. of the 1st DOSCAR.orbatomX.dat file
    set ia0 [lindex $atomList 0]
    
    foreach ia $atomList {
        set dat $head.orb${type}$ia.dat
        ::pwtk::fileMustExist $dat "${type}X datafied LOBSTER"
        
        set fid($ia) [open $dat r]        
        gets $fid($ia); # skip the 1st line
    }
    
    set oid [open $out w]
    
    if { $orb } {
        puts $oid "# Summed $type PDOS with orbital components for $types [join $atomList ,]"
        puts $oid "# orbital components: $info(proj1)"
        # other comment line are printed below
    } else {
        puts $oid "# Summed $type PDOS for $types [join $atomList ,]
# nspin = $nspin
# # iX indices of DOS components:
# i0 = sum"
        if { $nspin == 1 } {
            puts $oid "# energy sum"
        } else {
            puts $oid "# energy sum(up) sum(dw)"
        }
    }

    # read & write the comment lines
        
    foreach ia $atomList {
        gets $fid($ia) line
        
        while { [regexp {^ *#} $line] } {
            if { $ia == $ia0 && $orb } {
                puts $oid $line
            }
            set pos [tell $fid($ia)]
            gets $fid($ia) line
        }
        seek $fid($ia) $pos
        
        # at this point $line contains the 1st DOS-data line
        if { $ia == $ia0 } {
            set len [llength $line]
        }
    }

    if { ! $orb } {
        # set 'len' to only summed atomic DOS
        set len [expr { $nspin == 1 ? 2 : 3 }]
    }

    # read & process xDOS data

    for {set ip 1} {$ip <= $info(npoints)} {incr ip} {

        # read 

        foreach ia $atomList {
            set i 0
            foreach f [gets $fid($ia)] {
                set d($ia,$i) $f
                incr i
            }
        }    

        # energy
        puts -nonewline $oid [format %9.5f $d($ia0,0)]

        # xDOS        
        for {set i 1} {$i < $len} {incr i} {
            set pdos 0.0
            foreach ia $atomList {
                set pdos [expr $pdos + $d($ia,$i)]
            }
            puts -nonewline $oid [format " %.5f" $pdos]
        }
        puts $oid ""
    }

    # close files

    close $oid
    foreach ia $atomList { close $fid($ia) }

    return $out
}

TABLE OF CONTENTS

::pwtk::lobster::sumdos