#!/usr/bin/env python import sys, subprocess import re, os, glob, shutil, socket, time from os.path import getsize from scipy import * utime = '/usr/bin/time' # This time function is compatible with w2k timing import utils,indmffile import oca, ctqmc import numpy if numpy.__version__ == '1.0.1': loadtxt = io.read_array # older version of numpy incompatible! def savetxt(filename, data): io.write_array(filename, data, precision=16) import copy class Params(object): """Class to store various control paramers which control the main flow of the program The class is constructed such that it can handle two types of parameters: - starting parameters which are read only at the instantiation - on fly parameters which can be refreshed (reread from the file) at any time by a member function refresh Input: f_sparams -- filename of the parameters read only at the beginning sparams -- a dictionary of pairs {name of the variable: default value} which are looked for in the f_sparams file f_params -- filename of the parameters which can be refreshed several times by refresh member function params -- a dictionary of {variable: default value} pairs which can be updated at runtime Output: None Variables which become members of the class are those listed in params and sparams dictionaries """ def __init__(self, f_params, params, f_sparams=None): self.p ={} self.p.update(params) # default values from the __main__ # remember params and its data file self.f_params = f_params # values from the 'params.dat' file # executes the parameters at startup if f_sparams is not None: execfile(f_sparams) sp = locals() # stores parameters del sp['self'] # deletes automatic locals self.p.update(sp) # new or updates variables def refresh(self, fh_info): # executes the parameters file execfile(self.f_params) sp = locals() # stores locals del sp['self'] # deletes automatics # prints what was changed for c in sp.keys(): if (c not in self.p.keys()) or (sp[c] != self.p[c]): print >>fh_info, '********------ parameter change ------***********' print >>fh_info, c, "=", sp[c] # updates the dictonary self.p.update(sp) def __getitem__(self, x): return self.p[x] def __setitem__(self, i, value): self.p[i]=value def keys(self): return self.p.keys() class Struct: def __init__(self, filename, fh_info=sys.stdout): """ Reads Win2K structure file""" print >> fh_info, '***** Reading Structure file *****' fh = open(filename, 'r') self.title = fh.next() line = fh.next() lattic = line[:4] self.nat = int(line[27:30]) print >> fh_info, 'Number of sorts, nat=', self.nat fh.next() line = fh.next() self.angle=zeros(3,dtype=float) (self.a,self.b,self.c,self.angle[0],self.angle[1],self.angle[2]) = map(float, [line[:10],line[10:20],line[20:30],line[30:40],line[40:50],line[50:60]]) #(self.a,self.b,self.c,self.angle[0],self.angle[1],self.angle[2]) = map(float, line.split()) for i in range(3): if abs(self.angle[i])<1e-5: self.angle[i]=90. print >> fh_info, 'unit vectors (a,b,c) are ', self.a, self.b, self.c print >> fh_info, 'angles betwen unit vectors are', self.angle print >> fh_info, '---- over all types and atoms ----' self.iatnr=[] self.pos=[] self.mult=[] self.isplit=[] self.aname=[] self.jrj=[] self.r0=[] self.rmt=[] self.Znuc=[] self.RotLoc=[] for iat in range(self.nat): print >> fh_info, '---------- atom type', iat, '------------' line = fh.next() self.iatnr.append( int(line[4:8]) ) print >> fh_info, 'iatnr=', self.iatnr[-1] tpos = map(float, [line[12:22], line[25:35],line[38:48]]) self.pos.append(tpos) print >> fh_info, 'Position inside the unit cell, pos=', tpos line = fh.next() self.mult.append( int(line[15:17]) ) print >> fh_info, 'Number of atoms of this type, mult=', self.mult[-1] self.isplit.append( int(line[34:37]) ) print >> fh_info, 'Symmetry of the atom, isplit=', self.isplit[-1] if self.mult[-1]>1: print >> fh_info, ' *** Equivalent atoms ***' for mu in range(self.mult[-1]-1): line = fh.next() self.iatnr.append( int(line[4:8]) ) print >> fh_info, ' iatnr=', self.iatnr[-1], self.pos.append( map(float, [line[12:22], line[25:35],line[38:48]]) ) print >> fh_info, ' pos=', self.pos[-1] line = fh.next() self.aname.append( line[:10] ) self.jrj.append( int(line[15:20]) ) self.r0.append( float(line[25:35]) ) self.rmt.append( float(line[40:50]) ) self.Znuc.append( float(line[55:65]) ) print >> fh_info, '--------- Information about the current type -------------' print >> fh_info, 'aname=', self.aname[-1], 'Znuc=', self.Znuc[-1], 'Rmt=', self.rmt[-1] rt=[] for i in range(3): line = fh.next() rt.append( map(float, [line[20:30], line[30:40],line[40:50]]) ) self.RotLoc.append(rt) print >> fh_info, 'RotLoc=', self.RotLoc[-1] print >> fh_info, '-----------------------------------------------------------' self.aZnuc=[] for iat in range(self.nat): for mu in range(self.mult[iat]): self.aZnuc.append(self.Znuc[iat]) print >> fh_info, 'Z=', self.aZnuc print >> fh_info def dmft1(fday, case, fh_info, extn, MPI, ROOT, m_extn=''): inp_name = 'indmf1'+m_extn name = 'dmft1'+m_extn if not os.path.isfile(case+'.'+inp_name): print ' stop error: the required input file '+case+'.'+inp_name+' for the next step could not be found!' print >> fday, ' stop error: the required input file '+case+'.'+inp_name+' for the next step could not be found!' if not os.path.isfile(name+'.def'): print ' stop error: the required input file '+name+'.def for the next step could not be found!' print >> fday, ' stop error: the required input file '+name+'.def for the next step could not be found!' tim = time.strftime("%H:%M:%S") fl = open(':log', 'a') print >> fl, time.strftime("%a %b %d %H:%M:%S %Z %Y")+'> '+name fl.close() print >> fh_info, 'Running ---- ksum -----' cmd = utime+' '+MPI+ ' ./dmft '+name+'.def >> dmft1_info.out ' print >> fh_info, cmd fh_info.flush() print cmd out, err = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE).communicate() print >> fh_info, out, err fh_info.flush() fe = open('dmft1_log','w'); print >> fe, err; fe.close() print err.splitlines()[0] print >> fday, '>%-10s ( %s ) %s' % (name, tim, err.splitlines()[-2]) fday.flush() for fe in glob.glob('dmft1.error*'): if getsize(fe) !=0: print 'ERROR in dmft1 from file:', fe, open(fe,'r').read() sys.exit(1) if m_extn=='': shutil.copy2(case+'.cdos', case+'.cdos.'+extn) shutil.copy2(case+'.gc1', case+'.gc1.'+extn) shutil.copy2(case+'.dlt1', case+'.dlt1.'+extn) def recompute_mu(fday, case, fh_info, MPI, ROOT): inp_name = 'indmf0' name = 'dmft0' if not os.path.isfile(case+'.'+inp_name): print ' stop error: the required input file '+case+'.'+inp_name+' for the next step could not be found!' print >> fday, ' stop error: the required input file '+case+'.'+inp_name+' for the next step could not be found!' if not os.path.isfile(name+'.def'): print ' stop error: the required input file '+name+'.def for the next step could not be found!' print >> fday, ' stop error: the required input file '+name+'.def for the next step could not be found!' tim = time.strftime("%H:%M:%S") fl = open(':log', 'a') print >> fl, time.strftime("%a %b %d %H:%M:%S %Z %Y")+'> '+name fl.close() print >> fh_info, '--------- Computing eigenvalues ---------' print >> fh_info, '\nRunning ---- ksum ----' cmd = utime+' '+MPI+' ./dmft dmft0.def >> dmft0_info.out' print >> fh_info, cmd fh_info.flush() print cmd sout, serr = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE).communicate() print >> fh_info, sout, serr fh_info.flush() for fe in glob.glob('dmft0.error*'): if getsize(fe) !=0: print 'ERROR in dmft0', fe, ':', open(fe,'r').readlines() #sys.exit(1) print >> fday, '>%-10s ( %s ) %s' % (name, tim, serr.splitlines()[1]) fday.flush() print >> fh_info, '--------- Determine the chemical potential ---------' fh_info.flush() cmd = ROOT+'/x_dmft.py mu' print >> fh_info, '\nRunning ---- x_dmft.py ----' print >> fh_info, cmd fh_info.flush() print cmd sout, serr = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE).communicate() print >> fh_info, sout, serr fh_info.flush() ################################################ ## The charge self-consistency part ################################################ def AStep(fday, case, name, inp_name, WIEN, fh_info): if not os.path.isfile(case+'.'+inp_name): print ' stop error: the required input file '+case+'.'+inp_name+' for the next step could not be found!' print >> fday, ' stop error: the required input file '+case+'.'+inp_name+' for the next step could not be found!' tim = time.strftime("%H:%M:%S") cmd = WIEN+'/x -f '+case+' '+name print cmd print >> fh_info, (name+': '), cmd,'\n' fh_info.flush() stdin, stdout, stderr = os.popen3(cmd) err,out = stderr.read().strip(), stdout.read().strip() print err print >> fday, '>%-10s ( %s ) %s' % (name, tim, out) fday.flush() for fe in glob.glob(name+'.error*'): if getsize(fe) !=0: print 'ERROR in '+fe+' from file:', open(fe,'r').readlines() #sys.exit(1) #if os.path.exists(name+'.error.0') and os.path.getsize(name+'.error.0')>0: # ferr = open(name+'.error.0', 'r') # print ferr.readlines() def lapw0(fday, case, WIEN, fh_info): AStep(fday, case, 'lapw0', 'in0', WIEN, fh_info) def lapw1(fday, case, WIEN, fh_info): AStep(fday, case, 'lapw1 ', 'in1', WIEN, fh_info) def lapwso(fday, case, WIEN, fh_info): AStep(fday, case, 'lapwso', 'inso', WIEN, fh_info) def lcore(fday, case, WIEN, fh_info): AStep(fday, case, 'lcore', 'inc', WIEN, fh_info) def mixer(fday, case, WIEN, fh_info): AStep(fday, case, 'mixer', 'inm', WIEN, fh_info) def dmft2(fday, case, fh_info, WIEN, MPI, m_extn=''): inp_name = 'indmf2'+m_extn name = 'dmft2'+m_extn if not os.path.isfile(case+'.'+inp_name): print ' stop error: the required input file '+case+'.'+inp_name+' for the next step could not be found!' print >> fday, ' stop error: the required input file '+case+'.'+inp_name+' for the next step could not be found!' if not os.path.isfile(name+'.def'): print ' stop error: the required input file '+name+'.def for the next step could not be found!' print >> fday, ' stop error: the required input file '+name+'.def for the next step could not be found!' tim = time.strftime("%H:%M:%S") cmd = utime+' '+MPI + ' ./dmft2 '+name+'.def >> dmft2_info.out ' print >> fh_info, (name+': '), cmd fh_info.flush() fl = open(':log', 'a') print >> fl, time.strftime("%a %b %d %H:%M:%S %Z %Y")+'> '+name fl.close() print cmd stdin, stdout, stderr = os.popen3(cmd) err,out = stderr.read().strip(), stdout.read().strip() print err.splitlines()[0] print >> fday, '>%-10s ( %s ) %s' % (name, tim, err.splitlines()[-2]) fday.flush() for fe in glob.glob(name+'.error*'): if getsize(fe) !=0: print 'ERROR in '+fe+' from file:', open(fe,'r').readlines() #sys.exit(1) dat = open('dmft2_info.out').readlines() dat.reverse() dEF=0 for line in dat: if re.match('Difference in EF', line): dEF = float(line.split()[3]) break return dEF def scf(fday, case, WIEN): dat=[] for i in ['0','1','so','2']: if os.path.exists(case+'.scf'+i): dat += '\n------- '+case+'.scf'+i+' --------\n' dat += open(case+'.scf'+i,'r').readlines() if os.path.exists('Eorb.dat'): dat += '\n---------- Eorb.dat --------------\n' dat += open('Eorb.dat','r').readlines() for i in ['1s','2s','c']: if os.path.exists(case+'.scf'+i): dat += '\n------- '+case+'.scf'+i+' --------\n' dat += open(case+'.scf'+i,'r').readlines() fs = open(case+'.scf', 'a') fs.writelines(dat) fs.close() def scfm(fday, case, WIEN): dat=[] for i in ['m']: dat += '\n------- '+case+'.scf'+i+' --------\n' if os.path.exists(case+'.scf'+i): dat += open(case+'.scf'+i,'r').readlines() fs = open(case+'.scf', 'a') fs.writelines(dat) fs.close() def scf1(fday, case, WIEN): for i in ['clmsum','vsp','vns','vrespsum']: name = case+'.'+i if os.path.exists(name) and os.path.getsize(name)>0: shutil.copy2(name, name+'_old') #save last cycle def Diff(fday, case, dEF): # fs = open(case+'.scf', 'r') dat = fs.readlines() fs.close() DEne=[] DChr=[] for line in dat: if re.match(r':ENE', line): DEne.append(float(line.split()[-1])) #print line.split()[-1] if re.match(r':DIS', line): DChr.append(float(line.split()[-1])) #print line.split()[-1] drho = DChr[-1] if size(DEne)>1: dene = abs(DEne[-1]-DEne[-2]) else: dene = 1.0 if dene==0: dene = 1.0 print >> fday, ':ENERGY convergence: ', dene print >> fday, ':CHARGE convergence: ', drho print >> fday, ':EF convergence: ', dEF return (drho, dene) def SolveImpurity(EF, asolver, iat, extn, p, UpdateAtom, fh_info, fh_pinfo): print 'Running ----- impurity solver -----' print >> fh_info, 'UpdateAtom=', UpdateAtom print >> fh_info, 'Running ----- impurity solver -----' fh_info.flush() iprms = 'iparams'+str(iat) if p['solver'] in ['OCA','NCA','CTQMC']: asolver._exe(p[iprms], extn, UpdateAtom, p['wbroad']) else: print 'ERROR: Impurity solver not defined!' fh_info.flush() print >> fh_info, 'Taking care of high frequency' ntot = asolver.HighFrequency(p[iprms], p['DCs'], p['da'], extn, p['wbroad'], p['kbroad']) # Reading impurity levels for printing fEimp = open('imp.'+str(iat)+'/Eimp.inp', 'r') (Eimp,Olap,Edc) = [array([float(x) for x in line.split('#')[0].split()]) for line in fEimp.readlines()] # throw out comments print >> fh_info, 'Eimp=', Eimp-Edc print >> fh_info, 'Edc=', Edc print >> fh_info, 'ncorr=', ntot icycle, itt = extn.split('.') print >> fh_pinfo, '%3s.%3s %12.6f %12.6f %12.6f %12.6f %12.6f' % (icycle, itt, EF, Eimp[0], Eimp[-1], Edc[0], ntot) fh_pinfo.flush() def SimplifySiginds(siginds): " Takes dictionary of Sigind's and creates list or non-zero columns" def union(data): " Takes a union of array or list" c = [] for d in data: if d not in c: c.append(d) return c cols={} for icix in siginds.keys(): Sigind = siginds[icix] col = sort(union(array(Sigind).flatten())).tolist() if 0 in col: col.remove(0) cols[icix] = col return cols #def ImpurityLatticeConnection( cols, impcols ): # """Should be done better! # We have correlated index stored in Sigind and # impurity index stored in iSigind. # The translation between the two is necessary when # the impurity would like to know in which reference frame # is impurity problem defined. Local rotation should hence # be passed to impurity, and for that we need the translator. # """ # print "cols" # for k,v in cols.iteritems(): # print k,v # print "impcols" # for k,v in impcols.iteritems(): # print k,v # impcols_ind={} # for i,impcol in impcols.iteritems(): # impcols_ind[i] = [] # for j,col in cols.iteritems(): # if len([0 for x in col if x in impcol]) == len(col): # impcols_ind[i].append(j) # return impcols_ind def ImpurityLatticeConnection( cols, icols ): """Should be done better! We have correlated index stored in Sigind and impurity index stored in iSigind. The translation between the two is necessary when the impurity would like to know in which reference frame is impurity problem defined. Local rotation should hence be passed to impurity, and for that we need the translator. """ cols = SimplifySiginds(inl.siginds) icols = SimplifySiginds(iSigind) icols_ind={} for i in icols.keys(): icols_ind[i]=[] for j in cols.keys(): cimp = sort(icols[i]) catm = sort(cols[j]) # equivalent to "catm in cimp" -> if atom j is described by impurity i if reduce(lambda x,y: x and y, [atm in cimp for atm in catm]) : icols_ind[i].append(j) #if sum(sort(icols[i])-sort(cols[j]))==0: icols_ind[i].append(j) return icols_ind def YLM_convert_w2k_to_standard(CF, l, qsplit): """WIEN2K uses a different convention for spherical harmonics than the ED code in the case of d-orbitals. Input: T(i,m) = 'Transformation matrix' printed inside case.indmfl - i runs over real harmonics (z^2, x^2-y^2, yz, xz, xy) - m runs over complex harmonics (-2, -1, 0, 1, 2) Output: T(i,m') The matrix converting between conventions is 1 0 0 0 0 0 -i 0 0 0 R(m,m') = 0 0 1 0 0 0 0 0 i 0 0 0 0 0 1 This function just does: T(i,m') = T(i,m) * R(m,m') """ if l == 2: # only apply conversion if this is NOT a spin-orbit run if qsplit not in [1, -1, 4, 11]: nspin = len(CF)/(2*l+1) # double size of R if two spins R = diag([1, -1j, 1, 1j, 1]*nspin) CFnew = dot(CF, R) else: CFnew = CF return CFnew def combineud(case, ROOT, fh_info): cmd = ROOT+'/combineud '+case print >> fh_info, ('combineud: '), cmd fh_info.flush() stdin, stdout, stderr = os.popen3(cmd) err,out = stderr.read().strip(), stdout.read().strip() if err: print err if out: print out shutil.move(case+'.clmval', case+'.clmval_up') shutil.move(case+'.clmvaldn', case+'.clmval_dn') shutil.move(case+'.clmval_aver', case+'.clmval') def SSplit(ROOT, matsubara, iparams, ntail, fh_info, m_extn=''): cmd = ROOT+'/ssplit.py' if matsubara : cmd += ' -n '+str(iparams['nom'][0])+' -t '+str(ntail) if m_extn: cmd += ' -l '+m_extn print cmd print >> fh_info, cmd sout, serr = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE).communicate() print >> fh_info, sout, serr fh_info.flush() def SGather(ROOT, fh_info, m_extn=''): cmd = ROOT+'/sgather.py' if m_extn: cmd += ' -l '+m_extn print cmd print >> fh_info, cmd sout, serr = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE).communicate() print >> fh_info, sout, serr fh_info.flush() if __name__ == '__main__': # -------------- Parameters which are set through input files sparams.dat and params.dat ----------- # Default values of some parameters params = {'runIMP' : True, # Run impurity solver 'runGF' : True, # Run dmft1 for G and Delta 'UpdateAtom' : 0, # recompute cix-file for the impurity solver 'DCs' : 'fixn', # the double counting scheme, which fixes Edc through n0 #'DCs' : 'default', # the default double counting scheme #'DCs' : 'S0', # experimental #'DCs' : 'fixEimp', # experimental 'max_dmft_iterations': 1, # number of iteration of the dmft-loop only 'max_lda_iterations' : 1, # number of iteration of the LDA-loop only 'finish' : 100, # number of iterations of the charge loop 'da' : 0.0, # experimental 'wbroad' : 0.0, # broadening of sigma on the imaginary axis 'kbroad' : 0.0, # broadening of sigma on the imaginary axis 'recompute_mu' : False, # whether to run dmft0 to recompute mu inside dmft loop. When recomputing charge, this is not necessary 'solver' : 'OCA', # OCA solver #'solver' : 'CTQMC', # CTQMC solver 'com' : 0, # when computing density on imaginary axis, we can subtract Sigma_oo (com=0) or Sigma(0) (com=1) 'ntail' :30, # on imaginary axis, number of points in the tail of the logarithmic mesh 'mix_delta' :1.0, # whether to mix delta, or not #'broyd' : True, # Are we using broyden for charge mixing 'riter' : 100, # How often to restart broyden for charge mixing 'sleeptime' : 2, # If broyden file are present at the submit time, user has some time to kill the job 'cc' : 1e-5, # the charge density precision to stop the LDA+DMFT run 'ec' : 1e-5, # the energy precision to stop the LDA+DMFT run 'so' : False, # spin-orbit coupling 'rCF' : None, # Reduction of the crystal field splitting, if necessary 'recomputeEF' : 1, # Recompute EF in dmft2 step. If recomputeEF=2, it tries to find an insulating gap. 'mixEF' : 1.0, # Chemical potential can be mixed with mixEF<1.0 'saver' : 0.0, # Average over the last few DMFT self-energies is used for charge. Should be between [0,1] 'WL' : 1.0, # If recomputeEF=2, we look for the pole in the interval [-WL,WL] } p = Params('params.dat', params) if len(sys.argv)>1 and sys.argv[1] in ['-h', '--help']: help=""" The script runs LDA+DMFT. It is a wrapper, which calls other scripts and executables. Usuall execution goes through the following steps in a loop: x lapw0 -- computes LDA potential on current LDA+DMFT charge x lapw1 -- solves LDA eigenvalue problem [x lapwso] -- adds spin-orbit [x dmft0] -- computes dmft eigenvalues [x mu] -- computes dmft chemical potential x dmft1 -- computes local green's function and hybridization run impurity -- runs impurity problem x dmft2 -- computes LDA+DMFT valence charge, and the chemical potential x lcore -- computes LDA core charge x mixer -- mixes total charge with the previous charge The most common parameters, which should be given through 'params.dat' file, include: name possible values default help -------------------------------------------------------------------------------------- solver [OCA | CTQMC] OCA # impurity solver max_iterations [int] 1 # number of iteration of the dmft-loop only finish [int] 100 # number of iterations of the charge+dmft loop cc [float] 1e-5, # the charge density precision to stop the LDA+DMFT run ec [float] 1e-5, # the energy precision to stop the LDA+DMFT run broyd [True| False] True # Are we using broyden for charge mixing riter [int] 99 # How often to restart broyden for charge mixing sleeptime [int] 2 # If broyden file are present at the submit time, user # has some time to kill the job so [True| False] False # spin-orbit coupling DCs [fixn| default| S0| fixEimp] fixn # double-counting scheme runIMP [True| False] True # Run impurity solver or skip this step runGF [True| False] True # Run dmft1 for G and Delta, or skip this step ntail [int] 30 # sigma.inp on imaginary axis, number of points in the # tail of the logarithmic mesh mix_delta [float] 1.0 # linear mixing parameter for Delta. UpdateAtom [int] 0 # How often to recompute cix-file inside the # impurity solver wbroad [float] 0.0 # constant broadening of sigma on the imaginary axis kbroad [float] 0.0, # linear broadening of sigma on the imaginary axis recomputeEF [0|1|2] 1 # If we want to recompute the chemical potential in dmft2 step. 2->for Mott insulator recompute_mu [True|False] False # whether to run dmft0 to recompute mu inside dmft # loop. When recomputing charge, this is not necessary com [0| 1] 0 # when computing chemical potential in [dmft0|mu] step # on imaginary axis, we can subtract Sigma_oo (com=0) # or Sigma(0) (com=1) mixEF 0.0-1.0 # mixing of the chemical potential da [float] 0.0 # experimental - connected with DC scheme """ print help sys.exit(0) dmfe = utils.DmftEnvironment() # DMFT paths w2k = utils.W2kEnvironment() # W2k filenames and paths # Processing 'case.indmfl' file inl = indmffile.Indmfl(w2k.case) # case.indmfl file inl.read() # case.indmfl read toclean = glob.glob('*.scf')+glob.glob('Edc.dat')+glob.glob('dmft1.error*')+glob.glob('dmft2.error*')+glob.glob('*.cdos.*')+glob.glob('*.dlt.*')+glob.glob('*.gc1.*')+glob.glob('*.outputdmf1.*')+glob.glob('*.outputdmf2.*') for f in toclean: os.remove(f) # Info files fday = open(w2k.case+'.dayfile', 'w') print >> fday, 'Calculating '+w2k.case+' in '+os.getcwd()+'\n'+'on '+socket.gethostname()+' with PID '+str(os.getpid()) print >> fday, '\n\n' fh_info = open('dmft_info.out','w') fh_pinfo = open('info.iterate', 'a') fh_pinfo.write('%3s.%3s %12s %12s %12s %12s %12s\n' % ('#','#','mu','Eimp','Eimp[-1]','Edc','nf')) fh_pinfo.flush() if len(glob.glob('dmft0_info.out'))!=0 : os.remove('dmft0_info.out') if len(glob.glob('dmft1_info.out'))!=0 : os.remove('dmft1_info.out') if len(glob.glob('dmft2_info.out'))!=0 : os.remove('dmft2_info.out') m_extn = 'dn' if os.path.exists(w2k.case+'.indmfl'+'dn') else '' if m_extn: print >> fh_info, 'INFO: case.indmfldn present => magnetic calculation with two dmft2 steps' inldn = indmffile.Indmfl(w2k.case, 'indmfl'+m_extn) # For parallel execution of the impurity solver # mpi_prefix = dmfe.MPI # Reading parameters from params.dat p.refresh(fh_info) # Reads structure file struct = Struct(w2k.case+'.struct', fh_info) # Check spin-orbit coupling if os.path.isfile(w2k.case+".inso") and os.path.getsize(w2k.case+".inso")>0 : print 'Found '+w2k.case+'.inso file, hence assuming so-coupling exists. Switching -so switch!' p['so'] = True # corelated indexes cixs = inl.siginds.keys() # all columns for cix # Nuclear charges of all correlated blocks Znuc={} for icix in cixs: atm = inl.cps[icix][0][0] # correlated atoms Znuc[icix] = int(struct.aZnuc[atm-1]+0.5) # nuclear charge of each correlated block print >> fh_info, 'Znucs=', Znuc # Independent impurity problems iSigind = utils.ParsIndmfi(w2k.case) print >> fh_info, 'iSiginds=', iSigind cols = SimplifySiginds(inl.siginds) icols = SimplifySiginds(iSigind) icols_ind = ImpurityLatticeConnection( cols, icols ) print >> fh_info, 'Impurity-lattice connection: icols_ind=', icols_ind # Solvers initialized if p['solver'] in ['OCA', 'NCA']: if inl.matsubara: print 'ERROR: Can not run OCA/NCA on imaginary axis! Change case.indmfl:matsubara or params.dat:solver' solver=[] for i in iSigind.keys(): ic = icols_ind[i][0] iatom, l, qsplit = inl.cps[ic][0] # Doesn't take into account cluster problems correctly cftrans = YLM_convert_w2k_to_standard(inl.cftrans[ic], l, qsplit) solver.append( oca.IMP_OCA(Znuc[ic], 'imp.'+str(i), p['iparams'+str(i)], iSigind[i], cftrans) ) elif p['solver'] in ['CTQMC']: if not inl.matsubara: print 'ERROR: Can not run CTQMC on real axis! Change case.indmfl:matsubara or params.dat:solver' solver=[] for i in iSigind.keys(): ic = icols_ind[i][0] iatom, l, qsplit = inl.cps[ic][0] # Doesn't take into account cluster problems correctly cftrans = YLM_convert_w2k_to_standard(inl.cftrans[ic], l, qsplit) solver.append( ctqmc.IMP_CTQMC(Znuc[ic], 'imp.'+str(i), p['iparams'+str(i)], iSigind[i], cftrans) ) else: print 'ERROR: No solver specified! Bolling out.' shutil.copy2(dmfe.ROOT+'/dmft', '.') # For parallel execution, executable has to be on the current directory shutil.copy2(dmfe.ROOT+'/dmft2', '.') # For parallel execution, executable has to be on the current directory print >> fh_info, '--------- Preparing GF calculation ---------' cmd = dmfe.ROOT+'/x_dmft.py -d dmft0 >> dmft0_info.out 2>&1' print >> fh_info, cmd stdin, stdout, stderr = os.popen3(cmd) print >> fh_info, stdout.read(), stderr.read() print >> fh_info, '--------- Preparing GF calculation ---------' cmd = dmfe.ROOT+'/x_dmft.py -d dmft1 >> dmft1_info.out 2>&1' print >> fh_info, cmd stdin, stdout, stderr = os.popen3(cmd) print >> fh_info, stdout.read(), stderr.read() if m_extn: print >> fh_info, '--------- Preparing GC calculation '+m_extn+' ---------' cmd = dmfe.ROOT+'/x_dmft.py -d dmft1 -l dn >> dmft1_info.out 2>&1' print >> fh_info, cmd stdin, stdout, stderr = os.popen3(cmd) print >> fh_info, stdout.read(), stderr.read() print >> fh_info, '--------- Preparing CHR calculation ---------' cmd = dmfe.ROOT+'/x_dmft.py -d -m '+str(p['recomputeEF'])+' -x '+str(p['mixEF'])+' -w '+str(p['WL'])+' dmft2 >> ksum_info.out 2>&1' print >> fh_info, cmd stdin, stdout, stderr = os.popen3(cmd) print >> fh_info, stdout.read(), stderr.read() if m_extn: print >> fh_info, '--------- Preparing CHR calculation '+m_extn+'---------' cmd = dmfe.ROOT+'/x_dmft.py -d -m 0 dmft2 -l dn >> ksum_infodn.out 2>&1' print >> fh_info, cmd stdin, stdout, stderr = os.popen3(cmd) print >> fh_info, stdout.read(), stderr.read() shutil.copy2('sig.inp', 'sig.inp.0') if p['finish']>1: # Charge self-consistent if os.path.isfile(w2k.case+'.clmsum_old'): shutil.copy2(w2k.case+'.clmsum_old', w2k.case+'.clmsum') if not os.path.isfile(w2k.case+'.clmsum'): print 'no '+w2k.case+'.clmsum(_old) file found, which is necessary for lapw0 !' print >> fday, 'no '+w2k.case+'.clmsum(_old) file found, which is necessary for lapw0 !' sys.exit(1) if os.path.isfile(w2k.case+'.broyd1'): print w2k.case+'.broyd* files present! You did not save_lapw a previous clculation.' print 'You have '+str(p['sleeptime'])+' seconds to kill this job ( ^C or kill '+str(os.getpid())+' )' print 'or the script will rm *.broyd* and continue (use -NI to avoid automatic rm)' time.sleep(p['sleeptime']) cmd = 'rm -f *.broyd*' os.popen3(cmd) stdin, stdout, stderr = os.popen3(cmd) err,out = stderr.read().strip(), stdout.read().strip() print >> fday, w2k.case+'.broyd* files removed!' next='lapw0' print >> fday, '\n start'+(' '*8)+ time.asctime()+' with '+next+' ('+str(1)+'/'+str(p['riter'])+' to go)' # Starting with LAPW print >> fday, '\n cycle %s \t%s %s/%s to go\n' % (0,time.asctime(),p['finish'],(p['finish'])%p['riter']) lapw0(fday,w2k.case,w2k.WIENROOT, fh_info); fday.flush() lapw1(fday,w2k.case,w2k.WIENROOT, fh_info); fday.flush() if p['so']: lapwso(fday,w2k.case,w2k.WIENROOT, fh_info) fday.flush() ##################### # Major charge loop # ##################### icycle = nimp = 0 while icycle < p['finish']: ##################### # DMFT loop # ##################### sigmas=[] itt = 0 while itt < p['max_dmft_iterations']: extn = str(icycle+1)+'.'+str(itt+1) # Extension of output files p.refresh(fh_info) # Recomputing the chemical potential if p['recompute_mu']: SSplit(dmfe.ROOT, inl.matsubara, p['iparams0'], p['ntail'], fh_info) recompute_mu(fday, w2k.case, fh_info, dmfe.MPI, dmfe.ROOT) # Computing Green's function if p['runGF']: SSplit(dmfe.ROOT, inl.matsubara, p['iparams0'], p['ntail'], fh_info) dmft1(fday, w2k.case, fh_info, extn, dmfe.MPI, dmfe.ROOT) if m_extn: SSplit(dmfe.ROOT, inl.matsubara, p['iparams0'], p['ntail'], fh_info, m_extn) dmft1(fday, w2k.case, fh_info, extn, dmfe.MPI, dmfe.ROOT, m_extn) cmd = dmfe.ROOT+'/sjoin.py -m '+str(p['mix_delta']) if m_extn: cmd += ' -l '+m_extn if p['rCF'] is not None: cmd += ' -c "'+ p['rCF']+'"' print cmd print >> fh_info, cmd sout, serr = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE).communicate() print >> fh_info, sout, serr fh_info.flush() if os.path.isfile('EF.dat'): EF = float(open('EF.dat').read()) else: open(w2k.case+'.indmf1') EF =float(open(w2k.case+'.indmf1').readlines()[1].split()[1]) # ------------ Impurity solver ------------------- if p['runIMP']: nimp += 1 UpdateAtom = (p['UpdateAtom']!=0 and mod(nimp,p['UpdateAtom'])==0) for iat in iSigind.keys(): # Over all inequivalent impurity problems SolveImpurity(EF, solver[iat], iat, extn, p, UpdateAtom, fh_info, fh_pinfo) SGather(dmfe.ROOT, fh_info,m_extn) Eorb = 0.0 dat=[] for iat in iSigind.keys(): # Over all inequivalent impurity problems if os.path.exists(solver[iat].dir+'/Eorb.dat'): dat += ['--- impurity '+str(iat)+' with degeneracy '+str(len(icols_ind[i]))+' ---\n'] datw = open(solver[iat].dir+'/Eorb.dat','r').readlines() for line in datw: if line[:5]==':EORB': Eorb += float(line[5:]) * len(icols_ind[i]) # Takes into account the degeneracy of the impurity else: dat += [line] dat += [':EORB '+ str(Eorb)] open('Eorb.dat','w').writelines(dat) # Saves on the current directory the combination of all impurity problems shutil.copy2('sig.inp','sig.inp.'+extn) sigmas.append('sig.inp.'+extn) # save the filename print >> fh_info, '------- DMFT part of the step # ', itt, 'done ---------\n' fh_info.flush() itt += 1 # END DMFT INNER LOOP if len(sigmas)>0: sigmas.reverse() sigmas = sigmas[0:max(1,int(len(sigmas)*p['saver']))] print >> fh_info, 'Averaging over the following self-energies', sigmas cmd = 'saverage.py -o sig.inp ' + (' '.join(sigmas)) print >> fh_info, cmd sout, serr = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE).communicate() print >> fh_info, sout, serr fh_info.flush() for ita in range(p['max_lda_iterations']): SSplit(dmfe.ROOT, inl.matsubara, p['iparams0'], p['ntail'], fh_info) dEF = dmft2(fday,w2k.case,fh_info,w2k.WIENROOT,dmfe.MPI); fday.flush() if m_extn: SSplit(dmfe.ROOT, inl.matsubara, p['iparams0'], p['ntail'], fh_info, m_extn) dmft2(fday,w2k.case,fh_info,w2k.WIENROOT,dmfe.MPI,m_extn) fday.flush() combineud(w2k.case, dmfe.ROOT, fh_info) lcore(fday,w2k.case,w2k.WIENROOT,fh_info); fday.flush() scf (fday,w2k.case,w2k.WIENROOT); fday.flush() scf1 (fday,w2k.case,w2k.WIENROOT); fday.flush() mixer(fday,w2k.case,w2k.WIENROOT,fh_info); fday.flush() scfm (fday,w2k.case,w2k.WIENROOT); fday.flush() drho, dene = Diff(fday, w2k.case, dEF) print >> fh_info, '------- LDA(charge) step', ita, icycle, 'done ---------' istep1 = icycle+1 istep2 = ita+1 if (drho> fday, w2k.case+'.broyd* files removed!' print >> fday, '\n:ITT%-2d cycle %s \t%s %s/%s to go\n' % (istep2, istep1,time.asctime(),p['finish']-istep1,(p['finish']-istep1)%p['riter']) lapw0(fday,w2k.case,w2k.WIENROOT, fh_info); fday.flush() lapw1(fday,w2k.case,w2k.WIENROOT, fh_info); fday.flush() if p['so']: lapwso(fday,w2k.case,w2k.WIENROOT, fh_info) fday.flush() if drho