From 6fff56fe8d6550fe307221ef42adcb0b763a9411 Mon Sep 17 00:00:00 2001 From: aichhorn Date: Thu, 13 Aug 2015 16:50:48 +0200 Subject: [PATCH] new version of single-shot dft+dmft --- doc/guide/dftdmft_singleshot.rst | 131 ++++++++++----- doc/guide/images_scripts/dft_dmft_cthyb.py | 149 ++++++++++++++++++ .../images_scripts/dft_dmft_cthyb_slater.py | 149 ++++++++++++++++++ 3 files changed, 390 insertions(+), 39 deletions(-) create mode 100644 doc/guide/images_scripts/dft_dmft_cthyb.py create mode 100644 doc/guide/images_scripts/dft_dmft_cthyb_slater.py diff --git a/doc/guide/dftdmft_singleshot.rst b/doc/guide/dftdmft_singleshot.rst index 51330c3f..02d68466 100644 --- a/doc/guide/dftdmft_singleshot.rst +++ b/doc/guide/dftdmft_singleshot.rst @@ -5,11 +5,10 @@ Single-shot DFT+DMFT ==================== -.. warning:: - TO BE UPDATED! After having set up the hdf5 archive, we can now do our DFT+DMFT calculation. It consists of -initialisation steps, and the actual DMFT self consistency loop. +initialization steps, and the actual DMFT self-consistency loop, as is +discussed below. Initialisation of the calculation --------------------------------- @@ -25,15 +24,27 @@ to get the local quantities used in DMFT. It is initialized by:: Setting up the impurity solver ------------------------------ -The next step is to setup the impurity solver. -For more details here, see the `CTHYB `_ documentation. +The next step is to setup an impurity solver. There are different +solvers available within the TRIQS framework. Below, we will discuss +the example of the hybridisation +expansion :ref:`CTHYB solver `. Later on, we will +see also the example of the Hubbard-I solver. They all have in common, +that they are called by a uniform command:: + + S.solve(params) + +where `params` are the solver parameters and depend on the actual +solver that is used. Before going into the details of the solver, let +us discuss in the next section how to perform the DMFT loop using +the methods of :program:`dft_tools`, assuming that we have set up a +working solver instance. Doing the DMFT loop ------------------- -Having initialised the SumK class and the Solver, we can proceed with the DMFT -loop itself. As explained in the tutorial, we have to set up the loop over DMFT +Having initialized the SumK class and the Solver, we can proceed with the DMFT +loop itself. We have to set up the loop over DMFT iterations and the self-consistency condition:: n_loops = 5 @@ -47,13 +58,18 @@ iterations and the self-consistency condition:: S.solve(h_int=h_int, **p) # now solve the impurity problem dm = S.G_iw.density() # Density matrix of the impurity problem - SK.calc_dc(dm, U_interact=U, J_hund=J, orb=0, use_dc_formula=dc_type) # Set the double counting term + SK.calc_dc(dm, U_interact=U, J_hund=J, orb=0, use_dc_formula=1) # Set the double counting term SK.save(['chemical_potential','dc_imp','dc_energ']) # Save data in the hdf5 archive These basic steps are enough to set up the basic DMFT Loop. For a detailed -description of the :class:`SumkDFT` routines, see the reference manual. After -the self-consistency steps, the solution of the Anderson impurity problem is -calculation by CTQMC. Different to model calculations, we have to do a few +description of the :class:`SumkDFT` routines, see the reference +manual. + +After +the self-consistency steps (extracting a new :math:`G^0(i\omega)`), +the Anderson impurity problem is solved. + +Different to model calculations, we have to do a few more steps after this, because of the double-counting correction. We first calculate the density of the impurity problem. Then, the routine `calc_dc` takes as parameters this density matrix, the Coulomb interaction, Hund's rule @@ -80,11 +96,21 @@ For full charge-self consistent calculations, there are some more things to consider, which we will see later on. -A more advanced example ------------------------ +A full DFT+DMFT calculation +--------------------------- + +We will discuss now how to set up a full working calculation, +including setting up the CTHYB solver, and specifying 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. For the convenience of the user, we provide also two +working python scripts in this documentation. One for a calculation +using Kanamori definitions (:download:`dft_dmft_cthyb.py +`) and one with a +rotational-invariant Slater interaction hamiltonian (:download:`dft_dmft_cthyb_slater.py +`). The user has to adapt these +scripts to his own needs. -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 * @@ -94,16 +120,17 @@ First, we load the necessary modules:: from pytriqs.operators.util import * from pytriqs.applications.impurity_solvers.cthyb import * +The last two lines load the modules for the construction of the CTHYB +solver. Then we define some parameters:: - dft_filename='srvo3' - U = 2.7 + dft_filename='SrVO3' + U = 4.0 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 @@ -114,7 +141,11 @@ Then we define some parameters:: 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. +Most of these parameters are self-explanatory. The first, +`dft_filename`, gives the filename of the input files. For more +details on the solver parameters, we refer the user to +the :ref:`CTHYB solver ` documentation. + The next step, as described in the previous section, is to convert the input files:: Converter = Wien2kConverter(filename=dft_filename, repacking=True) @@ -140,24 +171,57 @@ from scratch:: previous_runs = mpi.bcast(previous_runs) previous_present = mpi.bcast(previous_present) + +You can see in this code snipet, that all results of this calculation +will be stored in a separate subgroup in the hdf file, called +`dmft_output`. Removing this subgroup allows you to reset your +calculation to the starting point easily. + 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 :class:`Solver` class:: +The next step is to initialise the :class:`Solver` class. It consist +of two steps + +#. Calculating the multi-band interaction matrix, and setting up the + interaction hamiltonian +#. Setting up the solver class + +The first step is done using methods of +the :ref:`TRIQS ` library:: 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 + # Use GF structure determined by DFT blocks: gf_struct = SK.gf_struct_solver[0] - # Construct U matrix for density-density calculations + # 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_int = h_int_density(spin_names, orb_names, map_operator_structure=SK.sumk_to_solver[0], U=Umat, Uprime=Upmat, H_dump="H.txt") + +We assumed here that we want to use an interaction matrix with +Kanamori definitions of :math:`U` and :math:`J`. For +other choices (Slater interaction matrix for instance), and other +parameters, we refer to the reference manual +of the :ref:`TRIQS ` library. + +Next, we construct the hamiltonian and the solver:: + + h_int = h_int_density(spin_names, orb_names, map_operator_structure=SK.sumk_to_solver[0], U=Umat, Uprime=Upmat) S = Solver(beta=beta, gf_struct=gf_struct) +As you see, we take only density-density interactions into +account. Other choices for the hamiltonian are + +* h_int_kanamori +* h_int_slater + +These two include full rotational invariant interactions. Again, +options can be found in the :ref:`TRIQS ` library +reference manual. + + If there are previous runs stored in the hdf5 archive, we can now load the self energy of the last iteration:: @@ -174,7 +238,7 @@ 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:: +refinements:: for iteration_number in range(1,loops+1): if mpi.is_master_node(): print "Iteration = ", iteration_number @@ -192,24 +256,13 @@ refinement:: 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) + S.G0_iw << S.Sigma_iw + inverse(S.G_iw) + S.G0_iw << inverse(S.G0_iw) # Solve the impurity problem: S.solve(h_int=h_int, **p) - # Solved. Now do post-processing: + # Solved. Now do post-solution stuff: mpi.report("Total charge of impurity problem : %.6f"%S.G_iw.total_density()) # Now mix Sigma and G with factor sigma_mix, if wanted: diff --git a/doc/guide/images_scripts/dft_dmft_cthyb.py b/doc/guide/images_scripts/dft_dmft_cthyb.py new file mode 100644 index 00000000..808a31c9 --- /dev/null +++ b/doc/guide/images_scripts/dft_dmft_cthyb.py @@ -0,0 +1,149 @@ +import pytriqs.utility.mpi as mpi +from pytriqs.operators.util import * +from pytriqs.archive import HDFArchive +from pytriqs.applications.impurity_solvers.cthyb import * +from pytriqs.gf.local import * +from pytriqs.applications.dft.sumk_dft import * +from pytriqs.applications.dft.converters.wien2k_converter import * + +dft_filename='SrVO3' +U = U.0 +J = 0.65 +beta = 40 +loops = 10 # Number of DMFT sc-loops +sigma_mix = 1.0 # 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 +h_field = 0.0 + +# Solver parameters +p = {} +p["max_time"] = -1 +p["length_cycle"] = 50 +p["n_warmup_cycles"] = 50 +p["n_cycles"] = 5000 + +Converter = Wien2kConverter(filename=dft_filename, repacking=True) +Converter.convert_dft_input() +mpi.barrier() + +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) + +SK=SumkDFT(hdf_file=dft_filename+'.h5',use_dft_blocks=use_blocks,h_field=h_field) + +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 density-density Hamiltonian and solver +h_int = h_int_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 previous_present: + chemical_potential = 0 + dc_imp = 0 + dc_energ = 0 + 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) + chemical_potential = mpi.bcast(chemical_potential) + dc_imp = mpi.bcast(dc_imp) + dc_energ = mpi.bcast(dc_energ) + SK.set_mu(chemical_potential) + SK.set_dc(dc_imp,dc_energ) + +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_int=h_int, **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_tau'] = S.G_tau + ar['G_iw'] = S.G_iw + ar['Sigma_iw'] = S.Sigma_iw + ar['G0-%s'%(iteration_number)] = S.G0_iw + ar['G-%s'%(iteration_number)] = S.G_iw + ar['Sigma-%s'%(iteration_number)] = 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']) + +if mpi.is_master_node(): + ar = HDFArchive("dftdmft.h5",'w') + ar["G_tau"] = S.G_tau + ar["G_iw"] = S.G_iw + ar["Sigma_iw"] = S.Sigma_iw diff --git a/doc/guide/images_scripts/dft_dmft_cthyb_slater.py b/doc/guide/images_scripts/dft_dmft_cthyb_slater.py new file mode 100644 index 00000000..8ed8ccc4 --- /dev/null +++ b/doc/guide/images_scripts/dft_dmft_cthyb_slater.py @@ -0,0 +1,149 @@ +import pytriqs.utility.mpi as mpi +from pytriqs.operators.util import * +from pytriqs.archive import HDFArchive +from pytriqs.applications.impurity_solvers.cthyb import * +from pytriqs.gf.local import * +from pytriqs.applications.dft.sumk_dft import * +from pytriqs.applications.dft.converters.wien2k_converter import * + +dft_filename='Gd_fcc' +U = 9.6 +J = 0.8 +beta = 40 +loops = 10 # Number of DMFT sc-loops +sigma_mix = 1.0 # Mixing factor of Sigma after solution of the AIM +delta_mix = 1.0 # Mixing factor of Delta as input for the AIM +dc_type = 0 # DC type: 0 FLL, 1 Held, 2 AMF +use_blocks = True # use bloc structure from DFT input +prec_mu = 0.0001 +h_field = 0.0 + +# Solver parameters +p = {} +p["max_time"] = -1 +p["length_cycle"] = 50 +p["n_warmup_cycles"] = 50 +p["n_cycles"] = 5000 + +Converter = Wien2kConverter(filename=dft_filename, repacking=True) +Converter.convert_dft_input() +mpi.barrier() + +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) + +SK=SumkDFT(hdf_file=dft_filename+'.h5',use_dft_blocks=use_blocks,h_field=h_field) + +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 Slater U matrix +Umat = U_matrix(n_orb=n_orb, U_int=U, J_hund=J, basis='cubic',) + +# Construct Hamiltonian and solver +h_int = h_int_slater(spin_names, orb_names, map_operator_structure=SK.sumk_to_solver[0], U_matrix=Umat) +S = Solver(beta=beta, gf_struct=gf_struct) + +if previous_present: + chemical_potential = 0 + dc_imp = 0 + dc_energ = 0 + 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) + chemical_potential = mpi.bcast(chemical_potential) + dc_imp = mpi.bcast(dc_imp) + dc_energ = mpi.bcast(dc_energ) + SK.set_mu(chemical_potential) + SK.set_dc(dc_imp,dc_energ) + +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_int=h_int, **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_tau'] = S.G_tau + ar['G_iw'] = S.G_iw + ar['Sigma_iw'] = S.Sigma_iw + ar['G0-%s'%(iteration_number)] = S.G0_iw + ar['G-%s'%(iteration_number)] = S.G_iw + ar['Sigma-%s'%(iteration_number)] = 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']) + +if mpi.is_master_node(): + ar = HDFArchive("dftdmft.h5",'w') + ar["G_tau"] = S.G_tau + ar["G_iw"] = S.G_iw + ar["Sigma_iw"] = S.Sigma_iw