TABLE OF CONTENTS
::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 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 "$ecutrho $ecutwfc" list.
SOURCE
variable E variable gpData ::pwtk::printTitle PSEUDO "Automatic search of \"ecutwfc\" & \"ecutrho\" for a pseudo-potential\n(based on poor man's total-energy criterium)" if { [auto_execok gnuplot] == {} } { ::pwtk::error "The [::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 { {e.arg {} "list of ecutwfc value for initial scanning (minimum of 4 value needed; default = 30 60 90 120)"} {d.arg -16 "value of dual for US and PAW potentials"} {thr.arg 1e-3 "ecutwfc convergence treshold for total energy (in Ry)"} {u.arg 4 "upscale factor for ecutrho threshold, i.e., rhothr = thr/upscale"} {n.arg 6 "number of pw.x calcs during the refinement stage"} {l.arg 0.8 "low factor of ecutwfc for the refinement stage"} {h.arg 2.0 "high factor of ecutwfc for the refinement stage"} {m.arg 10 "minimum allowed ecutwfc (if ecutwfc gets lower, drop error)"} {M.arg 1000 "maximum allowed ecutwfc (if ecutwfc gets larger, drop error)"} {t.arg eps "a Gnuplot file terminal to plot to"} {v "visualize the resulting plot"} } ::pwtk::parseOpt_ ::pwtk::checkOTypeStrict_ -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 {30 60 90 120} } } 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] }