diff --git a/dft_dmft_cthyb.py b/dft_dmft_cthyb.py deleted file mode 100644 index e5a26b32..00000000 --- a/dft_dmft_cthyb.py +++ /dev/null @@ -1,147 +0,0 @@ -import pytriqs.utility.mpi as mpi -from pytriqs.operators.util import * -from pytriqs.archive import HDFArchive -from triqs_cthyb import * -from pytriqs.gf import * -from triqs_dft_tools.sumk_dft import * -from triqs_dft_tools.converters.wien2k_converter import * - -dft_filename='Gd_fcc' -U = 9.6 -J = 0.8 -beta = 40 -loops = 10 # Number of DMFT sc-loops -sigma_mix = 1.0 # Mixing factor of Sigma after solution of the AIM -delta_mix = 1.0 # Mixing factor of Delta as input for the AIM -dc_type = 0 # DC type: 0 FLL, 1 Held, 2 AMF -use_blocks = True # use bloc structure from DFT input -prec_mu = 0.0001 -h_field = 0.0 - -# Solver parameters -p = {} -p["max_time"] = -1 -p["length_cycle"] = 50 -p["n_warmup_cycles"] = 50 -p["n_cycles"] = 5000 - -Converter = Wien2kConverter(filename=dft_filename, repacking=True) -Converter.convert_dft_input() -mpi.barrier() - -previous_runs = 0 -previous_present = False -if mpi.is_master_node(): - f = HDFArchive(dft_filename+'.h5','a') - if 'dmft_output' in f: - ar = f['dmft_output'] - if 'iterations' in ar: - previous_present = True - previous_runs = ar['iterations'] - else: - f.create_group('dmft_output') - del f -previous_runs = mpi.bcast(previous_runs) -previous_present = mpi.bcast(previous_present) - -SK=SumkDFT(hdf_file=dft_filename+'.h5',use_dft_blocks=use_blocks,h_field=h_field) - -n_orb = SK.corr_shells[0]['dim'] -l = SK.corr_shells[0]['l'] -spin_names = ["up","down"] -orb_names = [i for i in range(n_orb)] - -# Use GF structure determined by DFT blocks -gf_struct = [(block, indices) for block, indices in SK.gf_struct_solver[0].iteritems()] -# Construct U matrix for density-density calculations -Umat, Upmat = U_matrix_kanamori(n_orb=n_orb, U_int=U, J_hund=J) -# Construct Hamiltonian and solver -h_int = h_int_density(spin_names, orb_names, map_operator_structure=SK.sumk_to_solver[0], U=Umat, Uprime=Upmat, H_dump="H.txt") -S = Solver(beta=beta, gf_struct=gf_struct) - -if previous_present: - chemical_potential = 0 - dc_imp = 0 - dc_energ = 0 - if mpi.is_master_node(): - S.Sigma_iw << HDFArchive(dft_filename+'.h5','a')['dmft_output']['Sigma_iw'] - chemical_potential,dc_imp,dc_energ = SK.load(['chemical_potential','dc_imp','dc_energ']) - S.Sigma_iw << mpi.bcast(S.Sigma_iw) - chemical_potential = mpi.bcast(chemical_potential) - dc_imp = mpi.bcast(dc_imp) - dc_energ = mpi.bcast(dc_energ) - SK.set_mu(chemical_potential) - SK.set_dc(dc_imp,dc_energ) - -for iteration_number in range(1,loops+1): - if mpi.is_master_node(): print "Iteration = ", iteration_number - - SK.symm_deg_gf(S.Sigma_iw,orb=0) # symmetrise Sigma - SK.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()) - - # Init the DC term and the real part of Sigma, if no previous runs found: - if (iteration_number==1 and previous_present==False): - dm = S.G_iw.density() - SK.calc_dc(dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type) - S.Sigma_iw << SK.dc_imp[0]['up'][0,0] - - # Calculate new G0_iw to input into the solver: - if mpi.is_master_node(): - # We can do a mixing of Delta in order to stabilize the DMFT iterations: - S.G0_iw << S.Sigma_iw + inverse(S.G_iw) - ar = HDFArchive(dft_filename+'.h5','a')['dmft_output'] - if (iteration_number>1 or previous_present): - mpi.report("Mixing input Delta with factor %s"%delta_mix) - Delta = (delta_mix * delta(S.G0_iw)) + (1.0-delta_mix) * ar['Delta_iw'] - S.G0_iw << S.G0_iw + delta(S.G0_iw) - Delta - ar['Delta_iw'] = delta(S.G0_iw) - S.G0_iw << inverse(S.G0_iw) - del ar - - S.G0_iw << mpi.bcast(S.G0_iw) - - # Solve the impurity problem: - S.solve(h_int=h_int, **p) - - # Solved. Now do post-processing: - mpi.report("Total charge of impurity problem : %.6f"%S.G_iw.total_density()) - - # Now mix Sigma and G with factor sigma_mix, if wanted: - if (iteration_number>1 or previous_present): - if mpi.is_master_node(): - ar = HDFArchive(dft_filename+'.h5','a')['dmft_output'] - mpi.report("Mixing Sigma and G with factor %s"%sigma_mix) - S.Sigma_iw << sigma_mix * S.Sigma_iw + (1.0-sigma_mix) * ar['Sigma_iw'] - S.G_iw << sigma_mix * S.G_iw + (1.0-sigma_mix) * ar['G_iw'] - del ar - S.G_iw << mpi.bcast(S.G_iw) - S.Sigma_iw << mpi.bcast(S.Sigma_iw) - - # Write the final Sigma and G to the hdf5 archive: - if mpi.is_master_node(): - ar = HDFArchive(dft_filename+'.h5','a')['dmft_output'] - if previous_runs: iteration_number += previous_runs - ar['iterations'] = iteration_number - ar['G_tau'] = S.G_tau - ar['G_iw'] = S.G_iw - ar['Sigma_iw'] = S.Sigma_iw - ar['G0-%s'%(iteration_number)] = S.G0_iw - ar['G-%s'%(iteration_number)] = S.G_iw - ar['Sigma-%s'%(iteration_number)] = S.Sigma_iw - del ar - - # Set the new double counting: - dm = S.G_iw.density() # compute the density matrix of the impurity problem - SK.calc_dc(dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type) - - # Save stuff into the dft_output group of hdf5 archive in case of rerun: - SK.save(['chemical_potential','dc_imp','dc_energ']) - -if mpi.is_master_node(): - ar = HDFArchive("dftdmft.h5",'w') - ar["G_tau"] = S.G_tau - ar["G_iw"] = S.G_iw - ar["Sigma_iw"] = S.Sigma_iw diff --git a/doc/conf.py.in b/doc/conf.py.in index 32f60a1e..fb18a948 100644 --- a/doc/conf.py.in +++ b/doc/conf.py.in @@ -19,7 +19,7 @@ extensions = ['sphinx.ext.autodoc', source_suffix = '.rst' project = u'TRIQS DFTTools' -copyright = u'2011-2013, M. Aichhorn, L. Pourovskii, V. Vildosola, C. Martins' +copyright = u'2011-2019' version = '@DFT_TOOLS_VERSION@' mathjax_path = "@TRIQS_MATHJAX_PATH@/MathJax.js?config=default" @@ -32,6 +32,7 @@ html_context = {'header_title': 'dft tools', 'header_subtitle': 'connecting TRIQS to DFT packages', 'header_links': [['Install', 'install'], ['Documentation', 'documentation'], + ['Tutorials', 'tutorials'], ['Issues', 'issues'], ['About DFTTools', 'about']]} html_static_path = ['@CMAKE_SOURCE_DIR@/doc/_static'] diff --git a/doc/documentation.rst b/doc/documentation.rst index b000a6c8..7aa3fa79 100644 --- a/doc/documentation.rst +++ b/doc/documentation.rst @@ -16,18 +16,31 @@ Basic notions basicnotions/structure -User guide ----------- +Construction of local orbitals from DFT +--------------------------------------- .. toctree:: :maxdepth: 2 guide/conversion + + +DFT+DMFT +-------- + +.. toctree:: + :maxdepth: 2 + guide/dftdmft_singleshot - guide/SrVO3 guide/dftdmft_selfcons + +Postprocessing +-------------- + +.. toctree:: + :maxdepth: 2 + guide/analysis - guide/full_tutorial guide/transport diff --git a/doc/guide/analysis.rst b/doc/guide/analysis.rst index 004d22aa..90aca2bd 100644 --- a/doc/guide/analysis.rst +++ b/doc/guide/analysis.rst @@ -40,8 +40,8 @@ If required, we have to load and initialise the real-frequency self energy. Most you have your self energy already stored as a real-frequency :class:`BlockGf ` object in a hdf5 file:: - ar = HDFArchive('case.h5', 'a') - SigmaReFreq = ar['dmft_output']['Sigma_w'] + with HDFArchive('case.h5', 'r') as ar: + SigmaReFreq = ar['dmft_output']['Sigma_w'] You may also have your self energy stored in text files. For this case the :ref:`TRIQS ` library offers the function :meth:`read_gf_from_txt`, which is able to load the data from text files of one Green function block @@ -73,7 +73,6 @@ and additionally set the chemical potential and the double counting correction f 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) - del ar .. _dos_wannier: diff --git a/doc/guide/conv_W90.rst b/doc/guide/conv_W90.rst new file mode 100644 index 00000000..fdc9a5f1 --- /dev/null +++ b/doc/guide/conv_W90.rst @@ -0,0 +1,117 @@ +.. _convW90: + +Wannier90 Converter +=================== + +Using this converter it is possible to convert the output of +`wannier90 `_ +Maximally Localized Wannier Functions (MLWF) and create a HDF5 archive +suitable for one-shot DMFT calculations with the +:class:`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 triqs_dft_tools.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 (do not use the text/comments in your input file): + +.. 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. In this example we have four lines for the +four V atoms. The MLWFs were constructed for the t\ :sub:`2g` subspace, and thus +we set ``l`` to 2 and ``dim`` to 3 for all V atoms. Further the spin-orbit coupling (``SO``) +is set to 0 and ``irep`` to 0. +As in this example all 4 V atoms are equivalent we set ``sort`` to 0. We note +that, e.g., for a magnetic DMFT calculation the correlated atoms can be made +inequivalent at this point by using different values for ``sort``. + +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 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 ` + were not tested with this converter. +* ``proj_mat_all`` are not used, so there are no projectors onto the + uncorrelated orbitals for now. + + diff --git a/doc/guide/conv_generalhk.rst b/doc/guide/conv_generalhk.rst new file mode 100644 index 00000000..de6da36b --- /dev/null +++ b/doc/guide/conv_generalhk.rst @@ -0,0 +1,100 @@ +.. _convgeneralhk: + +A general H(k) +============== + +In addition to the more extensive Wien2k, VASP, and W90 converters, +:program:`DFTTools` contains also a light converter. It takes only +one inputfile, and creates the necessary hdf outputfile for +the DMFT calculation. The header of this input file has a defined +format, an example is the following (do not use the text/comments in your +input file): + +.. literalinclude:: images_scripts/case.hk + +The lines of this header define + +#. Number of :math:`\mathbf{k}`-points used in the calculation +#. Electron density for setting the chemical potential +#. Number of total atomic shells in the hamiltonian matrix. In short, + this gives the number of lines described in the following. IN the + example file give above this number is 2. +#. The next line(s) contain four numbers each: index of the atom, index + of the equivalent shell, :math:`l` quantum number, dimension + of this shell. Repeat this line for each atomic shell, the number + of the shells is given in the previous line. + + In the example input file given above, we have two inequivalent + atomic shells, one on atom number 1 with a full d-shell (dimension 5), + and one on atom number 2 with one p-shell (dimension 3). + + Other examples for these lines are: + + #. Full d-shell in a material with only one correlated atom in the + unit cell (e.g. SrVO3). One line is sufficient and the numbers + are `1 1 2 5`. + #. Full d-shell in a material with two equivalent atoms in the unit + cell (e.g. FeSe): You need two lines, one for each equivalent + atom. First line is `1 1 2 5`, and the second line is + `2 1 2 5`. The only difference is the first number, which tells on + which atom the shell is located. The second number is the + same in both lines, meaning that both atoms are equivalent. + #. t2g orbitals on two non-equivalent atoms in the unit cell: Two + lines again. First line is `1 1 2 3`, second line `2 2 2 3`. The + difference to the case above is that now also the second number + differs. Therefore, the two shells are treated independently in + the calculation. + #. d-p Hamiltonian in a system with two equivalent atoms each in + the unit cell (e.g. FeSe has two Fe and two Se in the unit + cell). You need for lines. First line `1 1 2 5`, second + line + `2 1 2 5`. These two lines specify Fe as in the case above. For the p + orbitals you need line three as `3 2 1 3` and line four + as `4 2 1 3`. We have 4 atoms, since the first number runs from 1 to 4, + but only two inequivalent atoms, since the second number runs + only form 1 to 2. + + Note that the total dimension of the hamiltonian matrices that are + read in is the sum of all shell dimensions that you specified. For + example number 4 given above we have a dimension of 5+5+3+3=16. It is important + that the order of the shells that you give here must be the same as + the order of the orbitals in the hamiltonian matrix. In the last + example case above the code assumes that matrix index 1 to 5 + belongs to the first d shell, 6 to 10 to the second, 11 to 13 to + the first p shell, and 14 to 16 the second p shell. + +#. Number of correlated shells in the hamiltonian matrix, in the same + spirit as line 3. + +#. The next line(s) contain six numbers: index of the atom, index + of the equivalent shell, :math:`l` quantum number, dimension + of the correlated shells, a spin-orbit parameter, and another + parameter defining interactions. Note that the latter two + parameters are not used at the moment in the code, and only kept + for compatibility reasons. In our example file we use only the + d-shell as correlated, that is why we have only one line here. + +#. The last line contains several numbers: the number of irreducible + representations, and then the dimensions of the irreps. One + possibility is as the example above, another one would be 2 + 2 3. This would mean, 2 irreps (eg and t2g), of dimension 2 and 3, + resp. + +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 +matrix of the imaginary part, and then move on to the next :math:`\mathbf{k}`-point. + +The converter itself is used as:: + + from triqs_dft_tools.converters.hk_converter import * + Converter = HkConverter(filename = hkinputfile) + Converter.convert_dft_input() + +where :file:`hkinputfile` is the name of the input file described +above. This produces the hdf file that you need for a DMFT calculation. + +For more options of this converter, have a look at the +:ref:`refconverters` section of the reference manual. + + diff --git a/doc/guide/conv_vasp.rst b/doc/guide/conv_vasp.rst new file mode 100644 index 00000000..fb80eac7 --- /dev/null +++ b/doc/guide/conv_vasp.rst @@ -0,0 +1,144 @@ +.. _convVASP: + + +Interface with VASP +=================== + +.. warning:: + The VASP interface is in the alpha-version and the VASP part of it is not + yet publicly released. The documentation may, thus, be subject to changes + before the final release. + +*Limitations of the alpha-version:* + + * The interface works correctly only if the k-point symmetries + are turned off during the VASP run (ISYM=-1). + + * Generation of projectors for k-point lines (option `Lines` in KPOINTS) + needed for Bloch spectral function calculations is not possible at the moment. + + * The interface currently supports only collinear-magnetism calculation + (this implis no spin-orbit coupling) and + spin-polarized projectors have not been tested. + +A detailed description of the VASP converter tool PLOVasp can be found +in the :ref:`PLOVasp User's Guide `. Here, a quick-start guide is presented. + +The VASP interface relies on new options introduced since version +5.4.x. In particular, a new INCAR-option `LOCPROJ` +and new `LORBIT` modes 13 and 14 have been added. + +Option `LOCPROJ` selects a set of localized projectors that will +be written to file `LOCPROJ` after a successful VASP run. +A projector set is specified by site indices, +labels of the target local states, and projector type: + + | `LOCPROJ = : : ` + +where `` represents a list of site indices separated by spaces, +with the indices corresponding to the site position in the POSCAR file; +`` specifies local states (see below); +`` chooses a particular type of the local basis function. +The recommended projector type is `Pr 2`. The formalism for this type +of projectors is presented in +`M. Schüler et al. 2018 J. Phys.: Condens. Matter 30 475901 `_. + +The allowed labels of the local states defined in terms of cubic +harmonics are: + + * Entire shells: `s`, `p`, `d`, `f` + + * `p`-states: `py`, `pz`, `px` + + * `d`-states: `dxy`, `dyz`, `dz2`, `dxz`, `dx2-y2` + + * `f`-states: `fy(3x2-y2)`, `fxyz`, `fyz2`, `fz3`, + `fxz2`, `fz(x2-y2)`, `fx(x2-3y2)`. + +For projector type `Pr 2`, one should also set `LORBIT = 14` in the INCAR file +and provide parameters `EMIN`, `EMAX`, defining, in this case, an +energy range (energy window) corresponding to the valence states. +Note that, as in the case +of a DOS calculation, the position of the valence states depends on the +Fermi level, which can usually be found at the end of the OUTCAR file. + +For example, in case of SrVO3 one may first want to perform a self-consistent +calculation, then set `ICHARGE = 1` and add the following additional +lines into INCAR (provided that V is the second ion in POSCAR): + + | `EMIN = 3.0` + | `EMAX = 8.0` + | `LORBIT = 14` + | `LOCPROJ = 2 : d : Pr 2` + +The energy range does not have to be precise. Important is that it has a large +overlap with valence bands and no overlap with semi-core or high unoccupied states. + +Conversion for the DMFT self-consistency cycle +---------------------------------------------- + +The projectors generated by VASP require certain post-processing before +they can be used for DMFT calculations. The most important step is to normalize +them within an energy window that selects band states relevant for the impurity +problem. Note that this energy window is different from the one described above +and it must be chosen independently of the energy +range given by `EMIN, EMAX` in INCAR. + +Post-processing of `LOCPROJ` data is generally done as follows: + +#. Prepare an input file `.cfg` (e.g., `plo.cfg`) that describes the definition + of your impurity problem (more details below). + +#. Extract the value of the Fermi level from OUTCAR and paste it at the end of + the first line of LOCPROJ. + +#. Run :program:`plovasp` with the input file as an argument, e.g.: + + | `plovasp plo.cfg` + + This requires that the TRIQS paths are set correctly (see Installation + of TRIQS). + +If everything goes right one gets files `.ctrl` and `.pg1`. +These files are needed for the converter that will be invoked in your +DMFT script. + +The format of input file `.cfg` is described in details in +the :ref:`User's Guide `. Here we just consider a simple example for the case +of SrVO3: + +.. literalinclude:: images_scripts/srvo3.cfg + +A projector shell is defined by a section `[Shell 1]` where the number +can be arbitrary and used only for user convenience. Several +parameters are required + +- **IONS**: list of site indices which must be a subset of indices + given earlier in `LOCPROJ`. +- **LSHELL**: :math:`l`-quantum number of the projector shell; the corresponding + orbitals must be present in `LOCPROJ`. +- **EWINDOW**: energy window in which the projectors are normalized; + note that the energies are defined with respect to the Fermi level. + +Option **TRANSFORM** is optional but here, it is specified to extract +only three :math:`t_{2g}` orbitals out of five `d` orbitals given by +:math:`l = 2`. + +The conversion to a h5-file is performed in the same way as for Wien2TRIQS:: + + from triqs_dft_tools.converters.vasp_converter import * + Converter = VaspConverter(filename = filename) + Converter.convert_dft_input() + +As usual, the resulting h5-file can then be used with the SumkDFT class. + +Note that the automatic detection of the correct block structure might +fail for VASP inputs. +This can be circumvented by setting a bigger value of the threshold in +:class:`SumkDFT `, e.g.:: + + SK.analyse_block_structure(threshold = 1e-4) + +However, do this only after a careful study of the density matrix and +the projected DOS in the localized basis. + diff --git a/doc/guide/conv_wien2k.rst b/doc/guide/conv_wien2k.rst new file mode 100644 index 00000000..5d888cc0 --- /dev/null +++ b/doc/guide/conv_wien2k.rst @@ -0,0 +1,174 @@ +.. _convWien2k: + +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 in/output files of Wien2k, and how to run +the DFT code. + +Conversion for the DMFT self-consistency cycle +---------------------------------------------- + +First, we have to write the necessary +quantities into a file that can be processed further by invoking in a +shell the command + + `x lapw2 -almd` + +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 +:program:`dmftproj`. For an extensive manual to this program see +:download:`TutorialDmftproj.pdf `. +Here we will only describe the basic steps. + +Let us take the compound SrVO3, a commonly used +example for DFT+DMFT calculations. The input file for +:program:`dmftproj` looks like + +.. literalinclude:: images_scripts/SrVO3.indmftpr + +The first three lines give the number of inequivalent sites, their +multiplicity (to be in accordance with the Wien2k *struct* file) and +the maximum orbital quantum number :math:`l_{max}`. In our case our +struct file contains the atoms in the order Sr, V, O. + +Next we have to +specify for each of the inequivalent sites, whether we want to treat +their orbitals as correlated or not. This information is given by the +following 3 to 5 lines: + +#. We specify which basis set is used (complex or cubic + harmonics). +#. The four numbers refer to *s*, *p*, *d*, and *f* electrons, + resp. Putting 0 means doing nothing, putting 1 will calculate + **unnormalized** projectors in compliance with the Wien2k + definition. The important flag is 2, this means to include these + electrons as correlated electrons, and calculate normalized Wannier + functions for them. In the example above, you see that only for the + vanadium *d* we set the flag to 2. If you want to do simply a DMFT + calculation, then set everything to 0, except one flag 2 for the + correlated electrons. +#. In case you have a irrep splitting of the correlated shell, you can + specify here how many irreps you have. You see that we put 2, since + eg and t2g symmetries are irreps in this cubic case. If you don't + want to use this splitting, just put 0. +#. (optional) If you specifies a number different from 0 in above line, you have + to tell now, which of the irreps you want to be treated + correlated. We want to t2g, and not the eg, so we set 0 for eg and + 1 for t2g. Note that the example above is what you need in 99% of + the cases when you want to treat only t2g electrons. For eg's only + (e.g. nickelates), you set 10 and 01 in this line. +#. (optional) If you have specified a correlated shell for this atom, + you have to tell if spin-orbit coupling should be taken into + account. 0 means no, 1 is yes. + +These lines have to be repeated for each inequivalent atom. + +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! + +After setting up this input file, you run: + + `dmftproj` + +Again, adding possible flags like -so for spin-orbit coupling. This +program produces the following files (in the following, take *case* as +the standard Wien2k place holder, to be replaced by the actual working +directory name): + + * :file:`case.ctqmcout` and :file:`case.symqmc` containing projector + operators and symmetry operations for orthonormalized Wannier + orbitals, respectively. + * :file:`case.parproj` and :file:`case.sympar` containing projector + operators and symmetry operations for uncorrelated states, + respectively. These files are needed for projected + density-of-states or spectral-function calculations in + post-processing only. + * :file:`case.oubwin` needed for the charge density recalculation in + the case of fully self-consistent DFT+DMFT run (see below). + +Now we convert these files into an hdf5 file that can be used for the +DMFT calculations. For this purpose we +use the python module :class:`Wien2kConverter `. It is initialized as:: + + from triqs_dft_tools.converters.wien2k_converter import * + Converter = Wien2kConverter(filename = case) + +The only necessary parameter to this construction is the parameter `filename`. +It has to be the root of the files produces by dmftproj. For our +example, the :program:`Wien2k` naming convention is that all files are +called the same, for instance +:file:`SrVO3.*`, so you would give `filename = "SrVO3"`. The constructor opens +an hdf5 archive, named :file:`case.h5`, where all the data is +stored. For other parameters of the constructor please visit the +:ref:`refconverters` section of the reference manual. + +After initializing the interface module, we can now convert the input +text files to the hdf5 archive by:: + + Converter.convert_dft_input() + +This reads all the data, and stores it in the file :file:`case.h5`. +In this step, the files :file:`case.ctqmcout` and +:file:`case.symqmc` +have to be present in the working directory. + +After this step, all the necessary information for the DMFT loop is +stored in the hdf5 archive, where the string variable +`Converter.hdf_filename` gives the file name of the archive. + +At this point you should use the method :meth:`dos_wannier_basis ` +contained in the module :class:`SumkDFTTools ` to check the density of +states of the Wannier orbitals (see :ref:`analysis`). + +You have now everything for performing a DMFT calculation, and you can +proceed with the section on :ref:`single-shot DFT+DMFT calculations `. + +Data for post-processing +------------------------ + +In case you want to do post-processing of your data using the module +:class:`SumkDFTTools `, some more files +have to be converted to the hdf5 archive. For instance, for +calculating the partial density of states or partial charges +consistent with the definition of :program:`Wien2k`, you have to invoke:: + + Converter.convert_parproj_input() + +This reads and converts the files :file:`case.parproj` and +:file:`case.sympar`. + +If you want to plot band structures, one has to do the +following. First, one has to do the Wien2k calculation on the given +:math:`\mathbf{k}`-path, and run :program:`dmftproj` on that path: + + | `x lapw1 -band` + | `x lapw2 -band -almd` + | `dmftproj -band` + + +Again, maybe with the optional additional extra flags according to +Wien2k. Now we use a routine of the converter module allows to read +and convert the input for :class:`SumkDFTTools `:: + + Converter.convert_bands_input() + +After having converted this input, you can further proceed with the +:ref:`analysis`. For more options on the converter module, please have +a look at the :ref:`refconverters` section of the reference manual. + +Data for transport calculations +------------------------------- + +For the transport calculations, the situation is a bit more involved, +since we need also the :program:`optics` package of Wien2k. Please +look at the section on :ref:`Transport` to see how to do the necessary +steps, including the conversion. + + diff --git a/doc/guide/conversion.rst b/doc/guide/conversion.rst index 40771090..8e4f2bd4 100644 --- a/doc/guide/conversion.rst +++ b/doc/guide/conversion.rst @@ -1,492 +1,27 @@ .. _conversion: -Orbital construction and conversion -=================================== +Supported interfaces +==================== The first step for a DMFT calculation is to provide the necessary input based on a DFT calculation. We will not review how to do the DFT calculation here in this documentation, but refer the user to the documentation and tutorials that come with the actual DFT -package. Here, we will describe how to use output created by Wien2k, -as well as how to use the light-weight general interface. +package. At the moment, there are two full charge self consistent interfaces, for the +Wien2k and the VASP DFT packages, resp. In addition, there is an interface to Wannier90, as well +as a light-weight general-purpose interface. In the following, we will describe the usage of these +conversion tools. -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 in/output files of Wien2k, and how to run -the DFT code. - -Conversion for the DMFT self-consistency cycle -"""""""""""""""""""""""""""""""""""""""""""""" - -First, we have to write the necessary -quantities into a file that can be processed further by invoking in a -shell the command - - `x lapw2 -almd` - -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 -:program:`dmftproj`. For an extensive manual to this program see -:download:`TutorialDmftproj.pdf `. -Here we will only describe the basic steps. - -Let us take the compound SrVO3, a commonly used -example for DFT+DMFT calculations. The input file for -:program:`dmftproj` looks like - -.. literalinclude:: images_scripts/SrVO3.indmftpr - -The first three lines give the number of inequivalent sites, their -multiplicity (to be in accordance with the Wien2k *struct* file) and -the maximum orbital quantum number :math:`l_{max}`. In our case our -struct file contains the atoms in the order Sr, V, O. - -Next we have to -specify for each of the inequivalent sites, whether we want to treat -their orbitals as correlated or not. This information is given by the -following 3 to 5 lines: - -#. We specify which basis set is used (complex or cubic - harmonics). -#. The four numbers refer to *s*, *p*, *d*, and *f* electrons, - resp. Putting 0 means doing nothing, putting 1 will calculate - **unnormalized** projectors in compliance with the Wien2k - definition. The important flag is 2, this means to include these - electrons as correlated electrons, and calculate normalized Wannier - functions for them. In the example above, you see that only for the - vanadium *d* we set the flag to 2. If you want to do simply a DMFT - calculation, then set everything to 0, except one flag 2 for the - correlated electrons. -#. In case you have a irrep splitting of the correlated shell, you can - specify here how many irreps you have. You see that we put 2, since - eg and t2g symmetries are irreps in this cubic case. If you don't - want to use this splitting, just put 0. -#. (optional) If you specifies a number different from 0 in above line, you have - to tell now, which of the irreps you want to be treated - correlated. We want to t2g, and not the eg, so we set 0 for eg and - 1 for t2g. Note that the example above is what you need in 99% of - the cases when you want to treat only t2g electrons. For eg's only - (e.g. nickelates), you set 10 and 01 in this line. -#. (optional) If you have specified a correlated shell for this atom, - you have to tell if spin-orbit coupling should be taken into - account. 0 means no, 1 is yes. - -These lines have to be repeated for each inequivalent atom. - -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! - -After setting up this input file, you run: - - `dmftproj` - -Again, adding possible flags like -so for spin-orbit coupling. This -program produces the following files (in the following, take *case* as -the standard Wien2k place holder, to be replaced by the actual working -directory name): - - * :file:`case.ctqmcout` and :file:`case.symqmc` containing projector - operators and symmetry operations for orthonormalized Wannier - orbitals, respectively. - * :file:`case.parproj` and :file:`case.sympar` containing projector - operators and symmetry operations for uncorrelated states, - respectively. These files are needed for projected - density-of-states or spectral-function calculations in - post-processing only. - * :file:`case.oubwin` needed for the charge density recalculation in - the case of fully self-consistent DFT+DMFT run (see below). - -Now we convert these files into an hdf5 file that can be used for the -DMFT calculations. For this purpose we -use the python module :class:`Wien2kConverter `. It is initialized as:: - - from triqs_dft_tools.converters.wien2k_converter import * - Converter = Wien2kConverter(filename = case) - -The only necessary parameter to this construction is the parameter `filename`. -It has to be the root of the files produces by dmftproj. For our -example, the Wien2k naming convention is that all files are -called the same, for instance -:file:`SrVO3.*`, so you would give `filename = "SrVO3"`. The constructor opens -an hdf5 archive, named :file:`case.h5`, where all the data is -stored. For other parameters of the constructor please visit the -:ref:`refconverters` section of the reference manual. - -After initializing the interface module, we can now convert the input -text files to the hdf5 archive by:: - - Converter.convert_dft_input() - -This reads all the data, and stores it in the file :file:`case.h5`. -In this step, the files :file:`case.ctqmcout` and -:file:`case.symqmc` -have to be present in the working directory. - -After this step, all the necessary information for the DMFT loop is -stored in the hdf5 archive, where the string variable -`Converter.hdf_filename` gives the file name of the archive. - -At this point you should use the method :meth:`dos_wannier_basis ` -contained in the module :class:`SumkDFTTools ` to check the density of -states of the Wannier orbitals (see :ref:`analysis`). - -You have now everything for performing a DMFT calculation, and you can -proceed with the section on :ref:`single-shot DFT+DMFT calculations `. - -Data for post-processing -"""""""""""""""""""""""" - -In case you want to do post-processing of your data using the module -:class:`SumkDFTTools `, some more files -have to be converted to the hdf5 archive. For instance, for -calculating the partial density of states or partial charges -consistent with the definition of Wien2k, you have to invoke:: - - Converter.convert_parproj_input() - -This reads and converts the files :file:`case.parproj` and -:file:`case.sympar`. - -If you want to plot band structures, one has to do the -following. First, one has to do the Wien2k calculation on the given -:math:`\mathbf{k}`-path, and run :program:`dmftproj` on that path: - - | `x lapw1 -band` - | `x lapw2 -band -almd` - | `dmftproj -band` - - -Again, maybe with the optional additional extra flags according to -Wien2k. Now we use a routine of the converter module allows to read -and convert the input for :class:`SumkDFTTools `:: - - Converter.convert_bands_input() - -After having converted this input, you can further proceed with the -:ref:`analysis`. For more options on the converter module, please have -a look at the :ref:`refconverters` section of the reference manual. - -Data for transport calculations -""""""""""""""""""""""""""""""" - -For the transport calculations, the situation is a bit more involved, -since we need also the :program:`optics` package of Wien2k. Please -look at the section on :ref:`Transport` to see how to do the necessary -steps, including the conversion. - -Interface with VASP ---------------------- - -.. warning:: - The VASP interface is in the alpha-version and the VASP part of it is not - yet publicly released. The documentation may, thus, be subject to changes - before the final release. - -Note that this VASP interface relies on new options introduced since version -5.4.x. - -Additionally, the interface only works correctly if the k-point symmetries -are turned off during the VASP run (ISYM=-1). - -The output of raw (non-normalized) projectors is controlled by an INCAR option -LOCPROJ whose complete syntax is described in the VASP documentaion. - -The definition of a projector set starts with specifying which sites -and which local states we are going to project onto. -This information is provided by option LOCPROJ - - | `LOCPROJ = : : ` - -where `` represents a list of site indices separated by spaces, -with the indices corresponding to the site position in the POSCAR file; -`` specifies local states (e.g. :math:`s`, :math:`p`, :math:`d`, -:math:`d_{x^2-y^2}`, etc.); -`` chooses a particular type of the local basis function. - -Some projector types also require parameters `EMIN`, `EMAX` in INCAR to -be set to define an (approximate) energy window corresponding to the -valence states. - -When either a self-consistent (`ICHARG < 10`) or a non-self-consistent -(`ICHARG >= 10`) calculation is done VASP produces file `LOCPROJ` which -will serve as the main input for the conversion routine. - - -Conversion for the DMFT self-consistency cycle -"""""""""""""""""""""""""""""""""""""""""""""" - -In order to use the projectors generated by VASP for defining an -impurity problem they must be processed, i.e. normalized, possibly -transformed, and then converted to a format suitable for DFT_tools scripts. - -Currently, it is necessary to add the Fermi energy by hand as the fifth value -in the first line of the LOCPROJ file before the next steps can be executed. - -The processing of projectors is performed by the program :program:`plovasp` -invoked as - - | `plovasp ` - -where `` is a input file controlling the conversion of projectors. - -The format of input file `` is described in details in -:ref:`plovasp`. Here we just give a simple example for the case -of SrVO3: - -.. literalinclude:: images_scripts/srvo3.cfg - -A projector shell is defined by a section `[Shell 1]` where the number -can be arbitrary and used only for user convenience. Several -parameters are required - -- **IONS**: list of site indices which must be a subset of indices - given earlier in `LOCPROJ`. -- **LSHELL**: :math:`l`-quantum number of the projector shell; the corresponding - orbitals must be present in `LOCPROJ`. -- **EWINDOW**: energy window in which the projectors are normalized; - note that the energies are defined with respect to the Fermi level. - -Option **TRANSFORM** is optional but here it is specified to extract -only three :math:`t_{2g}` orbitals out of five `d` orbitals given by -:math:`l = 2`. - -For the conversion to a h5 file we use on the python level (in analogy to the Wien2kConverter):: - - from triqs_dft_tools.converters.vasp_converter import * - Converter = VaspConverter(filename = filename) - Converter.convert_dft_input() - -As usual, the resulting h5-file can then be used with the SumkDFT class. - -Note that the automatic detection of the correct blockstructure might fail for VASP inputs. -This can be circumvented by increase the :class:`SumkDFT ` threshold to e.g.:: - - SK.analyse_block_structure(threshold = 1e-4) - -However, only do this after a careful study of the density matrix and the dos in the wannier basis. - -A general H(k) --------------- - -In addition to the more complicated Wien2k converter, -:program:`DFTTools` contains also a light converter. It takes only -one inputfile, and creates the necessary hdf outputfile for -the DMFT calculation. The header of this input file has a defined -format, an example is the following (do not use the text/comments in your -input file): - -.. literalinclude:: images_scripts/case.hk - -The lines of this header define - -#. Number of :math:`\mathbf{k}`-points used in the calculation -#. Electron density for setting the chemical potential -#. Number of total atomic shells in the hamiltonian matrix. In short, - this gives the number of lines described in the following. IN the - example file give above this number is 2. -#. The next line(s) contain four numbers each: index of the atom, index - of the equivalent shell, :math:`l` quantum number, dimension - of this shell. Repeat this line for each atomic shell, the number - of the shells is given in the previous line. - - In the example input file given above, we have two inequivalent - atomic shells, one on atom number 1 with a full d-shell (dimension 5), - and one on atom number 2 with one p-shell (dimension 3). - - Other examples for these lines are: - - #. Full d-shell in a material with only one correlated atom in the - unit cell (e.g. SrVO3). One line is sufficient and the numbers - are `1 1 2 5`. - #. Full d-shell in a material with two equivalent atoms in the unit - cell (e.g. FeSe): You need two lines, one for each equivalent - atom. First line is `1 1 2 5`, and the second line is - `2 1 2 5`. The only difference is the first number, which tells on - which atom the shell is located. The second number is the - same in both lines, meaning that both atoms are equivalent. - #. t2g orbitals on two non-equivalent atoms in the unit cell: Two - lines again. First line is `1 1 2 3`, second line `2 2 2 3`. The - difference to the case above is that now also the second number - differs. Therefore, the two shells are treated independently in - the calculation. - #. d-p Hamiltonian in a system with two equivalent atoms each in - the unit cell (e.g. FeSe has two Fe and two Se in the unit - cell). You need for lines. First line `1 1 2 5`, second - line - `2 1 2 5`. These two lines specify Fe as in the case above. For the p - orbitals you need line three as `3 2 1 3` and line four - as `4 2 1 3`. We have 4 atoms, since the first number runs from 1 to 4, - but only two inequivalent atoms, since the second number runs - only form 1 to 2. +.. toctree:: + :maxdepth: 2 - Note that the total dimension of the hamiltonian matrices that are - read in is the sum of all shell dimensions that you specified. For - example number 4 given above we have a dimension of 5+5+3+3=16. It is important - that the order of the shells that you give here must be the same as - the order of the orbitals in the hamiltonian matrix. In the last - example case above the code assumes that matrix index 1 to 5 - belongs to the first d shell, 6 to 10 to the second, 11 to 13 to - the first p shell, and 14 to 16 the second p shell. - -#. Number of correlated shells in the hamiltonian matrix, in the same - spirit as line 3. + conv_wien2k + conv_vasp + conv_W90 + conv_generalhk -#. The next line(s) contain six numbers: index of the atom, index - of the equivalent shell, :math:`l` quantum number, dimension - of the correlated shells, a spin-orbit parameter, and another - parameter defining interactions. Note that the latter two - parameters are not used at the moment in the code, and only kept - for compatibility reasons. In our example file we use only the - d-shell as correlated, that is why we have only one line here. - -#. The last line contains several numbers: the number of irreducible - representations, and then the dimensions of the irreps. One - possibility is as the example above, another one would be 2 - 2 3. This would mean, 2 irreps (eg and t2g), of dimension 2 and 3, - resp. - -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 -matrix of the imaginary part, and then move on to the next :math:`\mathbf{k}`-point. - -The converter itself is used as:: - - from triqs_dft_tools.converters.hk_converter import * - Converter = HkConverter(filename = hkinputfile) - Converter.convert_dft_input() - -where :file:`hkinputfile` is the name of the input file described -above. This produces the hdf file that you need for a DMFT calculation. - -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 -`wannier90 `_ -Maximally Localized Wannier Functions (MLWF) and create a HDF5 archive -suitable for one-shot DMFT calculations with the -:class:`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 triqs_dft_tools.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 (do not use the text/comments in your input file): - -.. 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. In this example we have four lines for the -four V atoms. The MLWFs were constructed for the t\ :sub:`2g` subspace, and thus -we set ``l`` to 2 and ``dim`` to 3 for all V atoms. Further the spin-orbit coupling (``SO``) -is set to 0 and ``irep`` to 0. -As in this example all 4 V atoms are equivalent we set ``sort`` to 0. We note -that, e.g., for a magnetic DMFT calculation the correlated atoms can be made -inequivalent at this point by using different values for ``sort``. - -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 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 ` - 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 ----------- +========== The interface packages are written such that all the file operations are done only on the master node. In general, the philosophy of the @@ -495,8 +30,9 @@ yourself, you have to *manually* broadcast it to the nodes. An exception to this rule is when you use routines from :class:`SumkDFT ` or :class:`SumkDFTTools `, where the broadcasting is done for you. + Interfaces to other packages ----------------------------- +============================ Because of the modular structure, it is straight forward to extend the :ref:`TRIQS ` package in order to work with other band-structure codes. The only necessary requirement is that diff --git a/doc/guide/dftdmft_singleshot.rst b/doc/guide/dftdmft_singleshot.rst index 16d9e75e..1ec43e0a 100644 --- a/doc/guide/dftdmft_singleshot.rst +++ b/doc/guide/dftdmft_singleshot.rst @@ -106,15 +106,15 @@ are present, or if the calculation should start from scratch:: previous_runs = 0 previous_present = False if mpi.is_master_node(): - f = HDFArchive(dft_filename+'.h5','a') - if 'dmft_output' in f: - ar = f['dmft_output'] - if 'iterations' in ar: - previous_present = True - previous_runs = ar['iterations'] + with HDFArchive(dft_filename+'.h5','a') as f: + if 'dmft_output' in f: + ar = f['dmft_output'] + if 'iterations' in ar: + previous_present = True + previous_runs = ar['iterations'] else: f.create_group('dmft_output') - del f + previous_runs = mpi.bcast(previous_runs) previous_present = mpi.bcast(previous_present) @@ -126,9 +126,8 @@ double counting values of the last iteration:: if previous_present: if mpi.is_master_node(): - ar = HDFArchive(dft_filename+'.h5','a') - S.Sigma_iw << ar['dmft_output']['Sigma_iw'] - del ar + with HDFArchive(dft_filename+'.h5','r') as ar: + S.Sigma_iw << ar['dmft_output']['Sigma_iw'] S.Sigma_iw << mpi.bcast(S.Sigma_iw) chemical_potential,dc_imp,dc_energ = SK.load(['chemical_potential','dc_imp','dc_energ']) @@ -153,11 +152,10 @@ functions) can be necessary in order to ensure convergence:: mix = 0.8 # mixing factor if (iteration_number>1 or previous_present): if mpi.is_master_node(): - ar = HDFArchive(dft_filename+'.h5','a') - mpi.report("Mixing Sigma and G with factor %s"%mix) - S.Sigma_iw << mix * S.Sigma_iw + (1.0-mix) * ar['dmft_output']['Sigma_iw'] - S.G_iw << mix * S.G_iw + (1.0-mix) * ar['dmft_output']['G_iw'] - del ar + with HDFArchive(dft_filename+'.h5','r') as ar: + mpi.report("Mixing Sigma and G with factor %s"%mix) + S.Sigma_iw << mix * S.Sigma_iw + (1.0-mix) * ar['dmft_output']['Sigma_iw'] + S.G_iw << mix * S.G_iw + (1.0-mix) * ar['dmft_output']['G_iw'] S.G_iw << mpi.bcast(S.G_iw) S.Sigma_iw << mpi.bcast(S.Sigma_iw) diff --git a/doc/guide/plovasp.rst b/doc/guide/plovasp.rst index a53b1ba3..99f5e1d9 100644 --- a/doc/guide/plovasp.rst +++ b/doc/guide/plovasp.rst @@ -1,37 +1,35 @@ .. _plovasp: -PLOVasp input file -================== +PLOVasp +======= -The general purpose of the PLOVasp tool is to transform -raw, non-normalized projectors generated by VASP into normalized -projectors corresponding to user-defined projected localized orbitals (PLO). -The PLOs can then be used for DFT+DMFT calculations with or without -charge self-consistency. PLOVasp also provides some utilities -for basic analysis of the generated projectors, such as outputting -density matrices, local Hamiltonians, and projected -density of states. +The general purpose of the PLOVasp tool is to transform raw, non-normalized +projectors generated by VASP into normalized projectors corresponding to +user-defined projected localized orbitals (PLO). The PLOs can then be used for +DFT+DMFT calculations with or without charge self-consistency. PLOVasp also +provides some utilities for basic analysis of the generated projectors, such as +outputting density matrices, local Hamiltonians, and projected density of +states. -PLOs are determined by the energy window in which the raw projectors -are normalized. This allows to define either atomic-like strongly -localized Wannier functions (large energy window) or extended -Wannier functions focusing on selected low-energy states (small -energy window). +PLOs are determined by the energy window in which the raw projectors are +normalized. This allows to define either atomic-like strongly localized Wannier +functions (large energy window) or extended Wannier functions focusing on +selected low-energy states (small energy window). -In PLOVasp all projectors sharing the same energy window are combined -into a `projector group`. Technically, this allows one to define -several groups with different energy windows for the same set of -raw projectors. Note, however, that DFTtools does not support projector -groups at the moment but this feature might appear in future releases. +In PLOVasp, all projectors sharing the same energy window are combined into a +`projector group`. Technically, this allows one to define several groups with +different energy windows for the same set of raw projectors. Note, however, +that DFTtools does not support projector groups at the moment but this feature +might appear in future releases. -A set of projectors defined on sites realted to each other either by symmetry -or by sort along with a set of :math:`l`, :math:`m` quantum numbers forms a -`projector shell`. There could be several projectors shells in a -projector group, implying that they will be normalized within -the same energy window. +A set of projectors defined on sites related to each other either by symmetry +or by an atomic sort, along with a set of :math:`l`, :math:`m` quantum numbers, +forms a `projector shell`. There could be several projectors shells in a +projector group, implying that they will be normalized within the same energy +window. -Projector shells and groups are specified by a user-defined input file -whose format is described below. +Projector shells and groups are specified by a user-defined input file whose +format is described below. Input file format ----------------- @@ -43,7 +41,7 @@ Parameters (or 'options') are grouped into sections specified as A PLOVasp input file can contain three types of sections: #. **[General]**: includes parameters that are independent - of a particular projector set, such as the Fermi level, additional + of a particular projector set, such as the Fermi level, additional output (e.g. the density of states), etc. #. **[Group ]**: describes projector groups, i.e. a set of projectors sharing the same energy window and normalization type. @@ -51,8 +49,8 @@ A PLOVasp input file can contain three types of sections: there should be no more than one projector group. #. **[Shell ]**: contains parameters of a projector shell labelled with ``. If there is only one group section and one shell section, - the group section can be omitted and its required parameters can be - given inside the single shell section. + the group section can be omitted but in this case, the group required + parameters must be provided inside the shell section. Section [General] """"""""""""""""" @@ -61,24 +59,24 @@ The entire section is optional and it contains three parameters: * **BASENAME** (string): provides a base name for output files. Default filenames are :file:`vasp.*`. -* **DOSMESH** ([float float] integer): if this parameter is given - projected density of states for each projected orbital will be +* **DOSMESH** ([float float] integer): if this parameter is given, + the projected density of states for each projected orbital will be evaluated and stored to files :file:`pdos_.dat`, where `n` is the - orbital number. The energy - mesh is defined by three numbers: `EMIN` `EMAX` `NPOINTS`. The first two + orbital index. The energy + mesh is defined by three numbers: `EMIN` `EMAX` `NPOINTS`. The first two can be omitted in which case they are taken to be equal to the projector energy window. **Important note**: at the moment this option works only if the tetrahedron integration method (`ISMEAR = -4` or `-5`) is used in VASP to produce `LOCPROJ`. * **EFERMI** (float): provides the Fermi level. This value overrides the one extracted from VASP output files. - + There are no required parameters in this section. Section [Shell] """"""""""""""" -This section specifies a projector shell. Each shell section must be +This section specifies a projector shell. Each `[Shell]` section must be labeled by an index, e.g. `[Shell 1]`. These indices can then be referenced in a `[Group]` section. @@ -87,17 +85,17 @@ In each `[Shell]` section two parameters are required: * **IONS** (list of integer): indices of sites included in the shell. The sites can be given either by a list of integers `IONS = 5 6 7 8` or by a range `IONS = 5..8`. The site indices must be compatible with - POSCAR file. + the POSCAR file. * **LSHELL** (integer): :math:`l` quantum number of the desired local states. It is important that a given combination of site indices and local states -given by `LSHELL` must be present in LOCPROJ file. +given by `LSHELL` must be present in the LOCPROJ file. There are additional optional parameters that allow one to transform the local states: * **TRANSFORM** (matrix): local transformation matrix applied to all states - in the projector shell. The matrix is defined by (multiline) block + in the projector shell. The matrix is defined by a (multiline) block of floats, with each line corresponding to a row. The number of columns must be equal to :math:`2 l + 1`, with :math:`l` given by `LSHELL`. Only real matrices are allowed. This parameter can be useful to select certain subset of @@ -105,14 +103,14 @@ the local states: * **TRANSFILE** (string): name of the file containing transformation matrices for each site. This option allows for a full-fledged functionality when it comes to local state transformations. The format of this file - is described in :ref:`_transformation_file`. + is described :ref:`below `. Section [Group] """"""""""""""" Each defined projector shell must be part of a projector group. In the current -implementation of DFTtools only a single group is supported which can be -labeled by any integer, e.g. `[Group 1]`. This implies that all projector shells +implementation of DFTtools only a single group (labelled by any integer, e.g. `[Group 1]`) +is supported. This implies that all projector shells must be included in this group. Required parameters for any group are the following: @@ -121,34 +119,49 @@ Required parameters for any group are the following: All defined shells must be grouped. * **EWINDOW** (float float): the energy window specified by two floats: bottom and top. All projectors in the current group are going to be normalized within - this window. + this window. *Note*: This option must be specified inside the `[Shell]` section + if only one shell is defined and the `[Group]` section is omitted. Optional group parameters: * **NORMALIZE** (True/False): specifies whether projectors in the group are - to be noramlized. The default value is **True**. + to be normalized. The default value is **True**. * **NORMION** (True/False): specifies whether projectors are normalized on - a per-site (per-ion) basis. That is, if `NORMION = True` the orthogonality + a per-site (per-ion) basis. That is, if `NORMION = True`, the orthogonality condition will be enforced on each site separately but the Wannier functions - on different sites will not be orthogonal. If `NORMION = False` Wannier functions + on different sites will not be orthogonal. If `NORMION = False`, the Wannier functions on different sites included in the group will be orthogonal to each other. - -.. _transformation_file + +.. _transformation_file: File of transformation matrices """"""""""""""""""""""""""""""" .. warning:: - The description below applies only to collinear cases (i.e. without spin-orbit - coupling). In this case the matrices are spin-independent. + The description below applies only to collinear cases (i.e., without spin-orbit + coupling). In this case, the matrices are spin-independent. The file specified by option `TRANSFILE` contains transformation matrices for each ion. Each line must contain a series of floats whose number is either equal to the number of orbitals :math:`N_{orb}` (in this case the transformation matrices are assumed to be real) or to :math:`2 N_{orb}` (for the complex transformation matrices). -The number of lines :math:`N` must be a multiple of the number of ions :math:`N_{ion}` +The total number of lines :math:`N` must be a multiple of the number of ions :math:`N_{ion}` and the ratio :math:`N / N_{ion}`, then, gives the dimension of the transformed orbital space. The lines with floats can be separated by any number of empty or -comment lines which are ignored. +comment lines (starting from `#`), which are ignored. + +A very simple example is a transformation matrix that selects the :math:`t_{2g}` manifold. +For two correlated sites, one can define the file as follows: +:: + + # Site 1 + 1.0 0.0 0.0 0.0 0.0 + 0.0 1.0 0.0 0.0 0.0 + 0.0 0.0 0.0 1.0 0.0 + + # Site 2 + 1.0 0.0 0.0 0.0 0.0 + 0.0 1.0 0.0 0.0 0.0 + 0.0 0.0 0.0 1.0 0.0 diff --git a/doc/guide/transport.rst b/doc/guide/transport.rst index 3f48bbbb..2aa1d1e5 100644 --- a/doc/guide/transport.rst +++ b/doc/guide/transport.rst @@ -96,12 +96,11 @@ The converter :meth:`convert_transport_input 1 or previous_present): if (mpi.is_master_node() and (mixing<1.0)): - ar = HDFArchive(dft_filename+'.h5','a') - 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 + with HDFArchive(dft_filename+'.h5','r') as ar: + 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'] S.G << mpi.bcast(S.G) S.Sigma << mpi.bcast(S.Sigma) @@ -106,11 +103,10 @@ for iteration_number in range(1,Loops+1): # store the impurity self-energy, GF as well as correlation energy in h5 if mpi.is_master_node(): - ar = HDFArchive(dft_filename+'.h5','a') - ar['dmft_output']['iterations'] = iteration_number + previous_runs - ar['dmft_output']['G'] = S.G - ar['dmft_output']['Sigma'] = S.Sigma - del ar + with HDFArchive(dft_filename+'.h5','a') as ar: + ar['dmft_output']['iterations'] = iteration_number + previous_runs + ar['dmft_output']['G'] = S.G + ar['dmft_output']['Sigma'] = S.Sigma #Save essential SumkDFT data: SK.save(['chemical_potential','dc_imp','dc_energ','correnerg']) diff --git a/doc/guide/images_scripts/Ce-gamma.struct b/doc/tutorials/images_scripts/Ce-gamma.struct similarity index 100% rename from doc/guide/images_scripts/Ce-gamma.struct rename to doc/tutorials/images_scripts/Ce-gamma.struct diff --git a/doc/guide/images_scripts/Ce-gamma_DOS.py b/doc/tutorials/images_scripts/Ce-gamma_DOS.py similarity index 100% rename from doc/guide/images_scripts/Ce-gamma_DOS.py rename to doc/tutorials/images_scripts/Ce-gamma_DOS.py diff --git a/doc/guide/images_scripts/Ce-gamma_DOS_script.rst b/doc/tutorials/images_scripts/Ce-gamma_DOS_script.rst similarity index 100% rename from doc/guide/images_scripts/Ce-gamma_DOS_script.rst rename to doc/tutorials/images_scripts/Ce-gamma_DOS_script.rst diff --git a/doc/guide/images_scripts/Ce-gamma_script.rst b/doc/tutorials/images_scripts/Ce-gamma_script.rst similarity index 100% rename from doc/guide/images_scripts/Ce-gamma_script.rst rename to doc/tutorials/images_scripts/Ce-gamma_script.rst diff --git a/doc/guide/images_scripts/Ce_DOS.png b/doc/tutorials/images_scripts/Ce_DOS.png similarity index 100% rename from doc/guide/images_scripts/Ce_DOS.png rename to doc/tutorials/images_scripts/Ce_DOS.png diff --git a/doc/tutorials/images_scripts/SrVO3.indmftpr b/doc/tutorials/images_scripts/SrVO3.indmftpr new file mode 100644 index 00000000..6f36362d --- /dev/null +++ b/doc/tutorials/images_scripts/SrVO3.indmftpr @@ -0,0 +1,15 @@ +3 ! Nsort +1 1 3 ! Mult(Nsort) +3 ! lmax +complex ! choice of angular harmonics +1 0 0 0 ! l included for each sort +0 0 0 0 ! If split into ireps, gives number of ireps. for a given orbital (otherwise 0) +cubic ! choice of angular harmonics +1 1 2 0 ! l included for each sort +0 0 2 0 ! If split into ireps, gives number of ireps. for a given orbital (otherwise 0) +01 ! +0 ! SO flag +complex ! choice of angular harmonics +1 1 0 0 ! l included for each sort +0 0 0 0 ! If split into ireps, gives number of ireps. for a given orbital (otherwise 0) +-0.11 0.14 \ No newline at end of file diff --git a/doc/tutorials/images_scripts/SrVO3.struct b/doc/tutorials/images_scripts/SrVO3.struct new file mode 100644 index 00000000..d0f2642f --- /dev/null +++ b/doc/tutorials/images_scripts/SrVO3.struct @@ -0,0 +1,25 @@ +SrVO3 +P LATTICE,NONEQUIV.ATOMS: 3221_Pm-3m +MODE OF CALC=RELA unit=bohr + 7.261300 7.261300 7.261300 90.000000 90.000000 90.000000 +ATOM 1: X=0.00000000 Y=0.00000000 Z=0.00000000 + MULT= 1 ISPLIT= 2 +Sr NPT= 781 R0=0.00001000 RMT= 2.50000 Z: 38.0 +LOCAL ROT MATRIX: 1.0000000 0.0000000 0.0000000 + 0.0000000 1.0000000 0.0000000 + 0.0000000 0.0000000 1.0000000 +ATOM 2: X=0.50000000 Y=0.50000000 Z=0.50000000 + MULT= 1 ISPLIT= 2 +V NPT= 781 R0=0.00005000 RMT= 1.91 Z: 23.0 +LOCAL ROT MATRIX: 1.0000000 0.0000000 0.0000000 + 0.0000000 1.0000000 0.0000000 + 0.0000000 0.0000000 1.0000000 +ATOM -3: X=0.00000000 Y=0.50000000 Z=0.50000000 + MULT= 3 ISPLIT=-2 + -3: X=0.50000000 Y=0.00000000 Z=0.50000000 + -3: X=0.50000000 Y=0.50000000 Z=0.00000000 +O NPT= 781 R0=0.00010000 RMT= 1.70 Z: 8.0 +LOCAL ROT MATRIX: 0.0000000 0.0000000 1.0000000 + 0.0000000 1.0000000 0.0000000 + -1.0000000 0.0000000 0.0000000 + 0 NUMBER OF SYMMETRY OPERATIONS diff --git a/doc/guide/images_scripts/SrVO3_Sigma_iw_it1.png b/doc/tutorials/images_scripts/SrVO3_Sigma_iw_it1.png similarity index 100% rename from doc/guide/images_scripts/SrVO3_Sigma_iw_it1.png rename to doc/tutorials/images_scripts/SrVO3_Sigma_iw_it1.png diff --git a/doc/guide/images_scripts/dft_dmft_cthyb.py b/doc/tutorials/images_scripts/dft_dmft_cthyb.py similarity index 73% rename from doc/guide/images_scripts/dft_dmft_cthyb.py rename to doc/tutorials/images_scripts/dft_dmft_cthyb.py index 3fcdceb6..e0b7dad1 100644 --- a/doc/guide/images_scripts/dft_dmft_cthyb.py +++ b/doc/tutorials/images_scripts/dft_dmft_cthyb.py @@ -38,15 +38,15 @@ p["fit_max_n"] = 60 previous_runs = 0 previous_present = False if mpi.is_master_node(): - f = HDFArchive(dft_filename+'.h5','a') - if 'dmft_output' in f: - ar = f['dmft_output'] - if 'iterations' in ar: - previous_present = True - previous_runs = ar['iterations'] - else: - f.create_group('dmft_output') - del f + with HDFArchive(dft_filename+'.h5','a') as f: + if 'dmft_output' in f: + ar = f['dmft_output'] + if 'iterations' in ar: + previous_present = True + previous_runs = ar['iterations'] + else: + f.create_group('dmft_output') + previous_runs = mpi.bcast(previous_runs) previous_present = mpi.bcast(previous_present) @@ -72,9 +72,8 @@ if previous_present: dc_imp = 0 dc_energ = 0 if mpi.is_master_node(): - ar = HDFArchive(dft_filename+'.h5','a') - S.Sigma_iw << ar['dmft_output']['Sigma_iw'] - del ar + with HDFArchive(dft_filename+'.h5','r') as ar: + S.Sigma_iw << ar['dmft_output']['Sigma_iw'] chemical_potential,dc_imp,dc_energ = SK.load(['chemical_potential','dc_imp','dc_energ']) S.Sigma_iw << mpi.bcast(S.Sigma_iw) chemical_potential = mpi.bcast(chemical_potential) @@ -103,14 +102,13 @@ for iteration_number in range(1,loops+1): # We can do a mixing of Delta in order to stabilize the DMFT iterations: S.G0_iw << S.Sigma_iw + inverse(S.G_iw) # The following lines are uncommented until issue #98 is fixed in TRIQS - # ar = HDFArchive(dft_filename+'.h5','a') - # if (iteration_number>1 or previous_present): - # mpi.report("Mixing input Delta with factor %s"%delta_mix) - # Delta = (delta_mix * delta(S.G0_iw)) + (1.0-delta_mix) * ar['dmft_output']['Delta_iw'] - # S.G0_iw << S.G0_iw + delta(S.G0_iw) - Delta - # ar['dmft_output']['Delta_iw'] = delta(S.G0_iw) + # with HDFArchive(dft_filename+'.h5','a') as ar: + # if (iteration_number>1 or previous_present): + # mpi.report("Mixing input Delta with factor %s"%delta_mix) + # Delta = (delta_mix * delta(S.G0_iw)) + (1.0-delta_mix) * ar['dmft_output']['Delta_iw'] + # S.G0_iw << S.G0_iw + delta(S.G0_iw) - Delta + # ar['dmft_output']['Delta_iw'] = delta(S.G0_iw) S.G0_iw << inverse(S.G0_iw) - # del ar S.G0_iw << mpi.bcast(S.G0_iw) @@ -123,25 +121,24 @@ for iteration_number in range(1,loops+1): # Now mix Sigma and G with factor sigma_mix, if wanted: if (iteration_number>1 or previous_present): if mpi.is_master_node(): - ar = HDFArchive(dft_filename+'.h5','a') - mpi.report("Mixing Sigma and G with factor %s"%sigma_mix) - S.Sigma_iw << sigma_mix * S.Sigma_iw + (1.0-sigma_mix) * ar['dmft_output']['Sigma_iw'] - S.G_iw << sigma_mix * S.G_iw + (1.0-sigma_mix) * ar['dmft_output']['G_iw'] - del ar + with HDFArchive(dft_filename+'.h5','r') as ar: + mpi.report("Mixing Sigma and G with factor %s"%sigma_mix) + S.Sigma_iw << sigma_mix * S.Sigma_iw + (1.0-sigma_mix) * ar['dmft_output']['Sigma_iw'] + S.G_iw << sigma_mix * S.G_iw + (1.0-sigma_mix) * ar['dmft_output']['G_iw'] + S.G_iw << mpi.bcast(S.G_iw) S.Sigma_iw << mpi.bcast(S.Sigma_iw) # Write the final Sigma and G to the hdf5 archive: if mpi.is_master_node(): - ar = HDFArchive(dft_filename+'.h5','a') - ar['dmft_output']['iterations'] = iteration_number + previous_runs - ar['dmft_output']['G_tau'] = S.G_tau - ar['dmft_output']['G_iw'] = S.G_iw - ar['dmft_output']['Sigma_iw'] = S.Sigma_iw - ar['dmft_output']['G0-%s'%(iteration_number)] = S.G0_iw - ar['dmft_output']['G-%s'%(iteration_number)] = S.G_iw - ar['dmft_output']['Sigma-%s'%(iteration_number)] = S.Sigma_iw - del ar + with ar = HDFArchive(dft_filename+'.h5','a') as ar: + ar['dmft_output']['iterations'] = iteration_number + previous_runs + ar['dmft_output']['G_tau'] = S.G_tau + ar['dmft_output']['G_iw'] = S.G_iw + ar['dmft_output']['Sigma_iw'] = S.Sigma_iw + ar['dmft_output']['G0-%s'%(iteration_number)] = S.G0_iw + ar['dmft_output']['G-%s'%(iteration_number)] = S.G_iw + ar['dmft_output']['Sigma-%s'%(iteration_number)] = S.Sigma_iw # Set the new double counting: dm = S.G_iw.density() # compute the density matrix of the impurity problem diff --git a/doc/guide/images_scripts/dft_dmft_cthyb_slater.py b/doc/tutorials/images_scripts/dft_dmft_cthyb_slater.py similarity index 73% rename from doc/guide/images_scripts/dft_dmft_cthyb_slater.py rename to doc/tutorials/images_scripts/dft_dmft_cthyb_slater.py index 5150f6e4..1bf40ba2 100644 --- a/doc/guide/images_scripts/dft_dmft_cthyb_slater.py +++ b/doc/tutorials/images_scripts/dft_dmft_cthyb_slater.py @@ -39,15 +39,14 @@ p["fit_max_n"] = 60 previous_runs = 0 previous_present = False if mpi.is_master_node(): - f = HDFArchive(dft_filename+'.h5','a') - if 'dmft_output' in f: - ar = f['dmft_output'] - if 'iterations' in ar: - previous_present = True - previous_runs = ar['iterations'] - else: - f.create_group('dmft_output') - del f + with HDFArchive(dft_filename+'.h5','a') as f: + if 'dmft_output' in f: + ar = f['dmft_output'] + if 'iterations' in ar: + previous_present = True + previous_runs = ar['iterations'] + else: + f.create_group('dmft_output') previous_runs = mpi.bcast(previous_runs) previous_present = mpi.bcast(previous_present) @@ -75,9 +74,8 @@ if previous_present: dc_imp = 0 dc_energ = 0 if mpi.is_master_node(): - ar = HDFArchive(dft_filename+'.h5','a') - S.Sigma_iw << ar['dmft_output']['Sigma_iw'] - del ar + with HDFArchive(dft_filename+'.h5','r') as ar: + S.Sigma_iw << ar['dmft_output']['Sigma_iw'] chemical_potential,dc_imp,dc_energ = SK.load(['chemical_potential','dc_imp','dc_energ']) S.Sigma_iw << mpi.bcast(S.Sigma_iw) chemical_potential = mpi.bcast(chemical_potential) @@ -106,14 +104,13 @@ for iteration_number in range(1,loops+1): # We can do a mixing of Delta in order to stabilize the DMFT iterations: S.G0_iw << S.Sigma_iw + inverse(S.G_iw) # The following lines are uncommented until issue #98 is fixed in TRIQS - # ar = HDFArchive(dft_filename+'.h5','a') - # if (iteration_number>1 or previous_present): - # mpi.report("Mixing input Delta with factor %s"%delta_mix) - # Delta = (delta_mix * delta(S.G0_iw)) + (1.0-delta_mix) * ar['dmft_output']['Delta_iw'] - # S.G0_iw << S.G0_iw + delta(S.G0_iw) - Delta - # ar['dmft_output']['Delta_iw'] = delta(S.G0_iw) + # with HDFArchive(dft_filename+'.h5','a') as ar: + # if (iteration_number>1 or previous_present): + # mpi.report("Mixing input Delta with factor %s"%delta_mix) + # Delta = (delta_mix * delta(S.G0_iw)) + (1.0-delta_mix) * ar['dmft_output']['Delta_iw'] + # S.G0_iw << S.G0_iw + delta(S.G0_iw) - Delta + # ar['dmft_output']['Delta_iw'] = delta(S.G0_iw) S.G0_iw << inverse(S.G0_iw) - # del ar S.G0_iw << mpi.bcast(S.G0_iw) @@ -126,25 +123,23 @@ for iteration_number in range(1,loops+1): # Now mix Sigma and G with factor sigma_mix, if wanted: if (iteration_number>1 or previous_present): if mpi.is_master_node(): - ar = HDFArchive(dft_filename+'.h5','a') - mpi.report("Mixing Sigma and G with factor %s"%sigma_mix) - S.Sigma_iw << sigma_mix * S.Sigma_iw + (1.0-sigma_mix) * ar['dmft_output']['Sigma_iw'] - S.G_iw << sigma_mix * S.G_iw + (1.0-sigma_mix) * ar['dmft_output']['G_iw'] - del ar + with HDFArchive(dft_filename+'.h5','r') as ar: + mpi.report("Mixing Sigma and G with factor %s"%sigma_mix) + S.Sigma_iw << sigma_mix * S.Sigma_iw + (1.0-sigma_mix) * ar['dmft_output']['Sigma_iw'] + S.G_iw << sigma_mix * S.G_iw + (1.0-sigma_mix) * ar['dmft_output']['G_iw'] S.G_iw << mpi.bcast(S.G_iw) S.Sigma_iw << mpi.bcast(S.Sigma_iw) # Write the final Sigma and G to the hdf5 archive: if mpi.is_master_node(): - ar = HDFArchive(dft_filename+'.h5','a') - ar['dmft_output']['iterations'] = iteration_number + previous_runs - ar['dmft_output']['G_tau'] = S.G_tau - ar['dmft_output']['G_iw'] = S.G_iw - ar['dmft_output']['Sigma_iw'] = S.Sigma_iw - ar['dmft_output']['G0-%s'%(iteration_number)] = S.G0_iw - ar['dmft_output']['G-%s'%(iteration_number)] = S.G_iw - ar['dmft_output']['Sigma-%s'%(iteration_number)] = S.Sigma_iw - del ar + with HDFArchive(dft_filename+'.h5','a') as ar: + ar['dmft_output']['iterations'] = iteration_number + previous_runs + ar['dmft_output']['G_tau'] = S.G_tau + ar['dmft_output']['G_iw'] = S.G_iw + ar['dmft_output']['Sigma_iw'] = S.Sigma_iw + ar['dmft_output']['G0-%s'%(iteration_number)] = S.G0_iw + ar['dmft_output']['G-%s'%(iteration_number)] = S.G_iw + ar['dmft_output']['Sigma-%s'%(iteration_number)] = S.Sigma_iw # Set the new double counting: dm = S.G_iw.density() # compute the density matrix of the impurity problem diff --git a/doc/guide/SrVO3.rst b/doc/tutorials/srvo3.rst similarity index 74% rename from doc/guide/SrVO3.rst rename to doc/tutorials/srvo3.rst index 64adb23b..3239da3f 100644 --- a/doc/guide/SrVO3.rst +++ b/doc/tutorials/srvo3.rst @@ -1,14 +1,10 @@ .. _SrVO3: -SrVO3 (single-shot) -=================== - -We will discuss now how to set up a full working calculation, +On the example of SrVO3 we will discuss now how to set up a full working calculation, including the initialization of the :ref:`CTHYB solver `. Some additional parameter are introduced to make the calculation more efficient. This is a more advanced example, which is -also suited for parallel execution. The conversion, which -we assume to be carried out already, is discussed :ref:`here `. +also suited for parallel execution. For the convenience of the user, we provide also two working python scripts in this documentation. One for a calculation @@ -18,6 +14,63 @@ rotational-invariant Slater interaction Hamiltonian (:download:`dft_dmft_cthyb_s `). The user has to adapt these scripts to his own needs. How to execute your script is described :ref:`here`. +The conversion will now be discussed in detail for the Wien2k and VASP packages. +For more details we refer to the :ref:`documentation `. + + +DFT (Wien2k) and Wannier orbitals +================================= + +DFT setup +--------- + +First, we do a DFT calculation, using the Wien2k package. As main input file we have to provide the so-called struct file :file:`SrVO3.struct`. We use the following: + +.. literalinclude:: images_scripts/SrVO3.struct + +Instead of going through the whole initialisation process, we can use :: + + init -b -vxc 5 -numk 5000 + +This is setting up a non-magnetic calculation, using the LDA and 5000 k-points in the full Brillouin zone. As usual, we start the DFT self consistent cycle by the Wien2k script :: + + run + +Wannier orbitals +---------------- + +As a next step, we calculate localised orbitals for the t2g orbitals. We use the same input file for :program:`dmftproj` as it was used in the :ref:`documentation`: + +.. literalinclude:: images_scripts/SrVO3.indmftpr + +To prepare the input data for :program:`dmftproj` we execute lapw2 with the `-almd` option :: + + x lapw2 -almd + +Then :program:`dmftproj` is executed in its default mode (i.e. without spin-polarization or spin-orbit included) :: + + dmftproj + +This program produces the necessary files for the conversion to the hdf5 file structure. This is done using +the python module :class:`Wien2kConverter `. A simple python script that initialises the converter is:: + + from triqs_dft_tools.converters.wien2k_converter import * + Converter = Wien2kConverter(filename = "SrVO3") + +After initializing the interface module, we can now convert the input +text files to the hdf5 archive by:: + + Converter.convert_dft_input() + +This reads all the data, and stores everything that is necessary for the DMFT calculation in the file :file:`SrVO3.h5`. + + +The DMFT calculation +==================== + +The DMFT script itself is, except very few details, independent of the DFT package that was used to calculate the local orbitals. +As soon as one has converted everything to the hdf5 format, the following procedure is practially the same. + Loading modules --------------- @@ -152,23 +205,21 @@ some additional refinements:: # Now mix Sigma and G with factor mix, if wanted: if (iteration_number>1 or previous_present): if mpi.is_master_node(): - ar = HDFArchive(dft_filename+'.h5','a') - mpi.report("Mixing Sigma and G with factor %s"%mix) - S.Sigma_iw << mix * S.Sigma_iw + (1.0-mix) * ar['dmft_output']['Sigma_iw'] - S.G_iw << mix * S.G_iw + (1.0-mix) * ar['dmft_output']['G_iw'] - del ar + with HDFArchive(dft_filename+'.h5','r') as ar: + mpi.report("Mixing Sigma and G with factor %s"%mix) + S.Sigma_iw << mix * S.Sigma_iw + (1.0-mix) * ar['dmft_output']['Sigma_iw'] + S.G_iw << mix * S.G_iw + (1.0-mix) * ar['dmft_output']['G_iw'] S.G_iw << mpi.bcast(S.G_iw) S.Sigma_iw << mpi.bcast(S.Sigma_iw) # Write the final Sigma and G to the hdf5 archive: if mpi.is_master_node(): - ar = HDFArchive(dft_filename+'.h5','a') - ar['dmft_output']['iterations'] = iteration_number - ar['dmft_output']['G_0'] = S.G0_iw - ar['dmft_output']['G_tau'] = S.G_tau - ar['dmft_output']['G_iw'] = S.G_iw - ar['dmft_output']['Sigma_iw'] = S.Sigma_iw - del ar + with HDFArchive(dft_filename+'.h5','a') as ar: + ar['dmft_output']['iterations'] = iteration_number + ar['dmft_output']['G_0'] = S.G0_iw + ar['dmft_output']['G_tau'] = S.G_tau + ar['dmft_output']['G_iw'] = S.G_iw + ar['dmft_output']['Sigma_iw'] = S.Sigma_iw # Set the new double counting: dm = S.G_iw.density() # compute the density matrix of the impurity problem diff --git a/doc/vasp/source/plotools.rst b/doc/vasp/source/plotools.rst index b9fe5d52..61d99513 100644 --- a/doc/vasp/source/plotools.rst +++ b/doc/vasp/source/plotools.rst @@ -46,7 +46,7 @@ electronic structure data. At this stage simple consistency checks are performed All electronic structure from VASP is stored in a class ElectronicStructure: -.. autoclass:: elstruct.ElectronicStructure +.. autoclass:: triqs_dft_tools.converters.plovasp.elstruct.ElectronicStructure :members: @@ -95,7 +95,7 @@ Order of operations: * distribute back the arrays assuming that the order is preserved -.. autoclass:: proj_shell.ProjectorShell +.. autoclass:: triqs_dft_tools.converters.plovasp.proj_shell.ProjectorShell :members: diff --git a/doc/vasp/source/vaspio.rst b/doc/vasp/source/vaspio.rst index 4cadd5c1..b8d568cf 100644 --- a/doc/vasp/source/vaspio.rst +++ b/doc/vasp/source/vaspio.rst @@ -1,4 +1,4 @@ -.. sec_vaspio +.. _vaspio: VASP input-output ################# diff --git a/python/converters/hk_converter.py b/python/converters/hk_converter.py index 31cccd47..510a7cca 100644 --- a/python/converters/hk_converter.py +++ b/python/converters/hk_converter.py @@ -260,13 +260,12 @@ class HkConverter(ConverterTools): R.close() # Save to the HDF5: - ar = HDFArchive(self.hdf_file, 'a') - if not (self.dft_subgrp in ar): - ar.create_group(self.dft_subgrp) - things_to_save = ['energy_unit', 'n_k', 'k_dep_projection', 'SP', 'SO', 'charge_below', 'density_required', + with HDFArchive(self.hdf_file, 'a') as ar: + if not (self.dft_subgrp in ar): + ar.create_group(self.dft_subgrp) + 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 + for it in things_to_save: + ar[self.dft_subgrp][it] = locals()[it] diff --git a/python/converters/plovasp/examples/ce/test_ham_hf.py b/python/converters/plovasp/examples/ce/test_ham_hf.py index 1243b85c..56fb20df 100644 --- a/python/converters/plovasp/examples/ce/test_ham_hf.py +++ b/python/converters/plovasp/examples/ce/test_ham_hf.py @@ -44,10 +44,9 @@ class TestSumkDFT(SumkDFT): fermi_weights = 0 band_window = 0 if mpi.is_master_node(): - ar = HDFArchive(self.hdf_file,'r') - fermi_weights = ar['dft_misc_input']['dft_fermi_weights'] - band_window = ar['dft_misc_input']['band_window'] - del ar + with HDFArchive(self.hdf_file,'r') as ar: + fermi_weights = ar['dft_misc_input']['dft_fermi_weights'] + band_window = ar['dft_misc_input']['band_window'] fermi_weights = mpi.bcast(fermi_weights) band_window = mpi.bcast(band_window) @@ -184,10 +183,9 @@ class TestSumkDFT(SumkDFT): fermi_weights = 0 band_window = 0 if mpi.is_master_node(): - ar = HDFArchive(self.hdf_file,'r') - fermi_weights = ar['dft_misc_input']['dft_fermi_weights'] - band_window = ar['dft_misc_input']['band_window'] - del ar + with HDFArchive(self.hdf_file,'r') as ar: + fermi_weights = ar['dft_misc_input']['dft_fermi_weights'] + band_window = ar['dft_misc_input']['band_window'] fermi_weights = mpi.bcast(fermi_weights) band_window = mpi.bcast(band_window) @@ -282,14 +280,13 @@ def dmft_cycle(): previous_present = False if mpi.is_master_node(): - ar = HDFArchive(HDFfilename,'a') - if 'iterations' in ar: - previous_present = True - previous_runs = ar['iterations'] - else: - previous_runs = 0 - previous_present = False - del ar + with HDFArchive(HDFfilename,'a') as ar: + if 'iterations' in ar: + previous_present = True + previous_runs = ar['iterations'] + else: + previous_runs = 0 + previous_present = False mpi.barrier() previous_runs = mpi.bcast(previous_runs) @@ -315,9 +312,8 @@ def dmft_cycle(): 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 + with HDFArchive(HDFfilename,'a') as ar: + S.Sigma <<= ar['SigmaF'] things_to_load=['chemical_potential','dc_imp'] old_data=SK.load(things_to_load) chemical_potential=old_data[0] @@ -365,13 +361,12 @@ def dmft_cycle(): # 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'] - del ar + with HDFArchive(HDFfilename,'r') as ar: + 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'] S.G = mpi.bcast(S.G) S.Sigma = mpi.bcast(S.Sigma) @@ -386,14 +381,13 @@ def dmft_cycle(): # 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 - del ar + with HDFArchive(HDFfilename,'a') as ar: + 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 #Save essential SumkDFT data: things_to_save=['chemical_potential','dc_energ','dc_imp'] @@ -428,11 +422,10 @@ def dmft_cycle(): # 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 + with HDFArchive(HDFfilename) as ar: + itn = ar['iterations'] + correnerg = ar['correnerg%s'%itn] + DCenerg = ar['DCenerg%s'%itn] correnerg -= DCenerg[0] f=open(lda_filename+'.qdmft','a') f.write("%.16f\n"%correnerg) diff --git a/python/converters/vasp_converter.py b/python/converters/vasp_converter.py index 8642644c..6be2d75e 100644 --- a/python/converters/vasp_converter.py +++ b/python/converters/vasp_converter.py @@ -269,22 +269,23 @@ class VaspConverter(ConverterTools): # Save it to the HDF: - 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', + with HDFArchive(self.hdf_file,'a') as ar: + 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] + for it in things_to_save: ar[self.dft_subgrp][it] = locals()[it] -# Store Fermi weights to 'dft_misc_input' - if not (self.misc_subgrp in ar): ar.create_group(self.misc_subgrp) - ar[self.misc_subgrp]['dft_fermi_weights'] = f_weights - ar[self.misc_subgrp]['band_window'] = band_window - del ar + # Store Fermi weights to 'dft_misc_input' + if not (self.misc_subgrp in ar): ar.create_group(self.misc_subgrp) + ar[self.misc_subgrp]['dft_fermi_weights'] = f_weights + ar[self.misc_subgrp]['band_window'] = band_window + # Symmetries are used, so now convert symmetry information for *correlated* orbitals: self.convert_symmetry_input(ctrl_head, orbits=self.corr_shells, symm_subgrp=self.symmcorr_subgrp) + # TODO: Implement misc_input # self.convert_misc_input(bandwin_file=self.bandwin_file,struct_file=self.struct_file,outputs_file=self.outputs_file, # misc_subgrp=self.misc_subgrp,SO=self.SO,SP=self.SP,n_k=self.n_k) @@ -381,10 +382,9 @@ class VaspConverter(ConverterTools): raise "convert_misc_input: reading file %s failed" %self.outputs_file # Save it to the HDF: - ar=HDFArchive(self.hdf_file,'a') - if not (misc_subgrp in ar): ar.create_group(misc_subgrp) - for it in things_to_save: ar[misc_subgrp][it] = locals()[it] - del ar + with HDFArchive(self.hdf_file,'a') as ar: + if not (misc_subgrp in ar): ar.create_group(misc_subgrp) + for it in things_to_save: ar[misc_subgrp][it] = locals()[it] def convert_symmetry_input(self, ctrl_head, orbits, symm_subgrp): @@ -405,10 +405,8 @@ class VaspConverter(ConverterTools): mat_tinv = [numpy.identity(1)] # Save it to the HDF: - ar=HDFArchive(self.hdf_file,'a') - if not (symm_subgrp in ar): ar.create_group(symm_subgrp) - things_to_save = ['n_symm','n_atoms','perm','orbits','SO','SP','time_inv','mat','mat_tinv'] - for it in things_to_save: -# print "%s:"%(it), locals()[it] - ar[symm_subgrp][it] = locals()[it] - del ar + with HDFArchive(self.hdf_file,'a') as ar: + if not (symm_subgrp in ar): ar.create_group(symm_subgrp) + things_to_save = ['n_symm','n_atoms','perm','orbits','SO','SP','time_inv','mat','mat_tinv'] + for it in things_to_save: + ar[symm_subgrp][it] = locals()[it] diff --git a/python/converters/wannier90_converter.py b/python/converters/wannier90_converter.py index 14f2f71b..925d14e3 100644 --- a/python/converters/wannier90_converter.py +++ b/python/converters/wannier90_converter.py @@ -345,18 +345,17 @@ class Wannier90Converter(ConverterTools): 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', + with HDFArchive(self.hdf_file, 'a') as ar: + 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 + for it in things_to_save: + ar[self.dft_subgrp][it] = locals()[it] def read_wannier90hr(self, hr_filename="wannier_hr.dat"): """ diff --git a/python/converters/wien2k_converter.py b/python/converters/wien2k_converter.py index e8419fa9..082b1cbd 100644 --- a/python/converters/wien2k_converter.py +++ b/python/converters/wien2k_converter.py @@ -258,18 +258,17 @@ class Wien2kConverter(ConverterTools): # Reading done! # Save it to the HDF: - 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', + with HDFArchive(self.hdf_file, 'a') as ar: + 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 + for it in things_to_save: + ar[self.dft_subgrp][it] = locals()[it] # Symmetries are used, so now convert symmetry information for # *correlated* orbitals: @@ -292,15 +291,14 @@ class Wien2kConverter(ConverterTools): return # get needed data from hdf file - ar = HDFArchive(self.hdf_file, 'a') - things_to_read = ['SP', 'SO', 'n_shells', + with HDFArchive(self.hdf_file, 'a') as ar: + things_to_read = ['SP', 'SO', 'n_shells', 'n_k', 'n_orbitals', 'shells'] - for it in things_to_read: - if not hasattr(self, it): - setattr(self, it, ar[self.dft_subgrp][it]) - self.n_spin_blocs = self.SP + 1 - self.SO - del ar + for it in things_to_read: + if not hasattr(self, it): + setattr(self, it, ar[self.dft_subgrp][it]) + self.n_spin_blocs = self.SP + 1 - self.SO mpi.report("Reading input from %s..." % self.parproj_file) @@ -368,16 +366,15 @@ class Wien2kConverter(ConverterTools): # Reading done! # Save it to the HDF: - ar = HDFArchive(self.hdf_file, 'a') - if not (self.parproj_subgrp in ar): - ar.create_group(self.parproj_subgrp) - # The subgroup containing the data. If it does not exist, it is - # created. If it exists, the data is overwritten! - things_to_save = ['dens_mat_below', 'n_parproj', + with HDFArchive(self.hdf_file, 'a') as ar: + if not (self.parproj_subgrp in ar): + ar.create_group(self.parproj_subgrp) + # The subgroup containing the data. If it does not exist, it is + # created. If it exists, the data is overwritten! + things_to_save = ['dens_mat_below', 'n_parproj', 'proj_mat_all', 'rot_mat_all', 'rot_mat_all_time_inv'] - for it in things_to_save: - ar[self.parproj_subgrp][it] = locals()[it] - del ar + for it in things_to_save: + ar[self.parproj_subgrp][it] = locals()[it] # Symmetries are used, so now convert symmetry information for *all* # orbitals: @@ -395,15 +392,14 @@ class Wien2kConverter(ConverterTools): try: # get needed data from hdf file - ar = HDFArchive(self.hdf_file, 'a') - things_to_read = ['SP', 'SO', 'n_corr_shells', + with HDFArchive(self.hdf_file, 'a') as ar: + things_to_read = ['SP', 'SO', 'n_corr_shells', 'n_shells', 'corr_shells', 'shells', 'energy_unit'] - for it in things_to_read: - if not hasattr(self, it): - setattr(self, it, ar[self.dft_subgrp][it]) - self.n_spin_blocs = self.SP + 1 - self.SO - del ar + for it in things_to_read: + if not hasattr(self, it): + setattr(self, it, ar[self.dft_subgrp][it]) + self.n_spin_blocs = self.SP + 1 - self.SO mpi.report("Reading input from %s..." % self.band_file) R = ConverterTools.read_fortran_file( @@ -482,16 +478,15 @@ class Wien2kConverter(ConverterTools): # Reading done! # Save it to the HDF: - ar = HDFArchive(self.hdf_file, 'a') - if not (self.bands_subgrp in ar): - ar.create_group(self.bands_subgrp) - # The subgroup containing the data. If it does not exist, it is - # created. If it exists, the data is overwritten! - things_to_save = ['n_k', 'n_orbitals', 'proj_mat', + with HDFArchive(self.hdf_file, 'a') as ar: + if not (self.bands_subgrp in ar): + ar.create_group(self.bands_subgrp) + # The subgroup containing the data. If it does not exist, it is + # created. If it exists, the data is overwritten! + things_to_save = ['n_k', 'n_orbitals', 'proj_mat', 'hopping', 'n_parproj', 'proj_mat_all'] - for it in things_to_save: - ar[self.bands_subgrp][it] = locals()[it] - del ar + for it in things_to_save: + ar[self.bands_subgrp][it] = locals()[it] def convert_misc_input(self): """ @@ -510,13 +505,12 @@ class Wien2kConverter(ConverterTools): return # Check if SP, SO and n_k are already in h5 - ar = HDFArchive(self.hdf_file, 'r') - if not (self.dft_subgrp in ar): - raise IOError, "convert_misc_input: No %s subgroup in hdf file found! Call convert_dft_input first." % self.dft_subgrp - SP = ar[self.dft_subgrp]['SP'] - SO = ar[self.dft_subgrp]['SO'] - n_k = ar[self.dft_subgrp]['n_k'] - del ar + with HDFArchive(self.hdf_file, 'r') as ar: + if not (self.dft_subgrp in ar): + raise IOError, "convert_misc_input: No %s subgroup in hdf file found! Call convert_dft_input first." % self.dft_subgrp + SP = ar[self.dft_subgrp]['SP'] + SO = ar[self.dft_subgrp]['SO'] + n_k = ar[self.dft_subgrp]['n_k'] things_to_save = [] @@ -612,12 +606,11 @@ class Wien2kConverter(ConverterTools): raise IOError, "convert_misc_input: reading file %s failed" % self.outputs_file # Save it to the HDF: - ar = HDFArchive(self.hdf_file, 'a') - if not (self.misc_subgrp in ar): - ar.create_group(self.misc_subgrp) - for it in things_to_save: - ar[self.misc_subgrp][it] = locals()[it] - del ar + with HDFArchive(self.hdf_file, 'a') as ar: + if not (self.misc_subgrp in ar): + ar.create_group(self.misc_subgrp) + for it in things_to_save: + ar[self.misc_subgrp][it] = locals()[it] def convert_transport_input(self): """ @@ -633,13 +626,12 @@ class Wien2kConverter(ConverterTools): return # Check if SP, SO and n_k are already in h5 - ar = HDFArchive(self.hdf_file, 'r') - if not (self.dft_subgrp in ar): - raise IOError, "convert_transport_input: No %s subgroup in hdf file found! Call convert_dft_input first." % self.dft_subgrp - SP = ar[self.dft_subgrp]['SP'] - SO = ar[self.dft_subgrp]['SO'] - n_k = ar[self.dft_subgrp]['n_k'] - del ar + with HDFArchive(self.hdf_file, 'r') as ar: + if not (self.dft_subgrp in ar): + raise IOError, "convert_transport_input: No %s subgroup in hdf file found! Call convert_dft_input first." % self.dft_subgrp + SP = ar[self.dft_subgrp]['SP'] + SO = ar[self.dft_subgrp]['SO'] + n_k = ar[self.dft_subgrp]['n_k'] # Read relevant data from .pmat/up/dn files ########################################### @@ -691,15 +683,14 @@ class Wien2kConverter(ConverterTools): R.close() # Reading done! # Put data to HDF5 file - ar = HDFArchive(self.hdf_file, 'a') - if not (self.transp_subgrp in ar): - ar.create_group(self.transp_subgrp) - # The subgroup containing the data. If it does not exist, it is - # created. If it exists, the data is overwritten!!! - things_to_save = ['band_window_optics', 'velocities_k'] - for it in things_to_save: - ar[self.transp_subgrp][it] = locals()[it] - del ar + with HDFArchive(self.hdf_file, 'a') as ar: + if not (self.transp_subgrp in ar): + ar.create_group(self.transp_subgrp) + # The subgroup containing the data. If it does not exist, it is + # created. If it exists, the data is overwritten!!! + things_to_save = ['band_window_optics', 'velocities_k'] + for it in things_to_save: + ar[self.transp_subgrp][it] = locals()[it] def convert_symmetry_input(self, orbits, symm_file, symm_subgrp, SO, SP): """ @@ -781,11 +772,10 @@ class Wien2kConverter(ConverterTools): # Reading done! # Save it to the HDF: - ar = HDFArchive(self.hdf_file, 'a') - if not (symm_subgrp in ar): - ar.create_group(symm_subgrp) - things_to_save = ['n_symm', 'n_atoms', 'perm', + with HDFArchive(self.hdf_file, 'a') as ar: + if not (symm_subgrp in ar): + ar.create_group(symm_subgrp) + things_to_save = ['n_symm', 'n_atoms', 'perm', 'orbits', 'SO', 'SP', 'time_inv', 'mat', 'mat_tinv'] - for it in things_to_save: - ar[symm_subgrp][it] = locals()[it] - del ar + for it in things_to_save: + ar[symm_subgrp][it] = locals()[it] diff --git a/python/sumk_dft.py b/python/sumk_dft.py index 825b6b5d..e358be3e 100644 --- a/python/sumk_dft.py +++ b/python/sumk_dft.py @@ -187,23 +187,22 @@ class SumkDFT(object): subgroup_present = 0 if mpi.is_master_node(): - ar = HDFArchive(self.hdf_file, 'r') - if subgrp in ar: - subgroup_present = True - # first read the necessary things: - for it in things_to_read: - if it in ar[subgrp]: - setattr(self, it, ar[subgrp][it]) - else: - mpi.report("Loading %s failed!" % it) - value_read = False - else: - if (len(things_to_read) != 0): - mpi.report( - "Loading failed: No %s subgroup in hdf5!" % subgrp) - subgroup_present = False - value_read = False - del ar + with HDFArchive(self.hdf_file, 'r') as ar: + if subgrp in ar: + subgroup_present = True + # first read the necessary things: + for it in things_to_read: + if it in ar[subgrp]: + setattr(self, it, ar[subgrp][it]) + else: + mpi.report("Loading %s failed!" % it) + value_read = False + else: + if (len(things_to_read) != 0): + mpi.report( + "Loading failed: No %s subgroup in hdf5!" % subgrp) + subgroup_present = False + value_read = False # now do the broadcasting: for it in things_to_read: setattr(self, it, mpi.bcast(getattr(self, it))) @@ -226,18 +225,16 @@ class SumkDFT(object): if not (mpi.is_master_node()): return # do nothing on nodes - ar = HDFArchive(self.hdf_file, 'a') - if not subgrp in ar: - ar.create_group(subgrp) - for it in things_to_save: - if it in [ "gf_struct_sumk", "gf_struct_solver", - "solver_to_sumk", "sumk_to_solver", "solver_to_sumk_block"]: - warn("It is not recommended to save '{}' individually. Save 'block_structure' instead.".format(it)) - try: - ar[subgrp][it] = getattr(self, it) - except: - mpi.report("%s not found, and so not saved." % it) - del ar + with HDFArchive(self.hdf_file, 'a') as ar: + if not subgrp in ar: ar.create_group(subgrp) + for it in things_to_save: + if it in [ "gf_struct_sumk", "gf_struct_solver", + "solver_to_sumk", "sumk_to_solver", "solver_to_sumk_block"]: + warn("It is not recommended to save '{}' individually. Save 'block_structure' instead.".format(it)) + try: + ar[subgrp][it] = getattr(self, it) + except: + mpi.report("%s not found, and so not saved." % it) def load(self, things_to_load, subgrp='user_data'): r""" @@ -258,16 +255,15 @@ class SumkDFT(object): if not (mpi.is_master_node()): return # do nothing on nodes - ar = HDFArchive(self.hdf_file, 'r') - if not subgrp in ar: - mpi.report("Loading %s failed!" % subgrp) - list_to_return = [] - for it in things_to_load: - try: - list_to_return.append(ar[subgrp][it]) - except: - raise ValueError, "load: %s not found, and so not loaded." % it - del ar + with HDFArchive(self.hdf_file, 'r') as ar: + if not subgrp in ar: + mpi.report("Loading %s failed!" % subgrp) + list_to_return = [] + for it in things_to_load: + try: + list_to_return.append(ar[subgrp][it]) + except: + raise ValueError, "load: %s not found, and so not loaded." % it return list_to_return ################ @@ -1822,10 +1818,9 @@ class SumkDFT(object): fermi_weights = 0 band_window = 0 if mpi.is_master_node(): - ar = HDFArchive(self.hdf_file,'r') - fermi_weights = ar['dft_misc_input']['dft_fermi_weights'] - band_window = ar['dft_misc_input']['band_window'] - del ar + with HDFArchive(self.hdf_file,'r') as ar: + fermi_weights = ar['dft_misc_input']['dft_fermi_weights'] + band_window = ar['dft_misc_input']['band_window'] fermi_weights = mpi.bcast(fermi_weights) band_window = mpi.bcast(band_window) diff --git a/python/symmetry.py b/python/symmetry.py index 461587e8..6e1ee825 100644 --- a/python/symmetry.py +++ b/python/symmetry.py @@ -58,16 +58,15 @@ class Symmetry: if mpi.is_master_node(): # Read the stuff on master: - ar = HDFArchive(hdf_file, 'r') - if subgroup is None: - ar2 = ar - else: - ar2 = ar[subgroup] + with HDFArchive(hdf_file, 'r') as ar: + if subgroup is None: + ar2 = ar + else: + ar2 = ar[subgroup] - for it in things_to_read: - setattr(self, it, ar2[it]) + for it in things_to_read: + setattr(self, it, ar2[it]) del ar2 - del ar # Broadcasting for it in things_to_read: diff --git a/test/plovasp/converter/test_converter_lunio3.py b/test/plovasp/converter/test_converter_lunio3.py index 95a87d50..1d9a2eb0 100644 --- a/test/plovasp/converter/test_converter_lunio3.py +++ b/test/plovasp/converter/test_converter_lunio3.py @@ -3,8 +3,8 @@ import os import rpath _rpath = os.path.dirname(rpath.__file__) + '/' -from pytriqs.applications.dft.converters.plovasp.converter import generate_and_output_as_text -from pytriqs.applications.dft.converters import VaspConverter +from triqs_dft_tools.converters.plovasp.converter import generate_and_output_as_text +from triqs_dft_tools.converters import VaspConverter import mytest ################################################################################ diff --git a/test/plovasp/converter/test_converter_one_site.py b/test/plovasp/converter/test_converter_one_site.py index 3b5cfc32..a53b9bce 100644 --- a/test/plovasp/converter/test_converter_one_site.py +++ b/test/plovasp/converter/test_converter_one_site.py @@ -3,8 +3,8 @@ import os import rpath _rpath = os.path.dirname(rpath.__file__) + '/' -from pytriqs.applications.dft.converters.plovasp.converter import generate_and_output_as_text -from pytriqs.applications.dft.converters import VaspConverter +from triqs_dft_tools.converters.plovasp.converter import generate_and_output_as_text +from triqs_dft_tools.converters import VaspConverter import mytest ################################################################################ @@ -29,7 +29,7 @@ class TestConverterOneSite(mytest.MyTestCase): generate_and_output_as_text(_rpath + 'example.cfg', _rpath + 'one_site/') test_file = _rpath + 'pg_output.test.h5' - converter = VaspConverter(filename=_rpath + 'one_site', + converter = VaspConverter(filename=_rpath + 'one_site', hdf_filename=test_file) converter.convert_dft_input() diff --git a/test/srvo3_transp.py b/test/srvo3_transp.py index 91b612b2..4d37523b 100644 --- a/test/srvo3_transp.py +++ b/test/srvo3_transp.py @@ -34,12 +34,11 @@ Converter.convert_transport_input() SK = SumkDFTTools(hdf_file='SrVO3.h5', use_dft_blocks=True) -ar = HDFArchive('SrVO3_Sigma.h5', 'a') -Sigma = ar['dmft_transp_input']['Sigma_w'] -SK.set_Sigma([Sigma]) -SK.chemical_potential = ar['dmft_transp_input']['chemical_potential'] -SK.dc_imp = ar['dmft_transp_input']['dc_imp'] -del ar +with HDFArchive('SrVO3_Sigma.h5', 'a') as ar: + Sigma = ar['dmft_transp_input']['Sigma_w'] + SK.set_Sigma([Sigma]) + SK.chemical_potential = ar['dmft_transp_input']['chemical_potential'] + SK.dc_imp = ar['dmft_transp_input']['dc_imp'] SK.transport_distribution(directions=['xx'], broadening=0.0, energy_window=[-0.3,0.3], Om_mesh=[0.00, 0.02] , beta=beta, with_Sigma=True) #SK.save(['Gamma_w','Om_meshr','omega','directions'])