Thanks for the report. Modified eplot_lapw and SRC_w2web/htdocs/exec/optimize.pl
attached. On 05/16/2018 04:20 PM, Fecher, Gerhard wrote:
Dear c/a fitters,
This concerns the latest Wien2k version
I receive only the content of
test_opt.analysis
when I try with w2web to plot E vs c/a
but neither the result of the fit nor the plot are shown,
this seems to be a problem with the present version of the
eplot
script
when I use the eplot script of version 14.2 it is nearly ok, however,
there are still two issues: instead of the result of the fit, the "correlation
matrix of the fit parameters" is shown
and the figure is missing.
Reason is that eplot and optimize.pl do not work well together:
optimize.pl
prints the last 5 lines of fit.log (but the result is before these lines) and
expects the graph as case.c_over_a.png (but has a different name .coa.)
this can be solved by changing the two lines
line 169 change "CASE.c_over_a.png"
$umps = qx(cp $DIR/$CASE.coa.png $tempdir/$SID-$$.png);
line 173 change "tail -5"
$OUT .= qx(cd $DIR;echo ' ';echo "Fit of: E = a1 + a2*x + a3*x^2 +
a4*x^3 + a5*x^4";tail -15 fit.log);
or indeed, by changing eplot (I just did not find fast how to supress the
output of the correlation matrix)
It would also nice to have the minimum (should be prepared by findmincoa called
in eplot) in the fit.log and w2web output.
Ciao
Gerhard
DEEP THOUGHT in D. Adams; Hitchhikers Guide to the Galaxy:
"I think the problem, to be quite honest with you,
is that you have never actually known what the question is."
====================================
Dr. Gerhard H. Fecher
Institut of Inorganic and Analytical Chemistry
Johannes Gutenberg - University
55099 Mainz
and
Max Planck Institute for Chemical Physics of Solids
01187 Dresden
________________________________________
Von: Wien [wien-boun...@zeus.theochem.tuwien.ac.at] im Auftrag von Gavin Abo
[gs...@crimson.ua.edu]
Gesendet: Dienstag, 15. Mai 2018 06:40
An: wien@zeus.theochem.tuwien.ac.at
Betreff: Re: [Wien] Finite temperature DFT: electronic free energy
I think you found a bug in init_lapw. The fix seems like it should be simple
though. On line 498 in the init_lapw script in the sed command, it looks like
$fermit just needs changed to $fermits.
I made a patch file called init_lapw.patch [
https://github.com/gsabo/WIEN2k-Patches/tree/master/17.1 ], which you should be
able to apply and try simply using:
username@computername:~/Desktop$ cd $WIENROOT
username@computername:~/WIEN2k$ wget
https://raw.githubusercontent.com/gsabo/WIEN2k-Patches/master/17.1/init_lapw.patch
username@computername:~/WIEN2k$ patch -b init_lapw init_lapw.patch
patching file init_lapw
Maybe F is the :ENE (TOTAL ENERGY) in case.scf.
I tried following what a previous user did [
https://www.mail-archive.com/wien@zeus.theochem.tuwien.ac.at/msg10319.html ]
but with a quick TiC calculation using run_lapw in WIEN2k 17.1:
T = 1000 K = 0.00633 Ry
init_lapw -b -fermits 0.00633
TiC.scf::ENE : ********** TOTAL ENERGY IN Ry = -1783.95041939
TiC.output2: Kb * T = 0.00633000
TiC.output2: -T*Entr = -0.00054181
T = 0 K = 0.0 Ry
init_lapw -b -fermits 0.0
TiC.scf::ENE : ********** TOTAL ENERGY IN Ry = -1783.94928338
TiC.output2: Kb * T = 0.00180000
TiC.output2: -T*Entr = -0.00005329
F_1000K - F_0K = (E - T*S) - (E - 0) = -T*S = -1783.95041939 - (-1783.94928338) =
-0.001136 => -T*S = -0.001136
If -T*S divided by 2:
-T*S = -0.001136/2 = - 0.000568 <- This seems reasonably close to the above
-0.00054181 for T = 1000 K.
As I recall, the -T*S term was doubled in WIEN2k versions after 2014 [
https://www.mail-archive.com/wien@zeus.theochem.tuwien.ac.at/msg10319.html ].
Though, there may be a flaw in my above calculation.
On 5/14/2018 5:35 PM, Sabry Moustafa wrote:
Hi;
I am interested in the electronic free energy (F=E-TS) using the finite
temperature DFT using Fermi-Dirac statistics (i.e., Mermin's extension to 0 K
DFT) -- i.e., I am not just looking for Fermi-Dirac as a broadening/smearing
technique. So,how this can be done in WIEN2k?
It looks like I have to define "-fermits X" in the init_lapw command (where X=
kT, in Ry). But when I did (T=1000K in my case):
init_lapw -b -fermits 0.00633
I got "fermit: Undefined variable." at the end. This is fixed though when
defined fermit as well:
init_lapw -b -fermit 0.00633 -fermits 0.00633
So, do I have to define both?
Also, where can I find this "free energy F". It does not seem to be the "TOTAL ENERGY" in case.scf.
Seems like I need to add this "TOTAL ENERGY" (E) to "-T*Entr" (-TS) in case.output2 file to get F?
Thanks ;
Sabry
:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:-:
Sabry G. Moustafa, Ph.D.
Department of Chemical and Biological Engineering,
University at Buffalo, The State University of New York.
511 Furnas Hall
Buffalo, NY 14260-4200
716-645-1186 (office)
716-239-8543 (cell)
E-mail: sabry...@buffalo.edu<mailto:sabry...@buffalo.edu>
_______________________________________________
Wien mailing list
Wien@zeus.theochem.tuwien.ac.at
http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien
SEARCH the MAILING-LIST at:
http://www.mail-archive.com/wien@zeus.theochem.tuwien.ac.at/index.html
--
P.Blaha
--------------------------------------------------------------------------
Peter BLAHA, Inst.f. Materials Chemistry, TU Vienna, A-1060 Vienna
Phone: +43-1-58801-165300 FAX: +43-1-58801-165982
Email: bl...@theochem.tuwien.ac.at WIEN2k: http://www.wien2k.at
WWW: http://www.imc.tuwien.ac.at/TC_Blaha
--------------------------------------------------------------------------
#!/bin/csh -f
# interface for plotting E vs. c/a curves
# data is generated with: optimize and "Analyze multiple SCF Files"
#
unalias rm
echo ""
echo ""
echo "####################################"
echo "# #"
echo "# E - PLOT #"
echo "# Extension by Morteza Jamal #"
echo "# (2014) #"
echo "####################################"
echo ""
echo ""
set tmp = tmp
set tmp2 = tmp2
set print = eplot.ps
unset type
unset analysis
unset terminal
unset savestring
set file = `pwd`
set file = $file:t
#set file=(`ls *.analysis`)
#set file = $file[1]:r
unset help
while ($#argv)
switch ($1)
case -h:
set help
shift; breaksw
case -t:
shift
set type = $1
shift; breaksw
case -p:
# shift
set terminal = png
shift; breaksw
case -f:
shift; set file = $1
shift; breaksw
case -a:
shift; set analysis;set savestring="$1"
if( "$savestring" == " ") set savestring
# grepline :ene "$string" 1 > $file.analysis
# grepline :vol "$string" 1 >> $file.analysis
shift; breaksw
default:
shift; breaksw
endsw
end
if ($?help) goto help
if !($?type) then
settype: #added
30|08|16, RUH
echo 'type "boa" for b/a or'
echo ' "coa" for c/a or'
echo ' "vol" for volume curve '
set type=$<
if ( $type == '' ) then #added
30|08|16, RUH
goto settype
else if ( $type != 'vol' && $type != 'coa' && $type != 'boa' ) then
goto settype
endif #end
addition
#if !($?type) then
# echo 'type "boa" for b/a or'
# echo ' "coa" for c/a or'
# echo ' "vol" for volume curve '
# set type=$<
endif
if($?analysis) then
set string="*${type}*${savestring}.scf";set savestring="${type}_$savestring"
grepline :ene "$string" 1 > $file.analysis
grepline :vol "$string" 1 >> $file.analysis
endif
if !(-e $file.analysis) goto error
if($?savestring) then
set print=$file.eplot_$savestring.ps
endif
if ( $type == 'vol' ) then
# set ene=`grepline_lapw :ene '*.scf' 1 | cut -f2 -d=`
# set vol=`grepline_lapw :vol '*.scf' 1 | cut -f2 -d=`
#echo $ene
#echo $vol
# set i=4
#touch $file.vol
#rm $file.vol
#loop:
# echo $vol[$i] $ene[$i] >>$file.vol
# @ i ++
#if ( $i <= $#ene ) goto loop
set ene=`grep :ENE $file.analysis | cut -f2 -d=`
set vol=`grep :VOL $file.analysis | cut -f2 -d=`
if (-e $file.vol) rm $file.vol
set i = 0
loop:
echo $vol[$i] $ene[$i] >>$file.vol
@ i ++
if ( $i <= $#ene ) goto loop
#bulk <$file.vol
x eosfit -f $file
cat $file.outputeos
echo " The above output is also in $file.outputeos"
echo " "
echo " Murnaghan-data are in $file.eosfit"
echo " Birch-Murnaghan-data are in $file.eosfitb"
echo " Vinet-Rose-data are in $file.eosfitv"
echo " Poirier-Tarantola-data are in $file.eosfitp"
echo " "
echo ' Plot Murnaghan,Birch-Murnaghan, Vinet-Rose or Poirier-Tarantola fit:
[M/b/v/p]" '
set fit=$<
echo "You may want to print $file.outputeos"
switch ($fit)
case [P,p]:
set murna=`grep V0, $file.outputeos | grep -v \* |tail -1`
set plotfile=$file.eosfitp
breaksw
case [V,v]:
set murna=`grep V0, $file.outputeos | grep -v \* |tail -2|head -1`
set plotfile=$file.eosfitv
breaksw
case [B,b]:
set murna=`grep V0, $file.outputeos | grep -v \* |tail -3|head -1`
set plotfile=$file.eosfitb
breaksw
default:
set murna=`grep V0, $file.outputeos | grep -v \* |tail -4|head -1`
set plotfile=$file.eosfit
breaksw
endsw
#echo $murna
#grep $type $file.analysis |sed "1d"| sed "s/.*$type//" |cut -c1-6,56- | tr "_"
" " |sort -n >$tmp2
#grep $type $file.analysis | grep -v "Analysis generated" |\
# sed "s/.*$type//" | tr ":a-z" " " | awk '{print $1 " " $NF}'|\
# cut -c1-6,8- | tr "_" " " |sort -n >$tmp2
set xmin
set xmax
set ymin
set ymax
replot:
rm -f $tmp
if ($?terminal) then
cat <<EOF >$tmp
set terminal png
set output 'eplot.png'
EOF
endif
cat <<EOF >>$tmp
set format y "%.4f"
set title "$file"
set xlabel "Volume (bohr^3)"
set ylabel "Energy (Ry)"
set xrange [${xmin}:${xmax}]
set yrange [${ymin}:${ymax}]
set xzeroaxis
#plot "$tmp2" title "$file"
plot "$file.vol" title "Murnaghan: $murna[1]" w p
replot "$plotfile" title "$murna[2-]" w l
EOF
if (! $?terminal) echo pause -1 >>$tmp
echo "press RETURN to continue"
gnuplot $tmp
echo -n "Do you want to set ranges? (y/N)"
set yn = ($<)
if ($yn == y) then
unset yn
echo -n " Limit x-range? (y/N)"
set yn = ($<)
if ($yn == y) then
unset yn
echo -n " xmin="
set xmin = ($<)
echo -n " xmax="
set xmax = ($<)
endif
echo -n " Limit y-range? (y/N)"
set yn = ($<)
if ($yn == y) then
unset yn
echo -n " ymin="
set ymin = ($<)
echo -n " ymax="
set ymax = ($<)
endif
goto replot
endif
echo -n "Do you want a hardcopy? (Y/n)"
set hardcopy = ($<)
if ($hardcopy != n) then
echo -n "Specify a filename (default is $print)"
set out = ($<)
echo "Printing hardcopy"
if ($out == "") set out = "$print"
cat <<EOF >$tmp
set terminal postscript
set output "$out"
#set data style linespoints
set format y "%.4f"
set title "$file"
set xlabel "Volume (bohr^3)"
set ylabel "Energy (Ry)"
set xrange [${xmin}:${xmax}]
set yrange [${ymin}:${ymax}]
plot "$file.vol" title "Murnaghan: $murna[1]" w p,"$plotfile" title
"$murna[2-]" w l
EOF
gnuplot $tmp
endif
else
########### by Morteza############
set namef=`echo "_initial.struct"`
cp $file$namef init.struct
if ( $type == coa ) then
set type2 = "c/a"
else
set type2 = "b/a"
endif
echo "$type" > tmp3
##################################
# here we do coa's
#grep $type $file.analysis |sed "1d"| sed "s/.*$type//" |cut -c1-6,56- | tr "_"
" " |sort -n >$tmp2
grep $type $file.analysis | grep -v "Analysis generated" | grep -v ":VOL " |\
sed "s/.*$type//" | tr ":a-zA-Z" " " | awk '{print $1 " " $NF}'|\
cut -c1-9,10- | tr "_" " " |sort -n >$tmp2
# Modified 24|04|2017 Sohaib
echo ' ' >$tmp
if ($?terminal) then
cat <<EOF >$tmp
set terminal png
set output '$file.$type.png'
EOF
endif
cat <<EOF >>$tmp
set xlabel "deviation from exp. $type2 ratio (%)"
set ylabel "Energy [Ry]"
set format y "%.4f"
f(x)=a1+a2*x+a3*x**2+a4*x**3+a5*x**4
fit f(x) '$tmp2' via '.fitparam'
plot "$tmp2" title "data" w p , f(x) title "polyfit\\\_4^{th}order"
save var 'varfit.dat'
#plot "$tmp2" title "$file"
EOF
if (! $?terminal) echo pause -1 >>$tmp
#############################################
#set e1 = `head -1 <tmp2 | awk '{print $2}'`
#############################################
cat <<EOF >.fitparam
a1=1
a2=1
a3=1
a4=1
a5=1
EOF
gnuplot $tmp
if ($?terminal) then
set out=eplot.ps
set hardcopy=y
else
echo -n "Do you want a hardcopy? (Y/n)"
set hardcopy = ($<)
if ($hardcopy != n) then
echo -n "Specify a filename (default is $print)"
set out = ($<)
echo "Printing hardcopy"
if ($out == "") set out = "$print"
endif
endif
#echo "press RETURN to continue" $hardcopy
if ($hardcopy != n) then
cat <<EOF >$tmp
set terminal postscript
set output "$out"
#set data style linespoints
set format y "%.4f"
f(x)=a1+a2*x+a3*x**2+a4*x**3+a5*x**4
fit f(x) '$tmp2' via '.fitparam'
plot "$tmp2" title "data" w p , f(x) title "polyfit\\\_4^{th}order"
#plot "$tmp2" title "$file"
EOF
gnuplot $tmp >& /dev/null
if ($?terminal) echo $file.$type.png generated
echo $out generated
endif
##################by Morteza#######################
#clear
echo ""
echo "******************************"
#************here we take varaible of fit*********************
if !(-e "varfit.dat") then
set filer = 'varfit.dat'
goto error
endif
grep "a[1-5]" varfit.dat > .var
head -5 < .var|awk '{print $3}' > varfit.dat
findMINcboa
echo "******************************"
echo ""
##################################################
endif
#rm $tmp $tmp2
mv $tmp eplot.gnuplot
if($?savestring) then
cp $file.outputeos $file.outputeos_$savestring
endif
exit 0
error:
echo ">>>"
echo ">>> ERROR: $file.analysis not found\!"
echo ">>> ERROR:"
echo '>>> ERROR: You should "Anaylze multiple SCF Files" first'
echo ">>>"
#exit(1)
help:
set type=vol
cat <<EOF
EPLOT is a plotting interface to plot E vs. Vol or c/a or b/a curves.
Once you have several scf calculations at different volumes (usually generated
with "optimize.job") you can generate the required "$file.analysis" using:
grepline :ENE "*$type*.scf" 1 > $file.analysis
grepline :VOL "*$type*.scf" 1 >> $file.analysis
Alternatively, you can also use:
"eplot -a searchstring",
which will generate $file.analysis for you and create
$file.outputeos_searchstring
Example: eplot -a " " will analyse all scf files '*$type*.scf'
Example: eplot -a pbe will analyse all scf files '*$type*pbe.scf'
(and $type will be vol/coa/boa)
Generates plots in X-window (default), png (-p) and ps format.
eplot [-t vol/coa/boa] [-p] [-f FILEHEAD] [-a search-string_in_scf-files]
EOF
optimize.pl
Description: Perl program
_______________________________________________ Wien mailing list Wien@zeus.theochem.tuwien.ac.at http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien SEARCH the MAILING-LIST at: http://www.mail-archive.com/wien@zeus.theochem.tuwien.ac.at/index.html
