mirror of
https://github.com/triqs/dft_tools
synced 2024-11-07 06:33:48 +01:00
self consistent part updated
This commit is contained in:
aichhorn
parent
5ebedd8d85
commit
365f77d623
@ -24,8 +24,8 @@ User guide
|
||||
guide/conversion
|
||||
guide/dftdmft_singleshot
|
||||
guide/dftdmft_selfcons
|
||||
guide/full_tutorial
|
||||
guide/analysis
|
||||
guide/full_tutorial
|
||||
guide/transport
|
||||
|
||||
|
||||
|
||||
@ -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<images_scripts/TutorialDmftproj.pdf>`.
|
||||
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<images_scripts/TutorialDmftproj.pdf>`.
|
||||
|
||||
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.
|
||||
|
||||
@ -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` ::
|
||||
|
||||
|
||||
@ -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']
|
||||
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:
|
||||
previous_runs = 0
|
||||
previous_present = False
|
||||
del ar
|
||||
|
||||
mpi.barrier()
|
||||
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()
|
||||
|
||||
@ -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,37 +21,16 @@ 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)
|
||||
|
||||
Loading…
Reference in New Issue
Block a user