### Imports run data from other routines from math import sqrt from variables import alphalist,hubbylist,molname,charge,spin,vis from atreader import * import re ### Checks for scipy and matplotlib try: from scipy import linalg scipy_loaded = True except ImportError: scipy_loaded = False try: import numpy as np import matplotlib.pyplot as plt import matplotlib.colors as colors import matplotlib.colorbar as colorbar import matplotlib.cm as cm matplotlib_loaded = True except ImportError: matplotlib_loaded = False ### Checks for json version try: import simplejson as json except ImportError: import json as json ### Basic linear regression routine especially needed if no scipy. count=0 def linreg(X, Y): """ Summary Linear regression of y = ax + b Usage real, real, real = linreg(list, list) Returns coefficients to the regression line "y=ax+b" from x[] and y[], and R^2 Value """ if len(X) != len(Y): raise ValueError, 'unequal length' N = len(X) Sx = Sy = Sxx = Syy = Sxy = 0.0 for x, y in map(None, X, Y): Sx = Sx + x Sy = Sy + y Sxx = Sxx + x*x Syy = Syy + y*y Sxy = Sxy + x*y det = Sxx * N - Sx * Sx a, b = (Sxy * N - Sy * Sx)/det, (Sxx * Sy - Sx * Sxy)/det meanerror = residual = 0.0 for x, y in map(None, X, Y): meanerror = meanerror + (y - Sy/N)**2 residual = residual + (y - a * x - b)**2 RR = 1 - residual/meanerror ss = residual / (N-2) Var_a, Var_b = ss * N / det, ss * Sxx / det return a, b, RR ### Start collecting occupations from output files if __name__=='__main__': bareresp=[] convresp=[] xtraresp=[] for num in range(0,len(typno)*len(typno)): bareresp.append(alphalist[1:1000]) convresp.append(alphalist[1:1000]) xtraresp.append(alphalist[1:1000]) alphalist2=alphalist[1:1000] hubbylist2=alphalist[1:1000] for lines in range(0,len(alphalist)-1): for line2 in range(0,len(typno)*len(typno)): bareresp[line2][lines]=0.0 convresp[line2][lines]=0.0 xtraresp[line2][lines]=0.0 alphalist2[lines]=float(alphalist[lines+1]) hubbylist2[lines]=0.0 hubbylist2.append(0.0) results=open('linrespmat.dat','w') if (len(hubbylist) > 1): uinl=hubbylist[0:1000] uoutl=hubbylist[0:1000] for line4 in range(0,len(hubbylist)): ### switch between Python and fortran numbering conventions if ( 'D' in hubbylist[line4]): hubbylist2[line4]=hubbylist[line4].split('D')[0]+'E'+hubbylist[line4].split('D')[1] else: hubbylist2[line4]=hubbylist[line4] count3=-1 for hubby in hubbylist: count3+=1 for hubtyp in typno: for alpha in alphalist[1:1000]: ### Output files from jobrun.py, same format. pre=molname+'.s'+str(spin+1)+'.n'+species[hubtyp-1]+'.u'+str(hubby) file=pre+'.a'+alpha+'.out' output=open(file,'r').readlines() count=0 for lines in range(0,len(output)): for hubnums in typno: ### We loop over the different Hubbard atom types. match=re.search('atom'+' '*(5-len(str(hubnums).split()))+str(hubnums)+' Tr',output[lines]) if match is not None: if len(typno)<=count=len(typno)*2: convresp[typno.index(hubtyp)*len(typno)+typno.index(hubnums)][alphalist.index(alpha)-1]=float(output[lines].split()[3]) count+=1 chi0,chif=[],[] if (len(typno) == 2): if scipy_loaded is False: ### This is only needed if we don't have scipy and we have a 2x2 chii0,chiif,umat=[[0.0,0.0],[0.0,0.0]],[[0.0,0.0],[0.0,0.0]],[[0.0,0.0],[0.0,0.0]] for rstyp1 in range(0,len(typno)): chi0.append([]) chif.append([]) for rstyp2 in range(0,len(typno)): chi0[rstyp1].append(linreg(alphalist2,bareresp[rstyp1*len(typno)+rstyp2])[0]) chif[rstyp1].append(linreg(alphalist2,convresp[rstyp1*len(typno)+rstyp2])[0]) if scipy_loaded is True: ### Scipy matrix inversion print "U matrix\n" results.write("U matrix\n") umat=linalg.inv(chi0)-linalg.inv(chif) print umat np.savetxt(results,umat) elif len(typno)==1: ### Manual 2x2 matrix inversion det0=chi0[1][1]*chi0[0][0]-chi0[0][1]*chi0[1][0] detf=chif[1][1]*chif[0][0]-chif[0][1]*chif[1][0] chii0[1][1]=chi0[0][0]/det0 chiif[1][1]=chif[0][0]/detf umat[1][1]=chii0[1][1]-chiif[1][1] chii0[0][0]=chi0[1][1]/det0 chiif[0][0]=chif[1][1]/detf umat[0][0]=chii0[0][0]-chiif[0][0] chii0[0][1]=-chi0[0][1]/det0 chii0[1][0]=-chi0[1][0]/det0 chiif[0][1]=-chif[0][1]/detf chiif[1][0]=-chif[1][0]/detf umat[0][1]=chii0[0][1]-chiif[0][1] umat[1][0]=chii0[1][0]-chiif[1][0] print "Umatrix\n" print ' '*(6+len(str(umat[0][0]))/2-len(str(typnum[0]))/2), typnum[0],' '*(len(str(umat[1][1]))-len(str(typnum[0]))),typnum[1] print typnum[0],' '*(5-len(str(typnum[0]))),umat[0][0],umat[0][1] print typnum[1],' '*(5-len(str(typnum[1]))),umat[1][0],umat[1][1] results.write("Umatrix\n") results.write("%s %s %s %s\n" %(' '*(6+len(str(umat[0][0]))/2-len(str(typnum[0]))/2), typnum[0],' '*(len(str(umat[1][1]))-len(str(typnum[0]))),typnum[1])) results.write("%s %s %s %s \n" %(typnum[0],' '*(5-len(str(typnum[0]))),umat[0][0],umat[0][1])) results.write("%s %s %s %s\n" %(typnum[1],' '*(5-len(str(typnum[1]))),umat[1][0],umat[1][1])) else: ### Larger than 2x2 and no scipy. results.write("Bare response matrix elements\n") json.dump(chi0,results) results.write("\nConverged response matrix elements\n") json.dump(chif,results) results.write("\nRows and columns ordered by atom type order\n") json.dump(typnum,results) if matplotlib_loaded is True: if scipy_loaded is True: if vis=='yes': ###additional plotting of matrix elements fig=plt.figure() fig.subplots_adjust(wspace=0.3) ticknums,typnums=[-0.5],[''] for num in range(0,len(typnum)): ticknums.append(num) typnums.append(typnum[num]) ticknums.append(len(typnum)-0.5) typnums.append('') ax = fig.add_subplot(131, title='Bare Resp') colors.Normalize(vmin=None, vmax=None, clip=False) ax.imshow(chi0, cmap=cm.RdYlBu, interpolation='nearest') ax.locs,ax.labels=plt.xticks(ticknums,typnums) ax.locs,ax.labels=plt.yticks(ticknums,typnums) for q in range(0,len(umat)): for s in range(0,len(umat)): ax.annotate('%1.2f' %(chi0[s][q]), xy=(q-0.25, s), xycoords='data', horizontalalignment='left', verticalalignment='bottom',color='white') ax2 = fig.add_subplot(132, title='Converged Resp') ax2.imshow(chif, cmap=cm.RdYlBu, interpolation='nearest') ax2.locs,ax.labels=plt.xticks(ticknums,typnums) ax2.locs,ax.labels=plt.yticks(ticknums,typnums) for q in range(0,len(umat)): for s in range(0,len(umat)): ax2.annotate('%1.2f' %(chif[s][q]), xy=(q-0.25, s), xycoords='data', horizontalalignment='left', verticalalignment='bottom',color='white') ax3 = fig.add_subplot(133,title='U matrix') ax3.imshow(umat, cmap=cm.RdYlBu, interpolation='nearest') ax3.locs,ax.labels=plt.xticks(ticknums,typnums) ax3.locs,ax.labels=plt.yticks(ticknums,typnums) for q in range(0,len(umat)): for s in range(0,len(umat)): ax3.annotate('%1.2f' %(umat[s][q]), xy=(q-0.25, s), xycoords='data', horizontalalignment='left', verticalalignment='bottom', color='white') plt.savefig('respmats')