[Wien] Problems when trying to plot E vs c/a
Peter Blaha
pblaha at theochem.tuwien.ac.at
Thu May 17 12:32:24 CEST 2018
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-bounces at zeus.theochem.tuwien.ac.at] im Auftrag von Gavin Abo [gsabo at crimson.ua.edu]
> Gesendet: Dienstag, 15. Mai 2018 06:40
> An: wien at 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 at computername:~/Desktop$ cd $WIENROOT
> username at computername:~/WIEN2k$ wget https://raw.githubusercontent.com/gsabo/WIEN2k-Patches/master/17.1/init_lapw.patch
> username at 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: sabrygad at buffalo.edu<mailto:sabrygad at buffalo.edu>
> _______________________________________________
> Wien mailing list
> Wien at 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: blaha at theochem.tuwien.ac.at WIEN2k: http://www.wien2k.at
WWW: http://www.imc.tuwien.ac.at/TC_Blaha
--------------------------------------------------------------------------
-------------- next part --------------
#!/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
-------------- next part --------------
A non-text attachment was scrubbed...
Name: optimize.pl
Type: application/x-perl
Size: 5390 bytes
Desc: not available
URL: <http://zeus.theochem.tuwien.ac.at/pipermail/wien/attachments/20180517/c193c62f/attachment.pl>
More information about the Wien
mailing list