#!/usr/bin/env python # Author: Kristjan Haule, March 2007 import os,sys,errno ROOT = os.environ.get('LMTO_DMFT_ROOT') if ROOT is not None: sys.path.append( ROOT ) sys.path.append( ROOT + '/fort') else: print >> sys.stderr, "Environment variable LMTO_DMFT_ROOT must be set!" print "Environment variable LMTO_DMFT_ROOT must be set!" sys.exit(1) #from ld_base import * from ldap import * # impurity solver from oca import * from ctqmc import * # specific imports to this module import shutil, copy import symeig # MPI start from mpi4py import MPI Master = 0 # for the version 0.5 MPI.rprint = MPI._rprint MPI.rank = MPI.WORLD_RANK MPI.size = MPI.WORLD_SIZE MPI.rprint('mpi.(root,size) = %d %d'%(Master,MPI.size), MPI.COMM_WORLD, Master) def omega_limits(iom, om): if (iom!=0): omega_a = 0.5*(om[iom]+om[iom-1]) else : omega_a = om[0] if (iomend): ii=end proc_limits[-1][1]=end return proc_limits def pr_limits1(start, end, MPI_size): """ Returns a 2D-list of points pnts[nprocessors][:] describing which points are to be computed on each processor""" pnts = [ [] for i in range(MPI_size)] n=start while n=end : break return pnts def MaxOmega(om, dmuLook): for i in range(len(om)): if om[i]>dmuLook: break return i def ComputeEigenvalues(lmt, om, Sigc, Sigma_oo, maxiom, fh_info, gamma=[0.0,0.0]): """ Routine computes eigenvalues of Hk+Sigma in parallel Input: lmt -- Lmtart class which contains all LMTO information om -- frequency Sigc[nom][nc] -- the columns of dynamical self-energy Sigma_oo[ndim,ndim] -- matrix of Sigma[infinity] maxiom -- integer smaller than len(om). Not all omega's are requeste, only [0,nom] gamma[2] -- broadening of the self-energy gamma[0] for non-correlated and gamma[1] for correlated Output: Et[nom,ndim,nkp] -- eigenvlaues """ print >> fh_info, 'Computing eigenvalues for all frequencies' proc_limits = pr_limits(-1,maxiom) (ndim,nkp) = (lmt.ndim, lmt.nkp) aEt=[] # Stores all eigenvalues for iom in range(proc_limits[MPI.rank][0],proc_limits[MPI.rank][1]): sig_DMFT_base = copy.deepcopy(Sigma_oo) if (iom>=0): # Adds Hartree-term and subtracts double-counting sig_DMFT_base += CreateSelfEnergyMatrix(Sigc[iom], lmt.Sigind, gamma) # Computes eigensystem only Et = array(zeros((ndim,nkp), dtype=complex), order='F') for ikp in range(nkp): # Rotates Sigma to this base sig2 = lmt.Embed(ikp, sig_DMFT_base) Et[:,ikp] = Eigenvaluesf(array(lmt.hamf[:,:,ikp]+sig2,order='F'), lmt.olap[:,:,ikp]) aEt.append(Et) print >> fh_info, '(iom,om,E0,En)= %3d %12.6f %12.6f %12.6f' % (iom, om[iom], Et[0,0], Et[-1,0]) fh_info.flush() aEt = array(aEt).flatten() # aEt[iom,ndim,nkp] data = MPI.WORLD.Allgather(aEt) return hstack(tuple(data)).reshape(maxiom+1,ndim,nkp) def ComputeChargeReal(mu_current, Egns, lmt, om, fh_info, noccb, max_metropolis_steps, use_tetra, LowerBound=-20): """ Routine computes weights for the expression 1/(om-Hk-Sigma) to be used for computing the chemical potential and electronic charge density. Input: Egns[nom,ndim,nkp] -- Eigenvalues of Hk+Sigma lmt_wk -- weights for all irreducible k-point om -- frequency mu_current -- the chemical potential LowerBound -- instead of -infinity. Where to start computing valence charge density Output: density -- density up to zero frequency """ om0=[] for im in range(len(om)): if (om[im]<0.0): om0.append(om[im]) else: break om0.append(-om[im-1]) # the last two frequencies should be at [-x,x] such that the intergation is up to zero. proc_limits = pr_limits(-1, len(om0)-1) # Stores all tetrahedron weights int_wghs=[] dens = 0 for iom in range(proc_limits[MPI.rank][0],proc_limits[MPI.rank][1]): if iom==-1 : omab = [LowerBound, om0[0]] elif iom==0: omab = [om0[0], 0.5*(om0[0]+om0[1])] else: omab = [0.5*(om0[iom-1]+om0[iom]), 0.5*(om0[iom]+om0[iom+1])] if (use_tetra): int_wgh = tetra.tetra_zweights_int_only(omab[0]+mu_current+1e-10j, omab[1]+mu_current+1e-10j, Egns[iom+1], lmt.itt, max_metropolis_steps) else: int_wgh = tetra.cmp_simple_integral_weights(omab[0]+mu_current, omab[1]+mu_current, Egns[iom+1], lmt.wk) #int_wgh = tetra.cmp_simple_integral_weights(omab[0]+mu_current, omab[1]+mu_current, Egns[iom+1], lmt.wk) # try Egns[iom+1] only int_wgh.real=0 int_wghs.append(int_wgh) dens += sum(int_wgh).imag*(-1./pi) # print >> fh_info, '(iom,om_a,om_b,TOS)= %3d %12.6f,%12.6f %20.12f' % (iom, omab[0], omab[1], dens) # fh_info.flush() data = MPI.WORLD.Allgather(dens) dens = sum(data) print >> fh_info, mu_current, dens print 'mu, dens, noccb=', mu_current, dens, noccb return dens-noccb def ComputeChargeImag(mu_current, Egns, lmt, om, lom, fh_info, noccb, com=0): """ Routine computes weights for the expression 1/(om-Hk-Sigma) to be used for computing the chemical potential and electronic charge density. Input: Egns[nom,ndim,nkp] -- Eigenvalues of Hk+Sigma lmt_wk -- weights for all irreducible k-point om -- frequency mu_current -- the chemical potential com -- if com=0, we subtract 1/(iom-Hk-Sigma_oo); if com=1, we subtract 1/(iom-Hk-Sigma(0)) Output: density -- density up to zero frequency First we define E0[p,ikp] = Eig[omega=inf,p,ikp]. We compute then 2*T*sum_{iom,k,p}[1/(iom+mu-Eig[om,p,ikp])-1/(iom+mu-E0[p,ikp])] + \sum_{k,p}ferm(E0[p,ikp]-mu) """ T = lom[0]/pi # temperature can be found from lowest Matsubara frequency # Energy which will be sutracted E0 = zeros(shape(Egns[0])) for p in range(len(E0)): for ikp in range(len(E0[p])): E0[p,ikp] = Egns[com,p,ikp] # Parallel summation proc_limits = pr_limits(1, len(om)+1) fdens=[] for im in range(proc_limits[MPI.rank][0],proc_limits[MPI.rank][1]): #for im in range(1,len(om)+1): z = om[im-1]*1j + mu_current sm=0 for p in range(lmt.ndim): for ikp in range(lmt.nkp): gc = 1/(z-Egns[im,p,ikp]) - 1/(z-E0[p,ikp]) sm += lmt.wk[ikp]*gc.real fdens.append(sm) data = MPI.WORLD.Allgather(fdens) fdens = reduce(lambda x, y: x + y, data) # interpolate on big mesh tckr = interpolate.splrep(om, fdens, s=0) lfdens = interpolate.splev(lom, tckr) # sum over the big mesh to get density dsum = 0 for im in range(len(lom)): dsum += lfdens[im] nt1 = 2*T*dsum # Add the fermi functions nt0=0 for ikp in range(lmt.nkp): for p in range(lmt.ndim): nt0 += ferm((E0[p,ikp]-mu_current)/T)*lmt.wk[ikp] dens = nt0+nt1 print 'mu, dens, noccb=', mu_current, dens, noccb return dens-noccb def FindSignChange(fComputeCharge, old_mu, sdmu, args): """ Routine brackets the chemical potential. Input: old_mu -- start of looking for mu sdmu -- the interval will be looked in the steps of sdmu in the following algorithm: 1) try [mu, mu+sdmu] 2) try [mu+sdmu, mu+2*sdmu] 3) try [mu+2*sdmu, mu+2*2*sdmu] .... args contains typically: Egns[nom,ndim,nkp] -- eigenvalues lmt -- the lmt class with information about the k-points and tetrahedra om -- frequency fh_info -- info-file noccb -- looking for a par of [mu0, mu1] such that N(mu0)> fh_info, 'Looking for the chemical potential. Starting from old_mu=', old_mu curr_mu = old_mu curr_dens = apply(fComputeCharge, (curr_mu,)+args) #Egns, lmt, om, fh_info, noccb, max_metropolis_steps, use_tetra, LowerBound) print >> fh_info, '(mu,dens)=', curr_mu, curr_dens if abs(curr_dens)<1e-6: dtmu = dtmu=1e-4 return (curr_mu-dtmu, curr_mu+dtmu) elif (curr_dens<0): tdmu = sdmu else: tdmu = -sdmu while (True): new_dens = apply(fComputeCharge, (curr_mu+tdmu,)+args) #Egns, lmt, om, fh_info, noccb, max_metropolis_steps, use_tetra, LowerBound) print >> fh_info, '(mu,dens)=', curr_mu+tdmu, new_dens if curr_dens*new_dens<0: break curr_dens = new_dens curr_mu += tdmu tdmu *= 2. print >> fh_info, '(mu0,mu1), (dens0,dens1)', curr_mu, curr_mu+tdmu, curr_dens, new_dens return (curr_mu, curr_mu+tdmu) def DensityMatrixCoeffInside_new(omab, lmt, sig_DMFT_base, use_tetra, max_metropolis_steps, symmetrize=True): """ Computes coefficients for charge density (nhh,nhj,njh,njj) and new density matrix at given frequency Input: lmt int_wgh -- tetrahedron weight for this frequency interval sig_DMFT_base -- self-energy symmetrize -- whether to symmetrize the gxx results Output: (nhh,nhj,njh,njj) -- integrated coefficients mocc -- density matrix """ (ndim,nkp) = (lmt.ndim, lmt.nkp) # Computes eigensystem Et = array(zeros((ndim,nkp), dtype=complex), order='F') Al = array(zeros((ndim,ndim,nkp), dtype=complex), order='F') Ar = array(zeros((ndim,ndim,nkp), dtype=complex), order='F') for ikp in range(nkp): # Rotates Sigma to this base sig2 = lmt.Embed(ikp, sig_DMFT_base) (Et[:,ikp], Al[:,:,ikp], Ar[:,:,ikp]) = Eigensystemf(array(lmt.hamf[:,:,ikp]+sig2,order='F'), lmt.olap[:,:,ikp]) if (use_tetra): int_wgh = tetra.tetra_zweights_int_only(omab[0], omab[1], Et, lmt.itt, max_metropolis_steps) else: int_wgh = tetra.cmp_simple_integral_weights(omab[0], omab[1], Et, lmt.wk) int_wgh.real=0 ddens = sum(int_wgh).imag*(-1./pi) # Computes matrix elements for the local green's function gxx = crho.cmp_rho_coefficients(lmt.inpdir, int_wgh, Al, Ar, lmt.bndind, lmt.lbndind) if (symmetrize): lmt.symmetrize(gxx) nxx = lmt.give_density_coeff(gxx) moccn = crho.cmp_occupancy(nxx[0], nxx[1], nxx[2], nxx[3], lmt.moo, lmt.bndind, lmt.isrt, lmt.llmtind) if (lmt.norbs>1 and lmt.add_sqrt2): for i in range(len(nxx)): n1 = shape(nxx[i])[0]/2 n2 = shape(nxx[i])[1]/2 nxx[i][0:n1, n2:2*n2] *= sqrt(2.) nxx[i][n1:2*n1, 0:n2] *= -sqrt(2.) return (nxx, moccn, ddens) def ComputeDensityMatrixCoeff_new(mu_current, lmt, om, Sigc, Sigma_oo, fh_info, noccb, gamma, max_metropolis_steps, use_tetra, LowerBound=-20): """ Computes coefficients for charge density (nhh,nhj,njh,njj) and new density matrix - integral over frequency Input: mu_current -- chemical potential lmt -- om -- frequency Sigc -- vector of self-energy Sigma_oo -- sigma_{infinity} fh_info -- some information noccb -- total valence charge gamma[2] -- broadening of the self-energy gamma[0] for non-correlated and gamma[1] for correlated max_metropolis_steps -- metropolis steps for tetrahedron method use_tetra -- using tetra or not Output: (nhh,nhj,njh,njj) -- integrated coefficients mocc -- density matrix """ #noccb = lmt.zvaltot # Need this charge print >> fh_info, 'Computing Density Matrix Coefficients' om0=[] for im in range(len(om)): if (om[im]<0.0): om0.append(om[im]) else: break om0.append(-om[im-1]) # the last two frequencies should be at [-x,x] such that the intergation is up to zero. proc_limits = pr_limits(-1, len(om0)-1) dens=0 mocc_new = zeros((lmt.ndim,lmt.ndim),dtype=complex,order='F') nxx=[zeros((lmt.ndim,lmt.ndim),dtype=complex,order='F'),zeros((lmt.ndim,lmt.lndim),dtype=complex,order='F'),zeros((lmt.lndim,lmt.ndim),dtype=complex,order='F'),zeros((lmt.lndim,lmt.lndim),dtype=complex,order='F')] for iom in range(proc_limits[MPI.rank][0],proc_limits[MPI.rank][1]): if iom==-1 : omab = array([LowerBound, om0[0]]) +mu_current+1e-10j elif iom==0: omab = array([om0[0], 0.5*(om0[0]+om0[1])]) +mu_current+1e-10j else: omab = array([0.5*(om0[iom-1]+om0[iom]), 0.5*(om0[iom]+om0[iom+1])]) +mu_current+1e-10j sig_DMFT_base = copy.deepcopy(Sigma_oo) if (iom>=0): # Adds Hartree-term and subtracts double-counting sig_DMFT_base += CreateSelfEnergyMatrix(Sigc[iom], lmt.Sigind, gamma) (nxxn, moccn, ddens) = DensityMatrixCoeffInside_new(omab, lmt, sig_DMFT_base, use_tetra, max_metropolis_steps, symmetrize=True) dens += ddens for i in range(len(nxx)): nxx[i] += nxxn[i] mocc_new += array(moccn) print >> fh_info, '(iom,om,TOS,TOS)= %3d %16.10f %20.12f %20.12f' %(iom, om[iom], dens, trace(mocc_new)) fh_info.flush() mocc_new = MPI.WORLD.Allreduce(mocc_new, MPI.SUM) for i in range(len(nxx)): nxx[i] = MPI.WORLD.Allreduce(nxx[i], MPI.SUM) tdenso = trace(mocc_new).real renorm = noccb/tdenso mocc_new *= renorm tdens = trace(mocc_new).real print '(Density1,Density_new,Density_old,renormalization)=', dens, tdens, tdenso, renorm.real print >> fh_info, '(Density1,Density_new,Density_old,renormalization)=', dens, tdens, tdenso, renorm.real return (nxx, mocc_new) def cmp_Dos_Gloc_LMTO(lmt, omega_mu, sig_DMFT_base, max_metropolis_steps): """ Computes partial DOS and gloc in LMTO base - hamf and olap need to be given in LMTO base Input: lmt -- lmtart class omega_mu -- om+mu olap[ndim,ndim,nkp] -- overlap matrix in DMFT base hamf[ndim,ndim,nkp] -- hamiltonian matrix in DMFT base Sigma[ndim,ndim] -- self-energy matrix max_metropolis_steps -- metropolis steps in sorting bands --------------------------------------------- The following members of class lmt are used: symmetrize -- symmetrization in spherics cmp_Gloc_coeff cmp_tgf2 give_density_coeff ndim -- size of the ham matrix Output: gloc[ndim,ndim] -- local green's function dDOS[ndim] -- partial DOS """ #sig_DMFT_base = Sigma_dyn + Sigma_oo # Total self-energy in DMFT base gxx = lmt.cmp_Gloc_coeff(omega_mu, sig_DMFT_base, tetra.tetra_weights, lmt.bndind, max_metropolis_steps) # coefficienst ghh,ghj,gjh,gjj lmt.symmetrize(gxx) # to simulate summation over the entire Brillouin zone gloc = lmt.cmp_tgf2(gxx) # the local green's function in DMFT base nxx = lmt.give_density_coeff(gxx) # coefficients for DOS dDOS2 = crho.cmp_occupancy(nxx[0], nxx[1], nxx[2], nxx[3], lmt.moo, lmt.bndind, lmt.isrt, lmt.lmtind) # partial dos dDOS = zeros(lmt.ndim,dtype=float) for p in range(lmt.ndim): dDOS[p] = dDOS2[p,p].real return (gloc, dDOS) def cmp_Dos_Gloc_DMFT(tetrahedronQ, lmt, omega_mu, olap, hamf, sig_DMFT_base, max_metropolis_steps, cut_ab=[-100,100]): """ Computes partial DOS and gloc in DMFT base - hamf and olap need to be given in DMFT base Input: lmt -- lmtart class omega_mu -- om+mu olap[ndim,ndim,nkp] -- overlap matrix in DMFT base hamf[ndim,ndim,nkp] -- hamiltonian matrix in DMFT base Sigma[ndim,ndim] -- self-energy matrix max_metropolis_steps -- metropolis steps in sorting bands --------------------------------------------- The following members of class lmt are used: itt[5,ntet] -- index of tetrahedron mesh T2C[ndim,ndim] -- matrix which rotates from spheric to cubis/relativistic coordinates ndim -- size of the ham matrix Output: gloc[ndim,ndim] -- local green's function dDOS[ndim] -- partial DOS """ if tetrahedronQ: # Actually computes local Green's function - sum over irreducible BZ #Sigma = Sigma_dyn + Sigma_oo gloc = cmp_Gloc(omega_mu, olap, hamf, sig_DMFT_base, lmt.itt, max_metropolis_steps, cut_ab[0], cut_ab[1]) # Array has to be converted to Fortran type array to symmetrize gloc = array(gloc, order='F') else: ############ gloc = zeros((lmt.ndim,lmt.ndim),dtype=complex,order='F') for ikp in range(lmt.nkp): gloc += linalg.inv(omega_mu*olap[:,:,ikp]-hamf[:,:,ikp]-sig_DMFT_base)*lmt.wk[ikp] ############ # Symmetrization replaces summation over full 1-st BZ sym.symmetrize_cubic(gloc, lmt.T2C) dDOS = zeros(lmt.ndim,dtype=float) for p in range(lmt.ndim): dDOS[p] = -gloc[p,p].imag/pi return (gloc, dDOS) def cmp_Eimp(mu, olap, hamf, Sigma_oo_lattice, Sigma_oo_IMP, lmt, base, fh_info, printQ=True): if (base=='DMFTbase'): oho = crho.cmp_oho(olap,hamf,lmt.wk) # O^{-1}*Hk*O^{-1} oso = crho.cmp_oso(olap,Sigma_oo_lattice,lmt.wk) # O^{-1}*Sigma_{oo}*O^{-1} elif (base=='LMTObase'): oho = crho.cmp_ohol(olap,hamf,lmt.Utk,lmt.wk) # O^{-1}*Hk*O^{-1} oso = crho.cmp_osol(olap,Sigma_oo_lattice,lmt.Utk,lmt.wk) # O^{-1}*Sigma_{oo}*O^{-1} !maybe bug! Sigma_oo is in DMFT base else: sys.stder('This base not implemented') sys.exit(1) sym.symmetrize_cubic(oho, lmt.T2C) sym.symmetrize_cubic(oso, lmt.T2C) Eimp_matrix = oho - mu*identity(lmt.ndim) + oso - Sigma_oo_IMP # Eimp of the correlated block is stored below # We need information about the nonzero entries in the matrix sind = lmt.Sigind.flatten() cols = utl.union(sind) deg = utl.repeat(sind,cols) ohov = zeros(max(cols),dtype=float) osov = zeros(max(cols),dtype=float) Eimp = zeros(max(cols),dtype=float) for p in range(lmt.ndim): for q in range(lmt.ndim): if (lmt.Sigind[p,q]>0): ii = lmt.Sigind[p,q]-1 Eimp[ii] += Eimp_matrix[p,q] ohov[ii] += oho[p,q] osov[ii] += oso[p,q] for p in range(len(Eimp)): Eimp[p] /= deg[p+1] ohov[p] /= deg[p+1] osov[p] /= deg[p+1] Eimps=[] for impr in lmt.Impr: Eimps.append([Eimp[r] for r in impr]) if (printQ): print 'oho=', ohov, 'mu=', mu if (printQ): print 'oso=', osov #if (printQ): print 'Eimps=', Eimps return (Eimp, Eimp_matrix, Eimps) def Give_Dos_Gloc_Delta(gloc, dDOS, lmt, omega, Eimp_matrix, Sigma_imp): """ Computes partial DOS, Delta and correlated gf Input: lmt -- lmtart class omega_mu -- om+mu Sigma[ndim,ndim] -- self-energy matrix gloc[ndim,ndim] -- matrix of local green's function dDOS[ndim] -- partial DOS --------------------------------------------- The following members of class lmt are used: itt[5,ntet] -- index of tetrahedron mesh T2C[ndim,ndim] -- matrix which rotates from spheric to cubis/relativistic coordinates bndindx[ndim,4] -- index from orbital to atom and l ndim -- size of the ham matrix Sigind[ndim,ndim] -- correlated index matrix Output: tDOS -- total DOS pDOS[natom,nmax] -- partial DOS gc -- columns of the Green's function of the correlated orbitals Delta -- hybridization function of the correlated orbitals """ #### VERY IMPORTANT: NEEDS TO ZERO non-f GF lmt.ZeroNonCorrelated(gloc) MDelta = omega*identity(lmt.ndim) -Eimp_matrix - linalg.inv(matrix(gloc)) - matrix(Sigma_imp) # Partial DOS is computed for each atom and each l # We need information about lmax and natom to do that # Information is obtained from lmt.bndindx lmax = int(max([x[1] for x in lmt.bndindx])) natom = int(max([x[0] for x in lmt.bndindx])) # Partial DOS is stored below dos_tot=0 pdos = zeros((natom,lmax+1),dtype=float) for p in range(lmt.ndim): iatom = int(lmt.bndindx[p][0]-1) il = int(lmt.bndindx[p][1]) pdos[iatom, il] += dDOS[p] dos_tot += dDOS[p] # Green's function of the correlated block is stored below # We need information about the nonzero entries in the matrix sind = lmt.Sigind.flatten() cols = utl.union(sind) deg = utl.repeat(sind,cols) gc = zeros(max(cols),dtype=complex) Delta = zeros(max(cols),dtype=complex) for p in range(lmt.ndim): for q in range(lmt.ndim): if (lmt.Sigind[p,q]>0): ii = lmt.Sigind[p,q]-1 gc[ii] += gloc[p,q] Delta[ii] += MDelta[p,q] for p in range(len(gc)): gc[p] /= deg[p+1] Delta[p] /= deg[p+1] return (gc, Delta, dos_tot, pdos) def Print_Dos_Gloc_Delta(outdir, om, gcDeDo): fh_dos = open(outdir+'/dos.out','w') fh_gf = open(outdir+'/gf.out', 'w') fh_delta = open(outdir+'/Delta.out', 'w') for iom in range(len(om)): fh_dos.write("%f \t%f " % (om[iom], gcDeDo[iom][2])) pdos = gcDeDo[iom][3] for iatom in range(shape(pdos)[0]): for il in range(shape(pdos)[1]): fh_dos.write("\t%f " % (pdos[iatom,il]) ) fh_dos.write("\n") fh_gf.write("%f " % (om[iom])) for g in gcDeDo[iom][0]: fh_gf.write("\t%f %f " % (g.real, g.imag)) fh_gf.write("\n") fh_delta.write("%f " % (om[iom])) for g in gcDeDo[iom][1]: fh_delta.write("\t%f %f " % (g.real, g.imag)) fh_delta.write("\n") fh_dos.close() fh_gf.close() fh_delta.close() #BRISI def Read_Dos_Gloc_Delta(outdir): fh_dos = open(outdir+'/dos.out','r') fh_gf = open(outdir+'/gf.out', 'r') fh_delta = open(outdir+'/Delta.out', 'r') om=[] tdos=[] pdos=[] for line in fh_dos: data = map(float,line.split()) om.append(data[0]) tdos.append(data[1]) pdos.append(array(data[2:]).reshape(2,4)) gc=[] for line in fh_gf: data = map(float,line.split()) tgc=[] for i in range((len(data)-1)/2): tgc.append(data[2*i+1]+data[2*i+2]*1j) gc.append(tgc) Delta=[] for line in fh_delta: data = map(float,line.split()) tdelta=[] for i in range((len(data)-1)/2): tdelta.append(data[2*i+1]+data[2*i+2]*1j) Delta.append(tdelta) gcDeDo=[] for iom in range(len(om)): gcDeDo.append([gc[iom], Delta[iom], tdos[iom], pdos[iom]]) return gcDeDo #BRISI def Cmp_Dos_Gloc_Delta(fh_info, tetrahedronOm, imag, om, mu, lmt, Sigc, Sigma_oo_lattice, Sigma_oo_IMP, Eimp_matrix, p_base, p_max_metropolis_steps, p_cut_ab, p_gamma): """ """ #proc_limits = pr_limits(0, len(om)) proc_limits = pr_limits1(0, len(om), MPI.size) print >> fh_info, 'proc_limits=', proc_limits #gcDeDo=[] # contains (gc, Delta, dos_tot, pdos) #for iom in range(proc_limits[MPI.rank][0],proc_limits[MPI.rank][1]): #gcDeDo = [[] for i in proc_limits[MPI.rank]] gdata=[] for im,iom in enumerate(proc_limits[MPI.rank]): if imag: omega = om[iom]*1j else: omega = om[iom] if iom> fh_info, 'DOS= ', om[iom], gcDeDo[-1][2] #gcDeDo[im] = (iom, Give_Dos_Gloc_Delta(gloc, dDOS, lmt, omega, Eimp_matrix, Sigma_imp)) gdata.append((iom, Give_Dos_Gloc_Delta(gloc, dDOS, lmt, omega, Eimp_matrix, Sigma_imp))) #print >> fh_info, 'DOS= ', om[iom], gcDeDo[im][1][2] print >> fh_info, 'DOS= ', om[iom], gdata[-1][1][2] fh_info.flush() #data = MPI.WORLD.Allgather(gcDeDo) data = MPI.WORLD.Allgather(gdata) # sort resulting data gcDeDo = [[] for i in range(len(om))] for dp in data: for d in dp: iom = d[0] gcDeDo[iom] = d[1] #reduce(lambda x, y: x + y, data) return gcDeDo def MakeDiagMatrix(a): A=zeros((len(a),len(a)),dtype=complex,order='F') for i in range(len(a)): A[i,i]=a[i] return A def first(x): if type(x)==int : return x else: return x[0] def WriteSigma(filename, om, Sig, sinf, Edc): fs = open(filename, 'w') print >> fs, '# s_oo=', sinf.tolist() print >> fs, '# Edc=', Edc.tolist() for i in range(len(om)): print >> fs, om[i], for b in range(len(Sig[i])): print >> fs, Sig[i][b].real, Sig[i][b].imag, print >> fs fs.close() def create_log_mesh(om, Sig, nom, ntail_): """Creates logarithmic mesh on Matsubara axis Takes first istart points from mesh om and the rest of om mesh is replaced by ntail poinst redistribued logarithmically. Input: om -- original long mesh istart -- first istart points unchanged ntail -- tail replaced by ntail points only Output: som -- smaller mesh created from big om mesh sSig[nom,nc] -- Sig on small mesh Also computed but not returned: ind_om -- index array which conatins index to kept Mastubara points """ istart = min(nom, len(om)) ntail = min(ntail_, len(om)-istart) istart = min(nom,len(om)) ntail = min(ntail, len(om)-istart) ind_om=[] alpha = log((len(om)-1.)/istart)/(ntail-1.) for i in range(istart): ind_om.append(i) for i in range(ntail): t = int(istart*exp(alpha*i)+0.5) if (t != ind_om[-1]): ind_om.append(t) som = zeros(len(ind_om), dtype='float') sSig = zeros((len(som), shape(Sig)[1]), dtype=complex) for i in range(len(ind_om)): som[i] = om[ind_om[i]] sSig[i,:] = Sig[ind_om[i]][:] return (som, sSig) def reorthogonalize(olap, hamf): nkp = shape(lmt.olap)[2] ndim = shape(lmt.olap)[0] print 'nkp=', nkp for ikp in range(nkp): eig1 = symeig.symeig(olap[:,:,ikp]) T = matrix(eig1[1]) D = matrix(zeros(shape(T),dtype=complex)) for i in range(len(D)): D[i,i] = 1/sqrt(abs(eig1[0][i])) Ut = T*D*conjugate(T.transpose()) # This is 1/Sqrt(Overlap) hamf[:,:,ikp] = Ut*hamf[:,:,ikp]*Ut # Orthogonalizing Hk ones more on imaginary axis # Finally, set overlap to unity for ip in range(ndim): for iq in range(ndim): if ip==iq : olap[ip-1,iq-1,ikp]=1.0 else: olap[ip-1,iq-1,ikp]=0.0 if __name__ == '__main__': # -------------- Parameters which are set through input files sparams.dat and params.dat ----------- # Default values of some parameters params = {'inpdir': '.', 'scratch': False, 'cix': None, 'start': 0, 'UpdateAtom': False, 'U': 0.0, 'T': 0.03, 'Niter': 1, 'max_iterations': 1000, 'finish': 100, 'admixDM': 0.5, 'admix_large':0.3, 'broyden': True, 'admix0': 0.3, 'admix_rho': 1., 'admix_mag': 1., 'restBrdn': 15, 'om_ab': [-10,10], 'Nom': 200, 'mix_ham': 1., 'admix_frozen': 1., 'max_metropolis_steps': 0, 'use_tetra' : 1, 'LowerBound' : -20, 'a' : 0.0, 'cut_ab' : [-100,100], 'base' : '"DMFTbase"', 'runLDA' : True, 'runIMP' : True, 'gamma' : [0.0, 0.0], 'gammag' : [0.0, 0.0], 'gbroad' : 0.01, 'recompute_mu': True, 'sdmu' : 0.1, 'mix_mu': 1., 'solver': 'OCA', #'solver': 'CTQMC', 'DoubleCounting': '"impurity"', 'com' : 0, # when computing density on imaginary axis, we can subtract Sigma_oo (com=0) or Sigma(0) (com=1) 'tetrahedronOm':1000, 'orthog':False, #'DCs' : 'default', #'DCs' : 'fixn', #'DCs' : 'fixEimp', } p = Params('sparams.dat', 'params.dat', params) fh_info = open(os.path.join(p['inpdir'], 'phi.info.'+str(MPI.rank)), 'w') print >> fh_info, '-'*80 print >> fh_info, '*'*20, 'Iteration', 1, '*'*20 print >> fh_info, '-'*80, '\n' p.refresh(fh_info) dEdc = 0.0 if (MPI.rank==Master): fh_pinfo = open(os.path.join(p['inpdir'], 'info.iterate'), 'a') #print >> fh_pinfo, params for k in p.keys(): print >> fh_pinfo, k,'=', p[k] fh_pinfo.write('%3s.%3s %12s %12s %12s %12s\n' % ('#','#','mu','Eimp','Edc','nf')) fh_pinfo.flush() dir_imp = 'imp' # LmtArt is executed dir_ksum = 'ksum' if MPI.rank==Master: # Prepares for lda execution ldaexe = LDApU(p['inpdir'], dir_ksum, fh_info, p['scratch']) if (p['runLDA']): # Runs the starting lda ldaexe.RunFirstLDA() dir_ksum = ldaexe.ksum MPI.WORLD.Barrier() print >> fh_info, "reading lmtart data" fh_info.flush() # allocated k-sum memory and reads the lmtart data lmt = Lmtart(os.path.join(p['inpdir'], dir_ksum), fh_info) if MPI.rank==Master: lmt.check_cix(p['cix']) MPI.WORLD.Barrier() lmt.read_cix() print >> fh_info, "lmtart data read." fh_info.flush() if p['solver'] in ['OCA', 'NCA']: imag = False elif p['solver'] in ['CTQMC']: imag = True else: print 'ERROR: No solver specified! Bolling out.' atoms = [first(x) for x in p['cix'].keys()] Znucs = [int(lmt.znuc[lmt.isrt[iatom-1]-1]) for iatom in atoms] print 'Znuc=', Znucs solver_J=None if MPI.rank==Master: mpi_prefix='' if len(glob.glob('mpi_prefix.dat'))!=0: mpi_prefix = open('mpi_prefix.dat', 'r').next().strip() if p['solver'] in ['OCA', 'NCA']: solver = [ IMP_OCA(Znucs[im], dir_imp+'.'+str(im), lmt.Impr[im], p['iparams'+str(im)], p['UpdateAtom']) for im in range(len(lmt.Impc))] elif p['solver'] in ['CTQMC']: solver = [ IMP_CTQMC(Znucs[im], dir_imp+'.'+str(im), lmt.Impr[im], p['iparams'+str(im)], p['UpdateAtom']) for im in range(len(lmt.Impc))] else: print 'ERROR: No solver specified! Bolling out.' solver_J = [isolver.J for isolver in solver] solver_J = MPI.WORLD[Master].Bcast(solver_J) # BUG !!! UC = zeros(lmt.ndim) JC = zeros(lmt.ndim) for im in lmt.Impw.keys(): for q in lmt.Impw[im]: UC[q]= p['iparams'+str(im)]['U'][0] JC[q]= solver_J[im] #print 'UC=', UC.tolist() #print 'JC=', JC.tolist() lmt_corr_index = lmt.Correlated_LS_index() # Checks if there is an old density matrix available / need for Sigma_oo (StartExist, mocc) = ReadMocc(os.path.join(p['inpdir'],"mocc.out"), lmt.ndim) # Checks if there is an old self-energy available fsiginp = p['inpdir']+"/sig.inp" if (os.path.exists(fsiginp)): if (MPI.rank==Master): print 'Found input self-energy ', fsiginp #(om,Sigc) = ReadSelfEnergyColumns(fsiginp) else: # Creates zero as start om = linspace(p['om_ab'][0], p['om_ab'][1], p['Nom']) Sigc = zeros((p['Nom'],lmt.Sigind.max())) if (MPI.rank==Master): print 'Could not find input self-energy ', fsiginp, '. Starting from default zero!' fs = open(fsiginp, 'w') for im in range(len(om)): print >> fs, om[im], for j in range(len(Sigc[im])): print >> fs, Sigc[im,j].real, Sigc[im,j].imag, print >> fs fs.close() MPI.WORLD.Barrier() old_mu = lmt.mu mu = old_mu print >> fh_info, "start iteration" fh_info.flush() noccb = int(lmt.zvaltot*lmt.nspin/2.+1e-5) # Need this charge for itt in range(p['max_iterations']): print >> fh_info, "iteration number", itt fh_info.flush() if (itt!=0 or StartExist): #p_a = 0 #if p['DCs']=='default' : p_a = p['a'] # BUG !!! p['U'] DOES NOT EXIST! IT IS PROPERTY OF THE CORRELATED ATOM, NOT THE SYSTEM. (Sigma_oo_LDA,Edc_LDAv) = ldau.hartree(UC, JC, mocc, lmt.bndind, lmt_corr_index, p['a'], subtractdc=0) Edc_LDA = MakeDiagMatrix(Edc_LDAv) la.ztransform(Sigma_oo_LDA, lmt.T2C,'N') # From spheric -> cubic la.ztransform(Edc_LDA, lmt.T2C,'N') # From spheric -> cubic if (MPI.rank==Master): WriteMatrix('occup.dat', mocc) WriteMatrix('sigm_oo.dat', Sigma_oo_LDA) WriteMatrix('Edc_lda.dat', Edc_LDA) else: Sigma_oo_LDA = zeros((lmt.ndim, lmt.ndim), dtype=complex, order='F') Edc_LDA = zeros((lmt.ndim, lmt.ndim), dtype=complex, order='F') for i in range(p['Niter']): MPI.WORLD.Barrier() p.refresh(fh_info) # update parameters (om,Sigc,Sigma_oo_IMPv,Edc_IMP) = ReadSelfEnergyColumns(fsiginp, contains_DC=True) # New self-energy Sigma_oo_IMP = CreateSelfEnergyMatrix(Sigma_oo_IMPv, lmt.Sigind) ##### if imag: lom = om[:] lSigc = copy.deepcopy(Sigc) (om, Sigc) = create_log_mesh(om, Sigc, p['iparams0']['nom'][0], p['ntail']) print >> fh_info, "Small mesh of size", len(om), "frequency points created" ##### if (p['DoubleCounting']=='LDA'): Sigma_oo_lattice = Sigma_oo_LDA-Edc_LDA else: Sigma_oo_lattice_v = Sigma_oo_IMPv-Edc_IMP Sigma_oo_lattice = CreateSelfEnergyMatrix(Sigma_oo_lattice_v, lmt.Sigind) # this corrects Edc if one wants to fix or slowly mix impurity levels # one needs to keep high-frequency expansion correct at each iteration if p['DCs']=='fixEimp': Sigma_oo_lattice += dEdc*identity(len(Sigma_oo_lattice)) if (MPI.rank==Master): WriteMatrix('Sigma_oo_lattice', Sigma_oo_lattice) if (p['base']=='DMFTbase'): lmt.read_ham() # Reads hamiltonian and overlap again if necessary Egns=None if p['recompute_mu']: # Eigenvalues for all frequencies if imag: maxomega = len(om) else: maxomega = MaxOmega(om, 0.0) Egns = ComputeEigenvalues(lmt, om, Sigc, Sigma_oo_lattice, maxomega, fh_info, p['gamma']) ####################### if imag: # bracketing zero (mu0, mu1) = FindSignChange(ComputeChargeImag, old_mu, p['sdmu'], args=(Egns, lmt, om, lom, fh_info, noccb, p['com'])) # find the chemical potential new_mu = optimize.brentq(ComputeChargeImag, mu0, mu1, args=(Egns, lmt, om, lom, fh_info, noccb, p['com'])) else: # bracketing zero (mu0, mu1) = FindSignChange(ComputeChargeReal, old_mu, p['sdmu'], args=(Egns, lmt, om, fh_info, noccb, p['max_metropolis_steps'], p['use_tetra'], p['LowerBound'])) # find the chemical potential new_mu = optimize.brentq(ComputeChargeReal, mu0, mu1, args=(Egns, lmt, om, fh_info, noccb, p['max_metropolis_steps'], p['use_tetra'], p['LowerBound'])) ################ dmu = new_mu-old_mu if itt==0 and i==0: mu = new_mu else: mu = old_mu*(1-p['mix_mu']) + new_mu*p['mix_mu'] if (MPI.rank==Master): print 'New chemical potential found at ', dmu, ' -> Chemical potential becomes ', mu print >> fh_info, 'New chemical potential found at ', dmu, ' -> Chemical potential becomes ', mu old_mu = mu if (p['runIMP']): if (MPI.rank==Master): print 'Using ',p['base'], 'to compute gf' # If one wants to compute Gloc in DMFT base, ham and olap need to be transformed to DMFT base if (p['base']=='DMFTbase'): trn.transform2(lmt.olap, lmt.hamf, lmt.Utk) ################ if p['orthog']: reorthogonalize(lmt.olap, lmt.hamf) ################ # Impurity levels (Eimp, Eimp_matrix, Eimps) = cmp_Eimp(mu, lmt.olap, lmt.hamf, Sigma_oo_lattice, Sigma_oo_IMP, lmt, p['base'], fh_info, MPI.rank==Master) if p['DCs']=='fixEimp': #BUGS!!! NEEDS SUBSTANTIAL REVISION DUE TO MULTIPLE ATOM POSSIBILITY AND CHANGE OF PARAMETERS for imp in range(len(solver)): iparam = p['iparams'+str(imp)] tEimp = Eimp[lmt.Impc[imp][0] : lmt.Impc[imp][1]] aEimp = sum(tEimp)/len(tEimp) Eimp0 = -(iparam['nf0']-0.5)*iparam['U'] dEdc = Eimp0 - aEimp # This is the discrepancy between the wanted and current Eimp -> needs to correct double counting tEimp += dEdc*ones(len(Eimp)) Eimp[lmt.Impc[imp][0] : lmt.Impc[imp][1]] = tEimp print >> fh_info, 'dEdc=', dEdc, ' a=', dEdc/iparam['U'] if (MPI.rank==Master): print 'Eimp=', Eimp.tolist() print >> fh_info, 'Eimp=', Eimp.tolist() # Below we compute gloc_{ff}, Delta, and partial DOS print >> fh_info, 'Computing gloc and DOS' gcDeDo = Cmp_Dos_Gloc_Delta(fh_info, p['tetrahedronOm'], imag, om, mu, lmt, Sigc, Sigma_oo_lattice, Sigma_oo_IMP, Eimp_matrix, p['base'], p['max_metropolis_steps'], p['cut_ab'], p['gammag']) print >> fh_info, 'Running impurity solver'; fh_info.flush() if (MPI.rank==Master): # Here we print dos.out, gf.out, Delta.out Print_Dos_Gloc_Delta('.', om, gcDeDo) shutil.copy2('dos.out', 'dos.out.'+str(itt)+'.'+str(i)) extn = str(itt)+'.'+str(i) if imag: omc = lom else: omc = om sinf_new = []; Edc_new=[]; nf_IMP=[] Sig_new = [[] for iom in range(len(omc))] for imp,isolver in enumerate(solver): # Here we run impurity solvers for each impurity problem cols = [lmt.Impc[imp][0],lmt.Impc[imp][1]] # columns for this impurity problem Delta = [gcDeDom[1][cols[0]:cols[1]] for gcDeDom in gcDeDo] theipars = 'iparams'+str(imp) #p[theipars]['U'] = [p['U'], "# Coulomb"] #p[theipars]['T'] = [p['T'], "# Temperature"] p[theipars]['Ed'] = [Eimp[cols[0]:cols[1]].tolist(), "# Impurity levels"] if imag: isolver._exe(p[theipars], mpi_prefix, lom, om, Delta, extn, Eimps[imp]) else: isolver._exe(p[theipars], mpi_prefix, om, Delta, p['gbroad'], extn, Eimps[imp]) #p_a = 0 #if p['DCs']=='default' : p_a = p['a'] (tSig_new, tsinf_new, tEdc_new, tnf_IMP) = isolver.HighFrequency(p[theipars], p['DCs'], p['a'], extn) sinf_new.append(tsinf_new) Edc_new.append(tEdc_new) nf_IMP.append(tnf_IMP) for iom in range(len(omc)): for b in range(len(tSig_new)): Sig_new[iom].append(tSig_new[b][iom]) sinf_new = hstack(tuple(sinf_new)) Edc_new = hstack(tuple(Edc_new)) nf_IMP = hstack(tuple(nf_IMP)) WriteSigma('sig.inp', omc, Sig_new, sinf_new, Edc_new) shutil.copy2('sig.inp', 'sig.inp.'+extn) MPI.WORLD.Barrier() if MPI.rank==Master and p['runIMP']: print >> fh_pinfo, '%3d.%3d %12.6f %12.6f %12.6f ' % (itt,i,mu,Eimp[0],Edc_new[0]), print >> fh_pinfo, '%12.6f '*len(nf_IMP) % tuple(nf_IMP) fh_pinfo.flush() if (itt+2 > p['finish']-p['start']): break if (p['base']=='DMFTbase'): lmt.read_ham() # Reads hamiltonian and overlap again if necessary # Current chemical potential is used to compute electron density. Here we compute LMTO coefficients for the density (nxx, mocc_new) = ComputeDensityMatrixCoeff_new(mu, lmt, om, Sigc, Sigma_oo_lattice, fh_info, noccb, p['gamma'], p['max_metropolis_steps'], p['use_tetra'], p['LowerBound']) # mixing of density matrix if (mocc != None): mocc = array(matrix(mocc)*(1-p['admixDM']) + matrix(mocc_new)*p['admixDM'], order='F') else: mocc = mocc_new if (MPI.rank==Master): # correlated occupancy nf = lmt.N_correlated(mocc_new) # New charge density is computed nrho = lmt.recompute_density(nxx) # Reading old scs file scs_f = LmtScsf() (rho_s, rho_cor_s, pot_frz_s, efmu_s, Enu_s) = scs_f.read(lmt.project_name+'.scs') # Computes difference between the old and new density diff = lmt.RhoDiff(rho_s, nrho) # Print occupancy in case we need to restart WriteMocc(os.path.join(p['inpdir'],"ksum/mocc.out"), mocc) # Computing energy on output density scs_f.write(lmt.project_name+'.scse', nrho, rho_cor_s, pot_frz_s, mu, Enu_s) Enes = ldaexe.RunLDAforEnergy(lmt.project_name+'.scse') print >> fh_info, '[EPot, ECoul, EExch]=', Enes print '[EPot, ECoul, EExch]=', Enes # Mixes total density if (p['broyden']): rho_mix = lmt.BroydenMix(rho_s, nrho, p['admix_large'], p['admix0'], p['admix_rho'], p['admix_mag']) else: rho_mix = smix(rho_s, nrho, p['admix_rho']) # Reads scf file to get new frozen potential scf_f = LmtScsf() (rho_f, rho_cor_f, pot_frz_f, efmu_f, Enu_f) = scf_f.read(lmt.project_name+'.scf') # Every restBrnd steps we restart Broyden mixing if (itt % p['restBrdn'] == 0): tmp_files = glob.glob(os.path.join(ldaexe.ksmdir, '*.tmp')) for fu in tmp_files: os.remove(fu) # Mixes the rest rho_cor_mix = smix(rho_cor_s, lmt.rho_cor, p['admix_frozen']) pot_frz_mix = smix(pot_frz_s, pot_frz_f, p['admix_frozen']) Enu_mix=[] for ispin in range(len(Enu_s)): tmp=[] for isort in range(len(Enu_s[ispin])): #print ispin, isort, shape(Enu_s[ispin][isort]), shape(Enu_f[ispin][isort]) #print Enu_s[ispin][isort] #print Enu_f[ispin][isort] elen = len(Enu_s[ispin][isort]) tmp.append((array(Enu_s[ispin][isort])*(1-p['admix_frozen']) + array(Enu_f[ispin][isort][:elen])*p['admix_frozen']).tolist()) Enu_mix.append(tmp) # Writes new density to file scf_f.write(lmt.project_name+'.scsn', rho_mix, rho_cor_f, pot_frz_f, mu, Enu_mix) shutil.move(lmt.project_name+'.scsn', lmt.project_name+'.scs') if MPI.rank==Master: diff = MPI.WORLD[Master].Bcast(diff) else: diff = MPI.WORLD[Master].Bcast() #diff = MPI.WORLD[Master].Bcast(diff) if (diff[0]<1e-5 and diff[1]<1e-5): break if (itt+2 > p['finish']-p['start']): break print >> fh_info, '-'*80 print >> fh_info, '*'*20, 'Iteration', itt+2, '*'*20 print >> fh_info, '-'*80, '\n' if (MPI.rank==Master): print '-'*80 print '*'*20, 'Iteration', itt+2, '*'*20 print '-'*80, '\n' if (MPI.rank==Master): if (p['runLDA']): # LDA computed new Hamiltonian ond Overlap ldaexe.RunLDA() MPI.WORLD.Barrier() # Reds new hamiltonian and computes new DMFT base lmt.read_ham_mix(p['mix_ham']) lmt.cmpTransformation() # Print occupancy in case we need to restart #WriteMocc(os.path.join(p['inpdir'],"ksum/mocc.out"), mocc)