from triqs_dft_tools.sumk_dft import * from triqs_dft_tools.converters.wien2k import * from triqs.applications.impurity_solvers.hubbard_I.hubbard_solver import Solver import os dft_filename = os.getcwd().rpartition('/')[2] beta = 40 U_int = 6.00 J_hund = 0.70 Loops = 5 # Number of DMFT sc-loops mixing = 0.7 # Mixing factor DC_type = 0 # 0...FLL, 1...Held, 2... AMF, 3...Lichtenstein chemical_potential_init=0.0 # initial chemical potential # Convert DMFT input: Converter = Wien2kConverter(filename=dft_filename) Converter.convert_dft_input() mpi.barrier() #check if there are previous runs: previous_runs = 0 previous_present = False if mpi.is_master_node(): with HDFArchive(dft_filename+'.h5','a') as f: if 'dmft_output' in f: ar = f['dmft_output'] if 'iterations' in ar: previous_present = True previous_runs = ar['iterations'] else: f.create_group('dmft_output') previous_runs = mpi.bcast(previous_runs) previous_present = mpi.bcast(previous_present) # Init the SumK class SK=SumkDFT(hdf_file=dft_filename+'.h5',use_dft_blocks=False) Norb = SK.corr_shells[0]['dim'] l = SK.corr_shells[0]['l'] # Init the Hubbard-I solver: S = Solver(beta = beta, l = l) chemical_potential=chemical_potential_init # load previous data: old self-energy, chemical potential, DC correction if previous_present: if mpi.is_master_node(): with HDFArchive(dft_filename+'.h5','r') as ar: S.Sigma << ar['dmft_output']['Sigma'] SK.chemical_potential,SK.dc_imp,SK.dc_energ = SK.load(['chemical_potential','dc_imp','dc_energ']) S.Sigma << mpi.bcast(S.Sigma) SK.chemical_potential = mpi.bcast(SK.chemical_potential) SK.dc_imp = mpi.bcast(SK.dc_imp) SK.dc_energ = mpi.bcast(SK.dc_energ) # DMFT loop: for iteration_number in range(1,Loops+1): itn = iteration_number + previous_runs # put Sigma into the SumK class: SK.set_Sigma([ S.Sigma ]) # Compute the SumK, possibly fixing mu by dichotomy chemical_potential = SK.calc_mu( precision = 0.000001 ) # Density: S.G <<= SK.extract_G_loc()[0] mpi.report("Total charge of Gloc : %.6f"%S.G.total_density()) # calculated DC at the first run to have reasonable initial non-interacting atomic level positions if ((iteration_number==1)and(previous_present==False)): dc_value_init=U_int/2.0 dm=S.G.density() SK.calc_dc( dm, U_interact = U_int, J_hund = J_hund, orb = 0, use_dc_formula = DC_type, use_dc_value=dc_value_init) # calculate non-interacting atomic level positions: eal = SK.eff_atomic_levels()[0] S.set_atomic_levels( eal = eal ) # solve it: S.solve(U_int = U_int, J_hund = J_hund, verbosity = 1) # Now mix Sigma and G with factor Mix, if wanted: if (iteration_number>1 or previous_present): if (mpi.is_master_node() and (mixing<1.0)): with HDFArchive(dft_filename+'.h5','r') as ar: mpi.report("Mixing Sigma and G with factor %s"%mixing) S.Sigma << mixing * S.Sigma + (1.0-mixing) * ar['dmft_output']['Sigma'] S.G << mixing * S.G + (1.0-mixing) * ar['dmft_output']['G'] S.G << mpi.bcast(S.G) S.Sigma << mpi.bcast(S.Sigma) # after the Solver has finished, set new double counting: dm = S.G.density() SK.calc_dc( dm, U_interact = U_int, J_hund = J_hund, orb = 0, use_dc_formula = DC_type ) # correlation energy calculations: SK.correnerg = 0.5 * (S.G * S.Sigma).total_density() mpi.report("Corr. energy = %s"%SK.correnerg) # store the impurity self-energy, GF as well as correlation energy in h5 if mpi.is_master_node(): with HDFArchive(dft_filename+'.h5','a') as ar: ar['dmft_output']['iterations'] = iteration_number + previous_runs ar['dmft_output']['G'] = S.G ar['dmft_output']['Sigma'] = S.Sigma #Save essential SumkDFT data: SK.save(['chemical_potential','dc_imp','dc_energ','correnerg']) if (mpi.is_master_node()): print('DC after solver: ',SK.dc_imp[0]) # print out occupancy matrix of Ce 4f mpi.report("Orbital densities of impurity Green function:") for s in dm: mpi.report("Block %s: "%s) for ii in range(len(dm[s])): str = '' for jj in range(len(dm[s])): if (dm[s][ii,jj].real>0): str += " %.4f"%(dm[s][ii,jj].real) else: str += " %.4f"%(dm[s][ii,jj].real) mpi.report(str) mpi.report("Total charge of impurity problem : %.6f"%S.G.total_density()) # find exact chemical potential SK.chemical_potential = SK.calc_mu( precision = 0.000001 ) # calculate and save occupancy matrix in the Bloch basis for Wien2k charge denity recalculation dN,d = SK.calc_density_correction(filename = dft_filename+'.qdmft') mpi.report("Trace of Density Matrix: %s"%d) # store correlation energy contribution to be read by Wien2ki and then included to DFT+DMFT total energy if (mpi.is_master_node()): SK.correnerg -= SK.dc_energ[0] f=open(dft_filename+'.qdmft','a') f.write("%.16f\n"%SK.correnerg) f.close()