pseudo.tcl

::pwtk::pseudo::etot_test

::pwtk::pseudo::etot_test

SYNOPSIS

proc ::pwtk::pseudo::etot_test {args} {

USAGE

   ::pwtk::pseudo::etot_test [options] pseudopotentialFile

DESCRIPTION

Test a pseudopotential for the minimum ecutwfc and, for US & PAW, also ecutrho, using the criterium based on total-energy. If the -d option is explicitly specified, then ecutrho is tested also for a NC pseudopotential.

OPTIONS

In most case no options are needed, hence use them only if the default setup fails.

    -e ECUTWFC_LIST --- (default = {30 60 90 120}) list of ecutwfc values for initial scanning,
                         minimum of 4 values needed
    -d DUAL         --- (default = 16) initial value of dual for US and PAW potentials
    -thr THRESHOLD  --- (default = 1e-3) ecutwfc convergence treshold for total energy (in Ry)
    -u UPSCALE      --- (default = 4) upscale factor for the ecutrho threshold, i.e.,
                        rhothr = wfcthr/upscale
    -n NCALCS       --- (default = 6) number of the pw.x SCF calcs during the refinement stage
    -l LOW_FACTOR   --- (default = 0.8) "low" factor for ecutwfc for the refinement stage, i.e.,
                        ecutwfc_min = low_factor*ecutwfc
    -h HIGH_FACTOR  --- (default = 2.0) "high" factor for ecutwfc for the refinement stage, i.e.,
                        ecutwfc_max = high_factor*ecutwfc
    -m MIN_ECUTWFC  --- (default = 10) minimum allowed ecutwfc (in Ry)
    -M MAX_ECUTWFC  --- (default = 1000) maximum allowed ecutwfc (in Ry)
    -t TERM         --- a Gnuplot file terminal to plot to (one of eps, pdf, png, svg)
    -v              --- visualize the resulting plot

ARGUMENTS

pseudopotentialFile -- name of the PP file

RETURN VALUE

The determined ecutwfc and ecutrho values as the "$ecutwfc $ecutrho" list.

SOURCE

    variable E
    variable gpData
    variable opts
    
    ::pwtk::printTitle PSEUDO "Automatic search of \"ecutwfc\" & \"ecutrho\" for a pseudo-potential\n(based on simple total-energy criterium)"

    if { [auto_execok gnuplot] == {} } {
        ::pwtk::error "[::pwtk::procName] requires gnuplot, but gnuplot was not found" 1
    }

    # just in case
    catch { unset E }
    
    ::pwtk::input::push; # do not move this push down the code, it must be called at the beginning

    # ------------------------------------------------------------------------
    #
    # 0. option-processing & init setup
    #
    # ------------------------------------------------------------------------
    
    set narg 1
    set usage "?-e ECUTWFC_LIST?  ?-d DUAL?  ?-t THRESHOLD?  ?-u UPSCALE?  ?-l LOW_FACTOR?  ?-h HIGH_FACTOR?  ?-m MIN_ECUTWFC?  ?-M MAX_ECUTWFC?  ?-v?  pseudopotentialFile"
    set options [subst {
        {e.arg   {}  "list of ecutwfc value for initial scanning (minimum of 4 value needed; default = 30 60 90 120)"}
        {d.arg   $opts(dual)    "value of dual for US and PAW potentials"}
        {thr.arg $opts(thr)     "ecutwfc convergence treshold for total energy (in Ry)"}
        {u.arg   $opts(upscale) "upscale factor for ecutrho threshold, i.e., rhothr = thr/upscale"}
        {n.arg   $opts(ncalc)   "number of pw.x calcs during the refinement stage"}
        {l.arg   $opts(low_factor)  "low factor of ecutwfc for the refinement stage"}
        {h.arg   $opts(high_factor) "high factor of ecutwfc for the refinement stage"}
        {m.arg   $opts(min_ecutwfc) "minimum allowed ecutwfc (if ecutwfc gets lower, drop error)"}
        {M.arg   $opts(max_ecutwfc) "maximum allowed ecutwfc (if ecutwfc gets larger, drop error)"}
        {t.arg   $::pwtk::gp::gp(workflow.terminal)  "a Gnuplot file terminal to plot to"}
        {v       "visualize the resulting plot"}
    }]
    
    ::pwtk::parseOpt_
    ::pwtk::checkOType_ -e $opt(e) {numberlist double}  "list of numbers"
    ::pwtk::checkOTypeStrict_ -d $opt(d) double number; # if dual < 0 ==> dual tests performed only for US & PAW pseudos
    ::pwtk::checkOTypeStrict_ -thr $opt(thr) {number posreal} "positive number"
    ::pwtk::checkOTypeStrict_ -u $opt(u) {number posreal} "positive number"
    ::pwtk::checkOTypeStrict_ -n $opt(n) {number posint}  "positive integer"
    ::pwtk::checkOTypeStrict_ -l $opt(l) {number posreal} "positive number"
    ::pwtk::checkOTypeStrict_ -h $opt(h) {number posreal} "positive number"
    ::pwtk::checkOTypeStrict_ -M $opt(M) {number posreal} "positive number"
    ::pwtk::checkOTypeStrict_ -m $opt(m) {number posreal} "positive number"
    ::pwtk::checkOTypeStrict_ -t $opt(t) {optionlist eps pdf png svg}

    set TERM   [string toupper $opt(t)]
    set ecutL  [lsort -real [lmap e [regsub -all , $opt(e) { }] { expr round($e) }]]
    set dual   $opt(d)
    set thr    $opt(thr)
    set rhothr [format %.1e [expr $opt(thr)/$opt(u)]]
    set thr_orig $thr

    if { $dual < 0 } {
        if { [needDual $args] } {
            # use the default daul 
            set dual [expr abs($dual)]        
        } else {
            # no dual tests for norm-conserving PP
            set dual 4
        }
    }

    # get info from the PP filename

    set ppFile $args
    set name [file rootname [file tail $ppFile]]
    set elem [elemFromPPName $name]
    set type [typeFromPPName $name]

    set suggested {}
    if { $ecutL eq {} } {
        # try to read suggested minimum cutoff for wavefunctions from PP file
        set e_ [ecutwfcFromPPFile $ppFile]
        if { $e_ ne {} } {
            set ecutL  [::pwtk::doubleSeq $opt(l)*$e_  $opt(h)*$e_  4 %.1f]
            set suggested "(the suggested minimum ecutwfc = $e_ Ry from the PP file was used)\n"
        } else {
            # use the defaults
            set ecutL $opts(ecutwfc_list)
        }
    }
    set eFirst [lindex $ecutL 0]
    set eLast  [lindex $ecutL end]

    # setup pw.x input
    
    minimalSetup $elem $ppFile
    
    # printout
    
    ::pwtk::print PSEUDO "Basic info:
   pseudopotential file     :  $ppFile
   pseudopotential type     :  $type
   chemical element         :  $elem
   total energy threshold   :  $thr Ry (for ecutwfc),  $rhothr Ry (for ecutrho)
   initial dual for ecutrho :  $dual
   initial ecutwfc scan     :  [join $ecutL {, }] Ry
   $suggested   
   Using ...          ibrav :  [::pwtk::pwi::ibrav]
                       alat :  [::pwtk::pwi::alat]
"

    # defs for running pw.x calcs and storing the results
    
    set runPW_ {
        SYSTEM [subst {
            ecutwfc = $ecut
            ecutrho = $dual*$ecut
        }]
        set scf scf.$name.ecut_$ecut.dual_$dual
        if { [info exists E($ecut,$dual)] } {
            ::pwtk::print "The pw.x calculation \"$scf.out\" already performed.\n"
        } else {          
            ::pwtk::runPW $scf
        }
        set E($ecut,$dual) [::pwtk::pwo::totene $scf.out]
    }
    
    set buildDat_ {
        set dat "# (dual = $dual)
# ecutwfc   E(total)  (in Ry)\n"
        foreach ecut $ecutL {
            append dat [format "   %6.1f %12.6f\n" $ecut $E($ecut,$dual)]
        }
        ::pwtk::print $dat
        ::pwtk::writeFile $head.dat $dat
    }

    # ------------------------------------------------------------------------
    #
    # 1. initial SCAN of the ecut estimate
    #
    # ------------------------------------------------------------------------

    set head ecutwfc.init.$name    
    set stage_ init
    set term {}
    set plot {}
    set param [::fileutil::tempfile param.]; # file where gnuplot stores the parameters
    
    foreach ecut $ecutL {
        eval $runPW_
    }
    
    # if total-energy difference is too low/high, perform one more calculation
    
    if { [expr abs($E($eLast,$dual) - $E($eFirst,$dual))] < 2*$thr } {
        # too-low
        set ecut [expr round(0.67*$eFirst)]        
        set ecutL [linsert $ecutL 0 [set eFirst $ecut]]

        ::pwtk::print PSEUDO "Total-energy difference between min and max tested ecutwfc is too low.
Performing additional calculation for ecutwfc = $ecut Ry
"        
        eval $runPW_
        
    } elseif { [expr abs($E($eLast,$dual) - $E([lindex $ecutL end-1],$dual))] >= $thr } {
        # too-high
        set ecut [expr $eLast + 30]
        set eFirst [lindex $ecutL 1]; # drop the 1st point (take the 2nd point as the first)
        set eLast  $ecut
        lappend ecutL $ecut

        ::pwtk::print PSEUDO "Total-energy difference between the last two points is too high.
Performing additional calculation for ecutwfc = $ecut Ry
"
        eval $runPW_
    }

    set Ntrials 3;  # maximum number of trials
    for {set i 1} {$i <= $Ntrials} {incr i} {

        eval $buildDat_
        ::pwtk::print "Initial energy datafile written to:   $head.dat"

        # initial Gnuplot fitting, f(ecut) = E0 + a*exp(-b*ecut)
        # parameters E0, a, b will be estimated and "loaded" into PWTK

        gpFit_
    
        # thr = a*exp(-b*ecutwfc) --> ecutwfc = -1/b ln(thr/a)
    
        set ecutwfc [expr round(-1.0/$b * log ($thr/$a))]
        ::pwtk::print PSEUDO "Initial estimate of ecutwfc = $ecutwfc Ry  (attempt #.$i)"

        # check that ecutwfc is within a reasonable range

        if { $ecutwfc > 0.0 && $ecutwfc <= $opt(M) } {
            # success !
            break
        } else {
            # perform another calculation or give-up

            if { $ecutwfc <= 0.0 } {
                set str "estimated ecutwfc <= 0"
            } elseif { $ecutwfc > $opt(M) } {
                set str "estimated ecutwfc is too high"
            }

            if { $i < $Ntrials } {
                ::pwtk::warning "The estimation of ecutwfc failed ($str)."
                ::pwtk::print PSEUDO "Expanding initial ecutwfc scan..."

                if { $eFirst >= $opt(m) && [expr abs($E($eLast,$dual) - $E($eFirst,$dual))] < 2*$thr } {
                    # super-soft && $eFirst high-enough
                    set ecut  [expr round(0.67*$eFirst)]
                    set ecutL [linsert $ecutL 0 [set eFirst $ecut]]
                } elseif { $ecutwfc < $opt(M) } {
                    # super-soft && $eFirst too-low
                    set ecut [expr $eLast + 30]
                    lappend ecutL [set eLast $ecut]
                } else {
                    # super-hard case
                    set ecut [expr $eLast + 100]
                    lappend ecutL [set eLast $ecut]
                }
                ::pwtk::print PSEUDO "Performing additional calculation for ecutwfc = $ecut Ry\n"
                eval $runPW_
            } elseif { $ecutwfc > $opt(M) } {
                # maximum number of trials exceeded & still $ecutwfc > $opt(M);
                # give-up
                clean_
                ::pwtk::error "The estimation of ecutwfc failed ($str)." 1
            }
        }
    }
    
    # ------------------------------------------------------------------------
    #
    # 2. Refined estimation of the ecut
    #
    # ------------------------------------------------------------------------

    set head ecutwfc.$name
    
    if { $ecutwfc < $opt(m) } {
        set thr [expr 0.5*$thr]; # new
        ::pwtk::warning "The estimated ecutwfc ($ecutwfc Ry) is too low. Using $opt(m) Ry instead.
The ecutwfc threshold was made tighter, $thr Ry."
        set ecutwfc $opt(m)
    }

    set stage_ refine        
    set eFirst [expr round($opt(l)*$ecutwfc)]    
    set dEmin 5; # dE should not be too low (min = $dEmin Ry)
    set dE [expr max($dEmin, round(($opt(h)*$ecutwfc - $eFirst)/($opt(n)-1)))]
    unset ecutL 
    for {set i 0} {$i < $opt(n)} {incr i} {
        set ecut [expr $eFirst + $i*$dE]
        lappend ecutL $ecut
    }
    set eLast $ecut

    ::pwtk::print PSEUDO "Performing refining calculations at ecutwfc = [join $ecutL {, }] Ry\n"

    foreach ecut $ecutL {
        eval $runPW_
    }
    eval $buildDat_
    ::pwtk::print "Refined energy datafile written to:   $head.dat"

    # refined gnuplot fitting & plotting
    
    gpPlot_
    gpFit_
    
    set ecutwfc [expr round(-1.0/$b * log ($thr/$a))]
    ::pwtk::print PSEUDO "Refined estimate of ecutwfc = $ecutwfc Ry\n"

    set super_soft 0
    if { $ecutwfc < $opt(m) } {        
        ::pwtk::warning "Pseudopotential \"$ppFile\" appears super soft, ecutwfc < $opt(m) Ry
Using $opt(m) Ry instead."
        set ecutwfc $opt(m)
        set super_soft 1

    } elseif { $ecutwfc > $opt(M) } {
        set posErr "The estimation of ecutwfc failed (predicted ecutwfc is too high).
Maximum allowed ecutwfc = $opt(M) Ry. It can be increased with the -M option"        
        ::pwtk::warning $posErr

        # make a new calc with ecutwfc += 100 Ry
        
        lappend ecutL [set ecut [incr eLast 100]]
        ::pwtk::print PSEUDO "Performing additional pw.x calculation at ecutwfc = $ecut Ry\n"
        eval $runPW_
        eval $buildDat_

        # refined gnuplot fitting & plotting
        
        gpPlot_
        gpFit_
        
        set ecutwfc [expr round(-1.0/$b * log ($thr/$a))]
        ::pwtk::print PSEUDO "2nd refined estimate of ecutwfc = $ecutwfc Ry\n"
        
        if { $ecutwfc > $opt(M) } {
            # still too-high, give-up
            clean_
            ::pwtk::error $posErr 1
        }
    }
    
    set txt "  ecutwfc  =  [format %3d $ecutwfc] Ry (for threshold = $thr Ry)"
    if { $super_soft } {
        ::pwtk::print $txt\n        
    } else {
        set E2 [expr round(-1.0/$b * log (2*$thr/$a))]
        set E3 [expr round(-1.0/$b * log (3*$thr/$a))]
        ::pwtk::print "$txt
         ( =  [format %3d $E2] Ry (for threshold = 2*$thr Ry) )
         ( =  [format %3d $E3] Ry (for threshold = 3*$thr Ry) )\n"
    }
    
    # for norm-conserving PP (dual = 4)
    set ecutrho [expr $dual*$ecutwfc]
    set ecutrho_dat N/A
    set ecutrho_thr N/A
    
    # ------------------------------------------------------------------------
    #
    # 3. check if dual can be lowered
    #
    # ------------------------------------------------------------------------
    
    if { $dual > 4 } {
        set dHigh $dual
        set dLow  4
        set ecut  $ecutwfc; # BEWARE: this is needed by $runPW_ !
        for {set i $dual} { $i >= $dLow } {incr i -2} {
            lappend dualL $i
        }
        ::pwtk::print PSEUDO "Checking if dual can be decreased...\n"
        ::pwtk::print "  Performing calculations for ecutwfc = $ecut Ry
                             and dual = [join $dualL {, }] 
"
        foreach dual $dualL {
            eval $runPW_
        }
        
        # write the ecutrho datafile
        
        set dat "# (data for ecutwfc = $ecut Ry)
# dual   E(total)      ecutrho  (in Ry)\n"
        foreach dual $dualL {
            append dat [format "   %3d %12.6f   %6.1f\n" $dual $E($ecut,$dual) [expr $dual*$ecut]]
        }
        ::pwtk::print $dat
        ::pwtk::writeFile ecutrho.$name.dat $dat

        # get a minimum dual: be $opt(u)-times more tight for the dual threshold

        set N [expr int($opt(u)) + 1]        
        for {set i 1} {$i <= $N} {incr i} {            
            set dualMin($i) $dLow      
            set t($i) [format %.1e [expr double($i*$thr_orig)/$opt(u)]]
            foreach dual $dualL {
                if { abs($E($ecut,$dual) - $E($ecut,$dHigh)) <= $t($i) } {
                    set dualMin($i) $dual
                } else {
                    break
                }
            }
        }

        # print-out
        ::pwtk::print PSEUDO "Estimated dual based on total-energy threshold:"
        ::pwtk::print "  dual = [format %2d $dualMin(1)], ecutrho = [format %3d [expr $ecut*$dualMin(1)]] Ry (for threshold = $t(1) Ry)"                
        for {set i 2} {$i <= $N} {incr i} {
            set i1 [expr $i-1]
            if { $dualMin($i1) != $dualMin($i) } {
                ::pwtk::print "     ( = [format %2d $dualMin($i)],         = [format %3d [expr $ecut*$dualMin($i)]] Ry (for threshold = $t($i) Ry) )"
            }
        }                
        puts ""
        
        set dual    $dualMin(1)
        set ecutrho [expr $dual*$ecutwfc]
        set ecutrho_dat ecutrho.$name.dat
        set ecutrho_thr [format {%.1e Ry} $t(1)]
    }
    
    ::pwtk::print PSEUDO "Final results for \"$ppFile\""
    ::pwtk::print "
  Generated files:
     ecutwfc data files    :  ecutwfc.$name.dat  ecutwfc.init.$name.dat
     Gnuplot plotting files:  ecutwfc.$name.gp   ecutwfc.init.$name.gp
     results plotted to $TERM:  ecutwfc.$name.$opt(t)

     ecutrho data file     :  $ecutrho_dat

  Used total-energy based thresholds:
     ecutwfc_threshold = [format %.1e $thr] Ry
     ecutrho_threshold = $ecutrho_thr

  Estimated energy cutoffs:
     ecutwfc = $ecutwfc Ry
     ecutrho = $ecutrho Ry  (dual = $dual)
"

    if { $opt(v) } {
        ::pwtk::display ecutwfc.$name.$opt(t)
    }

    # cleanup
    clean_

    return [list $ecutwfc $ecutrho]
}

TABLE OF CONTENTS

::pwtk::pseudo::etot_test