1#! /usr/bin/env python3
2
3import sys
4import os
5from numpy import array,polyfit,polyval,linspace,sqrt
6from qmcpack_analyzer import QmcpackAnalyzer
7from unit_converter import convert
8from plotting import *
9from debug import *
10
11
12args = sys.argv[1:]
13
14
15prefix = 'qmc'
16nfit   = 2
17Eatom  = None
18Eatom_err = 0
19if len(args)>0:
20    prefix = str(args[0])
21#end if
22if len(args)>1:
23    nfit = int(args[1])
24#end if
25if len(args)>2:
26    Eatom = float(args[2])
27#end if
28if len(args)>3:
29    Eatom_err = float(args[3])
30#end if
31if Eatom is None:
32    Elabel = 'total energy'
33    Eatom = 0
34else:
35    Elabel = 'binding energy'
36#end if
37
38
39dirs = os.listdir('./')
40
41scales = []
42energies = []
43errors = []
44for dir in dirs:
45    if dir.startswith('scale_'):
46        scale = float(dir.split('_')[1])
47        infile = os.path.join(dir,prefix+'.in.xml')
48        outfile = os.path.join(dir,prefix+'.s001.scalar.dat')
49        if os.path.exists(outfile):
50            qa = QmcpackAnalyzer(infile,analyze=True)
51            le = qa.qmc[1].scalars.LocalEnergy
52            scales.append(scale)
53            E = convert(le.mean,'Ha','eV')
54            Eerr = convert(le.error,'Ha','eV')
55            E -= 2*Eatom
56            Eerr = sqrt(Eerr**2+(2*Eatom_err)**2)
57            energies.append(E)
58            errors.append(Eerr)
59        #end if
60    #end if
61#end for
62
63
64
65se = 1.2074
66ee = -5.165
67
68seps = array(scales)*se
69order = seps.argsort()
70
71seps     = seps[order]
72energies = array(energies)[order]
73errors   = array(errors)[order]
74
75p = polyfit(seps,energies,nfit)
76sf = linspace(seps.min(),seps.max(),200)
77ef = polyval(p,sf)
78
79print()
80print('DMC '+Elabel+' vs separation distance')
81for i in range(len(seps)):
82    print('  {0:6.4f}  {1:6.4f} +/- {2:6.4f}'.format(seps[i],energies[i],errors[i]))
83#end for
84
85mloc = ef.argmin()
86print()
87print('DMC bond length, binding energy:',sf[mloc],ef[mloc])
88print('Exp bond length, binding energy:',se,ee)
89
90
91figure()
92errorbar(seps,energies,errors,fmt='r.')
93plot(sf,ef,'g-',lw=2,label='n={0} fit'.format(nfit))
94xlabel('bond length (A)')
95ylabel(Elabel+' (eV)')
96title('DMC '+Elabel+' vs. bond length')
97legend()
98
99show()
100