From 365f77d623c99b20308378863af4a5ff6e52f5df Mon Sep 17 00:00:00 2001 From: aichhorn Date: Fri, 14 Aug 2015 15:12:35 +0200 Subject: [PATCH] self consistent part updated --- doc/documentation.rst | 2 +- doc/guide/dftdmft_selfcons.rst | 90 +++++++++++------ doc/guide/full_tutorial.rst | 96 +++++++++++------- doc/guide/images_scripts/Ce-gamma.py | 123 ++++++++++------------- doc/guide/images_scripts/Ce-gamma_DOS.py | 34 +------ 5 files changed, 173 insertions(+), 172 deletions(-) diff --git a/doc/documentation.rst b/doc/documentation.rst index caea7da3..af0720eb 100644 --- a/doc/documentation.rst +++ b/doc/documentation.rst @@ -24,8 +24,8 @@ User guide guide/conversion guide/dftdmft_singleshot guide/dftdmft_selfcons - guide/full_tutorial guide/analysis + guide/full_tutorial guide/transport diff --git a/doc/guide/dftdmft_selfcons.rst b/doc/guide/dftdmft_selfcons.rst index 55042a29..08f20133 100644 --- a/doc/guide/dftdmft_selfcons.rst +++ b/doc/guide/dftdmft_selfcons.rst @@ -6,28 +6,30 @@ Full charge self consistency Wien2k + dmftproj ----------------- -.. warning:: - TO BE UPDATED! .. warning:: Before using this tool, you should be familiar with the band-structure package :program:`Wien2k`, since - the calculation is controlled by the :program:`Wien2k` scripts! See also the :download:`dmftproj tutorial`. + the calculation is controlled by the :program:`Wien2k` scripts! Be + sure that you also understand how :program:`dmftproj` is used to + construct the Wannier functions. For this step, see either sections + :ref:`conversion`, or the extensive :download:`dmftproj manual`. In order to do charge self-consistent calculations, we have to tell the band structure program about the changes in the charge density due to correlation effects. In the following, we discuss how to use the -:program:`TRIQS` tools in combination with the :program:`Wien2k` program, although an extension to other -codes is also possible. +:program:`TRIQS` tools in combination with the :program:`Wien2k` program. -We can use the DMFT script as introduced in sections :ref:`DFTDMFTmain` and :ref:`advanced`, with a few simple +We can use the DMFT script as introduced in section :ref:`singleshot`, +with just a few simple modifications. First, in order to be compatible with the :program:`Wien2k` standards, the DMFT script has to be -named ``case.py``, where `case` is the name of the :program:`Wien2k` calculation, see the section -:ref:`interfacetowien` for details. We can then set the variable `dft_filename` dynamically:: +named :file:`case.py`, where `case` is the place holder name of the :program:`Wien2k` calculation, see the section +:ref:`conversion` for details. We can then set the variable `dft_filename` dynamically:: import os dft_filename = os.getcwd().rpartition('/')[2] -This sets the `dft_filename` to the name of the current directory. The remainder of the script is identical to -that for one-shot calculations. Only at the very end do we have to calculate the modified charge density, +This sets the `dft_filename` to the name of the current directory. The +remaining part of the script is identical to +that for one-shot calculations. Only at the very end we have to calculate the modified charge density, and store it in a format such that :program:`Wien2k` can read it. Therefore, after the DMFT loop that we saw in the previous section, we symmetrise the self energy, and recalculate the impurity Green function:: @@ -48,15 +50,20 @@ First we find the chemical potential with high precision, and after that the rou ``SK.calc_density_correction(filename)`` calculates the density matrix including correlation effects. The result is stored in the file `dft_filename.qdmft`, which is later read by the :program:`Wien2k` program. The last statement saves the chemical potential into the hdf5 archive. + We need also the correlation energy, which we evaluate by the Migdal formula:: correnerg = 0.5 * (S.G_iw * S.Sigma_iw).total_density() +Other ways of calculating the correlation energy are possible, for +instance a direct measurment of the expectation value of the +interacting hamiltonian. However, the Migdal formula works always, +independent of the solver that is used to solve the impurity problem. From this value, we substract the double counting energy:: correnerg -= SK.dc_energ[0] -and save this value too:: +and save this value in the file, too:: if (mpi.is_master_node()): f=open(dft_filename+'.qdmft','a') @@ -68,34 +75,55 @@ where you will apply this method. That is the reason why we give the index `0` i If you have more than one correlated atom in the unit cell, but all of them are equivalent atoms, you have to multiply the `correnerg` by their multiplicity before writing it to the file. The multiplicity is easily found in the main input file of the :program:`Wien2k` package, i.e. `case.struct`. In case of -non-equivalent atoms, the correlation energy has to be calculated for all of them separately (FOR EXPERTS ONLY). +non-equivalent atoms, the correlation energy has to be calculated for +all of them separately and summed up. As mentioned above, the calculation is controlled by the :program:`Wien2k` scripts and not by :program:`python` -routines. Therefore, at the command line, you start your calculation for instance by:: +routines. You should think of replacing the lapw2 part of the +:program:`Wien2k` self-consistency cycle by - me@home $ run -qdmft -i 10 + | `lapw2 -almd` + | `dmftproj` + | `pytriqs case.py` + | `lapw2 -qdmft` -The flag `-qdmft` tells the script that the density matrix including correlation effects is to be read in from the `case.qdmft` -file and that 10 self-consistency iterations are to be done. If you run the code on a parallel machine, you can specify the number of -nodes to be used with the `-np` flag:: +In other words, for the calculation of the density matrix in lapw2, we +add the DMFT corrections through our python scripts. +Therefore, at the command line, you start your calculation for instance by: - me@home $ run -qdmft -np 64 -i 10 + `me@home $ run -qdmft 1 -i 10` -In that case, you have to give the proper `MPI` execution statement, e.g. `mpiexec`, in the `run_lapw` script (see the -corresponding :program:`Wien2k` documentation). In many cases it is advisable to start from a converged one-shot -calculation. +The flag `-qdmft` tells the :program:`Wien2k` script that the density +matrix including correlation effects is to be read in from the +`case.qdmft` file, and that you want the code to run on one computing +core only. Moreover, we ask for 10 self-consistency iterations are to be +done. If you run the code on a parallel machine, you can specify the +number of nodes to be used: -For practical purposes, you keep the number of DMFT loops within one DFT cycle low, or even to `loops=1`. If you encouter + `me@home $ run -qdmft 64 -i 10` + +In that case, you will run on 64 computing cores. As standard setting, +we use `mpirun` as the proper MPI execution statement. If you happen +to have a differnet, non-standard MPI setup, you have to give the +proper MPI execution statement, in the `run_lapw` script (see the +corresponding :program:`Wien2k` documentation). + +In many cases it is advisable to start from a converged one-shot +calculation. For practical purposes, you keep the number of DMFT loops +within one DFT cycle low, or even to `loops=1`. If you encouter unstable convergence, you have to adjust the parameters such as -`loops`, `mix`, or `Delta_mix` to improve the convergence. +the number of DMFT loops, or some mixing of the self energy to improve +the convergence. -In the next section, :ref:`DFTDMFTtutorial`, we will see in a detailed -example how such a self consistent calculation is performed. - -VASP + wannier90 ----------------- - -.. warning:: - IN PROGRESS! +In the section :ref:`DFTDMFTtutorial` we will see in a detailed +example how such a self-consistent calculation is performed from scratch. +Other DFT codes +--------------- + +The extension to other DFT codes is straight forward. As described +here, one needs to implement the correlated density matrix to be used +for the calculation of the charge density. This implementation will of +course depend on the DFT package, and might be easy to do or a quite +involved project. The formalism, however, is straight forward. diff --git a/doc/guide/full_tutorial.rst b/doc/guide/full_tutorial.rst index e82b3e4e..21ace0d5 100644 --- a/doc/guide/full_tutorial.rst +++ b/doc/guide/full_tutorial.rst @@ -3,8 +3,11 @@ DFT+DMFT tutorial: Ce with Hubbard-I approximation ================================================== -In this tutorial we will perform DFT+DMFT :program:`Wien2k` calculations of the high-temperature :math:`\gamma`-phase of Ce employing the -Hubbard-I approximation for its localized *4f* shell. +In this tutorial we will perform DFT+DMFT :program:`Wien2k` +calculations from scratch, including all steps described in the +previous sections. As example, we take the high-temperature +:math:`\gamma`-phase of Ce employing the Hubbard-I approximation for +its localized *4f* shell. Wien2k setup ------------ @@ -15,7 +18,7 @@ for the :math:`\gamma`-Ce fcc structure with lattice parameter of 9.75 a.u. .. literalinclude:: Ce-gamma.struct We initalize non-magnetic :program:`Wien2k` calculations using the :program:`init` script as described in the same manual. -For this example we specify 3000 *k*-points in the full Brillouin zone +For this example we specify 3000 :math:`\mathbf{k}`-points in the full Brillouin zone and LDA exchange-correlation potential (*vxc=5*), other parameters are defaults. The Ce *4f* electrons are treated as valence states. Hence, the initialization script is executed as follows :: @@ -25,15 +28,18 @@ Hence, the initialization script is executed as follows :: and then LDA calculations of non-magnetic :math:`\gamma`-Ce are performed by launching the :program:`Wien2k` :program:`run` script. These self-consistent LDA calculations will typically take a couple of minutes. -DMFTPROJ --------- +Wannier orbitals: dmftproj +-------------------------- -Then we create :file:`Ce-gamma.indmftpr` file specifying parameters for construction of Wannier orbitals representing *4f* states: +Then we create the :file:`Ce-gamma.indmftpr` file specifying parameters for construction of Wannier orbitals representing *4f* states: .. literalinclude:: Ce-gamma.indmftpr -First three lines give the number of inequivalent sites, their multiplicity (to be in accordance with the *struct* file) and the maximum orbital quantum number :math:`l_{max}` -The following four lines describe the treatment of Ce *spdf* orbitals by the :program:`dmftproj` program :: +As we learned in the section :ref:`conversion`, the first three lines +give the number of inequivalent sites, their multiplicity (to be in +accordance with the *struct* file) and the maximum orbital quantum +number :math:`l_{max}`. +The following four lines describe the treatment of Ce *spdf* orbitals by the :program:`dmftproj` program:: complex 1 1 1 2 ! l included for each sort @@ -41,19 +47,21 @@ The following four lines describe the treatment of Ce *spdf* orbitals by the :pr 0 where `complex` is the choice for the angular basis to be used (spherical complex harmonics), in the next line we specify, for each orbital -quantum number, whether it is treated as correlated ('2') and, hence, the corresponding Wannier orbitals will be generated or uncorrelated ('1'). +quantum number, whether it is treated as correlated ('2') and, hence, the corresponding Wannier orbitals will be generated, or uncorrelated ('1'). In the latter case the :program:`dmftproj` program will generate projectors to be used in calculations of corresponding partial densities of states (see below). In the present case we choose the fourth (i. e. *f*) orbitals as correlated. The next line specify the number of irreducible representations into which a given correlated shell should be split (or '0' if no splitting is desired, as in the present case). The fourth line specifies whether the spin-orbit interaction should be switched on ('1') or off ('0', as in the present case). -Finally, the last line iof the file :: +Finally, the last line of the file :: -.40 0.40 ! Energy window relative to E_f -specify the energy window for Wannier functions' construction. For a more complete description of :program:`dmftproj` options see its manual. +specifies the energy window for Wannier functions' construction. For a +more complete description of :program:`dmftproj` options see its +manual. -To prepaire input data for :program:`dmftproj` we execute :program:`lapw2` with the `-almd` option :: +To prepare input data for :program:`dmftproj` we execute lapw2 with the `-almd` option :: x lapw2 -almd @@ -74,18 +82,15 @@ DMFT setup: Hubbard-I calculations in TRIQS -------------------------------------------- In order to run DFT+DMFT calculations within Hubbard-I we need the corresponding python script, :ref:`Ce-gamma_script`. -It is generally similar to the script for the case of DMFT calculations with the CT-QMC solver (see :ref:`advanced`), -however there are also some differences. First, instead of *pytriqs.applications.dft.solver_multiband* we import Hubbard-I solver :: +It is generally similar to the script for the case of DMFT calculations with the CT-QMC solver (see :ref:`dftdmft_singleshot`), +however there are also some differences. First difference is that we import the Hubbard-I solver by:: from pytriqs.applications.impurity_solvers.hubbard_I.hubbard_solver import Solver -The Hubbard-I solver is very fast and we do not need to take into account the DFT blocking structure or use any approximation for the *U*-matrix :: - - use_blocks = False # use bloc structure from DFT input - use_matrix = True # use the U matrix calculated from Slater coefficients instead of (U+2J, U, U-J) - -We load and convert the :program:`dmftproj` output and initialize the *SumkDFT* class as described in :ref:`DFTDMFTmain` and :ref:`advanced` and then set up the Hubbard-I solver :: - +The Hubbard-I solver is very fast and we do not need to take into account the DFT block structure or use any approximation for the *U*-matrix. +We load and convert the :program:`dmftproj` output and initialize the +:class:`SumkDFT` class as described in :ref:`conversion` and +:ref:`singleshot` and then set up the Hubbard-I solver :: S = Solver(beta = beta, l = l) @@ -106,10 +111,10 @@ The `Solver.solve(U_int, J_hund)` statement has two necessary parameters, the Hu * `Iteration_Number`: the iteration number of the DMFT loop. Used only for printing. By default `Iteration_Number=1` * `Test_Convergence`: convergence criterion. Once the self-energy is converged below `Test_Convergence` the Hubbard-I solver is not called anymore. By default `Test_Convergence=0.0001`. -We need also to introduce some changes in the DMFT loop with respect that used for CT-QMC calculations in :ref:`advanced`. +We need also to introduce some changes in the DMFT loop with respect that used for CT-QMC calculations in :ref:`singleshot`. The hybridization function is neglected in the Hubbard-I approximation, and only non-interacting level positions (:math:`\hat{\epsilon}=-\mu+\langle H^{ff} \rangle - \Sigma_{DC}`) are required. -Hence, instead of computing `S.G0` as in :ref:`advanced` we set the level positions:: +Hence, instead of computing `S.G0` as in :ref:`singleshot` we set the level positions:: # set atomic levels: eal = SK.eff_atomic_levels()[0] @@ -120,36 +125,38 @@ Green's function and then save them in the hdf5 file . Then the double counting is recalculated and the correlation energy is computed with the Migdal formula and stored in hdf5. Finally, we compute the modified charge density and save it as well as correlational correction to the total energy in -:file:`Ce-gamma.qdmft` file, which is then read by :program:`lapw2` in the case of self-consistent DFT+DMFT calculations. +:file:`Ce-gamma.qdmft` file, which is then read by lapw2 in the case of self-consistent DFT+DMFT calculations. Running single-shot DFT+DMFT calculations ------------------------------------------ After having prepared the script one may run one-shot DMFT calculations by -executing :ref:`Ce-gamma_script` with :program:`pytriqs` on a single processor:: +executing :ref:`Ce-gamma_script` with :program:`pytriqs` on a single processor: - pytriqs Ce-gamma.py + `pytriqs Ce-gamma.py` -or in parallel mode:: +or in parallel mode: - mpirun pytriqs Ce-gamma.py + `mpirun -np 64 pytriqs Ce-gamma.py` where :program:`mpirun` launches these calculations in parallel mode and enables MPI. The exact form of this command will, of course, depend on -mpi-launcher installed in your system. +mpi-launcher installed in your system, but the form above applies to +99% of the system setups. Running self-consistent DFT+DMFT calculations --------------------------------------------- -Instead of doing one-shot run one may also perform fully self-consistent -DFT+DMFT calculations, as we will do in this tutorial. We launch these -calculations as follows :: +Instead of doing a one-shot run one may also perform fully self-consistent +DFT+DMFT calculations, as we will do now. We launch these +calculations as follows : - run_triqs -qdmft + `run -qdmft 1` -where `-qdmft` flag turns on DFT+DMFT calculations with :program:`Wien2k`. We +where `-qdmft` flag turns on DFT+DMFT calculations with +:program:`Wien2k`, and one computing core. We use here the default convergence criterion in :program:`Wien2k` (convergence to 0.1 mRy in energy). @@ -180,8 +187,9 @@ Post-processing and data analysis Within Hubbard-I one may also easily obtain the angle-resolved spectral function (band structure) and integrated spectral function (density of states or DOS). In -difference with the CT-QMC approach one does not need to provide the -real-frequency self-energy (see :ref:`analysis`) as it can be calculated directly +difference with the CT-QMC approach one does not need to do an +analytic continuations to get the +real-frequency self-energy, as it can be calculated directly in the Hubbard-I solver. The corresponding script :ref:`Ce-gamma_DOS_script` contains several new parameters :: @@ -191,18 +199,28 @@ The corresponding script :ref:`Ce-gamma_DOS_script` contains several new paramet N_om=2001 # number of points on the real-energy axis mesh broadening = 0.02 # broadening (the imaginary shift of the real-energy mesh) -Then one needs to load projectors needed for calculations of corresponding projected densities of states, as well as corresponding symmetries. To get access to analysing tools we initialize the `SumkDFTTools` class :: +Then one needs to load projectors needed for calculations of +corresponding projected densities of states, as well as corresponding +symmetries:: + + Converter.convert_parpoj_input() + +To get access to analysing tools we initialize the +:class:`SumkDFTTools` class :: SK = SumkDFTTools(hdf_file=dft_filename+'.h5', use_dft_blocks=False) -Then after the solver initialization we load the previously calculated chemical potential and double-counting correction. Having set up atomic levels we then compute the atomic Green's function and self-energy on the real axis:: +After the solver initialization, we load the previously calculated +chemical potential and double-counting correction. Having set up +atomic levels we then compute the atomic Green's function and +self-energy on the real axis:: S.set_atomic_levels( eal = eal ) S.GF_realomega(ommin=ommin, ommax = ommax, N_om=N_om,U_int=U_int,J_hund=J_hund) put it into SK class and then calculated the actual DOS:: - SK.dos_partial(broadening=broadening) + SK.dos_parproj_basis(broadening=broadening) We may first increase the number of **k**-points in BZ to 10000 by executing :program:`Wien2k` program :program:`kgen` :: diff --git a/doc/guide/images_scripts/Ce-gamma.py b/doc/guide/images_scripts/Ce-gamma.py index dfd3be4b..b347a699 100644 --- a/doc/guide/images_scripts/Ce-gamma.py +++ b/doc/guide/images_scripts/Ce-gamma.py @@ -2,45 +2,43 @@ from pytriqs.applications.dft.sumk_dft import * from pytriqs.applications.dft.converters.wien2k_converter import * from pytriqs.applications.impurity_solvers.hubbard_I.hubbard_solver import Solver -lda_filename = 'Ce-gamma' + +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 Mix = 0.7 # Mixing factor in QMC DC_type = 0 # 0...FLL, 1...Held, 2... AMF, 3...Lichtenstein -DC_Mix = 1.0 # 1.0 ... all from imp; 0.0 ... all from Gloc -useBlocs = False # use bloc structure from LDA input -useMatrix = True # use the U matrix calculated from Slater coefficients instead of (U+2J, U, U-J) chemical_potential_init=0.0 # initial chemical potential -HDFfilename = lda_filename+'.h5' +HDFfilename = dft_filename+'.h5' # Convert DMFT input: -# Can be commented after the first run -Converter = Wien2kConverter(filename=lda_filename) +Converter = Wien2kConverter(filename=filename) Converter.convert_dft_input() +mpi.barrier() #check if there are previous runs: previous_runs = 0 previous_present = False - if mpi.is_master_node(): - ar = HDFArchive(HDFfilename,'a') - if 'iterations' in ar: - previous_present = True - previous_runs = ar['iterations'] - else: - previous_runs = 0 - previous_present = False - del ar - -mpi.barrier() + f = HDFArchive(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) # Init the SumK class -SK=SumkDFT(hdf_file=lda_filename+'.h5',use_dft_blocks=False) +SK=SumkDFT(hdf_file=dft_filename+'.h5',use_dft_blocks=False) Norb = SK.corr_shells[0]['dim'] l = SK.corr_shells[0]['l'] @@ -50,44 +48,35 @@ S = Solver(beta = beta, l = l) chemical_potential=chemical_potential_init # load previous data: old self-energy, chemical potential, DC correction -if (previous_present): - mpi.report("Using stored data for initialisation") - if (mpi.is_master_node()): - ar = HDFArchive(HDFfilename,'a') - S.Sigma <<= ar['SigmaF'] - del ar - things_to_load=['chemical_potential','dc_imp'] - old_data=SK.load(things_to_load) - chemical_potential=old_data[0] - SK.dc_imp=old_data[1] - S.Sigma = mpi.bcast(S.Sigma) - chemical_potential=mpi.bcast(chemical_potential) - SK.dc_imp=mpi.bcast(SK.dc_imp) - +if previous_present: + if mpi.is_master_node(): + S.Sigma << HDFArchive(dft_filename+'.h5','a')['dmft_output']['Sigma'] + chemical_potential,dc_imp,dc_energ = SK.load(['chemical_potential','dc_imp','dc_energ']) + S.Sigma << mpi.bcast(S.Sigma) + SK.set_mu(chemical_potential) + SK.set_dc(dc_imp,dc_energ) + # DMFT loop: -for Iteration_Number in range(1,Loops+1): +for iteration_number in range(1,Loops+1): - itn = Iteration_Number + previous_runs + itn = iteration_number + previous_runs # put Sigma into the SumK class: SK.put_Sigma(Sigma_imp = [ S.Sigma ]) # Compute the SumK, possibly fixing mu by dichotomy - if SK.density_required and (Iteration_Number > 1): - chemical_potential = SK.calc_mu( precision = 0.000001 ) - else: - mpi.report("No adjustment of chemical potential\nTotal density = %.3f"%SK.total_density(mu=chemical_potential)) - + 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)): + 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) + 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] @@ -96,18 +85,17 @@ for Iteration_Number in range(1,Loops+1): # solve it: S.solve(U_int = U_int, J_hund = J_hund, verbosity = 1) - # Now mix Sigma and G: - if ((itn>1)or(previous_present)): - if (mpi.is_master_node()and (Mix<1.0)): - ar = HDFArchive(HDFfilename,'r') - mpi.report("Mixing Sigma and G with factor %s"%Mix) - if ('SigmaF' in ar): - S.Sigma <<= Mix * S.Sigma + (1.0-Mix) * ar['SigmaF'] - if ('GF' in ar): - S.G <<= Mix * S.G + (1.0-Mix) * ar['GF'] + # Now mix Sigma and G with factor Mix, if wanted: + if (iteration_number>1 or previous_present): + if (mpi.is_master_node() and (sigma_mix<1.0)): + ar = HDFArchive(dft_filename+'.h5','a')['dmft_output'] + mpi.report("Mixing Sigma and G with factor %s"%sigma_mix) + S.Sigma << sigma_mix * S.Sigma + (1.0-sigma_mix) * ar['Sigma'] + S.G << sigma_mix * S.G + (1.0-sigma_mix) * ar['G'] del ar - S.G = mpi.bcast(S.G) - S.Sigma = mpi.bcast(S.Sigma) + 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() @@ -118,19 +106,16 @@ for Iteration_Number in range(1,Loops+1): mpi.report("Corr. energy = %s"%correnerg) # store the impurity self-energy, GF as well as correlation energy in h5 - if (mpi.is_master_node()): - ar = HDFArchive(HDFfilename,'a') - ar['iterations'] = itn - ar['chemical_cotential%s'%itn] = chemical_potential - ar['SigmaF'] = S.Sigma - ar['GF'] = S.G - ar['correnerg%s'%itn] = correnerg - ar['DCenerg%s'%itn] = SK.dc_energ + 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'] = S.G + ar['Sigma'] = S.Sigma del ar #Save essential SumkDFT data: - things_to_save=['chemical_potential','dc_energ','dc_imp'] - SK.save(things_to_save) + SK.save(['chemical_potential','dc_imp','dc_energ','correnerg']) if (mpi.is_master_node()): print 'DC after solver: ',SK.dc_imp[0] @@ -150,22 +135,16 @@ for Iteration_Number in range(1,Loops+1): # find exact chemical potential -if (SK.density_required): - SK.chemical_potential = SK.calc_mu( precision = 0.000001 ) +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 = lda_filename+'.qdmft') +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()): - ar = HDFArchive(HDFfilename) - itn = ar['iterations'] - correnerg = ar['correnerg%s'%itn] - DCenerg = ar['DCenerg%s'%itn] - del ar correnerg -= DCenerg[0] - f=open(lda_filename+'.qdmft','a') + f=open(dft_filename+'.qdmft','a') f.write("%.16f\n"%correnerg) f.close() diff --git a/doc/guide/images_scripts/Ce-gamma_DOS.py b/doc/guide/images_scripts/Ce-gamma_DOS.py index 0f7a2d4a..b9b92aae 100644 --- a/doc/guide/images_scripts/Ce-gamma_DOS.py +++ b/doc/guide/images_scripts/Ce-gamma_DOS.py @@ -9,9 +9,6 @@ beta = 40 U_int = 6.00 J_hund = 0.70 DC_type = 0 # 0...FLL, 1...Held, 2... AMF, 3...Lichtenstein -load_previous = True # load previous results -useBlocs = False # use bloc structure from DFT input -useMatrix = True # use the U matrix calculated from Slater coefficients instead of (U+2J, U, U-J) ommin=-4.0 ommax=6.0 N_om=2001 @@ -24,38 +21,17 @@ Converter = Wien2kConverter(filename=dft_filename,repacking=True) Converter.convert_dft_input() Converter.convert_parproj_input() -#check if there are previous runs: -previous_runs = 0 -previous_present = False -if mpi.is_master_node(): - ar = HDFArchive(HDFfilename) - if 'iterations' in ar: - previous_present = True - previous_runs = ar['iterations'] - else: - previous_runs = 0 - previous_present = False - del ar - -mpi.barrier() -previous_runs = mpi.bcast(previous_runs) -previous_present = mpi.bcast(previous_present) # Init the SumK class SK = SumkDFTTools(hdf_file=dft_filename+'.h5',use_dft_blocks=False) # load old chemical potential and DC -chemical_potential=0.0 if mpi.is_master_node(): - ar = HDFArchive(HDFfilename) - things_to_load=['chemical_potential','dc_imp'] - old_data=SK.load(things_to_load) - chemical_potential=old_data[0] - SK.dc_imp=old_data[1] -SK.chemical_potential=mpi.bcast(chemical_potential) -SK.dc_imp=mpi.bcast(SK.dc_imp) - + chemical_potential,dc_imp,dc_energ = SK.load(['chemical_potential','dc_imp','dc_energ']) +SK.set_mu(chemical_potential) +SK.set_dc(dc_imp,dc_energ) + if (mpi.is_master_node()): print 'DC after reading SK: ',SK.dc_imp[0] @@ -74,4 +50,4 @@ S.GF_realomega(ommin=ommin, ommax = ommax, N_om=N_om,U_int=U_int,J_hund=J_hund) SK.put_Sigma(Sigma_imp = [S.Sigma]) # compute DOS -SK.dos_partial(broadening=broadening) +SK.dos_parproj_basis(broadening=broadening)