[Wien] Problems when trying to plot E vs c/a

Fecher, Gerhard fecher at uni-mainz.de
Thu May 17 13:30:48 CEST 2018


Hallo Peter,
thanks for the files.
unforunately, the otimize.pl still doesn't show the result of the fit (plot is there)
output is in a shortened version:

Fit of:  E = a1 + a2*x + a3*x^2 + a4*x^3 + a5*x^4
a1              1.000 
a2              0.000  1.000 
a3             -0.725 -0.000  1.000 
a4             -0.000 -0.930  0.000  1.000 
a5              0.648  0.000 -0.985 -0.000  1.000 

the line 174 should contain at least   tail -15    (instead of -5)    what results in the output of the parameters and the correlation matrix 

Fit of:  E = a1 + a2*x + a3*x^2 + a4*x^3 + a5*x^4
Final set of parameters            Asymptotic Standard Error
=======================            ==========================
a1              = -5573.9          +/- 3.634e-06    (6.519e-08%)
a2              = 4.23124e-06      +/- 9.205e-06    (217.5%)
a3              = 0.000137795      +/- 2.93e-05     (21.26%)
a4              = 7.61902e-06      +/- 1.037e-05    (136.1%)
a5              = -1.43164e-05     +/- 2.725e-05    (190.3%)

correlation matrix of the fit parameters:
                a1     a2     a3     a4     a5     
a1              1.000 
a2              0.000  1.000 
a3             -0.725 -0.000  1.000 
a4             -0.000 -0.930  0.000  1.000 
a5              0.648  0.000 -0.985 -0.000  1.000

or shorter versuion is to use  tail -15 fit.log  | head -7   because I don't think that the correlation matrix is needed in the w2web output (it's found in fit.log anyway)
the result is then only
 
Fit of:  E = a1 + a2*x + a3*x^2 + a4*x^3 + a5*x^4
Final set of parameters            Asymptotic Standard Error
=======================            ==========================
a1              = -5573.9          +/- 3.634e-06    (6.519e-08%)
a2              = 4.23124e-06      +/- 9.205e-06    (217.5%)
a3              = 0.000137795      +/- 2.93e-05     (21.26%)
a4              = 7.61902e-06      +/- 1.037e-05    (136.1%)
a5              = -1.43164e-05     +/- 2.725e-05    (190.3%)

the optimize.pl file changed in the latter way is attached


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 Peter Blaha [pblaha at theochem.tuwien.ac.at]
Gesendet: Donnerstag, 17. Mai 2018 12:32
An: wien at zeus.theochem.tuwien.ac.at
Betreff: Re: [Wien] Problems when trying to plot E vs c/a

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 --------------
A non-text attachment was scrubbed...
Name: optimize.pl
Type: application/octet-stream
Size: 5607 bytes
Desc: optimize.pl
URL: <http://zeus.theochem.tuwien.ac.at/pipermail/wien/attachments/20180517/bab90ede/attachment.obj>


More information about the Wien mailing list