TABLE OF CONTENTS

::pwtk::eos::findMinimum

SYNOPSIS

proc ::pwtk::eos::findMinimum {head datafile} {

PURPOSE

Read a-vs-E datafile and find an optimum a0, E0, and B0 (i.e. lattice-parameter, total-energy, and bulk-modulus). The datafile should consits of two-columns, the first is the lattice-parameter [in bohr], the second is the correspponding total-energy [in Ryd]. Comment lines starting with "#" are allowed.

ARGUMENTS

SOURCE

    # read points from datafile

    set ibrav [::pwtk::input::namelistGetVarValue SYSTEM ibrav]
    
    set npoints 0
    ::pwtk::lineread line $datafile {   

        if { [llength $line] >= 2 && ![regexp "^ *#" $line] } {     

            incr npoints

            if { abs($ibrav) > 0 && abs($ibrav) < 4 } {

                set a($npoints) [lindex $line 0]
                if { ![string is double $a($npoints)] } {
                    ::pwtk::error "syntax error reading $datafile"
                }
                set V($npoints) [getV0fromA0_ $a($npoints)]

            } elseif { $ibrav == 4 } {

                set V($npoints) [lindex $line 0]
                if { ![string is double $V($npoints)] } {
                    ::pwtk::error "syntax error reading $datafile"
                }
                set a($npoints) [getA0fromV0_ $V($npoints)]
                
            } else {
                ::pwtk::error "only ibrav in range of [1,4] supported"
            }
            set E($npoints) [lindex $line 1]
            if { ![string is double $E($npoints)] } {
                ::pwtk::error "syntax error reading $datafile"
            }
        }
    }

    if { $npoints < 3 } {
        ::pwtk::error "can't find minimum for less than three data points"
    }
    if { $npoints < 7 } {
        ::pwtk::warning "number of data points < 7"
    }

    # map between ibrav & ev.x lattice

    set lattice [::pwtk::eos::latticeType_]

    # map between id & equation-of-state for ev.x

    array set emap {
        1 birch1
        2 birch2
        3 keane
        4 murnaghan
    }
    
    # find the lowest-energy point

    set mini 1
    set E(mini) $E(1)
    for {set i 1} {$i <= $npoints} {incr i} {
        if { $E($i) < $E(mini) } { 
            set mini $i 
            set E(mini) $E($i)
        }
    }

    # calculate a0 and B0 for different number of points ...
    # TODO: this should be adopted and symmetric range of points (wrt minimum) determined !

    set done 0
    set first 1
    set last  $npoints
    set ii 1

    while { ! $done } {

        set npoi [expr {$last - $first + 1}]

        # write datafile with $npoi points
        
        set dataf ${head}.$first-$last.dat
        set dat [open $dataf w]     
        for {set i $first} {$i <= $last} {incr i} {         
            if { abs($ibrav) < 4 } {
                puts $dat " $a($i)   $E($i) "
            } else {
                puts $dat " $V($i)   $E($i) "
            }
        }
        close $dat

        # calculate a0, E0, and B0 for the current number of points,
        # $npoi, using all four equation of states supported by ev.x

        for {set eqs 1} {$eqs <= 4} {incr eqs} {
            
            set prefix ${head}.$emap($eqs).$first-$last
            set inp [open $prefix.ev.in w]
            
            if { [file exists $prefix.ev_out] } {
                file delete -- $prefix.ev_out
            }

            puts $inp "au\n$lattice\n$eqs\n$dataf\n$prefix.ev_out"
            close $inp

            # BEWARE: ev.x seems not to understand the "-in" command line option, use stdin redirection instead
            ::pwtk::run -serial -ihandle "<" ev.x $prefix.ev.in $prefix.ev.out
        
            # extract the a0 && E0 && B0

            set out [open $prefix.ev_out r]
            gets $out line
            gets $out line
            set line [regsub -all = $line {}]

            if { $ibrav != 4 } {
                # cubic systems
                set a0($npoi,$eqs) [lindex $line 2]
            } elseif { $ibrav == 4 } {
                # for hcp ev.x reports V0
                set a0($npoi,$eqs) [:::pwtk::eos::hexGetA0fromV0_ [lindex $line 2]]
            }
            set B0($npoi,$eqs) [lindex $line 5]
            set E0($npoi,$eqs) [lindex $line end]
            close $out
        }

        set im($npoi) $first
        set iM($npoi) $last
        
        # decrease the number of points 

        if { $ii % 2 } {
            incr last -1
        } else {
            incr first 
        }
        
        if { ! ($first < $mini && $last > $mini) || $npoi <= 7 } {
            set done 1
        }

        incr ii
    }  

    # print the result

    set res [open $head.RESULTS w]

    puts $res "\nFINAL RESULTS::"
    puts $res "---------------\n"
    puts $res "  Lattice type:  $lattice\n"
    puts $res "  OPTIMISED LATTICE PARAMETER, a0 (in Bohr units)::\n"
    puts $res "  ------------+------------+------------+------------+-----------"
    puts $res "   # points   | MURNAGHAN  | BIRCH O(1) | BIRCH O(2) |    KEANE  "
    puts $res "  ------------+------------+------------+------------+-----------"     
    for {set np $npoints} {$np >= $npoi} {incr np -1} { 
        puts $res [format "   %2i: %2i--%2i |  %8.3f  |  %8.3f  |  %8.3f  |  %8.3f " \
                       $np $im($np) $iM($np) \
                       $a0($np,4) $a0($np,1) $a0($np,2) $a0($np,3)]
    }
    puts $res "  ------------+------------+------------+------------+-----------\n" 


    puts $res "  OPTIMISED Total Energy, E0 (in Ryd units)::\n"
    puts $res "  ------------+------------+------------+------------+-----------"
    puts $res "   # points   | MURNAGHAN  | BIRCH O(1) | BIRCH O(2) |    KEANE  "
    puts $res "  ------------+------------+------------+------------+-----------"     
    for {set np $npoints} {$np >= $npoi} {incr np -1} { 
        puts $res [format "   %2i: %2i--%2i | %10.4f | %10.4f | %10.4f | %10.4f " \
                       $np $im($np) $iM($np) \
                       $E0($np,4) $E0($np,1) $E0($np,2) $E0($np,3)]
    }
    puts $res "  ------------+------------+------------+------------+-----------\n" 


    puts $res "  OPTIMISED Bulk Modulus, B0 (in kbar units)::\n"
    puts $res "  ------------+------------+------------+------------+-----------"
    puts $res "   # points   | MURNAGHAN  | BIRCH O(1) | BIRCH O(2) |    KEANE  "
    puts $res "  ------------+------------+------------+------------+-----------"     
    for {set np $npoints} {$np >= $npoi} {incr np -1} { 
        puts $res [format "   %2i: %2i--%2i |   %5i    |   %5i    |   %5i    |   %5i " \
                       $np $im($np) $iM($np) \
                       $B0($np,4) $B0($np,1) $B0($np,2) $B0($np,3)]
    }
    puts $res "  ------------+------------+------------+------------+-----------\n" 
    close $res

    # print also to stdout

    set res [open $head.RESULTS r]
    ::pwtk::print [read $res]
    ::pwtk::print "\nResult written to file $head.RESULTS\n"
    close $res    
}