#!/usr/bin/env python import sys, re, os, glob, shutil from scipy import * import ldau import krams import strns as sts import linalg as la ROOT = os.environ.get('LMTO_DMFT_ROOT') class IMP_OCA(object): """One crossing approximation Impurity solver """ def __init__(self, Znuc, dire, Impq, params, UpdateAtom=False): # impurity directory self.dir = dire+'/' self.UpdateAtom = UpdateAtom if len(glob.glob(dire))==0 : os.mkdir( dire ) # log-file is opened self.fh_info = open(self.dir + 'info.out', 'w') self.fh_infoj = open(self.dir + 'jinfo.out', 'w') self.sparams={'Sig':'Sigma.000', 'Ac':'Ac.inp', 'cix':'out.cix', 'AlocOut':'Aloc.imp', 'SigOut':'Sigma.000', 'sig': 'sig.inp', 'gloc':'gloc.out', 'pcore':0, 'acore':0} shutil.copy2('Sigma.000', self.dir+self.sparams['Sig']) # file from database which contains atomic information for this particular atom self.atom_arg = ROOT + '/impurity/database/' + 'atom.'+str(Znuc)+'.py' # executable which broadens self.broad = ROOT + '/impurity/bin/broad' self.PARAMS = 'PARAMS.oca' execfile(self.atom_arg) lcs = locals() # it seems python does not know about variables obtained by execfile # stores some of just defined variables self.l = lcs['l'] self.spr={} if params.has_key('nc'): self.spr['n'] = params['nc'][0] else: self.spr['n'] = lcs['n'] if params.has_key('Ncentral'): self.spr['Ncentral'] = params['Ncentral'][0] else: self.spr['Ncentral'] = self.spr['n'][1:len(self.spr['n'])-1] if params.has_key('para'): self.spr['para'] = params['para'][0] else: self.spr['para'] = lcs['para'] if params.has_key('kOCA'): self.spr['kOCA'] = params['kOCA'][0] else: self.spr['kOCA'] = -1 if params.has_key('mOCA'): self.spr['mOCA'] = params['mOCA'][0] else: self.spr['mOCA'] = 1e-3 if params.has_key('Ex'): self.spr['Ex'] = params['Ex'][0] else: self.spr['Ex'] = lcs['Ex'] if params.has_key('Ep'): self.spr['Ep'] = params['Ep'][0] else: self.spr['Ep'] = lcs['Ep'] if params.has_key('J'): self.spr['J'] = params['J'][0] else: self.spr['J'] = lcs['J'] if params.has_key('cx'): self.spr['cx'] = params['cx'][0] else: self.spr['cx'] = lcs['cx'] if params.has_key('Eoca'): self.spr['Eoca'] = params['Eoca'][0] else: self.spr['Eoca'] = lcs['Eoca'] if params.has_key('qatom'): self.spr['qatom'] = params['qatom'][0] else: self.spr['qatom'] = lcs['qatom'] #self.J = lcs['J'] #self.cx = lcs['cx'] self.J = self.spr['J'] self.cx = self.spr['cx'] if type(Impq)!= list : Impq = Impq.tolist() self.Impq = Impq print >> self.fh_info, 'Impq=', Impq # Below exact diagonalization of the atom is run to produce 'cix.out' file # information is read from the file in database if not self.UpdateAtom: # executable for ED of the atom self.atom_exe = ROOT + '/impurity/bin/atom' # argument list is created to later run the executable for exact diagonalization of the atom atom_arg = 'n="%s" l=%d J=%f cx=%f para=%s Ex="%s" Ep="%s" qOCA=%d Eoca=%f qatom=%d Ncentral=%s mOCA=%s kOCA="%s"' % (self.spr['n'], lcs['l'], self.spr['J'], self.spr['cx'], self.spr['para'], self.spr['Ex'], self.spr['Ep'], lcs['qOCA'], self.spr['Eoca'], self.spr['qatom'], self.spr['Ncentral'], self.spr['mOCA'], self.spr['kOCA']) # runs ED for the atom stdin, stdout, stderr = os.popen3('cd '+ self.dir + '; ' + self.atom_exe +' ' +atom_arg+' > nohup.out 2>&1') print >> self.fh_info, stdout.read() print >> self.fh_info, stderr.read() else: self.atom_exe = ROOT + '/impurity/bin/atom_d.py' def _exe(self, params, mpi_prefix, om, Delta, gbroad, extn, Eimps, kbroad=0.0, maxAc=20.): """ """ if self.UpdateAtom: #Eaver = sum(Eimps)/len(Eimps) # One should not double-count impurity levels, this is average #Eimps = array(Eimps)-Eaver # the cix file will contain only the deviations fromaverage #params['Ed'] = Eaver*ones(len(params['Ed'])) # the average will enter by PARAMS.oca file print 'Eimps=', Eimps print 'Ed=', params['Ed'] #atom_arg = self.atom_arg + ' \"Eimp=['+(('%f,'*len(Eimps))%tuple(Eimps))+']\"' atom_arg = self.atom_arg + ' \"Eimp='+str(Eimps)+'\" \"n='+str(self.spr['n'])+'\" \"Ncentral='+str(self.spr['Ncentral'])+'\"' atom_arg += ' "Impq='+str(self.Impq)+'"' # runs ED for the atom stdin, stdout, stderr = os.popen3('cd '+ self.dir + '; ' + self.atom_exe +' ' +atom_arg+' > nohup.out 2>&1') print >> self.fh_info, stdout.read() print >> self.fh_info, stderr.read() # OCA executable shutil.copy2(ROOT + '/impurity/bin/' + params['exe'][0], self.dir+params['exe'][0]) self.oca_exe = './'+params['exe'][0] # Preparing PARAMS.oca file fp = open(self.dir+self.PARAMS, 'w') for p in self.sparams.keys(): fp.write(p + '=' + str(self.sparams[p]) + '\n') for p in params.keys(): print >> fp, "%s=%s " % (p, params[p][0]), '\t', params[p][1] fp.close() # Preparing Ac.inp file fp = open(self.dir+self.sparams['Ac']+'b', 'w') QcutAc = False if 'cutAc' in params.keys(): QcutAc = True cutAc = params['cutAc'][0] for im in range(len(om)): print >> fp, om[im], for b in range(len(Delta[im])): Ac = -Delta[im][b].imag/pi if (Ac<0.0): Ac=0 if (Ac>maxAc): Ac=0 if (QcutAc): Ac *= ferm((cutAc[0]-om[im])/params['Th'][0])*ferm((om[im]-cutAc[1])/params['Th'][0]) print >> fp, Ac, print >> fp fp.close() if gbroad>0: broad_cmd = 'cd '+self.dir+'; '+self.broad+' -w '+str(gbroad)+' -k '+str(kbroad)+' '+self.sparams['Ac']+'b'+' > '+self.sparams['Ac'] print broad_cmd stdin, stdout, stderr = os.popen3(broad_cmd) print stderr.read() else: shutil.move(self.dir+self.sparams['Ac']+'b', self.dir+self.sparams['Ac']) shutil.copy2(self.dir+self.sparams['Ac'], self.dir+self.sparams['Ac']+'.'+extn) # Below we execute oca oca_cmd = 'cd '+self.dir+'; '+mpi_prefix+' '+self.oca_exe+' '+self.PARAMS+' > nohup_imp.out 2>&1 ' print oca_cmd print 'Running ---- OCA -----' stdin, stdout, stderr = os.popen3(oca_cmd) print >> self.fh_infoj, stdout.read(), stderr.read() def HighFrequency(self, params, DCs, da, extn): """ Corrects high-frequency of OCA spetral functions such that it gives correct nf, normalization and sigma_oo This is achieved by adding two Lorentzians at high frequency. One below EF in the interval (par['epsilon'][0][0],par[epsilon][0][1]) and one above EF in the interval (par['epsilon'][1][0],par[epsilon][1][1]) Input: params -- parameters must contain T Ed U epsilon - prefered positions of additional lorentzians Th - when high-frequency spectral function is cut, it is cut by fermi function with temperature Th Gh - the width of the two lorentzians added Output: Sig[b][iom] -- DMFT dynamic self-energy which vanishes at infinity sinf[b] -- DMFT self-energy at infinity Edc[b] -- double counting using DMFT occupancies The 'corrected' Green's function will be of the form: G(om) = int(A(x)/(om-x)) + a1/(om-eps1+i*Gh) + a2/(om-eps2+i*Gh) The equations to be satisfied are: normalization: m0 + a1 + a2 = 1 density: nc + a1 = nc_{exact} where nc=int(A(x)f(x)) and nc_{exact} is computed from pseudo spectral functions Sigma_oo: m1 + a1*eps1 + a2*eps2 = Eimp+Sigma_oo == dsinf The equations can be brought to the form x*u + y*v = w with u and v positive numbers and x and y unknown numbers in the interval [0,1] In terms of above quantities, we have u = a1*(eps[0][1]-eps[0][0]) v = a2*(eps[1][1]-eps[1][0]) w = Eimp+Sigma_oo - m1 - a1*eps[0][0] - a2*eps[1][0] x = (eps1-eps[0][0])/(eps[0][1]-eps[0][0]) y = (eps2-eps[1][0])/(eps[1][1]-eps[1][0]) The solution exists if 0u): x=0.5*(x_min+x_max) else : y=0.5*(y_min+y_max) """ ################################################ # Reading of parameters from impurity cix file # ################################################ # cix file is used to find the following variables: J, l, Ns, baths fc = open(self.dir + self.sparams['cix']) first_line = fc.next() # next line of cix contains Ns and baths cixdat = fc.next().split() baths = int(cixdat[0]) Ns = map(int, cixdat[1:1+baths]) # Reading Aloc.imp which was produced by OCA fA = open(self.dir + self.sparams['AlocOut']) # Aloc.imp contains nf, moment, ntot,... first_line = fA.next().strip() adat = first_line.split() for par in adat: m = re.search('[ntot|nf|moment]=', par) if m is not None : exec(par) om=[] Af=[] for p in fA: dat = p.split() om.append(float(dat[0])) Af.append(map(float,dat[1:1+baths])) Af=array(Af) om=array(om) # reading of bath Weiss field fAc = open(self.dir + self.sparams['Ac']) Acd=[] ii=0 for p in fAc: dat = p.split() omega = float(dat[0]) if (abs(omega-om[ii])>1e-5): print 'Seems that %s and %s are not compatible. Exiting!'%(self.sparams['AlocOut'],self.sparams['Ac']) sys.exit(1) ii += 1 Acd.append(map(float,dat[1:1+baths])) Acd=array(Acd) # define some common variables in this functions T = params['T'][0] Ed = params['Ed'][0] epsilon = params['epsilon'][0] # Here we construct few functions for faster computation # of moments and occupancies # # dh - mesh weights according to trapezoid rule # omdh - omega*dh for calculation of first moment # fedh - f(omega)*dh for calculation of occupancy # dh=[0.5*(om[1]-om[0])] omdh=[om[0]*dh[-1]] fedh=[ferm(om[0]/T)*dh[-1]] for im in range(1,len(om)-1): dh.append(0.5*(om[im+1]-om[im-1])) omdh.append(om[im]*dh[-1]) fedh.append(ferm(om[im]/T)*dh[-1]) dh.append(0.5*(om[-1]-om[-2])) omdh.append(om[-1]*dh[-1]) fedh.append(ferm(om[-1]/T)*dh[-1]) dh=array(dh) omdh=array(omdh) fedh=array(fedh) #da = 0 #if DCs['scheme']=='default' : da = DCs['a'] # Here we compute self-energy(infinity) and # double-counting using impurity occupancy (sinf, Edc) = Sinftyv(params['U'][0], self.J, self.l, baths, Ns, nf, da) # If user fixes double counting by certain predescribed nf0 # (fixn=True), we need to normalized actual nf[:] such that # sum(nf[:]) is equal to predefined nf0 # Here we do not change Sigma_{infinity} but only Edc. # Sigma_{infinity} is determined by actual nf, while Edc is determined by modified nf. if DCs=='fixn': nfdc = nf[:] nfdc = array(nfdc)*(params['nf0'][0]/sum(nfdc)) da = 0.0 (sinf_dummy, Edcn) = Sinftyv(params['U'][0], self.J, self.l, baths, Ns, nfdc, da) ta = (Edc-Edcn)/params['U'][0] print >>self.fh_info, ('### a=%f' % (sum(ta)/len(ta))) Edc = Edcn[:] # Some info up to now print >>self.fh_info, '# l=%d T=%f U=%f J=%f' % (self.l, T, params['U'][0], self.J) print >>self.fh_info, '# Eimp=', Ed print >>self.fh_info, '# baths=%d'%baths, 'Ns=', Ns print >>self.fh_info, '# ntot=%f'%ntot print >>self.fh_info, '# nf=', nf print >>self.fh_info, '# moment=', moment print >>self.fh_info, '# sinfty=', sinf.tolist() self.fh_info.flush() _Af=[] _Gf=[] _Sig=[] _Delta=[] _eps1=[] _eps2=[] _a1=[] _a2=[] for b in range(baths): # over all components of spectral functions (over baths) nf_exact = nf[b]/Ns[b] dsinf = sinf[b]+Ed[b] # the limit of vanishing spectral weight, i.e., when we have only two poles # this formula should always be obeyed if (dsinf<0 and epsilon[0][0]*nf_exact+epsilon[1][0]*(1-nf_exact)>dsinf): print >> self.fh_info, 'epsilon was not choosen correctly and the solution can not be found!' epsilon[1][0] = (dsinf - epsilon[0][0]*nf_exact)/(1-nf_exact) - 0.5 print >> self.fh_info, 'setting epsilon[1][0] to ', epsilon[1][0] if (dsinf>0 and epsilon[0][1]*nf_exact+epsilon[1][1]*(1-nf_exact)> self.fh_info, 'epsilon was not choosen correctly and the solution can not be found!' epsilon[0][1] = (dsinf - epsilon[1][1]*(1-nf_exact))/nf_exact + 0.5 print >> self.fh_info, 'setting epsilon[0][1] to ', epsilon[0][1] Afo = array(Af[:,b]) m0 = dot(Afo,dh) m1 = dot(Afo,omdh) nc = dot(Afo,fedh) print >>self.fh_info, '#start [%d]: m0=%f m1=%f nc=%f nf_exact=%f' % (b, m0, m1, nc, nf_exact) a1 = nf_exact - nc if (a1<0): # weight below Ef is to large -> need to cut it a bit for im in range(len(om)): (a1n, a2n, m1n) = trycutAf(om[im], om, Afo, dh, omdh, fedh, nf_exact, params['Th'][0], self.fh_info) print >>self.fh_info, '#ca1 [%d]: cat L=%f a1=%f a2=%f m1=%f' % (b, om[im], a1n, a2n, m1n) if a1n>0: L1 = om[im] break (a1, a2, m0, m1, nc) = cutAf(L1, om, Afo, dh, omdh, fedh, nf_exact, params['Th'][0], self.fh_info) a2 = 1-m0-a1 if (a2<0): # lorentzian weight above Ef is to large -> need to cut it a bit for im in range(len(om)-1,0,-1): (a1n, a2n, m1n) = trycutAf(om[im], om, Afo, dh, omdh, fedh, nf_exact, params['Th'][0], self.fh_info) print >>self.fh_info, '#ca2 [%d]: cat L=%f a1=%f a2=%f m1=%f' % (b, om[im], a1n, a2n, m1n) if a2n>0: L2 = om[im] break (a1, a2, m0, m1, nc) = cutAf(L2, om, Afo, dh, omdh, fedh, nf_exact, params['Th'][0], self.fh_info) # Find u,v,w for this set of parameters (a1, a2, u, v, w) = uvw(m0, m1, nf_exact, nc, dsinf, epsilon) print >>self.fh_info, '# [%d]: miss-nf=%f miss-weight=%f s_oo=%f a1=%f a2=%f u=%f v=%f w=%f' % (b, nf_exact-nc, 1-m0, sinf[b], a1, a2, u, v, w) self.fh_info.flush() # Here we compute the positions of the two lorentzian # which will allow the exact density, normalization and sigma_oo (success, x, y) = Soluvw(u, v, w) # If is not possible to get exact density, normalization and sigma_oo by adding # two lorentzians # Will try to cut high-frequency spectral function Lb=None if (not success): # Here we cut the spectral function to make it more symmetric # start cutting at large frequency # cuts until the condition is met L0 = om[-1] # cut at positive frequency (success, a1, a2, x, y) = ww(L0, params['Th'][0], om, Afo, dh, omdh, fedh, nf_exact, epsilon, dsinf, self.fh_info) L0 /= 1.05 for j in range(100): (success, a1, a2, x, y) = ww(L0, params['Th'][0], om, Afo, dh, omdh, fedh, nf_exact, epsilon, dsinf, self.fh_info) if success: break L0 /= 1.05 (successn, a1n, a2n, xn, yn) = ww_cut(L0, params['Th'][0], om, Afo, dh, omdh, fedh, nf_exact, epsilon, dsinf, self.fh_info) if (not success): print "Can't determin a way to get exact nf, norm and sigma_oo. You have to figure out the way to do that!" sys.exit(1) print >> self.fh_info, "# [%d] a1=%f a2=%f x=%f y=%f" %(b, a1, a2, x, y) eps1_ = epsilon[0][0] + x*(epsilon[0][1]-epsilon[0][0]) eps2_ = epsilon[1][0] + y*(epsilon[1][1]-epsilon[1][0]) print >>self.fh_info, '# [%d]: a1=%f a2=%f eps1=%f eps2=%f' % (b, a1, a2, eps1_, eps2_) print >> self.fh_info # Actual cutting of the functions and adding the Lorentzians Afn = Afo # use cutted function for i in range(len(om)): # Adding the Lorentzians Afn[i] += -(a1/pi/(om[i]-eps1_+params['Gh'][0]*1j) + a2/pi/(om[i]-eps2_+params['Gh'][0]*1j)).imag Ac = array(Acd[:,b]) gf=[] sig=[] Delta=[] for i in range(len(om)): gr = -pi*krams.kramarskronig(Afn, om, i) gc = gr - pi*Afn[i]*1j deltar = -pi*krams.kramarskronig(Ac, om, i) delta = deltar - pi*Ac[i]*1j gf.append(gc) Delta.append(delta) sigm = om[i] - Ed[b] - 1/gc - delta - sinf[b] if sigm.imag > 0: sigm = sigm.real -1e-12j sig.append( sigm ) _Af.append(Afn) _Gf.append(gf) _Sig.append(sig) _Delta.append(Delta) _eps1.append(eps1_) _eps2.append(eps2_) _a1.append(a1) _a2.append(a2) print >>self.fh_info for b in range(baths): print >>self.fh_info, '## [%d]: eps1=%f eps2=%f a1=%f a2=%f' % (b, _eps1[b], _eps2[b], _a1[b], _a2[b]) self.fh_info.flush() shutil.move(self.dir+self.sparams['gloc'], self.dir+self.sparams['gloc']+'_o') shutil.move(self.dir+self.sparams['sig'], self.dir+self.sparams['sig']+'_o') shutil.move(self.dir+self.sparams['AlocOut'], self.dir+self.sparams['AlocOut']+'_o') fg = open(self.dir+self.sparams['gloc'], 'w') fs = open(self.dir+self.sparams['sig'], 'w') fa = open(self.dir+self.sparams['AlocOut'], 'w') fd = open(self.dir+'Delta.out', 'w') print >> fs, '# s_oo=', sinf.tolist() print >> fs, '# Edc=', Edc.tolist() print >> fa, first_line for i in range(len(om)): print >> fg, om[i], for b in range(baths): print >> fg, _Gf[b][i].real, _Gf[b][i].imag, print >> fg print >> fs, om[i], for b in range(baths): print >> fs, _Sig[b][i].real, _Sig[b][i].imag, print >> fs print >> fa, om[i], for b in range(baths): print >> fa, _Af[b][i], print >> fa print >> fd, om[i], for b in range(baths): print >> fd, _Delta[b][i].real, _Delta[b][i].imag, print >> fd fg.close() fs.close() fa.close() fd.close() shutil.copy2(self.dir+self.sparams['AlocOut'], self.dir+self.sparams['AlocOut']+'.'+extn) shutil.copy2(self.dir+self.sparams['AlocOut']+'_o', self.dir+self.sparams['AlocOut']+'_o.'+extn) #shutil.copy2(self.dir+self.sparams['sig'], self.dir+self.sparams['sig']+'.'+extn) shutil.copy2(self.dir+'nohup_imp.out', self.dir+'nohup_imp.out'+'.'+extn) return (_Sig, sinf, Edc, ntot) def printm(a): """ Prints matrix in readable form""" for i in range(shape(a)[0]): for j in range(shape(a)[1]): print "%8.4f " % a[i,j].real, print def ferm(x): """Fermi function for T=1""" if (x>100): return 0 if (x<-100): return 1 return 1/(exp(x)+1) def Sinftyv(U, J, l, baths, Ns, nf, da_): """ Computes Sigma_oo using OCA occupancies (density matrix given as a vector of numbers) Input: l -- quantum number l baths -- number of nonequivalent baths for OCA solver Ns[baths] -- degeneracy of each bath nf[baths] -- occupancy of each bath U, J -- Coulomb interaction Output: sig_oo[baths] -- hartree-fock value of self-energy using impurity occupancies Edc[baths] -- double counting using impurity occupancies """ # occupation vector from OCA solver # OCA gives only few numbers and we need to create vector - matrix with this size = sum(Ns) occ0=[] for b in range(baths): aocc = nf[b]/Ns[b] for j in range(Ns[b]): occ0.append(aocc) # becomes occupation matrix occj occj = zeros((size,size), dtype=complex,order='F') for mj in range(size): occj[mj,mj] = occ0[mj] # correlated index - all are correlated corind=ones(size,order='F',dtype=int) # index (atom,l,m,s) in spheric bndind=[] for s in range(1,3): for m in range(-l,l+1): bndind.append([1,l,m,s]) bndind = array(bndind,order='F',dtype=int) # index (atom,l,m,s) in relativistic bndindJ = [] for j in arange(l-0.5,l+0.5+1): if (j<0): continue for mj in arange(-j,j+1): bndindJ.append([1,l,j,mj]) bndindJ = array(bndindJ,order='F',dtype=float) # transformation from spheric to relativistic T2C = sts.cmp_t2j(bndindJ, bndind) # occupation - From relativistic -> spheric la.ztransform(occj, T2C,'C') # The code for Sigma_oo works in spheric harmonics Uc = ones(len(bndindJ))*U Jc = ones(len(bndindJ))*J (sinfty, Edc) = ldau.hartree(Uc,Jc,occj,bndind,corind,da=da_, subtractdc=0) la.ztransform(sinfty, T2C,'N') # From spheric -> relativistic # matrix of Sigma(infinity) is collapsed into few numbers # just like OCA occupation was given sinftyv=[] Edcv=[] ii=0 for b in range(baths): averS=0 averE=0 for j in range(Ns[b]): averS += sinfty[ii,ii].real averE += Edc[ii].real ii += 1 averS/=Ns[b] averE/=Ns[b] sinftyv.append(averS) Edcv.append(averE) return (array(sinftyv), array(Edcv)) def trycutAf(L, om, Af, dh, omdh, fedh, nf_exact, Th, fh_info): m0 = 0 m1 = 0 nc = 0 for i in range(len(om)): if (L<0): wcut = ferm(-(om[i]-L)/Th) else : wcut = ferm((om[i]-L)/Th) m0 += Af[i]*dh[i]*wcut m1 += Af[i]*omdh[i]*wcut nc += Af[i]*fedh[i]*wcut miss_dp = nf_exact-nc a1 = miss_dp a2 = 1-m0-a1 return (a1, a2, m1) def cutAf(L, om, Af, dh, omdh, fedh, nf_exact, Th, fh_info): m0 = 0 m1 = 0 nc = 0 for i in range(len(om)): if (L<0): wcut = ferm(-(om[i]-L)/Th) else : wcut = ferm((om[i]-L)/Th) Af[i] *= wcut m0 += Af[i]*dh[i] m1 += Af[i]*omdh[i] nc += Af[i]*fedh[i] miss_dp = nf_exact-nc a1 = miss_dp a2 = 1-m0-a1 return (a1, a2, m0, m1, nc) def Soluvw(u, v, w): """ Given positive numbers u and v, gives one solution to the equation x*u + y*v = w where x and y are numbers in the interval [0,1]. If there is no solution, gives False, if solution exists, gives the 'average' solution x = (x_min+x_max)/2 or y = (y_min+y_max)/2. """ if (w<0 or w>u+v): return (False, 0.0, 0.0) # Solution does not exists because u,v,x,y are postive if (w==0): return (True, 0.0, 0.0) x_max = (w0 and (w-v)0 and (w-u)u : x = 0.5*(x_max+x_min) y = (w-x*u)/v else: y = 0.5*(y_max+y_min) x = (w-y*v)/u return (True, x, y) def uvw(m0, m1, nf_exact, nc, dsinf, epsilon): a1 = nf_exact-nc if (a1<0): a1=0 a2 = 1-m0-a1 if (a2<0): a2=0 w = dsinf-m1-a1*epsilon[0][0]-a2*epsilon[1][0] u = a1*(epsilon[0][1]-epsilon[0][0]) v = a2*(epsilon[1][1]-epsilon[1][0]) return (a1, a2, u, v, w) def ww(L, Th, om, Af, dh, omdh, fedh, nf_exact, epsilon, dsinf, fh_info): m0 = 0 m1 = 0 nc = 0 for i in range(len(om)): wcut = ferm(-(L+om[i])/Th)*ferm((om[i]-L)/Th) m0 += Af[i]*dh[i]*wcut m1 += Af[i]*omdh[i]*wcut nc += Af[i]*fedh[i]*wcut (a1, a2, u, v, w) = uvw(m0, m1, nf_exact, nc, dsinf, epsilon) (succ, x, y) = Soluvw(u, v, w) print >> fh_info, '# L=%9.6f u=%9.6f v=%9.6f w=%9.6f x=%9.6f y=%9.6f m0=%9.6f m1=%9.6f a1=%9.6f a2=%9.6f dsinf=%9.6f' % (L, u, v, w, x, y, m0, m1, a1, a2, dsinf) return (succ, a1, a2, x, y) def ww_cut(L, Th, om, Af, dh, omdh, fedh, nf_exact, epsilon, dsinf, fh_info): m0 = 0 m1 = 0 nc = 0 for i in range(len(om)): wcut = ferm(-(L+om[i])/Th)*ferm((om[i]-L)/Th) Af[i] *= wcut m0 += Af[i]*dh[i] m1 += Af[i]*omdh[i] nc += Af[i]*fedh[i] (a1, a2, u, v, w) = uvw(m0, m1, nf_exact, nc, dsinf, epsilon) (succ, x, y) = Soluvw(u, v, w) print >> fh_info, '# L=%9.6f u=%9.6f v=%9.6f w=%9.6f x=%9.6f y=%9.6f m0=%9.6f m1=%9.6f a1=%9.6f a2=%9.6f dsinf=%9.6f' % (L, u, v, w, x, y, m0, m1, a1, a2, dsinf) return (succ, a1, a2, x, y)