3
0
mirror of https://github.com/triqs/dft_tools synced 2024-12-30 16:15:47 +01:00

Merge branch 'master' into vasp

Conflicts:
	CMakeLists.txt
	python/CMakeLists.txt
	python/converters/__init__.py
	test/CMakeLists.txt
This commit is contained in:
Oleg E. Peil 2016-03-11 09:53:24 +01:00
commit f142db96fb
44 changed files with 4870 additions and 112 deletions

View File

@ -14,14 +14,19 @@ enable_testing()
# Load TRIQS, including all predefined variables from TRIQS installation
find_package(TRIQS REQUIRED)
if (NOT ${TRIQS_WITH_PYTHON_SUPPORT})
MESSAGE(FATAL_ERROR "dft_tools require Python support in TRIQS")
endif()
# Check that versions are compatible
if(NOT DFT_TOOLS_VERSION EQUAL TRIQS_VERSION)
message(FATAL_ERROR "The application version is not compatible with the TRIQS library (TRIQS library version: ${TRIQS_VERSION} while this application version: ${DFT_TOOLS_VERSION})")
endif()
if (NOT ${TRIQS_WITH_PYTHON_SUPPORT})
MESSAGE(FATAL_ERROR "dft_tools require Python support in TRIQS")
endif()
# Get hash
triqs_get_git_hash(${CMAKE_SOURCE_DIR} "DFT_TOOLS")
if(${GIT_RESULT} EQUAL 0)
message(STATUS "Hash: ${DFT_TOOLS_GIT_HASH}")
endif(${GIT_RESULT} EQUAL 0)
# We want to be installed in the TRIQS tree
set(CMAKE_INSTALL_PREFIX ${TRIQS_PATH})

8
cmake/sitecustomize.py Normal file
View File

@ -0,0 +1,8 @@
def application_pytriqs_import(name,*args,**kwargs):
if name.startswith('@package_name@'):
name = name[len('@package_name@')+1:]
return builtin_import(name,*args,**kwargs)
import __builtin__
__builtin__.__import__, builtin_import = application_pytriqs_import, __builtin__.__import__

View File

@ -77,7 +77,7 @@ 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
SK.set_Sigma([ S.Sigma_iw ]) # set 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())

View File

@ -2,6 +2,9 @@
#
# TRIQS documentation build configuration file
import sys
sys.path.insert(0, "@TRIQS_SPHINXEXT_PATH@/numpydoc")
extensions = ['sphinx.ext.autodoc',
'sphinx.ext.mathjax',
'sphinx.ext.intersphinx',

View File

@ -3,13 +3,14 @@ Frequently-Asked Questions
wien2k: FERMI ERROR when running `x lapw2 -almd -band`
------------------------------------------------------
In some versions of Wien2k, there is a problem in running `x lapw2 -almd -band`.
A hack solution is as follows:
1) `x lapw1 -band`
2) edit in2 file: replace 'TOT' with 'QTL', 'TETRA' with 'ROOT'
3) `x lapw2 -almd -band`
4) `dmftproj -band` (add the fermi energy to file, it can be found by running `grep :FER *.scf`)
4) `dmftproj -band` (add the Fermi energy to file, it can be found by running `grep :FER *.scf`)
How do I plot the output of `spaghettis`?
-----------------------------------------
@ -20,6 +21,20 @@ the parameters as desired.
.. literalinclude:: plotting_spaghettis.py
x optic does not write a case.pmat file
---------------------------------------
Make sure that you set line 6 to "ON" and put a "1" to the following line.
The "1" is undocumented in Wien2k, but needed to have `case.pmat` written.
However, we are working on reading directly the `case.mommat2` file.
No module named pytriqs.*** error when running a script
-------------------------------------------------------
Make sure that have properly build, tested and installed TRIQS and DFTTools
using, make, make test and make install. Additionally, you should always
use pytriqs to call your scripts, e.g. pytriqs yourscript.py
Why is my calculation not working?
----------------------------------

View File

@ -64,9 +64,9 @@ where:
It is important that each data file has to contain three columns: the real frequency mesh, the real part and the imaginary part
of the self energy - exactly in this order! The mesh should be the same for all files read in and non-uniform meshes are not supported.
Finally, we put the self energy into the `SK` object::
Finally, we set the self energy into the `SK` object::
SK.put_Sigma(Sigma_imp = [SigmaReFreq])
SK.set_Sigma([SigmaReFreq])
and additionally set the chemical potential and the double counting correction from the DMFT calculation::
@ -96,7 +96,7 @@ the output is printed into the files
* `DOS_wannier_(sp)_proj(i)_(m)_(n).dat`: As above, but printed as orbitally-resolved matrix in indices
`(m)` and `(n)`. For `d` orbitals, it gives the DOS separately for, e.g., :math:`d_{xy}`, :math:`d_{x^2-y^2}`, and so on,
otherwise, the ouptput is returned by the function for a further usage in :program:`python`.
otherwise, the output is returned by the function for a further usage in :program:`python`.
Partial charges
---------------
@ -104,7 +104,7 @@ Partial charges
Since we can calculate the partial charges directly from the Matsubara Green's functions, we also do not need a
real frequency self energy for this purpose. The calculation is done by::
SK.put_Sigma(Sigma_imp = SigmaImFreq)
SK.set_Sigma(SigmaImFreq)
dm = SK.partial_charges(beta=40.0, with_Sigma=True, with_dc=True)
which calculates the partial charges using the self energy, double counting, and chemical potential as set in the

View File

@ -15,7 +15,7 @@ Interface with Wien2k
We assume that the user has obtained a self-consistent solution of the
Kohn-Sham equations. We further have to require that the user is
familiar with the main inout/output files of Wien2k, and how to run
familiar with the main in/output files of Wien2k, and how to run
the DFT code.
Conversion for the DMFT self-consistency cycle
@ -31,7 +31,7 @@ We note that any other flag for lapw2, such as -c or -so (for
spin-orbit coupling) has to be added also to this line. This creates
some files that we need for the Wannier orbital construction.
The orbital construction itself is done by the fortran program
The orbital construction itself is done by the Fortran program
:program:`dmftproj`. For an extensive manual to this program see
:download:`TutorialDmftproj.pdf <images_scripts/TutorialDmftproj.pdf>`.
Here we will only describe only the basic steps.
@ -79,7 +79,7 @@ following 3 to 5 lines:
These lines have to be repeated for each inequivalent atom.
The last line gives the energy window, relativ to the Fermi energy,
The last line gives the energy window, relative to the Fermi energy,
that is used for the projective Wannier functions. Note that, in
accordance with Wien2k, we give energies in Rydberg units!
@ -207,7 +207,7 @@ The lines of this header define
2 3. Thiw would mean, 2 irreps (eg and t2g), of dimension 2 and 3,
resp.
After these header lines, the file has to contain the hamiltonian
After these header lines, the file has to contain the Hamiltonian
matrix in orbital space. The standard convention is that you give for
each
:math:`\mathbf{k}`-point first the matrix of the real part, then the
@ -226,6 +226,116 @@ with the
For more options of this converter, have a look at the
:ref:`refconverters` section of the reference manual.
Wannier90 Converter
-------------------
Using this converter it is possible to convert the output of
:program:`Wannier90` (http://wannier.org) calculations of
Maximally Localized Wannier Functions (MLWF) and create a HDF5 archive
suitable for one-shot DMFT calculations with the
:class:`SumkDFT <pytriqs.applications.dft.sumk_dft.SumkDFT>` class.
The user must supply two files in order to run the Wannier90 Converter:
#. The file :file:`seedname_hr.dat`, which contains the DFT Hamiltonian
in the MLWF basis calculated through :program:`wannier90` with ``hr_plot = true``
(please refer to the :program:`wannier90` documentation).
#. A file named :file:`seedname.inp`, which contains the required
information about the :math:`\mathbf{k}`-point mesh, the electron density,
the correlated shell structure, ... (see below).
Here and in the following, the keyword ``seedname`` should always be intended
as a placeholder for the actual prefix chosen by the user when creating the
input for :program:`wannier90`.
Once these two files are available, one can use the converter as follows::
from pytriqs.applications.dft.converters import Wannier90Converter
Converter = Wannier90Converter(seedname='seedname')
Converter.convert_dft_input()
The converter input :file:`seedname.inp` is a simple text file with
the following format:
.. literalinclude:: images_scripts/LaVO3_w90.inp
The example shows the input for the perovskite crystal of LaVO\ :sub:`3`
in the room-temperature `Pnma` symmetry. The unit cell contains four
symmetry-equivalent correlated sites (the V atoms) and the total number
of electrons per unit cell is 8 (see second line).
The first line specifies how to generate the :math:`\mathbf{k}`-point
mesh that will be used to obtain :math:`H(\mathbf{k})`
by Fourier transforming :math:`H(\mathbf{R})`.
Currently implemented options are:
* :math:`\Gamma`-centered uniform grid with dimensions
:math:`n_{k_x} \times n_{k_y} \times n_{k_z}`;
specify ``0`` followed by the three grid dimensions,
like in the example above
* :math:`\Gamma`-centered uniform grid with dimensions
automatically determined by the converter (from the number of
:math:`\mathbf{R}` vectors found in :file:`seedname_hr.dat`);
just specify ``-1``
Inside :file:`seedname.inp`, it is crucial to correctly specify the
correlated shell structure, which depends on the contents of the
:program:`wannier90` output :file:`seedname_hr.dat` and on the order
of the MLWFs contained in it.
The number of MLWFs must be equal to, or greater than the total number
of correlated orbitals (i.e., the sum of all ``dim`` in :file:`seedname.inp`).
If the converter finds fewer MLWFs inside :file:`seedname_hr.dat`, then it
stops with an error; if it finds more MLWFs, then it assumes that the
additional MLWFs correspond to uncorrelated orbitals (e.g., the O-\ `2p` shells).
When reading the hoppings :math:`\langle w_i | H(\mathbf{R}) | w_j \rangle`
(where :math:`w_i` is the :math:`i`-th MLWF), the converter also assumes that
the first indices correspond to the correlated shells (in our example,
the V-t\ :sub:`2g` shells). Therefore, the MLWFs corresponding to the
uncorrelated shells (if present) must be listed **after** those of the
correlated shells.
With the :program:`wannier90` code, this can be achieved this by listing the
projections for the uncorrelated shells after those for the correlated shells.
In our `Pnma`-LaVO\ :sub:`3` example, for instance, we could use::
Begin Projections
V:l=2,mr=2,3,5:z=0,0,1:x=-1,1,0
O:l=1:mr=1,2,3:z=0,0,1:x=-1,1,0
End Projections
where the ``x=-1,1,0`` option indicates that the V--O bonds in the octahedra are
rotated by (approximatively) 45 degrees with respect to the axes of the `Pbnm` cell.
The converter will analyse the matrix elements of the local hamiltonian
to find the symmetry matrices `rot_mat` needed for the global-to-local
transformation of the basis set for correlated orbitals
(see section :ref:`hdfstructure`).
The matrices are obtained by finding the unitary transformations that diagonalize
:math:`\langle w_i | H_I(\mathbf{R}=0,0,0) | w_j \rangle`, where :math:`I` runs
over the correlated shells and `i,j` belong to the same shell (more details elsewhere...).
If two correlated shells are defined as equivalent in :file:`seedname.inp`,
then the corresponding eigenvalues have to match within a threshold of 10\ :sup:`-5`,
otherwise the converter will produce an error/warning.
If this happens, please carefully check your data in :file:`seedname_hr.dat`.
This method might fail in non-trivial cases (i.e., more than one correlated
shell is present) when there are some degenerate eigenvalues:
so far tests have not shown any issue, but one must be careful in those cases
(the converter will print a warning message).
The current implementation of the Wannier90 Converter has some limitations:
* Since :program:`wannier90` does not make use of symmetries (symmetry-reduction
of the :math:`\mathbf{k}`-point grid is not possible), the converter always
sets ``symm_op=0`` (see the :ref:`hdfstructure` section).
* No charge self-consistency possible at the moment.
* Calculations with spin-orbit (``SO=1``) are not supported.
* The spin-polarized case (``SP=1``) is not yet tested.
* The post-processing routines in the module
:class:`SumkDFTTools <pytriqs.applications.dft.sumk_dft_tools.SumkDFTTools>`
were not tested with this converter.
* ``proj_mat_all`` are not used, so there are no projectors onto the
uncorrelated orbitals for now.
MPI issues
----------

View File

@ -41,7 +41,7 @@ These steps are not necessary, but can help to reduce fluctuations in the total
Now we calculate the modified charge density::
# find exact chemical potential
SK.put_Sigma(Sigma_imp = [ S.Sigma_iw ])
SK.set_Sigma([ S.Sigma_iw ])
chemical_potential = SK.calc_mu( precision = 0.000001 )
dN, d = SK.calc_density_correction(filename = dft_filename+'.qdmft')
SK.save(['chemical_potential','dc_imp','dc_energ'])
@ -56,10 +56,10 @@ 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,
instance a direct measurement 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::
From this value, we subtract the double counting energy::
correnerg -= SK.dc_energ[0]
@ -104,13 +104,13 @@ number of nodes to be used:
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
to have a different, 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
within one DFT cycle low, or even to `loops=1`. If you encounter
unstable convergence, you have to adjust the parameters such as
the number of DMFT loops, or some mixing of the self energy to improve
the convergence.

View File

@ -50,7 +50,7 @@ iterations and the self-consistency condition::
n_loops = 5
for iteration_number in range(n_loops) : # start the DMFT loop
SK.put_Sigma(Sigma_imp = [ S.Sigma ]) # Put self energy to the SumK class
SK.set_Sigma([ S.Sigma ]) # Put self energy to the SumK class
chemical_potential = SK.calc_mu() # calculate the chemical potential for the given density
S.G_iw << SK.extract_G_loc()[0] # extract the local Green function
S.G0_iw << inverse(S.Sigma_iw + inverse(S.G_iw)) # finally get G0, the input for the Solver
@ -107,7 +107,7 @@ 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
<images_scripts/dft_dmft_cthyb.py>`) and one with a
rotational-invariant Slater interaction hamiltonian (:download:`dft_dmft_cthyb_slater.py
rotational-invariant Slater interaction Hamiltonian (:download:`dft_dmft_cthyb_slater.py
<images_scripts/dft_dmft_cthyb.py>`). The user has to adapt these
scripts to his own needs.
@ -145,7 +145,7 @@ Most of these parameters are self-explanatory. The first,
details on the solver parameters, we refer the user to
the :ref:`CTHYB solver <triqscthyb:welcome>` documentation.
We assume that the conversion to the hdf5 archive is alreadz done. We
We assume that the conversion to the hdf5 archive is already done. We
can check now in this archive, if previous runs are present, or if we have to start
from scratch::
@ -165,7 +165,7 @@ from scratch::
previous_present = mpi.bcast(previous_present)
You can see in this code snipet, that all results of this calculation
You can see in this code snippet, that all results of this calculation
will be stored in a separate subgroup in the hdf5 file, called
`dmft_output`. Removing this subgroup allows you to reset your
calculation to the starting point easily.
@ -178,7 +178,7 @@ The next step is to initialise the :class:`Solver <pytriqs.applications.impurit
of two steps
#. Calculating the multi-band interaction matrix, and setting up the
interaction hamiltonian
interaction Hamiltonian
#. Setting up the solver class
The first step is done using methods of
@ -199,13 +199,13 @@ other choices (Slater interaction matrix for instance), and other
parameters, we refer to the reference manual
of the :ref:`TRIQS <triqslibs:welcome>` library.
Next, we construct the hamiltonian and the solver::
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
account. Other choices for the Hamiltonian are
* h_int_kanamori
* h_int_slater
@ -239,7 +239,7 @@ refinements::
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
SK.set_Sigma([ 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())

View File

@ -73,7 +73,7 @@ This program produces the following files:
* :file:`Ce-gamma.ctqmcout` and :file:`Ce-gamma.symqmc` containing projector operators and symmetry operations for orthonormalized Wannier orbitals, respectively.
* :file:`Ce-gamma.parproj` and :file:`Ce-gamma.sympar` containing projector operators and symmetry operations for uncorrelated states, respectively. These files are needed for projected density-of-states or spectral-function calculations.
* :file:`Ce-gamma.oubwin` needed for the charge desity recalculation in the case of fully self-consistent DFT+DMFT run (see below).
* :file:`Ce-gamma.oubwin` needed for the charge density recalculation in the case of fully self-consistent DFT+DMFT run (see below).
Now we have all necessary input from :program:`Wien2k` for running DMFT calculations.
@ -101,9 +101,9 @@ The Hubbard-I initialization `Solver` has also optional parameters one may use:
* `n_msb`: the number of Matsubara frequencies used. The default is `n_msb=1025`.
* `use_spin_orbit`: if set 'True' the solver is run with spin-orbit coupling included. To perform actual DFT+DMFT calculations with spin-orbit one should also run :program:`Wien2k` and :program:`dmftproj` in spin-polarized mode and with spin-orbit included. By default, `use_spin_orbit=False`.
* `Nmoments`: the number of moments used to describe high-ferquency tails of the Hubbard-I Green's function and self-energy. By default `Nmoments = 5`
* `Nmoments`: the number of moments used to describe high-frequency tails of the Hubbard-I Green's function and self-energy. By default `Nmoments = 5`
The `Solver.solve(U_int, J_hund)` statement has two necessary parameters, the Hubbard U parameter `U_int` and Hund's rule coupling `J_hund`. Notice that the solver constructs the full 4-index `U`-matrix by default, and the `U_int` parameter is in fact the Slatter `F0` integral. Other optional parameters are:
The `Solver.solve(U_int, J_hund)` statement has two necessary parameters, the Hubbard U parameter `U_int` and Hund's rule coupling `J_hund`. Notice that the solver constructs the full 4-index `U`-matrix by default, and the `U_int` parameter is in fact the Slater `F0` integral. Other optional parameters are:
* `T`: matrix that transforms the interaction matrix from complex spherical harmonics to a symmetry adapted basis. By default, the complex spherical harmonics basis is used and `T=None`.
* `verbosity`: tunes output from the solver. If `verbosity=0` only basic information is printed, if `verbosity=1` the ground state atomic occupancy and its energy are printed, if `verbosity=2` additional information is printed for all occupancies that were diagonalized. By default, `verbosity=0`.

View File

@ -2,7 +2,6 @@ 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
import os
dft_filename = os.getcwd().rpartition('/')[2]
@ -10,14 +9,12 @@ beta = 40
U_int = 6.00
J_hund = 0.70
Loops = 5 # Number of DMFT sc-loops
Mix = 0.7 # Mixing factor in QMC
mixing = 0.7 # Mixing factor
DC_type = 0 # 0...FLL, 1...Held, 2... AMF, 3...Lichtenstein
chemical_potential_init=0.0 # initial chemical potential
HDFfilename = dft_filename+'.h5'
# Convert DMFT input:
Converter = Wien2kConverter(filename=filename)
Converter = Wien2kConverter(filename=dft_filename)
Converter.convert_dft_input()
mpi.barrier()
@ -25,7 +22,7 @@ mpi.barrier()
previous_runs = 0
previous_present = False
if mpi.is_master_node():
f = HDFArchive(filename+'.h5','a')
f = HDFArchive(dft_filename+'.h5','a')
if 'dmft_output' in f:
ar = f['dmft_output']
if 'iterations' in ar:
@ -53,11 +50,11 @@ if previous_present:
ar = HDFArchive(dft_filename+'.h5','a')
S.Sigma << ar['dmft_output']['Sigma']
del ar
chemical_potential,dc_imp,dc_energ = SK.load(['chemical_potential','dc_imp','dc_energ'])
SK.chemical_potential,SK.dc_imp,SK.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)
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):
@ -65,7 +62,7 @@ for iteration_number in range(1,Loops+1):
itn = iteration_number + previous_runs
# put Sigma into the SumK class:
SK.put_Sigma(Sigma_imp = [ S.Sigma ])
SK.set_Sigma([ S.Sigma ])
# Compute the SumK, possibly fixing mu by dichotomy
chemical_potential = SK.calc_mu( precision = 0.000001 )
@ -78,7 +75,7 @@ for iteration_number in range(1,Loops+1):
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]
@ -89,11 +86,11 @@ for iteration_number in range(1,Loops+1):
# 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)):
if (mpi.is_master_node() and (mixing<1.0)):
ar = HDFArchive(dft_filename+'.h5','a')
mpi.report("Mixing Sigma and G with factor %s"%sigma_mix)
S.Sigma << sigma_mix * S.Sigma + (1.0-sigma_mix) * ar['dmft_output']['Sigma']
S.G << sigma_mix * S.G + (1.0-sigma_mix) * ar['dmft_output']['G']
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']
del ar
S.G << mpi.bcast(S.G)
S.Sigma << mpi.bcast(S.Sigma)
@ -104,8 +101,8 @@ for iteration_number in range(1,Loops+1):
SK.calc_dc( dm, U_interact = U_int, J_hund = J_hund, orb = 0, use_dc_formula = DC_type )
# correlation energy calculations:
correnerg = 0.5 * (S.G * S.Sigma).total_density()
mpi.report("Corr. energy = %s"%correnerg)
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():
@ -145,7 +142,7 @@ 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()):
correnerg -= DCenerg[0]
SK.correnerg -= SK.dc_energ[0]
f=open(dft_filename+'.qdmft','a')
f.write("%.16f\n"%correnerg)
f.write("%.16f\n"%SK.correnerg)
f.close()

View File

@ -14,24 +14,22 @@ ommax=6.0
N_om=2001
broadening = 0.02
HDFfilename = dft_filename+'.h5'
# Convert DMFT input:
Converter = Wien2kConverter(filename=dft_filename,repacking=True)
Converter.convert_dft_input()
Converter.convert_parproj_input()
# Init the SumK class
SK = SumkDFTTools(hdf_file=dft_filename+'.h5',use_dft_blocks=False)
# load old chemical potential and DC
if mpi.is_master_node():
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)
SK.chemical_potential,SK.dc_imp,SK.dc_energ = SK.load(['chemical_potential','dc_imp','dc_energ'])
SK.chemical_potential = mpi.bcast(SK.chemical_potential)
SK.dc_imp = mpi.bcast(SK.dc_imp)
SK.dc_energ = mpi.bcast(SK.dc_energ)
if (mpi.is_master_node()):
print 'DC after reading SK: ',SK.dc_imp[0]
@ -47,7 +45,7 @@ S.set_atomic_levels( eal = eal )
# Run the solver to get GF and self-energy on the real axis
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])
SK.set_Sigma([S.Sigma])
# compute DOS
SK.dos_parproj_basis(broadening=broadening)

View File

@ -0,0 +1,7 @@
0 6 4 6
8.0
4
0 0 2 3 0 0
1 0 2 3 0 0
2 0 2 3 0 0
3 0 2 3 0 0

View File

@ -82,7 +82,7 @@ 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
SK.set_Sigma([ S.Sigma_iw ]) # set 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())

View File

@ -83,7 +83,7 @@ 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
SK.set_Sigma([ S.Sigma_iw ]) # set 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())

View File

@ -1,6 +1,6 @@
.. _Transport:
Transport calculations
Transport calculations test
======================
Formalism
@ -65,7 +65,10 @@ The basics steps to calculate the matrix elements of the momentum operator with
6) Run `x optic`.
Additionally the input file :file:`case.inop` is required. A detail description on how to setup this file can be found in the Wien2k user guide [#userguide]_ on page 166.
Here the energy window can be chosen according to the window used for :program:`dmftproj`. However, keep in mind that energies have to be specified in absolute values! Furthermore it is important to set line 6 to ON for printing the matrix elements to the :file:`.pmat` file.
The optics energy window should be chosen according to the window used for :program:`dmftproj`. Note that the current version of the transport code uses only the smaller
of those two windows. However, keep in mind that the optics energy window has to be specified in absolute values and NOT relative to the Fermi energy!
You can read off the Fermi energy from the :file:`case.scf2` file. Please do not set the optional parameter NBvalMAX in :file:`case.inop`.
Furthermore it is necessary to set line 6 to "ON" and put a "1" in the following line to enable the printing of the matrix elements to :file:`case.pmat`.
Using the transport code
@ -86,7 +89,7 @@ reads the required data of the Wien2k output and stores it in the `dft_transp_in
Additionally we need to read and set the self energy, the chemical potential and the double counting::
ar = HDFArchive('case.h5', 'a')
SK.put_Sigma(Sigma_imp = [ar['dmft_output']['Sigma_w']])
SK.set_Sigma([ar['dmft_output']['Sigma_w']])
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)

View File

@ -5,19 +5,25 @@ Converters
Wien2k Converter
----------------
.. autoclass:: pytriqs.applications.dft.converters.wien2k_converter.Wien2kConverter
.. autoclass:: dft.converters.wien2k_converter.Wien2kConverter
:members:
:special-members:
:show-inheritance:
H(k) Converter
--------------
.. autoclass:: pytriqs.applications.dft.converters.hk_converter.HkConverter
.. autoclass:: dft.converters.hk_converter.HkConverter
:members:
:special-members:
Wannier90 Converter
--------------
.. autoclass:: dft.converters.wannier90_converter.Wannier90Converter
:members:
:special-members:
Converter Tools
---------------
.. autoclass:: pytriqs.applications.dft.converters.converter_tools.ConverterTools
.. autoclass:: dft.converters.converter_tools.ConverterTools
:members:
:special-members:

View File

@ -2,7 +2,7 @@ SumK DFT
========
.. autoclass:: pytriqs.applications.dft.sumk_dft.SumkDFT
.. autoclass:: sumk_dft.SumkDFT
:members:
:special-members:
:show-inheritance:

View File

@ -2,7 +2,7 @@ SumK DFT Tools
==============
.. autoclass:: pytriqs.applications.dft.sumk_dft_tools.SumkDFTTools
.. autoclass:: sumk_dft_tools.SumkDFTTools
:members:
:special-members:
:show-inheritance:

View File

@ -1,6 +1,6 @@
Symmetry
========
.. autoclass:: pytriqs.applications.dft.Symmetry
.. autoclass:: Symmetry
:members:
:special-members:

View File

@ -1,6 +1,6 @@
TransBasis
==========
.. autoclass:: pytriqs.applications.dft.trans_basis.TransBasis
.. autoclass:: trans_basis.TransBasis
:members:
:special-members:

View File

@ -1,8 +1,11 @@
# where will the python end up in triqs?
set(python_destination applications/dft)
set(python_destination pytriqs/applications/dft)
# site_customize for build
set(package_name "pytriqs.applications")
configure_file(${CMAKE_SOURCE_DIR}/cmake/sitecustomize.py ${CMAKE_CURRENT_BINARY_DIR}/sitecustomize.py @ONLY)
# make a local pytriqs copy
#triqs_prepare_local_pytriqs_merged_with_my_python(${python_destination})
triqs_prepare_local_pytriqs(${python_destination})
# to be able to run from toplevel

View File

@ -23,7 +23,8 @@
from wien2k_converter import Wien2kConverter
from hk_converter import HkConverter
from vasp_converter import VaspConverter
from wannier90_converter import Wannier90Converter
__all__ =['Wien2kConverter','HkConverter','VaspConverter']
__all__ =['Wien2kConverter','HkConverter','Wannier90Converter','VaspConverter']

View File

@ -0,0 +1,578 @@
##########################################################################
#
# TRIQS: a Toolbox for Research in Interacting Quantum Systems
#
# Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola
#
# TRIQS is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.
#
# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# TRIQS. If not, see <http://www.gnu.org/licenses/>.
#
##########################################################################
###
# Wannier90 to HDF5 converter for the SumkDFT class of dfttools/TRIQS;
#
# written by Gabriele Sclauzero (Materials Theory, ETH Zurich), Dec 2015 -- Jan 2016,
# under the supervision of Claude Ederer (Materials Theory).
# Partially based on previous work by K. Dymkovski and the DFT_tools/TRIQS team.
#
# Limitations of the current implementation:
# - the case with SO=1 is not considered at the moment
# - the T rotation matrices are not used in this implementation
# - projectors for uncorrelated shells (proj_mat_all) cannot be set
#
# Things to be improved/checked:
# - the case with SP=1 might work, but was never tested (do we need to define
# rot_mat_time_inv also if symm_op = 0?)
# - the calculation of rot_mat in find_rot_mat() relies on the eigenvalues of H(0);
# this might fail in presence of degenerate eigenvalues (now just prints warning)
# - the FFT is always done in serial mode (because all converters run serially);
# this can become very slow with a large number of R-vectors/k-points
# - make the code more MPI safe (error handling): if we run with more than one process
# and an error occurs on the masternode, the calculation does not abort
###
from types import *
import numpy
import math
from pytriqs.archive import *
from converter_tools import *
from itertools import product
import os.path
class Wannier90Converter(ConverterTools):
"""
Conversion from Wannier90 output to an hdf5 file that can be used as input for the SumkDFT class.
"""
def __init__(self, seedname, hdf_filename=None, dft_subgrp='dft_input',
symmcorr_subgrp='dft_symmcorr_input', repacking=False):
"""
Initialise the class.
Parameters
----------
seedname : string
Base name of Wannier90 files
hdf_filename : string, optional
Name of hdf5 archive to be created
dft_subgrp : string, optional
Name of subgroup storing necessary DFT data
symmcorr_subgrp : string, optional
Name of subgroup storing correlated-shell symmetry data
repacking : boolean, optional
Does the hdf5 archive need to be repacked to save space?
"""
self._name = "Wannier90Converter"
assert type(seedname) == StringType, self._name + \
": Please provide the DFT files' base name as a string."
if hdf_filename is None:
hdf_filename = seedname + '.h5'
self.hdf_file = hdf_filename
# if the w90 output is seedname_hr.dat, the input file for the
# converter must be called seedname.inp
self.inp_file = seedname + '.inp'
self.w90_seed = seedname
self.dft_subgrp = dft_subgrp
self.symmcorr_subgrp = symmcorr_subgrp
self.fortran_to_replace = {'D': 'E'}
# threshold below which matrix elements from wannier90 should be considered equal
self._w90zero = 2.e-6
# Checks if h5 file is there and repacks it if wanted:
if (os.path.exists(self.hdf_file) and repacking):
ConverterTools.repack(self)
def convert_dft_input(self):
"""
Reads the appropriate files and stores the data for the
- dft_subgrp
- symmcorr_subgrp
in the hdf5 archive.
"""
# Read and write only on the master node
if not (mpi.is_master_node()):
return
mpi.report("Reading input from %s..." % self.inp_file)
# R is a generator : each R.Next() will return the next number in the file
R = ConverterTools.read_fortran_file(
self, self.inp_file, self.fortran_to_replace)
shell_entries = ['atom', 'sort', 'l', 'dim']
corr_shell_entries = ['atom', 'sort', 'l', 'dim', 'SO', 'irep']
# First, let's read the input file with the parameters needed for the conversion
try:
# read k - point mesh generation option
kmesh_mode = int(R.next())
if kmesh_mode >= 0:
# read k-point mesh size from input
nki = [int(R.next()) for idir in range(3)]
else:
# some default grid, if everything else fails...
nki = [8, 8, 8]
# read the total number of electrons per cell
density_required = float(R.next())
# we do not read shells, because we have no additional shells beyond correlated ones,
# and the data will be copied from corr_shells into shells (see below)
# number of corr. shells (e.g. Fe d, Ce f) in the unit cell,
n_corr_shells = int(R.next())
# now read the information about the correlated shells (atom, sort, l, dim, SO flag, irep):
corr_shells = [{name: int(val) for name, val in zip(
corr_shell_entries, R)} for icrsh in range(n_corr_shells)]
except StopIteration: # a more explicit error if the file is corrupted.
mpi.report(self._name + ": reading input file %s failed!" %
self.inp_file)
# close the input file
R.close()
# Set or derive some quantities
# Wannier90 does not use symmetries to reduce the k-points
# the following might change in future versions
symm_op = 0
# copy corr_shells into shells (see above)
n_shells = n_corr_shells
shells = []
for ish in range(n_shells):
shells.append({key: corr_shells[ish].get(
key, None) for key in shell_entries})
###
SP = 0 # NO spin-polarised calculations for now
SO = 0 # NO spin-orbit calculation for now
charge_below = 0 # total charge below energy window NOT used for now
energy_unit = 1.0 # should be understood as eV units
###
# this is more general
n_spin = SP + 1 - SO
dim_corr_shells = sum([sh['dim'] for sh in corr_shells])
mpi.report(
"Total number of WFs expected in the correlated shells: %d" % dim_corr_shells)
# determine the number of inequivalent correlated shells and maps, needed for further processing
n_inequiv_shells, corr_to_inequiv, inequiv_to_corr = ConverterTools.det_shell_equivalence(
self, corr_shells)
mpi.report("Number of inequivalent shells: %d" % n_inequiv_shells)
mpi.report("Shell representatives: " + format(inequiv_to_corr))
shells_map = [inequiv_to_corr[corr_to_inequiv[ish]]
for ish in range(n_corr_shells)]
mpi.report("Mapping: " + format(shells_map))
# build the k-point mesh, if its size was given on input (kmesh_mode >= 0),
# otherwise it is built according to the data in the hr file (see below)
if kmesh_mode >= 0:
n_k, k_mesh, bz_weights = self.kmesh_build(nki, kmesh_mode)
self.n_k = n_k
self.k_mesh = k_mesh
# not used in this version: reset to dummy values?
n_reps = [1 for i in range(n_inequiv_shells)]
dim_reps = [0 for i in range(n_inequiv_shells)]
T = []
for ish in range(n_inequiv_shells):
ll = 2 * corr_shells[inequiv_to_corr[ish]]['l'] + 1
lmax = ll * (corr_shells[inequiv_to_corr[ish]]['SO'] + 1)
T.append(numpy.zeros([lmax, lmax], numpy.complex_))
spin_w90name = ['_up', '_down']
hamr_full = []
# TODO: generalise to SP=1 (only partially done)
rot_mat_time_inv = [0 for i in range(n_corr_shells)]
# Second, let's read the file containing the Hamiltonian in WF basis produced by Wannier90
for isp in range(n_spin):
# begin loop on isp
# build filename according to wannier90 conventions
if SP == 1:
mpi.report(
"Reading information for spin component n. %d" % isp)
hr_file = self.w90_seed + spin_w90name[isp] + '_hr.dat'
else:
hr_file = self.w90_seed + '_hr.dat'
# now grab the data from the H(R) file
mpi.report(
"The Hamiltonian in MLWF basis is extracted from %s ..." % hr_file)
nr, rvec, rdeg, nw, hamr = self.read_wannier90hr(hr_file)
# number of R vectors, their indices, their degeneracy, number of WFs, H(R)
mpi.report("... done: %d R vectors, %d WFs found" % (nr, nw))
if isp == 0:
# set or check some quantities that must be the same for both spins
self.nrpt = nr
# k-point grid: (if not defined before)
if kmesh_mode == -1:
# the size of the k-point mesh is determined from the largest R vector
nki = [2 * rvec[:, idir].max() + 1 for idir in range(3)]
# it will be the same as in the win only when nki is odd, because of the
# wannier90 convention: if we have nki k-points along the i-th direction,
# then we should get 2*(nki/2)+nki%2 R points along that direction
n_k, k_mesh, bz_weights = self.kmesh_build(nki)
self.n_k = n_k
self.k_mesh = k_mesh
# set the R vectors and their degeneracy
self.rvec = rvec
self.rdeg = rdeg
self.nwfs = nw
# check that the total number of WFs makes sense
if self.nwfs < dim_corr_shells:
mpi.report("ERROR: number of WFs in the file smaller than number of correlated orbitals!")
elif self.nwfs > dim_corr_shells:
# NOTE: correlated shells must appear before uncorrelated ones inside the file
mpi.report("Number of WFs larger than correlated orbitals:\n" +
"WFs from %d to %d treated as uncorrelated" % (dim_corr_shells + 1, self.nwfs))
else:
mpi.report("Number of WFs equal to number of correlated orbitals")
# we assume spin up and spin down always have same total number of WFs
n_orbitals = numpy.ones(
[self.n_k, n_spin], numpy.int) * self.nwfs
else:
# consistency check between the _up and _down file contents
if nr != self.nrpt:
mpi.report("Different number of R vectors for spin-up/spin-down!")
if nw != self.nwfs:
mpi.report("Different number of WFs for spin-up/spin-down!")
hamr_full.append(hamr)
# FIXME: when do we actually need deepcopy()?
# hamr_full.append(deepcopy(hamr))
for ir in range(nr):
# checks if the Hamiltonian is real (it should, if wannierisation worked fine)
if numpy.abs((hamr[ir].imag.max()).max()) > self._w90zero:
mpi.report("H(R) has large complex components at R %d" % ir)
# copy the R=0 block corresponding to the correlated shells
# into another variable (needed later for finding rot_mat)
if rvec[ir, 0] == 0 and rvec[ir, 1] == 0 and rvec[ir, 2] == 0:
ham_corr0 = hamr[ir][0:dim_corr_shells, 0:dim_corr_shells]
# checks if ham0 is Hermitian
if not numpy.allclose(ham_corr0.transpose().conjugate(), ham_corr0, atol=self._w90zero, rtol=1.e-9):
raise ValueError("H(R=0) matrix is not Hermitian!")
# find rot_mat symmetries by diagonalising the on-site Hamiltonian of the first spin
if isp == 0:
use_rotations, rot_mat = self.find_rot_mat(n_corr_shells, corr_shells, shells_map, ham_corr0)
else:
# consistency check
use_rotations_, rot_mat_ = self.find_rot_mat(n_corr_shells, corr_shells, shells_map, ham_corr0)
if (use_rotations and not use_rotations_):
mpi.report("Rotations cannot be used for spin component n. %d" % isp)
for icrsh in range(n_corr_shells):
if not numpy.allclose(rot_mat_[icrsh], rot_mat[icrsh], atol=self._w90zero, rtol=1.e-15):
mpi.report("Rotations for spin component n. %d do not match!" % isp)
# end loop on isp
mpi.report("The k-point grid has dimensions: %d, %d, %d" % tuple(nki))
# if calculations are spin-polarized, then renormalize k-point weights
if SP == 1:
bz_weights = 0.5 * bz_weights
# Third, compute the hoppings in reciprocal space
hopping = numpy.zeros([self.n_k, n_spin, numpy.max(n_orbitals), numpy.max(n_orbitals)], numpy.complex_)
for isp in range(n_spin):
# make Fourier transform H(R) -> H(k) : it can be done one spin at a time
hamk = self.fourier_ham(self.nwfs, hamr_full[isp])
# copy the H(k) in the right place of hoppings... is there a better way to do this??
for ik in range(self.n_k):
#hopping[ik,isp,:,:] = deepcopy(hamk[ik][:,:])*energy_unit
hopping[ik, isp, :, :] = hamk[ik][:, :] * energy_unit
# Then, initialise the projectors
k_dep_projection = 0 # we always have the same number of WFs at each k-point
proj_mat = numpy.zeros([self.n_k, n_spin, n_corr_shells, max(
[crsh['dim'] for crsh in corr_shells]), numpy.max(n_orbitals)], numpy.complex_)
iorb = 0
# Projectors simply consist in identity matrix blocks selecting those MLWFs that
# correspond to the specific correlated shell indexed by icrsh.
# NOTE: we assume that the correlated orbitals appear at the beginning of the H(R)
# file and that the ordering of MLWFs matches the corr_shell info from the input.
for icrsh in range(n_corr_shells):
norb = corr_shells[icrsh]['dim']
proj_mat[:, :, icrsh, 0:norb, iorb:iorb +
norb] = numpy.identity(norb, numpy.complex_)
iorb += norb
# Finally, save all required data into the HDF archive:
ar = HDFArchive(self.hdf_file, 'a')
if not (self.dft_subgrp in ar):
ar.create_group(self.dft_subgrp)
# The subgroup containing the data. If it does not exist, it is created. If it exists, the data is overwritten!
things_to_save = ['energy_unit', 'n_k', 'k_dep_projection', 'SP', 'SO', 'charge_below', 'density_required',
'symm_op', 'n_shells', 'shells', 'n_corr_shells', 'corr_shells', 'use_rotations', 'rot_mat',
'rot_mat_time_inv', 'n_reps', 'dim_reps', 'T', 'n_orbitals', 'proj_mat', 'bz_weights', 'hopping',
'n_inequiv_shells', 'corr_to_inequiv', 'inequiv_to_corr']
for it in things_to_save:
ar[self.dft_subgrp][it] = locals()[it]
del ar
def read_wannier90hr(self, hr_filename="wannier_hr.dat"):
"""
Method for reading the seedname_hr.dat file produced by Wannier90 (http://wannier.org)
Parameters
----------
hr_filename : string
full name of the H(R) file produced by Wannier90 (usually seedname_hr.dat)
Returns
-------
nrpt : integer
number of R vectors found in the file
rvec_idx : numpy.array of integers
Miller indices of the R vectors
rvec_deg : numpy.array of floats
weight of the R vectors
num_wf : integer
number of Wannier functions found
h_of_r : list of numpy.array
<w_i|H(R)|w_j> = Hamilonian matrix elements in the Wannier basis
"""
# Read only from the master node
if not (mpi.is_master_node()):
return
try:
with open(hr_filename, "r") as hr_filedesc:
hr_data = hr_filedesc.readlines()
hr_filedesc.close()
except IOError:
mpi.report("The file %s could not be read!" % hr_filename)
mpi.report("Reading %s..." % hr_filename + hr_data[0])
try:
# reads number of Wannier functions per spin
num_wf = int(hr_data[1])
nrpt = int(hr_data[2])
except ValueError:
mpi.report("Could not read number of WFs or R vectors")
# allocate arrays to save the R vector indexes and degeneracies and the Hamiltonian
rvec_idx = numpy.zeros((nrpt, 3), dtype=int)
rvec_deg = numpy.zeros(nrpt, dtype=int)
h_of_r = [numpy.zeros((num_wf, num_wf), dtype=numpy.complex_)
for n in range(nrpt)]
# variable currpos points to the current line in the file
currpos = 2
try:
ir = 0
# read the degeneracy of the R vectors (needed for the Fourier transform)
while ir < nrpt:
currpos += 1
for x in hr_data[currpos].split():
if ir >= nrpt:
raise IndexError("wrong number of R vectors??")
rvec_deg[ir] = int(x)
ir += 1
# for each direct lattice vector R read the block of the
# Hamiltonian H(R)
for ir, jj, ii in product(range(nrpt), range(num_wf), range(num_wf)):
# advance one line, split the line into tokens
currpos += 1
cline = hr_data[currpos].split()
# check if the orbital indexes in the file make sense
if int(cline[3]) != ii + 1 or int(cline[4]) != jj + 1:
mpi.report(
"Inconsistent indices at %s%s of R n. %s" % (ii, jj, ir))
rcurr = numpy.array(
[int(cline[0]), int(cline[1]), int(cline[2])])
if ii == 0 and jj == 0:
rvec_idx[ir] = rcurr
rprec = rcurr
else:
# check if the vector indices are consistent
if not numpy.array_equal(rcurr, rprec):
mpi.report(
"Inconsistent indices for R vector n. %s" % ir)
# fill h_of_r with the matrix elements of the Hamiltonian
h_of_r[ir][ii, jj] = complex(float(cline[5]), float(cline[6]))
except ValueError:
mpi.report("Wrong data or structure in file %s" % hr_filename)
# return the data into variables
return nrpt, rvec_idx, rvec_deg, num_wf, h_of_r
def find_rot_mat(self, n_sh, sh_lst, sh_map, ham0):
"""
Method for finding the matrices that bring from local to global coordinate systems
(and viceversa), based on the eigenvalues of H(R=0)
Parameters
----------
n_sh : integer
number of shells
sh_lst : list of shells-type dictionaries
contains the shells (could be correlated or not)
sh_map : list of integers
mapping between shells
ham0 : numpy.array of floats
local Hamiltonian matrix elements
Returns
-------
istatus : integer
if 0, something failed in the construction of the matrices
rot_mat : list of numpy.array
rotation matrix for each of the shell
"""
# initialize the rotation matrices to identities
rot_mat = [numpy.identity(sh_lst[ish]['dim'], dtype=complex)
for ish in range(n_sh)]
istatus = 0
hs = ham0.shape
if hs[0] != hs[1] or hs[0] != sum([sh['dim'] for sh in sh_lst]):
mpi.report(
"find_rot_mat: wrong block structure of input Hamiltonian!")
istatus = 0
# this error will lead into troubles later... early return
return istatus, rot_mat
# TODO: better handling of degenerate eigenvalue case
eigval_lst = []
eigvec_lst = []
iwf = 0
# loop over shells
for ish in range(n_sh):
# nw = number of orbitals in this shell
nw = sh_lst[ish]["dim"]
# diagonalize the sub-block of H(0) corresponding to this shell
eigval, eigvec = numpy.linalg.eigh(
ham0[iwf:iwf + nw, iwf:iwf + nw])
# find the indices sorting the eigenvalues in ascending order
eigsrt = eigval[0:nw].argsort()
# order eigenvalues and eigenvectors and save in a list
eigval_lst.append(eigval[eigsrt])
eigvec_lst.append(eigvec[eigsrt])
iwf += nw
# TODO: better handling of degenerate eigenvalue case
if sh_map[ish] != ish: # issue warning only when there are equivalent shells
for i in range(nw):
for j in range(i + 1, nw):
if (abs(eigval[j] - eigval[i]) < self._w90zero):
mpi.report("WARNING: degenerate eigenvalue of H(0) detected for shell %d: " % (ish) +
"global-to-local transformation might not work!")
for ish in range(n_sh):
try:
# build rotation matrices by combining the unitary
# transformations that diagonalize H(0)
rot_mat[ish] = numpy.dot(eigvec_lst[ish], eigvec_lst[
sh_map[ish]].conjugate().transpose())
except ValueError:
mpi.report(
"Global-to-local rotation matrices cannot be constructed!")
istatus = 1
# check that eigenvalues are the same (within accuracy) for
# equivalent shells
if not numpy.allclose(eigval_lst[ish], eigval_lst[sh_map[ish]], atol=self._w90zero, rtol=1.e-15):
mpi.report(
"ERROR: eigenvalue mismatch between equivalent shells! %d" % ish)
eigval_diff = eigval_lst[ish] - eigval_lst[sh_map[ish]]
mpi.report("Eigenvalue difference: " + format(eigval_diff))
istatus = 0
# TODO: add additional consistency check on rot_mat matrices?
return istatus, rot_mat
def kmesh_build(self, msize=None, mmode=0):
"""
Method for the generation of the k-point mesh.
Right now it only supports the option for generating a full grid containing k=0,0,0.
Parameters
----------
msize : list of 3 integers
the dimensions of the mesh
mmode : integer
mesh generation mode (right now, only full grid available)
Returns
-------
nkpt : integer
total number of k-points in the mesh
k_mesh : numpy.array[nkpt,3] of floats
the coordinates of all k-points
wk : numpy.array[nkpt] of floats
the weight of each k-point
"""
if mmode != 0:
raise ValueError("Mesh generation mode not supported: %s" % mmode)
# a regular mesh including Gamma point
# total number of k-points
nkpt = msize[0] * msize[1] * msize[2]
kmesh = numpy.zeros((nkpt, 3), dtype=float)
ii = 0
for ix, iy, iz in product(range(msize[0]), range(msize[1]), range(msize[2])):
kmesh[ii, :] = [float(ix) / msize[0], float(iy) / msize[1], float(iz) / msize[2]]
ii += 1
# weight is equal for all k-points because wannier90 uses uniform grid on whole BZ
# (normalization is always 1 and takes into account spin degeneracy)
wk = numpy.ones([nkpt], dtype=float) / float(nkpt)
return nkpt, kmesh, wk
def fourier_ham(self, norb, h_of_r):
"""
Method for obtaining H(k) from H(R) via Fourier transform
The R vectors and k-point mesh are read from global module variables
Parameters
----------
norb : integer
number of orbitals
h_of_r : list of numpy.array[norb,norb]
Hamiltonian H(R) in Wannier basis
Returns
-------
h_of_k : list of numpy.array[norb,norb]
transformed Hamiltonian H(k) in Wannier basis
"""
twopi = 2 * numpy.pi
h_of_k = [numpy.zeros((norb, norb), dtype=numpy.complex_) for ik in range(self.n_k)]
ridx = numpy.array(range(self.nrpt))
for ik, ir in product(range(self.n_k), ridx):
rdotk = twopi * numpy.dot(self.k_mesh[ik], self.rvec[ir])
factor = (math.cos(rdotk) + 1j * math.sin(rdotk)) / float(self.rdeg[ir])
h_of_k[ik][:, :] += factor * h_of_r[ir][:, :]
return h_of_k

View File

@ -448,7 +448,6 @@ class SumkDFT:
broadening = 0.01
else: # broadening = 2 * \Delta omega, where \Delta omega is the spacing of omega points
broadening = 2.0 * ( (mesh[1]-mesh[0])/(mesh[2]-1) )
n_iw = 1025 # Default number of Matsubara frequencies
# Are we including Sigma?
if with_Sigma:
@ -457,10 +456,16 @@ class SumkDFT:
if with_dc: sigma_minus_dc = self.add_dc(iw_or_w)
if iw_or_w == "iw":
beta = Sigma_imp[0].mesh.beta # override beta if Sigma_iw is present
n_iw = len(Sigma_imp[0].mesh)
mesh = Sigma_imp[0].mesh
elif iw_or_w == "w":
mesh = Sigma_imp[0].mesh
else:
if (iw_or_w == "w") and (mesh is None):
raise ValueError, "lattice_gf: Give the mesh=(om_min,om_max,n_points) for the lattice GfReFreq."
if iw_or_w == "iw":
if beta is None: raise ValueError, "lattice_gf: Give the beta for the lattice GfReFreq."
mesh = MeshImFreq(beta=beta, S='Fermion', n_max=1025) # Default number of Matsubara frequencies
elif iw_or_w == "w":
if mesh is None: raise ValueError, "lattice_gf: Give the mesh=(om_min,om_max,n_points) for the lattice GfReFreq."
mesh = MeshReFreq(mesh[0],mesh[1],mesh[2])
# Check if G_latt is present
set_up_G_latt = False # Assume not
@ -479,12 +484,9 @@ class SumkDFT:
gf_struct = [ (spn[isp], block_structure[isp]) for isp in range(self.n_spin_blocks[self.SO]) ]
block_ind_list = [block for block,inner in gf_struct]
if iw_or_w == "iw":
glist = lambda : [ GfImFreq(indices=inner,beta=beta,n_points=n_iw) for block,inner in gf_struct]
glist = lambda : [ GfImFreq(indices=inner,mesh=mesh) for block,inner in gf_struct ]
elif iw_or_w == "w":
if with_Sigma:
glist = lambda : [ GfReFreq(indices=inner,mesh=Sigma_imp[0].mesh) for block,inner in gf_struct]
else:
glist = lambda : [ GfReFreq(indices=inner,window=(mesh[0],mesh[1]),n_points=mesh[2]) for block,inner in gf_struct]
glist = lambda : [ GfReFreq(indices=inner,mesh=mesh) for block,inner in gf_struct ]
G_latt = BlockGf(name_list = block_ind_list, block_list = glist(), make_copies = False)
G_latt.zero()
@ -510,6 +512,8 @@ class SumkDFT:
return G_latt
def set_Sigma(self,Sigma_imp):
self.put_Sigma(Sigma_imp)
def put_Sigma(self, Sigma_imp):
r"""
@ -557,7 +561,6 @@ class SumkDFT:
for icrsh in range(self.n_corr_shells):
for bname,gf in SK_Sigma_imp[icrsh]: gf << self.rotloc(icrsh,gf,direction='toGlobal')
def extract_G_loc(self, mu=None, with_Sigma=True, with_dc=True):
r"""
Extracts the local downfolded Green function by the Brillouin-zone integration of the lattice Green's function.

View File

@ -47,6 +47,7 @@ class SumkDFTTools(SumkDFT):
misc_data=misc_data)
# Uses .data of only GfReFreq objects.
def dos_wannier_basis(self, mu=None, broadening=None, mesh=None, with_Sigma=True, with_dc=True, save_to_file=True):
"""
Calculates the density of states in the basis of the Wannier functions.
@ -163,6 +164,7 @@ class SumkDFTTools(SumkDFT):
return DOS, DOSproj, DOSproj_orb
# Uses .data of only GfReFreq objects.
def dos_parproj_basis(self, mu=None, broadening=None, mesh=None, with_Sigma=True, with_dc=True, save_to_file=True):
"""
Calculates the orbitally-resolved DOS.
@ -290,6 +292,7 @@ class SumkDFTTools(SumkDFT):
return DOS, DOSproj, DOSproj_orb
# Uses .data of only GfReFreq objects.
def spaghettis(self,broadening=None,plot_shift=0.0,plot_range=None,ishell=None,mu=None,save_to_file='Akw_'):
"""
Calculates the correlated band structure using a real-frequency self energy.
@ -377,6 +380,11 @@ class SumkDFTTools(SumkDFT):
for sp in spn:
Akw[sp][ish,ik,iom] = G_loc[sp].data[iom,ish,ish].imag/(-1.0*numpy.pi)
# Collect data from mpi
for sp in spn:
Akw[sp] = mpi.all_reduce(mpi.world, Akw[sp], lambda x,y : x+y)
mpi.barrier()
if save_to_file and mpi.is_master_node():
if ishell is None:
for sp in spn: # loop over GF blocs:
@ -394,7 +402,7 @@ class SumkDFTTools(SumkDFT):
else: # ishell is not None
for sp in spn:
for ish in range(self.shells[ishell]['dim']):
f = open(save_to_file+sp+'_proj'+str(ish)+'.dat','w') # Open file for storage:
f = open(save_to_file+str(ishell)+'_'+sp+'_proj'+str(ish)+'.dat','w') # Open file for storage:
for ik in range(self.n_k):
for iom in range(n_om):
if (mesh[iom] > om_minplot) and (mesh[iom] < om_maxplot):
@ -565,6 +573,7 @@ class SumkDFTTools(SumkDFT):
return vol_c, vol_p
# Uses .data of only GfReFreq objects.
def transport_distribution(self, beta, directions=['xx'], energy_window=None, Om_mesh=[0.0], with_Sigma=False, n_om=None, broadening=0.0):
r"""
Calculates the transport distribution

30
python/version.py.in Normal file
View File

@ -0,0 +1,30 @@
################################################################################
#
# TRIQS: a Toolbox for Research in Interacting Quantum Systems
#
# Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola
#
# TRIQS is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.
#
# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# TRIQS. If not, see <http://www.gnu.org/licenses/>.
#
################################################################################
version = "@DFT_TOOLS_VERSION@"
triqs_hash = "@TRIQS_GIT_HASH@"
cthyb_hash = "@CTHYB_GIT_HASH@"
def show_version():
print "\nYou are using the dft_tools version %s\n"%version
def show_git_hash():
print "\nYou are using the dft_tools git hash %s based on triqs git hash %s\n"%(cthyb_hash, triqs_hash)

View File

@ -1,12 +1,18 @@
# load triqs helper to set up tests
find_package(TriqsTest)
FILE(COPY SrVO3.h5 SrVO3_Sigma.h5 SrVO3.pmat SrVO3.struct SrVO3.outputs SrVO3.oubwin SrVO3.ctqmcout SrVO3.symqmc SrVO3.sympar SrVO3.parproj hk_convert_hamiltonian.hk DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
#triqs_add_test_hdf(wien2k_convert " -p 1.e-6" )
#triqs_add_test_hdf(hk_convert " -p 1.e-6" )
#triqs_add_test_hdf(sumkdft_basic " -d 1.e-6" )
#triqs_add_test_hdf(srvo3_Gloc " -d 1.e-6" )
#triqs_add_test_hdf(srvo3_transp " -d 1.e-6" )
#triqs_add_test_hdf(sigma_from_file " -d 1.e-6" )
# Copy h5 files to binary dir
FILE(GLOB all_h5_files RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h5)
file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/${all_h5_files} DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
# Copy other files
FILE(COPY SrVO3.pmat SrVO3.struct SrVO3.outputs SrVO3.oubwin SrVO3.ctqmcout SrVO3.symqmc SrVO3.sympar SrVO3.parproj hk_convert_hamiltonian.hk LaVO3-Pnma_hr.dat LaVO3-Pnma.inp DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
triqs_add_python_test(wien2k_convert)
triqs_add_python_test(hk_convert)
triqs_add_python_test(w90_convert)
triqs_add_python_test(sumkdft_basic)
triqs_add_python_test(srvo3_Gloc)
triqs_add_python_test(srvo3_transp)
triqs_add_python_test(sigma_from_file)
add_subdirectory(plovasp)

7
test/LaVO3-Pnma.inp Normal file
View File

@ -0,0 +1,7 @@
0 3 2 3
8.0
4
0 0 2 3 0 0
1 0 2 3 0 0
2 0 2 3 0 0
3 0 2 3 0 0

3893
test/LaVO3-Pnma_hr.dat Normal file

File diff suppressed because it is too large Load Diff

View File

@ -23,7 +23,12 @@
from pytriqs.applications.dft.converters import *
from pytriqs.archive import *
from pytriqs.utility.h5diff import h5diff
import pytriqs.utility.mpi as mpi
Converter = HkConverter(filename='hk_convert_hamiltonian.hk',hdf_filename='hk_convert.output.h5')
Converter = HkConverter(filename='hk_convert_hamiltonian.hk',hdf_filename='hk_convert.out.h5')
Converter.convert_dft_input()
if mpi.is_master_node():
h5diff("hk_convert.out.h5","hk_convert.ref.h5")

View File

@ -1,7 +1,29 @@
################################################################################
#
# TRIQS: a Toolbox for Research in Interacting Quantum Systems
#
# Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola
#
# TRIQS is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.
#
# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# TRIQS. If not, see <http://www.gnu.org/licenses/>.
#
################################################################################
from pytriqs.archive import *
from pytriqs.gf.local import *
from pytriqs.gf.local.tools import *
from pytriqs.applications.dft.sumk_dft_tools import *
from pytriqs.utility.comparison_tests import *
import numpy as np
# Read self energy from hdf file
@ -24,9 +46,9 @@ a_list = [a for a,al in SK.gf_struct_solver[0].iteritems()]
g_list = [read_gf_from_txt([['Sigma_' + a + '.dat']], a) for a in a_list]
Sigma_txt = BlockGf(name_list = a_list, block_list = g_list, make_copies=False)
SK.put_Sigma(Sigma_imp = [Sigma_txt])
SK.set_Sigma([Sigma_txt])
SK.hdf_file = 'sigma_from_file.output.h5'
SK.hdf_file = 'sigma_from_file.out.h5'
SK.save(['Sigma_imp_w'])
if ((Sigma_txt - Sigma_hdf).real < 1e-6) & ((Sigma_txt - Sigma_hdf).imag < 1e-6):

View File

@ -24,6 +24,8 @@ from pytriqs.gf.local import *
from pytriqs.applications.dft.sumk_dft import *
from pytriqs.applications.dft.converters.wien2k_converter import *
from pytriqs.operators.util import set_operator_structure
from pytriqs.utility.comparison_tests import *
from pytriqs.utility.h5diff import h5diff
# Basic input parameters
beta = 40
@ -41,9 +43,12 @@ gf_struct = set_operator_structure(spin_names,orb_names,orb_hybridized)
glist = [ GfImFreq(indices=inner,beta=beta) for block,inner in gf_struct.iteritems()]
Sigma_iw = BlockGf(name_list = gf_struct.keys(), block_list = glist, make_copies = False)
SK.put_Sigma([Sigma_iw])
Gloc=SK.extract_G_loc()
SK.set_Sigma([Sigma_iw])
Gloc = SK.extract_G_loc()
ar = HDFArchive('srvo3_Gloc.output.h5','w')
ar['Gloc'] = Gloc[0]
del ar
if mpi.is_master_node():
with HDFArchive('srvo3_Gloc.out.h5','w') as ar:
ar['Gloc'] = Gloc[0]
if mpi.is_master_node():
h5diff("srvo3_Gloc.out.h5","srvo3_Gloc.ref.h5")

View File

@ -23,6 +23,8 @@ from numpy import *
from pytriqs.applications.dft.converters.wien2k_converter import *
from pytriqs.applications.dft.sumk_dft import *
from pytriqs.applications.dft.sumk_dft_tools import *
from pytriqs.utility.comparison_tests import *
from pytriqs.utility.h5diff import h5diff
beta = 40
@ -34,7 +36,7 @@ SK = SumkDFTTools(hdf_file='SrVO3.h5', use_dft_blocks=True)
ar = HDFArchive('SrVO3_Sigma.h5', 'a')
Sigma = ar['dmft_transp_input']['Sigma_w']
SK.put_Sigma(Sigma_imp = [Sigma])
SK.set_Sigma([Sigma])
SK.chemical_potential = ar['dmft_transp_input']['chemical_potential']
SK.dc_imp = ar['dmft_transp_input']['dc_imp']
del ar
@ -43,6 +45,8 @@ SK.transport_distribution(directions=['xx'], broadening=0.0, energy_window=[-0.3
#SK.save(['Gamma_w','Om_meshr','omega','directions'])
#SK.load(['Gamma_w','Om_meshr','omega','directions'])
SK.conductivity_and_seebeck(beta=beta)
SK.hdf_file = 'srvo3_transp.output.h5'
SK.hdf_file = 'srvo3_transp.out.h5'
SK.save(['seebeck','optic_cond'])
if mpi.is_master_node():
h5diff("srvo3_transp.out.h5","srvo3_transp.ref.h5")

View File

@ -22,14 +22,18 @@
from pytriqs.archive import *
from pytriqs.applications.dft.sumk_dft_tools import SumkDFTTools
import pytriqs.utility.mpi as mpi
from pytriqs.utility.comparison_tests import *
from pytriqs.utility.h5diff import h5diff
SK = SumkDFTTools(hdf_file = 'SrVO3.h5')
dm = SK.density_matrix(method = 'using_gf', beta = 40)
dm_pc = SK.partial_charges(beta=40,with_Sigma=False,with_dc=False)
ar = HDFArchive('sumkdft_basic.output.h5','w')
ar['dm'] = dm
ar['dm_pc'] = dm_pc
del ar
with HDFArchive('sumkdft_basic.out.h5','w') as ar:
ar['dm'] = dm
ar['dm_pc'] = dm_pc
if mpi.is_master_node():
h5diff('sumkdft_basic.out.h5','sumkdft_basic.ref.h5')

34
test/w90_convert.py Normal file
View File

@ -0,0 +1,34 @@
################################################################################
#
# TRIQS: a Toolbox for Research in Interacting Quantum Systems
#
# Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola
#
# TRIQS is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.
#
# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# TRIQS. If not, see <http://www.gnu.org/licenses/>.
#
################################################################################
from pytriqs.applications.dft.converters import *
from pytriqs.archive import *
from pytriqs.utility.h5diff import h5diff
import pytriqs.utility.mpi as mpi
Converter = Wannier90Converter(seedname='LaVO3-Pnma',hdf_filename='w90_convert.out.h5')
Converter.convert_dft_input()
if mpi.is_master_node():
h5diff("w90_convert.out.h5","w90_convert.ref.h5")

BIN
test/w90_convert.ref.h5 Normal file

Binary file not shown.

View File

@ -22,13 +22,15 @@
from pytriqs.archive import *
from pytriqs.applications.dft.converters import Wien2kConverter
from pytriqs.utility.comparison_tests import *
from pytriqs.utility.h5diff import h5diff
import pytriqs.utility.mpi as mpi
Converter = Wien2kConverter(filename='SrVO3')
Converter.hdf_file = 'wien2k_convert.output.h5'
Converter.hdf_file = 'wien2k_convert.out.h5'
Converter.convert_dft_input()
Converter.convert_parproj_input()
if mpi.is_master_node():
h5diff('wien2k_convert.out.h5','wien2k_convert.ref.h5')