.. _advanced: A more advanced example ======================= Normally, one wants to adjust some more parameters in order to make the calculation more efficient. Here, we will see a more advanced example, which is also suited for parallel execution. First, we load the necessary modules:: from pytriqs.applications.dft.sumk_dft import * from pytriqs.applications.dft.converters.wien2k_converter import * from pytriqs.gf.local import * from pytriqs.archive import HDFArchive from pytriqs.operators.hamiltonians import * from pytriqs.applications.impurity_solvers.cthyb import * Then we define some parameters:: dft_filename='srvo3' U = 2.7 J = 0.65 beta = 40 loops = 10 # Number of DMFT sc-loops sigma_mix = 0.8 # Mixing factor of Sigma after solution of the AIM delta_mix = 1.0 # Mixing factor of Delta as input for the AIM dc_type = 1 # DC type: 0 FLL, 1 Held, 2 AMF use_blocks = True # use bloc structure from DFT input prec_mu = 0.0001 # Solver parameters p = {} p["length_cycle"] = 200 p["n_warmup_cycles"] = 2000 p["n_cycles"] = 20000 Most of these parameters are self-explanatory. The first, `dft_filename`, gives the filename of the input files. The next step, as described in the previous section, is to convert the input files:: Converter = Wien2kConverter(filename=dft_filename, repacking=True) Converter.convert_dft_input() mpi.barrier() The command ``mpi.barrier()`` ensures that all nodes wait until the conversion of the input is finished on the master node. After the conversion, we can check in the hdf5 archive, if previous runs are present, or if we have to start from scratch:: previous_runs = 0 previous_present = False if mpi.is_master_node(): f = HDFArchive(dft_filename+'.h5','a') 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') del f previous_runs = mpi.bcast(previous_runs) previous_present = mpi.bcast(previous_present) Now we can use all this information to initialise the :class:`SumkDFT` class:: SK = SumkDFT(hdf_file=dft_filename+'.h5',use_dft_blocks=use_blocks) The next step is to initialise the Solver:: n_orb = SK.corr_shells[0]['dim'] l = SK.corr_shells[0]['l'] spin_names = ["up","down"] orb_names = [i for i in range(n_orb)] # Use GF structure determined by DFT blocks gf_struct = SK.gf_struct_solver[0] # Construct U matrix for density-density calculations Umat, Upmat = U_matrix_kanamori(n_orb=n_orb, U_int=U, J_hund=J) # Construct Hamiltonian and solver h_loc = h_loc_density(spin_names, orb_names, map_operator_structure=SK.sumk_to_solver[0], U=Umat, Uprime=Upmat, H_dump="H.txt") S = Solver(beta=beta, gf_struct=gf_struct) If there are previous runs stored in the hdf5 archive, we can now load the self energy of the last iteration:: if previous_present: if mpi.is_master_node(): S.Sigma_iw << HDFArchive(dft_filename+'.h5','a')['dmft_output']['Sigma_iw'] chemical_potential,dc_imp,dc_energ = SK.load(['chemical_potential','dc_imp','dc_energ']) S.Sigma_iw << mpi.bcast(S.Sigma_iw) SK.set_mu(chemical_potential) SK.set_dc(dc_imp,dc_energ) The self-energy is broadcast from the master node to the slave nodes. Also, the last saved chemical potential and double counting values are read in and set. Now we can go to the definition of the self-consistency step. It consists again of the basic steps discussed in the previous section, with some additional refinement:: for iteration_number in range(1,loops+1): if mpi.is_master_node(): print "Iteration = ", iteration_number SK.symm_deg_gf(S.Sigma_iw,orb=0) # symmetrise Sigma SK.put_Sigma(Sigma_imp = [ S.Sigma_iw ]) # put Sigma into the SumK class chemical_potential = SK.calc_mu( precision = prec_mu ) # find the chemical potential for given density S.G_iw << SK.extract_G_loc()[0] # calc the local Green function mpi.report("Total charge of Gloc : %.6f"%S.G_iw.total_density()) # Init the DC term and the real part of Sigma, if no previous runs found: if (iteration_number==1 and previous_present==False): dm = S.G_iw.density() SK.calc_dc(dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type) S.Sigma_iw << SK.dc_imp[0]['up'][0,0] # Calculate new G0_iw to input into the solver: if mpi.is_master_node(): # We can do a mixing of Delta in order to stabilize the DMFT iterations: S.G0_iw << S.Sigma_iw + inverse(S.G_iw) ar = HDFArchive(dft_filename+'.h5','a')['dmft_output'] if (iteration_number>1 or previous_present): mpi.report("Mixing input Delta with factor %s"%delta_mix) Delta = (delta_mix * delta(S.G0_iw)) + (1.0-delta_mix) * ar['Delta_iw'] S.G0_iw << S.G0_iw + delta(S.G0_iw) - Delta ar['Delta_iw'] = delta(S.G0_iw) S.G0_iw << inverse(S.G0_iw) del ar S.G0_iw << mpi.bcast(S.G0_iw) # Solve the impurity problem: S.solve(h_loc=h_loc, **p) # Solved. Now do post-processing: mpi.report("Total charge of impurity problem : %.6f"%S.G_iw.total_density()) # Now mix Sigma and G with factor sigma_mix, if wanted: if (iteration_number>1 or previous_present): if mpi.is_master_node(): ar = HDFArchive(dft_filename+'.h5','a')['dmft_output'] mpi.report("Mixing Sigma and G with factor %s"%sigma_mix) S.Sigma_iw << sigma_mix * S.Sigma_iw + (1.0-sigma_mix) * ar['Sigma_iw'] S.G_iw << sigma_mix * S.G_iw + (1.0-sigma_mix) * ar['G_iw'] del ar S.G_iw << mpi.bcast(S.G_iw) S.Sigma_iw << mpi.bcast(S.Sigma_iw) # Write the final Sigma and G to the hdf5 archive: if mpi.is_master_node(): ar = HDFArchive(dft_filename+'.h5','a')['dmft_output'] if previous_runs: iteration_number += previous_runs ar['iterations'] = iteration_number ar['G_0'] = S.G0_iw ar['G_tau'] = S.G_tau ar['G_iw'] = S.G_iw ar['Sigma_iw'] = S.Sigma_iw del ar # Set the new double counting: dm = S.G_iw.density() # compute the density matrix of the impurity problem SK.calc_dc(dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type) # Save stuff into the dft_output group of hdf5 archive in case of rerun: SK.save(['chemical_potential','dc_imp','dc_energ']) This is all we need for the DFT+DMFT calculation. At the end, all results are stored in the hdf5 output file.