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