TABLE OF CONTENTS

::pwtk::eos::findMinimum

SYNOPSIS

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

PURPOSE

Depending on the latticeType (sc, bcc, fcc, hex, or noncubic), read either a-vs-E (cubic systems) or V-vs-E (noncubic systems) datafile and find optimum a0, E0, and B0 (i.e. lattice parameter, total energy, and bulk modulus). The datafile should consist of two-columns, the first is either the lattice parameter [in bohr] or the cell volume [in bohr^3], the second is the correspponding total energy [in Ryd].

Comment lines starting with "#" are allowed in the datafile.

This command is used by the EOS workflow.

ARGUMENTS

SOURCE

    # read points from datafile

    set ibrav [::pwtk::pwi::ibrav]
    
    set npoints 0
    ::pwtk::lineread line $datafile {   

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

            incr npoints

            if { $latticeType eq "sc" || $latticeType eq "bcc" || $latticeType eq "fcc" } {

                set a($npoints) [lindex $line 0]
                if { ![string is double $a($npoints)] } {
                    ::pwtk::error "syntax error reading $datafile: expected a real-number for A, but got $a($npoints)" 1
                }
                set V($npoints) [getV0fromA0_ $a($npoints)]

            } else {

                set V($npoints) [lindex $line 0]
                if { ![string is double $V($npoints)] } {
                    ::pwtk::error "syntax error reading $datafile: expected a real-number for V, but got $V($npoints)" 1                }
                set a($npoints) [getA0fromV0_ $V($npoints)]
                
            }
            
            set E($npoints) [lindex $line 1]
            if { ![string is double $E($npoints)] } {
                ::pwtk::error "syntax error reading $datafile: expected a real-number for E, but got $E($npoints)" 1
            }
        }
    }

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

    # 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) > 0 && 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
            
            if { [file exists $prefix.ev_out] } {
                file delete -- $prefix.ev_out
            }

            ::pwtk::writeFile $prefix.ev.in "au\n$latticeType\n$eqs\n$dataf\n$prefix.ev_out"

            # 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
            close $out
            set line [regsub -all = $line {}]

            if { abs($ibrav) > 0 && abs($ibrav) < 4 } {
                # cubic systems
                set a0($npoi,$eqs) [lindex $line 2]
            } else {
                # for hcp & noncubic ev.x reports V0
                set a0($npoi,$eqs) [getA0fromV0_ [lindex $line 2]]
            }
            set B0($npoi,$eqs) [lindex $line 5]
            set E0($npoi,$eqs) [lindex $line end]
        }

        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 "
FINAL RESULTS::
---------------

  Lattice type:  $latticeType

  OPTIMISED LATTICE PARAMETER, a0 (in Bohr units)::

  ------------+------------+------------+------------+-----------
   # points   | MURNAGHAN  | BIRCH O(1) | BIRCH O(2) |    KEANE  
  ------------+------------+------------+------------+-----------"
    
    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)::

  ------------+------------+------------+------------+-----------
   # points   | MURNAGHAN  | BIRCH O(1) | BIRCH O(2) |    KEANE  
  ------------+------------+------------+------------+-----------"
    
    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)::

  ------------+------------+------------+------------+-----------
   # points   | MURNAGHAN  | BIRCH O(1) | BIRCH O(2) |    KEANE  
  ------------+------------+------------+------------+-----------"
    
    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

    ::pwtk::print [pwtk::readFile -nonewline $head.RESULTS]
    ::pwtk::print "Result written to file $head.RESULTS\n\n"
}