LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-12-23 04:43:42 +01:00
Code Issues Releases Wiki Activity

Merge branch 'merge_vasp2dmft' into unstable

Browse Source
This commit is contained in:
aichhorn 2020-02-17 15:19:58 +01:00
commit b647762349
123 changed files with 4118 additions and 27357 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
doc/CMakeLists.txt
View File

@ -18,6 +18,6 @@ add_custom_target(doc_sphinx ALL DEPENDS ${sphinx_top} ${CMAKE_CURRENT_BINARY_DI
# ---------------------------------
install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/html/ COMPONENT documentation DESTINATION share/doc/triqs_dft_tools
FILES_MATCHING
REGEX "\\.(html|pdf|png|gif|jpg|js|xsl|css|py|txt|inv|bib)$"
REGEX "\\.(html|pdf|png|gif|jpg|js|xsl|css|py|txt|inv|bib|cfg)$"
PATTERN "_*"
)

6
doc/conf.py.in
View File

@ -14,7 +14,11 @@ extensions = ['sphinx.ext.autodoc',
'sphinx.ext.todo',
'sphinx.ext.viewcode',
'sphinx.ext.autosummary',
'numpydoc']
'nbsphinx',
'numpydoc',
'sphinx.ext.githubpages',
'IPython.sphinxext.ipython_console_highlighting'
]
source_suffix = '.rst'

437
doc/guide/conv_vasp.rst
View File

@ -1,50 +1,67 @@
.. _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.
The VASP interface relies on new options introduced since version 5.4.x In
particular, a new INCAR-option `LOCPROJ
<https://cms.mpi.univie.ac.at/wiki/index.php/LOCPROJ>`_, the new `LORBIT` modes
13 and 14 have been added, and the new `ICHARG` mode 5 for charge
self-consistent DFT+DMFT calculations have been added.
*Limitations of the alpha-version:*
The VASP interface methodologically builds on the so called projection on
localized orbitals (PLO) scheme, where the resulting KS states from DFT are
projected on localized orbitals, which defines a basis for setting up a
Hubbard-like model Hamiltonian. Resulting in lattice object stored in `SumkDFT`.
The implementation is presented in `M. Schüler et al. 2018 J. Phys.: Condens.
Matter 30 475901 <https://doi.org/10.1088/1361-648X/aae80a>`_.
* The interface works correctly only if the k-point symmetries
are turned off during the VASP run (ISYM=-1).
The interface consists of two parts, :ref:`PLOVASP<refPLOVASP>`, a collection of
python classes and functions converting the raw VASP output to proper projector
functions, and the python based :ref:`VaspConverter<refVASPconverter>`, which
creates a h5 archive from the :ref:`PLOVASP<refPLOVASP>` output readable by
`SumkDFT`. Therefore, the conversion consist always of two steps.
* Generation of projectors for k-point lines (option `Lines` in KPOINTS)
needed for Bloch spectral function calculations is not possible at the moment.
Here, we will present a guide how the interface `can` be used to create input for a DMFT calculation, using SrVO3 as an example. Full examples can be found in the :ref:`tutorial section of DFTTools<tutorials>`.
* The interface currently supports only collinear-magnetism calculation
(this implis no spin-orbit coupling) and
spin-polarized projectors have not been tested.
Limitations of the interface
============================
A detailed description of the VASP converter tool PLOVasp can be found
in the :ref:`PLOVasp User's Guide <plovasp>`. Here, a quick-start guide is presented.
* 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 implies no spin-orbit coupling) and spin-polarized projectors have not
been tested.
* The converter needs the correct Fermi energy from VASP, which is read from
the LOCPROJ file. However, VASP by default does not output this information.
Please see `Remarks on the VASP version`_.
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.
VASP: generating raw projectors
===============================
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:
The VASP **INCAR** option `LOCPROJ` selects a set of localized projectors that
will be written to the 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 = <sites> : <shells> : <projector type>`
where `<sites>` represents a list of site indices separated by spaces,
with the indices corresponding to the site position in the POSCAR file;
`<shells>` specifies local states (see below);
`<projector type>` 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 <https://doi.org/10.1088/1361-648X/aae80a>`_.
where `<sites>` represents a list of site indices separated by spaces, with the
indices corresponding to the site position in the **POSCAR** file; `<shells>`
specifies local states (see below); `<projector type>` chooses a particular type
of the local basis function. The recommended projector type is `Pr 2`. This will
perform a projection of the Kohn-Sham states onto the VASP PAW projector
functions. The number specified behind `Pr` is selecting a specific PAW channel,
see the `VASP wiki page <https://cms.mpi.univie.ac.at/wiki/index.php/LOCPROJ>`_
for more information. The formalism for this type of projectors is presented in
`M. Schüler et al. 2018 J. Phys.: Condens. Matter 30 475901
<https://doi.org/10.1088/1361-648X/aae80a>`_. For further details on the
`LOCPROJ` flag also have a look in the `VASP wiki
<https://cms.mpi.univie.ac.at/wiki/index.php/LOCPROJ>`_.
The allowed labels of the local states defined in terms of cubic
harmonics are:
harmonics are (mind the order):
* Entire shells: `s`, `p`, `d`, `f`
@ -55,16 +72,25 @@ harmonics are:
* `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 projector type `Pr`, one should ideally 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. Setting `LORBIT=14` will perform an automatic optimization of
the PAW projector channel as described in `M. Schüler et al. 2018 J. Phys.:
Condens. Matter 30 475901 <https://doi.org/10.1088/1361-648X/aae80a>`_, by
using a linear combination of the PAW channels, to maximize the overlap in
the chosen energy window between the projector and the Kohn-Sham state.
Therefore, setting `LORBIT=14` will let VASP ignore the channel specified
after `Pr`. This optimization is only performed for the projector type `Pr`,
not for `Ps` and obviously not for `Hy`. We recommend to specify the PAW
channel anyway, in case one forgets to set `LORBIT=14`.
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):
In case of SrVO3 one may first want to perform a self-consistent
calculation to know the Fermi level and the rough position of the target states.
In the next step one sets `ICHARG = 1` and adds the following additional lines
into INCAR (provided that V is the second ion in POSCAR):
| `EMIN = 3.0`
| `EMAX = 8.0`
@ -72,73 +98,320 @@ lines into INCAR (provided that V is the second ion in POSCAR):
| `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.
overlap with valence bands and no overlap with semi-core or high unoccupied
states. This **INCAR** will calculate and write-out projections for all five d-orbitals.
VASP input-output
-----------------
The calculated projections :math:`\langle \chi_L | \Psi_\mu \rangle` are written
into files **PROJCAR** and **LOCPROJ**. The difference between these two files
is that **LOCPROJ** contains raw matrices without any reference to
sites/orbitals, while **PROJCAR** is more detailed. In particular, the
information that can be obtained for each projector from **PROJCAR** is the
following:
* site (and species) index
* for each `k`-point and band: a set of complex numbers for labeled orbitals
At the same time, **LOCPROJ** contains the total number of projectors (as well
as the number of `k`-points, bands, and spin channels) in the first line, which
can be used to allocate the arrays before parsing.
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
The projectors generated by VASP require certain post-processing before they can
be used for DMFT calculations. The most important step is to (ortho-)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.
problem. This will create proper Wannier functions of the initial projections
produced by VASP. 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 the **INCAR** VASP input file. This part is done in `PLOVASP`.
Post-processing of `LOCPROJ` data is generally done as follows:
#. Prepare an input file `<name>.cfg` (e.g., `plo.cfg`) that describes the definition
of your impurity problem (more details below).
PLOVASP: converting VASP output
--------------------------------
#. Extract the value of the Fermi level from OUTCAR and paste it at the end of
the first line of LOCPROJ.
:ref:`PLOVASP<refPLOVASP>` is a collection of python functions and classes, post-processing the raw VASP `LOCPROJ` output creating proper projector functions.
#. Run :program:`plovasp` with the input file as an argument, e.g.:
The following VASP files are used by PLOVASP:
* PROJCAR, LOCPROJ: raw projectors generated by VASP-PLO interface
* EIGENVAL: Kohn-Sham eigenvalues as well as `k`-points with weights and Fermi weights
* IBZKPT: `k`-point data (:math:`\Gamma`)
* POSCAR: crystal structure data
| `plovasp plo.cfg`
To run `PLOVASP`, one first prepares an input file `<name>.cfg` (default name `plo.cfg`) that describes the definition of the correlated subspace. For SrVO3 this input file would look like this:
This requires that the TRIQS paths are set correctly (see Installation
of TRIQS).
.. literalinclude:: ../tutorials/svo_vasp/plo.cfg
If everything goes right one gets files `<name>.ctrl` and `<name>.pg1`.
These files are needed for the converter that will be invoked in your
DMFT script.
In the [section] general, the `DOSMESH` defines an energy window and number of
data points, which lets the converter calculate the density of states of the
created projector functions in a given energy window. Each 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
The format of input file `<name>.cfg` is described in details in
the :ref:`User's Guide <plovasp>`. 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`.
- **IONS**: list of site indices which must be a subset of indices given earlier
in the VASP INCAR `LOCPROJ` flag. Note: If projections are performed for
multiple sites one can specify symmetry equivalent sites with brackets: `[2
3]`. Here the projector are generated for ions 2 and 3, but they will be
marked as symmetry equivalent later in 'SumkDFT'.
- **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 Option **TRANSFORM** is optional here, and it is specified to extract
only the three :math:`t_{2g}` orbitals out of the five `d` orbitals given by
:math:`l = 2`. A detailed explanation of all input parameters can be found
further below `PLOVASP detailed guide`_.
The conversion to a h5-file is performed in the same way as for Wien2TRIQS::
Next, the converter is executed. This can be done by calling :program:`PLOVASP` directly in the command line with the input file as an argument, e.g.:
| `plovasp plo.cfg`
or embedded in a python script as::
import triqs_dft_tools.converters.plovasp.converter as plo_converter
# Generate and store PLOs
plo_converter.generate_and_output_as_text('plo.cfg', vasp_dir='./')
This will create the xml files `vasp.ctrl` and `vasp.pg1` containing the orthonormalized projector functions readable by the :ref:`VaspConverter<refVASPconverter>`. Moreover, `PLOVASP` will output important information of the orthonormalization process, such as the density matrix of the correlated shell and the local Hamiltonian.
Running the VASP converter
-------------------------------------
The actual conversion to a h5-file is performed with the orthonormalized projector functions readable by the :ref:`VaspConverter<refVASPconverter>` in the same fashion as with the other `DFTTools` converters::
from triqs_dft_tools.converters.vasp_converter import *
Converter = VaspConverter(filename = filename)
Converter = VaspConverter(filename = 'vasp')
Converter.convert_dft_input()
As usual, the resulting h5-file can then be used with the SumkDFT class.
As usual, the resulting h5-file can then be used with the SumkDFT class::
sk = SumkDFTTools(hdf_file='vasp.h5')
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 <dft.sumk_dft.SumkDFT>`, e.g.::
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 <dft.sumk_dft.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.
However, this should only be done after a careful study of the density matrix and the projected DOS in the localized basis. For the complete process for SrVO3 see the tutorial for the VASP interface `here <../tutorials/svo_vasp/svo_notebook.html>`_.
PLOVASP detailed guide
======================
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). To enhance the performance
parsing the raw VASP output files, the parser is implemented in plain C. The
idea is that the python part of the parser first reads the first line of
**LOCPROJ** and then calls the C-routine with necessary parameters to parse
**PROJCAR**. The resulting 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).
In PLOVASP, all projectors sharing the same energy window are combined into a
`projector group`. This allows one in principal to define several groups with
different energy windows for the same set of raw projectors. Note: multiple groups are not yet implemented.
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. Additionally, each shell can be marked correlated or non-correlated, to tell `SumkDFT` whether or not these should be treated in the DMFT impurity problem.
Input file format
-----------------
The input file is written in the standard config-file format.
Parameters (or 'options') are grouped into sections specified as
`[Section name]`. All parameters must be defined inside some section.
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
output (e.g. the density of states), etc.
#. **[Group <Ng>]**: describes projector groups, i.e. a set of
projectors sharing the same energy window and normalization type.
At the moment, DFTtools support only one projector group, therefore
there should be no more than one projector group.
#. **[Shell <Ns>]**: contains parameters of a projector shell labelled
with `<Ns>`. If there is only one group section and one 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]
"""""""""""""""""
The entire section is optional and it contains four parameters:
* **BASENAME** (string): provides a base name for output files.
Default filenames are :file:`vasp.*`.
* **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_<s>_<n>.dat`, where `s` is the
shell index and `n` the ion 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.
* **HK** (True/False): If True, the projectors are applied the the Kohn-Sham
eigenvalues which results in a Hamitlonian H(k) in orbital basis. The H(k)
is written for each group to a file :file:`Basename.hk<Ng>`. It is recommended
to also set `COMPLEMENT = True` (see below). Default is False.
There are no required parameters in this section.
Section [Shell]
"""""""""""""""
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.
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
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 the LOCPROJ file.
There are additional optional parameters that allow one to transform
the local states:
* **CORR** (True/False): Determines if shell is correlated or not. At least one
shell has to be correlated. Default is True.
* **SORT** (integer): Overrides the default detection of ion sorts by supplying
an integer. Default is `None`, for which the default behavior is retained.
* **TRANSFORM** (matrix): local transformation matrix applied to all states
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
orbitals or perform a simple global rotation.
* **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 :ref:`below <transformation_file>`.
Section [Group]
"""""""""""""""
Each defined projector shell must be part of a projector group. In the current
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:
* **SHELLS** (list of integers): indices of projector shells included in the group.
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. *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 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
condition will be enforced on each site separately but the 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. The default value is **True**
* **BANDS** (int int): the energy window specified by two ints: band index of
lowest band and band index of highest band. Using this overrides the selection
in `EWINDOW`.
* **COMPLEMENT** (True/False). If True, the orthogonal complement is calculated
resulting in unitary (quadratic) projectors, i.e., the same number of orbitals
as bands. It is required to have an equal number of bands in the energy window
at each k-point. Default is False.
.. _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 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 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 (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
Remarks on the VASP version
===============================
In the current version of the interface the Fermi energy is extracted from the
`LOCPROJ` file. The file should contain the Fermi energy in the header. One can
either copy the Fermi energy manually there after a successful VASP run, or
modify the VASP source code slightly, by replacing the following line in
`locproj.F` (around line 695):
::
WRITE(99,'(4I6," # of spin, # of k-points, # of bands, # of proj" )') NS,NK,NB,NF
with:
::
WRITE(99,'(4I6,F12.7," # of spin, # of k-points, # of bands, # of proj, Efermi" )') W%WDES%NCDIJ,NK,NB,NF,EFERMI
Another critical point for CSC calculations is the function call of
`LPRJ_LDApU` in VASP. This function is not needed, and was left there for debug
purposes, but is called every iteration. Removing the call to this function in `electron.F` in line 644 speeds up the calculation significantly in the `ICHARG=5` mode. Moreover, this prevents VASP from generating the `GAMMA` file, which should ideally only be done by the DMFT code after a successful DMFT step, and then be read by VASP.
Furthermore, there is a bug in `fileio.F` around line 1710 where VASP tries to
print "reading the density matrix from Gamma". This should be done only by the
master node, and VASP gets stuck sometimes. Adding a
::
IF (IO%IU0>=0) THEN
...
ENDIF
statement resolves this issue. A similar problem occurs, when VASP writes the
`OSZICAR` file and a buffer is stuck. Adding a `flush` to the buffer in
`electron.F` around line 580 after
::
CALL STOP_TIMING("G",IO%IU6,"DOS")
flush(17)
print *, ' '
resolves this issue. Otherwise the OSZICAR file is not written properly.

2
doc/guide/conversion.rst
View File

@ -13,7 +13,7 @@ as a light-weight general-purpose interface. In the following, we will describe
conversion tools.
.. toctree::
:maxdepth: 2
:maxdepth: 3
conv_wien2k
conv_vasp

68
doc/guide/dftdmft_selfcons.rst
View File

@ -124,28 +124,66 @@ example how such a self-consistent calculation is performed from scratch.
VASP + PLOVasp
--------------
.. warning::
This is a preliminary documentation valid for the alpha-version of the interface.
Modifications to the implementation might be introduced before the final release.
Unlike Wien2k implementation the charge self-consistent DMFT cycle in the
framework of PLOVasp interface is controlled by an external script. Because of
the specific way the DFT self-consistency is implemented in VASP the latter has
to run parallel to the DMFT script, with the synchronisation being ensured by a
lock file.
Unlike Wien2k implementation the charge self-consistent DMFT cycle
in the framework of PLOVasp interface is controlled by an external script.
Because of the specific way the DFT self-consistency is implemented in VASP
the latter has to run parallel to the DMFT script, with the synchronisation being
ensured by a lock file. PLOVasp interface provides a shell-script :program:`vasp_dmft.sh`
which takes care of the process management. The user must, however, specify a path
to VASP code and provide the DMFT Python-script.
Once VASP reaches the point where the projectors are generated
it creates a lock file `vasp.lock` and waits until the lock file is
removed. The shell script, in turn, waits for the VASP process and once
the lock file is created it starts a DMFT iteration. The DMFT iteration
must finish by generating a Kohn-Sham (KS) density matrix (file `GAMMA`)
and removing the lock file. The VASP process then reads in `GAMMA`
and proceeds with the next iteration. PLOVasp interface provides a shell-script :program:`vasp_dmft` (in the triqs bin directory):
::
vasp_dmft [-n <number of cores>] -i <number of iterations> -j <number of VASP iterations with fixed charge density> [-v <VASP version>] [-p <path to VASP directory>] [<dmft_script.py>]
If the number of cores is not specified it is set to 1 by default.
Set the number of times the dmft solver is called with -i <number of iterations>
Set the number of VASP iteration with a fixed charge density update
inbetween the dmft runs with -j <number of VASP iterations with fixed charge density>
Set the version of VASP by -v standard(default)/no_gamma_write to
specify if VASP writes the GAMMA file or not.
If the path to VASP directory is not specified it must be provided by a
variable VASP_DIR.
<dmft_script.py> must provide an importable function 'dmft_cycle()'
which is invoked once per DFT+DMFT iteration. If the script name is
omitted the default name 'csc_dmft.py' is used.
which takes care of the process management. The user must, however, specify a path to VASP code and provide the DMFT Python-script. See for an example :ref:`NiO CSC tutorial<nio_csc>`.
The user-provided script is almost the same as for Wien2k charge self-consistent
calculations with the main difference that its functionality (apart from
the lines importing other modules) should be placed inside a function `dmft_cycle()`
which will be called every DMFT cycle. Another difference is the way
function `calc_density_correction()` works.
calculations with the main difference that its functionality (apart from the
lines importing other modules) should be placed inside a function `dmft_cycle()`
which will be called every DMFT cycle.
VASP has a special INCAR `ICHARG=5` mode, that has to be switched on to make VASP wait for the `vasp.lock` file, and read the updated charge density after each step. One should add the following lines to the `INCAR` file::
ICHARG = 5
NELM = 1000
NELMIN = 1000
Technically, VASP runs with `ICHARG=5` in a SCF mode, and adding the DMFT
changes to the DFT density in each step, so that the full DFT+DMFT charge
density is constructed in every step. This is only done in VASP because only the
changes to the DFT density are read by VASP not the full DFT+DMFT density.
Moreover, one should always start with a converged `WAVECAR` file, or make sure,
that the KS states are well converged before the first projectors are created!
To understand the difference please make sure to read `ISTART flag VASP wiki
<https://www.vasp.at/wiki/index.php/ISTART>`_. Furthermore, the flags `NELM` and
`NELMIN` ensure that VASP does not terminate after the default number of
iterations of 60.
Other DFT codes
---------------
The extension to other DFT codes is straight forward. As described
The extension to other DFT codes is straightforward. As described
here, one needs to implement the correlated density matrix to be used
for the calculation of the charge density. This implementation will of
course depend on the DFT package, and might be easy to do or a quite

167
doc/guide/plovasp.rst
View File

@ -1,167 +0,0 @@
.. _plovasp:
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.
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.
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.
Input file format
-----------------
The input file is written in the standard config-file format.
Parameters (or 'options') are grouped into sections specified as
`[Section name]`. All parameters must be defined inside some section.
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
output (e.g. the density of states), etc.
#. **[Group <Ng>]**: describes projector groups, i.e. a set of
projectors sharing the same energy window and normalization type.
At the moment, DFTtools support only one projector group, therefore
there should be no more than one projector group.
#. **[Shell <Ns>]**: contains parameters of a projector shell labelled
with `<Ns>`. If there is only one group section and one 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]
"""""""""""""""""
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,
the projected density of states for each projected orbital will be
evaluated and stored to files :file:`pdos_<n>.dat`, where `n` is the
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
labeled by an index, e.g. `[Shell 1]`. These indices can then be referenced
in a `[Group]` section.
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
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 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 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
orbitals or perform a simple global rotation.
* **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 :ref:`below <transformation_file>`.
Section [Group]
"""""""""""""""
Each defined projector shell must be part of a projector group. In the current
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:
* **SHELLS** (list of integers): indices of projector shells included in the group.
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. *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 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
condition will be enforced on each site separately but the 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:
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 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 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 (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

36
doc/reference/converters.rst
View File

@ -22,6 +22,42 @@ Wannier90 Converter
:members:
:special-members:
PLOVASP
----------
.. _refPLOVASP:
PLOVASP reference, the classes / functions are sorted the way the converter uses them.
.. automodule:: triqs_dft_tools.converters.plovasp.converter
:members: generate_and_output_as_text
.. automodule:: triqs_dft_tools.converters.plovasp.inpconf
:members: ConfigParameters
.. automodule:: triqs_dft_tools.converters.plovasp.vaspio
:members: VaspData, Plocar, Poscar, Kpoints, Eigenval, Doscar, read_symmcar
.. automodule:: triqs_dft_tools.converters.plovasp.elstruct
:members: ElectronicStructure
.. automodule:: triqs_dft_tools.converters.plovasp.plotools
:members:
.. automodule:: triqs_dft_tools.converters.plovasp.proj_shell
:members:
.. automodule:: triqs_dft_tools.converters.plovasp.proj_group
:members:
VASP Converter
-------------------
.. _refVASPconverter:
.. autoclass:: triqs_dft_tools.converters.vasp_converter.VaspConverter
:members:
:special-members:
Converter Tools
---------------
.. autoclass:: triqs_dft_tools.converters.converter_tools.ConverterTools

19
doc/tutorials.rst
View File

@ -23,3 +23,22 @@ Full charge self consistency with Wien2k: :math:`\gamma`-Ce
tutorials/ce-gamma-fscs_wien2k
VASP interface examples
-----------------------
Simple Converter example: SrVO3
"""""""""""""""""""""""""""""""
.. toctree::
:maxdepth: 2
tutorials/svo_vasp/svo_notebook
Complex example: NiO
"""""""""""""""""""""""""""""
.. toctree::
:maxdepth: 2
tutorials/nio_csc

18
doc/tutorials/images_scripts/INCAR Normal file
View File

@ -0,0 +1,18 @@
System = NiO
ISMEAR = -5
# the energy window to optimize projector channels
EMIN = -3
EMAX = 7
NBANDS = 16
LMAXMIX = 6
# switch off all symmetries
ISYM = -1
# project to Ni d and O p states
LORBIT = 14
LOCPROJ = 1 : d : Pr 1
LOCPROJ = 2 : p : Pr 1

7
doc/tutorials/images_scripts/INCAR.rst Normal file
View File

@ -0,0 +1,7 @@
.. _INCAR:
INCAR
-----------
.. literalinclude:: INCAR
:language: bash

4
doc/tutorials/images_scripts/KPOINTS Normal file
View File

@ -0,0 +1,4 @@
Automatically generated mesh
0
Gamma
6 6 6

7
doc/tutorials/images_scripts/KPOINTS.rst Normal file
View File

@ -0,0 +1,7 @@
.. _KPOINTS:
KPOINTS
-----------
.. literalinclude:: KPOINTS
:language: bash

90
doc/tutorials/images_scripts/NiO_local_lattice_GF.py Normal file
View File

@ -0,0 +1,90 @@
from itertools import *
import numpy as np
import pytriqs.utility.mpi as mpi
from pytriqs.archive import *
from pytriqs.gf import *
from triqs_dft_tools.sumk_dft import *
from triqs_dft_tools.sumk_dft_tools import *
from pytriqs.operators.util.hamiltonians import *
from pytriqs.operators.util.U_matrix import *
from triqs_cthyb import *
import warnings
warnings.filterwarnings("ignore", category=FutureWarning)
filename = 'nio'
SK = SumkDFT(hdf_file = filename+'.h5', use_dft_blocks = False)
beta = 5.0
# We analyze the block structure of the Hamiltonian
Sigma = SK.block_structure.create_gf(beta=beta)
SK.put_Sigma([Sigma])
G = SK.extract_G_loc()
SK.analyse_block_structure_from_gf(G, threshold = 1e-3)
# Setup CTQMC Solver
n_orb = SK.corr_shells[0]['dim']
spin_names = ['up','down']
orb_names = [i for i in range(0,n_orb)]
# Print some information on the master node
iteration_offset = 0
Sigma_iw = 0
if mpi.is_master_node():
ar = HDFArchive(filename+'.h5','a')
if not 'DMFT_results' in ar: ar.create_group('DMFT_results')
if not 'Iterations' in ar['DMFT_results']: ar['DMFT_results'].create_group('Iterations')
if 'iteration_count' in ar['DMFT_results']:
iteration_offset = ar['DMFT_results']['iteration_count']+1
print('offset',iteration_offset)
Sigma_iw = ar['DMFT_results']['Iterations']['Sigma_it'+str(iteration_offset-1)]
SK.dc_imp = ar['DMFT_results']['Iterations']['dc_imp'+str(iteration_offset-1)]
SK.dc_energ = ar['DMFT_results']['Iterations']['dc_energ'+str(iteration_offset-1)]
SK.chemical_potential = ar['DMFT_results']['Iterations']['chemical_potential'+str(iteration_offset-1)].real
iteration_offset = mpi.bcast(iteration_offset)
Sigma_iw = mpi.bcast(Sigma_iw)
SK.dc_imp = mpi.bcast(SK.dc_imp)
SK.dc_energ = mpi.bcast(SK.dc_energ)
SK.chemical_potential = mpi.bcast(SK.chemical_potential)
SK.put_Sigma(Sigma_imp = [Sigma_iw])
ikarray = numpy.array(range(SK.n_k))
# set up the orbitally resolved local lattice greens function:
n_orbs = SK.proj_mat_csc.shape[2]
spn = SK.spin_block_names[SK.SO]
mesh = Sigma_iw.mesh
block_structure = [range(n_orbs) for sp in spn]
gf_struct = [(spn[isp], block_structure[isp])
for isp in range(SK.n_spin_blocks[SK.SO])]
block_ind_list = [block for block, inner in gf_struct]
glist = lambda: [GfImFreq(indices=inner, mesh=mesh)
for block, inner in gf_struct]
G_latt_orb = BlockGf(name_list=block_ind_list,
block_list=glist(), make_copies=False)
G_latt_orb.zero()
for ik in mpi.slice_array(ikarray):
G_latt_KS = SK.lattice_gf(ik=ik, beta=beta)
G_latt_KS *= SK.bz_weights[ik]
for bname, gf in G_latt_orb:
add_g_ik = gf.copy()
add_g_ik.zero()
add_g_ik << SK.downfold(ik, 0, bname, G_latt_KS[bname], gf, shells='csc', ir=None)
gf << gf + add_g_ik
G_latt_orb << mpi.all_reduce(
mpi.world, G_latt_orb, lambda x, y: x + y)
mpi.barrier()
if mpi.is_master_node():
ar['DMFT_results']['Iterations']['G_latt_orb_it'+str(iteration_offset-1)] = G_latt_orb
if mpi.is_master_node(): del ar

9
doc/tutorials/images_scripts/NiO_local_lattice_GF.py.rst Normal file
View File

@ -0,0 +1,9 @@
.. _NiO_local_lattice_GF.py:
NiO_local_lattice_GF.py
-----------
Download :download:`NiO_local_lattice_GF.py <./NiO_local_lattice_GF.py>`.
.. literalinclude:: NiO_local_lattice_GF.py
:language: python

10
doc/tutorials/images_scripts/POSCAR Normal file
View File

@ -0,0 +1,10 @@
NiO
4.17 #taken from 9x9x9 with sigma=0.2 ismear=2
+0.0000000000 +0.5000000000 +0.5000000000
+0.5000000000 +0.0000000000 +0.5000000000
+0.5000000000 +0.5000000000 +0.0000000000
Ni O
1 1
Direct
+0.0000000000 +0.0000000000 +0.0000000000
+0.5000000000 +0.5000000000 +0.5000000000

7
doc/tutorials/images_scripts/POSCAR.rst Normal file
View File

@ -0,0 +1,7 @@
.. _POSCAR:
POSCAR
-----------
.. literalinclude:: POSCAR
:language: bash

3
doc/tutorials/images_scripts/converter.py Normal file
View File

@ -0,0 +1,3 @@
from triqs_dft_tools.converters.vasp_converter import *
Converter = VaspConverter(filename = 'nio')
Converter.convert_dft_input()

10
doc/tutorials/images_scripts/converter.py.rst Normal file
View File

@ -0,0 +1,10 @@
.. _converter.py:
converter.py
-------------------
Download :download:`converter.py <./converter.py>`.
.. literalinclude:: converter.py
:language: python

51
doc/tutorials/images_scripts/maxent.py Normal file
View File

@ -0,0 +1,51 @@
from pytriqs.gf import *
from pytriqs.archive import *
from triqs_maxent import *
filename = 'nio'
ar = HDFArchive(filename+'.h5','a')
if 'iteration_count' in ar['DMFT_results']:
iteration_offset = ar['DMFT_results']['iteration_count']+1
G_latt = ar['DMFT_results']['Iterations']['G_latt_orb_it'+str(iteration_offset-1)]
tm = TauMaxEnt(cost_function='bryan', probability='normal')
print(G_latt['up'][0,0])
t2g_orbs = [0,1,3]
eg_orbs = [2,4]
op_orbs = [5,6,7]
orbs = [t2g_orbs, eg_orbs, op_orbs]
#orbs = [t2g_orbs]
for orb in orbs:
print '\n'+str(orb[0])+'\n'
gf = 0*G_latt['up'][0,0]
for iO in orb:
gf = gf + G_latt['up'][iO,iO]
tm.set_G_iw(gf)
tm.omega =LinearOmegaMesh(omega_min=-20, omega_max=20, n_points=201)
tm.alpha_mesh = LogAlphaMesh(alpha_min=0.01, alpha_max=20000, n_points=60)
tm.set_error(1.e-3)
result=tm.run()
result.get_A_out('LineFitAnalyzer')
if 'iteration_count' in ar['DMFT_results']:
iteration_offset = ar['DMFT_results']['iteration_count']+1
for oo in orb:
ar['DMFT_results']['Iterations']['G_latt_orb_w_o'+str(oo)+'_it'+str(iteration_offset-1)] = result.analyzer_results['LineFitAnalyzer']['A_out']
ar['DMFT_results']['Iterations']['w_it'+str(iteration_offset-1)] = result.omega
# you may be interested in the details of the line analyzer:
# from pytriqs.plot.mpl_interface import oplot
#plt.figure(2)
#result.analyzer_results['LineFitAnalyzer'].plot_linefit()
#plt.savefig('ana'+str(orb[0])+'.pdf',fmt='pdf')
del ar

9
doc/tutorials/images_scripts/maxent.py.rst Normal file
View File

@ -0,0 +1,9 @@
.. _maxent.py:
maxent.py
-----------
Download :download:`maxent.py <./maxent.py>`.
.. literalinclude:: maxent.py
:language: python

167
doc/tutorials/images_scripts/nio.py Normal file
View File

@ -0,0 +1,167 @@
from itertools import *
import numpy as np
import pytriqs.utility.mpi as mpi
from pytriqs.archive import *
from pytriqs.gf import *
import sys, pytriqs.version as triqs_version
from triqs_dft_tools.sumk_dft import *
from triqs_dft_tools.sumk_dft_tools import *
from pytriqs.operators.util.hamiltonians import *
from pytriqs.operators.util.U_matrix import *
from triqs_cthyb import *
import triqs_cthyb.version as cthyb_version
import triqs_dft_tools.version as dft_tools_version
import warnings
warnings.filterwarnings("ignore", category=FutureWarning)
filename = 'nio'
SK = SumkDFT(hdf_file = filename+'.h5', use_dft_blocks = False)
beta = 5.0
Sigma = SK.block_structure.create_gf(beta=beta)
SK.put_Sigma([Sigma])
G = SK.extract_G_loc()
SK.analyse_block_structure_from_gf(G, threshold = 1e-3)
for i_sh in range(len(SK.deg_shells)):
num_block_deg_orbs = len(SK.deg_shells[i_sh])
mpi.report('found {0:d} blocks of degenerate orbitals in shell {1:d}'.format(num_block_deg_orbs, i_sh))
for iblock in range(num_block_deg_orbs):
mpi.report('block {0:d} consists of orbitals:'.format(iblock))
for keys in SK.deg_shells[i_sh][iblock].keys():
mpi.report(' '+keys)
# Setup CTQMC Solver
n_orb = SK.corr_shells[0]['dim']
spin_names = ['up','down']
orb_names = [i for i in range(0,n_orb)]
#gf_struct = set_operator_structure(spin_names, orb_names, orb_hyb)
gf_struct = SK.gf_struct_solver[0]
mpi.report('Sumk to Solver: %s'%SK.sumk_to_solver)
mpi.report('GF struct sumk: %s'%SK.gf_struct_sumk)
mpi.report('GF struct solver: %s'%SK.gf_struct_solver)
S = Solver(beta=beta, gf_struct=gf_struct)
# Construct the Hamiltonian and save it in Hamiltonian_store.txt
H = Operator()
U = 8.0
J = 1.0
U_sph = U_matrix(l=2, U_int=U, J_hund=J)
U_cubic = transform_U_matrix(U_sph, spherical_to_cubic(l=2, convention=''))
Umat, Upmat = reduce_4index_to_2index(U_cubic)
H = h_int_density(spin_names, orb_names, map_operator_structure=SK.sumk_to_solver[0], U=Umat, Uprime=Upmat)
# Print some information on the master node
mpi.report('Greens function structure is: %s '%gf_struct)
mpi.report('U Matrix set to:\n%s'%Umat)
mpi.report('Up Matrix set to:\n%s'%Upmat)
# Parameters for the CTQMC Solver
p = {}
p["max_time"] = -1
p["random_name"] = ""
p["random_seed"] = 123 * mpi.rank + 567
p["length_cycle"] = 100
p["n_warmup_cycles"] = 20000
p["n_cycles"] = 200000
p["fit_max_moment"] = 4
p["fit_min_n"] = 30
p["fit_max_n"] = 50
p["perform_tail_fit"] = True
# Double Counting: 0 FLL, 1 Held, 2 AMF
DC_type = 0
DC_value = 59.0
# Prepare hdf file and and check for previous iterations
n_iterations = 10
iteration_offset = 0
if mpi.is_master_node():
ar = HDFArchive(filename+'.h5','a')
if not 'DMFT_results' in ar: ar.create_group('DMFT_results')
if not 'Iterations' in ar['DMFT_results']: ar['DMFT_results'].create_group('Iterations')
if not 'DMFT_input' in ar: ar.create_group('DMFT_input')
if not 'Iterations' in ar['DMFT_input']: ar['DMFT_input'].create_group('Iterations')
if not 'code_versions' in ar['DMFT_input']: ar['DMFT_input'].create_group('code_versio\
ns')
ar['DMFT_input']['code_versions']["triqs_version"] = triqs_version.version
ar['DMFT_input']['code_versions']["triqs_git"] = triqs_version.git_hash
ar['DMFT_input']['code_versions']["cthyb_version"] = cthyb_version.version
ar['DMFT_input']['code_versions']["cthyb_git"] = cthyb_version.cthyb_hash
ar['DMFT_input']['code_versions']["dft_tools_version"] = dft_tools_version.version
ar['DMFT_input']['code_versions']["dft_tools_version"] = dft_tools_version.dft_tools_hash
if 'iteration_count' in ar['DMFT_results']:
iteration_offset = ar['DMFT_results']['iteration_count']+1
S.Sigma_iw = ar['DMFT_results']['Iterations']['Sigma_it'+str(iteration_offset-1)]
SK.dc_imp = ar['DMFT_results']['Iterations']['dc_imp'+str(iteration_offset-1)]
SK.dc_energ = ar['DMFT_results']['Iterations']['dc_energ'+str(iteration_offset-1)]
SK.chemical_potential = ar['DMFT_results']['Iterations']['chemical_potential'+str(iteration_offset-1)].real
ar['DMFT_input']["dmft_script_it"+str(iteration_offset)] = open(sys.argv[0]).read()
iteration_offset = mpi.bcast(iteration_offset)
S.Sigma_iw = mpi.bcast(S.Sigma_iw)
SK.dc_imp = mpi.bcast(SK.dc_imp)
SK.dc_energ = mpi.bcast(SK.dc_energ)
SK.chemical_potential = mpi.bcast(SK.chemical_potential)
# Calc the first G0
SK.symm_deg_gf(S.Sigma_iw,orb=0)
SK.put_Sigma(Sigma_imp = [S.Sigma_iw])
SK.calc_mu(precision=0.01)
S.G_iw << SK.extract_G_loc()[0]
SK.symm_deg_gf(S.G_iw, orb=0)
#Init the DC term and the self-energy if no previous iteration was found
if iteration_offset == 0:
dm = S.G_iw.density()
SK.calc_dc(dm, U_interact=U, J_hund=J, orb=0, use_dc_formula=DC_type,use_dc_value=DC_value)
S.Sigma_iw << SK.dc_imp[0]['up'][0,0]
mpi.report('%s DMFT cycles requested. Starting with iteration %s.'%(n_iterations,iteration_offset))
# The infamous DMFT self consistency cycle
for it in range(iteration_offset, iteration_offset + n_iterations):
mpi.report('Doing iteration: %s'%it)
# Get G0
S.G0_iw << inverse(S.Sigma_iw + inverse(S.G_iw))
# Solve the impurity problem
S.solve(h_int = H, **p)
if mpi.is_master_node():
ar['DMFT_input']['Iterations']['solver_dict_it'+str(it)] = p
ar['DMFT_results']['Iterations']['Gimp_it'+str(it)] = S.G_iw
ar['DMFT_results']['Iterations']['Gtau_it'+str(it)] = S.G_tau
ar['DMFT_results']['Iterations']['Sigma_uns_it'+str(it)] = S.Sigma_iw
# Calculate double counting
dm = S.G_iw.density()
SK.calc_dc(dm, U_interact=U, J_hund=J, orb=0, use_dc_formula=DC_type,use_dc_value=DC_value)
# Get new G
SK.symm_deg_gf(S.Sigma_iw,orb=0)
SK.put_Sigma(Sigma_imp=[S.Sigma_iw])
SK.calc_mu(precision=0.01)
S.G_iw << SK.extract_G_loc()[0]
# print densities
for sig,gf in S.G_iw:
mpi.report("Orbital %s density: %.6f"%(sig,dm[sig][0,0]))
mpi.report('Total charge of Gloc : %.6f'%S.G_iw.total_density())
if mpi.is_master_node():
ar['DMFT_results']['iteration_count'] = it
ar['DMFT_results']['Iterations']['Sigma_it'+str(it)] = S.Sigma_iw
ar['DMFT_results']['Iterations']['Gloc_it'+str(it)] = S.G_iw
ar['DMFT_results']['Iterations']['G0loc_it'+str(it)] = S.G0_iw
ar['DMFT_results']['Iterations']['dc_imp'+str(it)] = SK.dc_imp
ar['DMFT_results']['Iterations']['dc_energ'+str(it)] = SK.dc_energ
ar['DMFT_results']['Iterations']['chemical_potential'+str(it)] = SK.chemical_potential
if mpi.is_master_node(): del ar

9
doc/tutorials/images_scripts/nio.py.rst Normal file
View File

@ -0,0 +1,9 @@
.. _nio.py:
nio.py
-----------
Download :download:`nio.py <./nio.py>`.
.. literalinclude:: nio.py
:language: python

BIN
doc/tutorials/images_scripts/nio_Aw.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 25 KiB

201
doc/tutorials/images_scripts/nio_csc.py Normal file
View File

@ -0,0 +1,201 @@
from itertools import *
import numpy as np
import pytriqs.utility.mpi as mpi
from pytriqs.archive import *
from pytriqs.gf import *
import sys, pytriqs.version as triqs_version
from triqs_dft_tools.sumk_dft import *
from triqs_dft_tools.sumk_dft_tools import *
from pytriqs.operators.util.hamiltonians import *
from pytriqs.operators.util.U_matrix import *
from triqs_cthyb import *
import triqs_cthyb.version as cthyb_version
import triqs_dft_tools.version as dft_tools_version
from triqs_dft_tools.converters.vasp_converter import *
import warnings
warnings.filterwarnings("ignore", category=FutureWarning)
def dmft_cycle():
filename = 'nio'
Converter = VaspConverter(filename=filename)
Converter.convert_dft_input()
SK = SumkDFT(hdf_file = filename+'.h5', use_dft_blocks = False)
beta = 5.0
Sigma = SK.block_structure.create_gf(beta=beta)
SK.put_Sigma([Sigma])
G = SK.extract_G_loc()
SK.analyse_block_structure_from_gf(G, threshold = 1e-2)
for i_sh in range(len(SK.deg_shells)):
num_block_deg_orbs = len(SK.deg_shells[i_sh])
mpi.report('found {0:d} blocks of degenerate orbitals in shell {1:d}'.format(num_block_deg_orbs, i_sh))
for iblock in range(num_block_deg_orbs):
mpi.report('block {0:d} consists of orbitals:'.format(iblock))
for keys in SK.deg_shells[i_sh][iblock].keys():
mpi.report(' '+keys)
# Setup CTQMC Solver
n_orb = SK.corr_shells[0]['dim']
spin_names = ['up','down']
orb_names = [i for i in range(0,n_orb)]
#gf_struct = set_operator_structure(spin_names, orb_names, orb_hyb)
gf_struct = SK.gf_struct_solver[0]
mpi.report('Sumk to Solver: %s'%SK.sumk_to_solver)
mpi.report('GF struct sumk: %s'%SK.gf_struct_sumk)
mpi.report('GF struct solver: %s'%SK.gf_struct_solver)
S = Solver(beta=beta, gf_struct=gf_struct)
# Construct the Hamiltonian and save it in Hamiltonian_store.txt
H = Operator()
U = 8.0
J = 1.0
U_sph = U_matrix(l=2, U_int=U, J_hund=J)
U_cubic = transform_U_matrix(U_sph, spherical_to_cubic(l=2, convention=''))
Umat, Upmat = reduce_4index_to_2index(U_cubic)
H = h_int_density(spin_names, orb_names, map_operator_structure=SK.sumk_to_solver[0], U=Umat, Uprime=Upmat)
# Print some information on the master node
mpi.report('Greens function structure is: %s '%gf_struct)
mpi.report('U Matrix set to:\n%s'%Umat)
mpi.report('Up Matrix set to:\n%s'%Upmat)
# Parameters for the CTQMC Solver
p = {}
p["max_time"] = -1
p["random_name"] = ""
p["random_seed"] = 123 * mpi.rank + 567
p["length_cycle"] = 100
p["n_warmup_cycles"] = 2000
p["n_cycles"] = 20000
p["fit_max_moment"] = 4
p["fit_min_n"] = 30
p["fit_max_n"] = 50
p["perform_tail_fit"] = True
# Double Counting: 0 FLL, 1 Held, 2 AMF
DC_type = 0
DC_value = 59.0
# Prepare hdf file and and check for previous iterations
n_iterations = 1
iteration_offset = 0
if mpi.is_master_node():
ar = HDFArchive(filename+'.h5','a')
if not 'DMFT_results' in ar: ar.create_group('DMFT_results')
if not 'Iterations' in ar['DMFT_results']: ar['DMFT_results'].create_group('Iterations')
if not 'DMFT_input' in ar: ar.create_group('DMFT_input')
if not 'Iterations' in ar['DMFT_input']: ar['DMFT_input'].create_group('Iterations')
if not 'code_versions' in ar['DMFT_input']: ar['DMFT_input'].create_group('code_versio\
ns')
ar['DMFT_input']['code_versions']["triqs_version"] = triqs_version.version
ar['DMFT_input']['code_versions']["triqs_git"] = triqs_version.git_hash
ar['DMFT_input']['code_versions']["cthyb_version"] = cthyb_version.version
ar['DMFT_input']['code_versions']["cthyb_git"] = cthyb_version.cthyb_hash
ar['DMFT_input']['code_versions']["dft_tools_version"] = dft_tools_version.version
ar['DMFT_input']['code_versions']["dft_tools_git"] = dft_tools_version.dft_tools_hash
if 'iteration_count' in ar['DMFT_results']:
iteration_offset = ar['DMFT_results']['iteration_count']+1
S.Sigma_iw = ar['DMFT_results']['Iterations']['Sigma_it'+str(iteration_offset-1)]
SK.dc_imp = ar['DMFT_results']['Iterations']['dc_imp'+str(iteration_offset-1)]
SK.dc_energ = ar['DMFT_results']['Iterations']['dc_energ'+str(iteration_offset-1)]
SK.chemical_potential = ar['DMFT_results']['Iterations']['chemical_potential'+str(iteration_offset-1)].real
ar['DMFT_input']["dmft_script_it"+str(iteration_offset)] = open(sys.argv[0]).read()
iteration_offset = mpi.bcast(iteration_offset)
S.Sigma_iw = mpi.bcast(S.Sigma_iw)
SK.dc_imp = mpi.bcast(SK.dc_imp)
SK.dc_energ = mpi.bcast(SK.dc_energ)
SK.chemical_potential = mpi.bcast(SK.chemical_potential)
# Calc the first G0
SK.symm_deg_gf(S.Sigma_iw,orb=0)
SK.put_Sigma(Sigma_imp = [S.Sigma_iw])
SK.calc_mu(precision=0.01)
S.G_iw << SK.extract_G_loc()[0]
SK.symm_deg_gf(S.G_iw, orb=0)
#Init the DC term and the self-energy if no previous iteration was found
if iteration_offset == 0:
dm = S.G_iw.density()
SK.calc_dc(dm, U_interact=U, J_hund=J, orb=0, use_dc_formula=DC_type,use_dc_value=DC_value)
S.Sigma_iw << SK.dc_imp[0]['up'][0,0]
mpi.report('%s DMFT cycles requested. Starting with iteration %s.'%(n_iterations,iteration_offset))
# The infamous DMFT self consistency cycle
for it in range(iteration_offset, iteration_offset + n_iterations):
mpi.report('Doing iteration: %s'%it)
# Get G0
S.G0_iw << inverse(S.Sigma_iw + inverse(S.G_iw))
# Solve the impurity problem
S.solve(h_int = H, **p)
if mpi.is_master_node():
ar['DMFT_input']['Iterations']['solver_dict_it'+str(it)] = p
ar['DMFT_results']['Iterations']['Gimp_it'+str(it)] = S.G_iw
ar['DMFT_results']['Iterations']['Gtau_it'+str(it)] = S.G_tau
ar['DMFT_results']['Iterations']['Sigma_uns_it'+str(it)] = S.Sigma_iw
# Calculate double counting
dm = S.G_iw.density()
SK.calc_dc(dm, U_interact=U, J_hund=J, orb=0, use_dc_formula=DC_type,use_dc_value=DC_value)
# Get new G
SK.symm_deg_gf(S.Sigma_iw,orb=0)
SK.put_Sigma(Sigma_imp=[S.Sigma_iw])
SK.calc_mu(precision=0.01)
S.G_iw << SK.extract_G_loc()[0]
# print densities
for sig,gf in S.G_iw:
mpi.report("Orbital %s density: %.6f"%(sig,dm[sig][0,0]))
mpi.report('Total charge of Gloc : %.6f'%S.G_iw.total_density())
if mpi.is_master_node():
ar['DMFT_results']['iteration_count'] = it
ar['DMFT_results']['Iterations']['Sigma_it'+str(it)] = S.Sigma_iw
ar['DMFT_results']['Iterations']['Gloc_it'+str(it)] = S.G_iw
ar['DMFT_results']['Iterations']['G0loc_it'+str(it)] = S.G0_iw
ar['DMFT_results']['Iterations']['dc_imp'+str(it)] = SK.dc_imp
ar['DMFT_results']['Iterations']['dc_energ'+str(it)] = SK.dc_energ
ar['DMFT_results']['Iterations']['chemical_potential'+str(it)] = SK.chemical_potential
if mpi.is_master_node():
print 'calculating mu...'
SK.chemical_potential = SK.calc_mu( precision = 0.000001 )
if mpi.is_master_node():
print 'calculating GAMMA'
SK.calc_density_correction(dm_type='vasp')
if mpi.is_master_node():
print 'calculating energy corrections'
correnerg = 0.5 * (S.G_iw * S.Sigma_iw).total_density()
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,use_dc_value=DC_value)
dc_energ = SK.dc_energ[0]
if mpi.is_master_node():
ar['DMFT_results']['Iterations']['corr_energy_it'+str(it)] = correnerg
ar['DMFT_results']['Iterations']['dc_energy_it'+str(it)] = dc_energ
if mpi.is_master_node(): del ar
return correnerg, dc_energ

9
doc/tutorials/images_scripts/nio_csc.py.rst Normal file
View File

@ -0,0 +1,9 @@
.. _nio_csc.py:
nio_csc.py
-------------
Download :download:`nio_csc.py <./nio_csc.py>`.
.. literalinclude:: nio_csc.py
:language: python

28
doc/tutorials/images_scripts/plo.cfg Normal file
View File

@ -0,0 +1,28 @@
[General]
BASENAME = nio
DOSMESH = -21 55 400
[Group 1]
SHELLS = 1 2
EWINDOW = -9 2
NORMION = False
NORMALIZE = True
BANDS = 2 10
[Shell 1] # Ni d shell
LSHELL = 2
IONS = 1
CORR = True
TRANSFORM = 1.0 0.0 0.0 0.0 0.0
0.0 1.0 0.0 0.0 0.0
0.0 0.0 1.0 0.0 0.0
0.0 0.0 0.0 1.0 0.0
0.0 0.0 0.0 0.0 1.0
[Shell 2] # O p shell
LSHELL = 1
IONS = 2
CORR = False
TRANSFORM = 1.0 0.0 0.0
0.0 1.0 0.0
0.0 0.0 1.0

9
doc/tutorials/images_scripts/plo.cfg.rst Normal file
View File

@ -0,0 +1,9 @@
.. _plo.cfg:
plo.cfg
-----------
Download :download:`plo.cfg<./plo.cfg>`.
.. literalinclude:: plo.cfg
:language: bash

102
doc/tutorials/nio_csc.rst Normal file
View File

@ -0,0 +1,102 @@
.. _nio_csc:
DFT and projections
==================================================
We will perform DFT+DMFT calcluations for the charge-transfer insulator NiO. We start from scratch and provide all necessary input files to do the calculations: First for doing a single-shot calculation.
.. (and then for charge-selfconsistency).
VASP setup
-------------------------------
We start by running a simple VASP calculation to converge the charge density initially.
Use the :ref:`INCAR`, :ref:`POSCAR`, and :ref:`KPOINTS` provided and use your
own :file:`POTCAR` file.
Let us take a look in the :file:`INCAR`, where we have to specify local orbitals as basis
for our many-body calculation.
.. literalinclude:: images_scripts/INCAR
`LORBIT = 14` takes care of optimizing the projectors in the energy window defined
by `EMIN` and `EMAX`. We switch off all symmetries with `ISYM=-1` since symmetries
are not implemented in the later DMFT scripts. Finally, we select the relevant orbitals
for atom 1 (Ni, d-orbitals) and 2 (O, p-orbitals) by the two `LOCPROJ` lines.
For details refer to the VASP wiki on the `LOCPROJ <https://cms.mpi.univi
e.ac.at/wiki/index.php/LOCPROJ>`_ flag. The projectors are stored in the file `LOCPROJ`.
PLOVASP
------------------------------
Next, we postprocess the projectors, which VASP stored in the file `LOCPROJ`.
We do this by invoking :program:`plovasp plo.cfg` which is configured by an input file, e.g., named :ref:`plo.cfg`.
.. literalinclude:: images_scripts/plo.cfg
Here, in `[General]` we set the basename and the grid for calculating the density of
states. In `[Group 1]` we define a group of two shells which are orthonormalized with
respect to states in an energy window from `-9` to `2` for all ions simultanously
(`NORMION = False`). We define the two shells, which correspond to the Ni d states
and the O p states. Only the Ni shell is treated as correlated (`CORR = True`), i.e.,
is supplemented with a Coulomb interaction later in the DMFT calculation.
Converting to hdf5 file
-------------------------------
We gather the output generated by :program:`plovasp` into a hdf5 archive which :program:`dft_tools` is able to read. We do this by running :program:`python converter.py` on the script :ref:`converter.py`:
.. literalinclude:: images_scripts/converter.py
Now we are all set to perform a dmft calculation.
DMFT
==================================================
dmft script
------------------------------
Since the python script for performing the dmft loop pretty much resembles that presented in the tutorial on :ref:`SrVO3 <srvo3>`, we will not go into detail here but simply provide the script :ref:`nio.py`. Following Kunes et al. in `PRB 75 165115 (2007) <https://journals.aps.org/prb/abstract/10.1103/PhysRevB.75.165115>`_ we use :math:`U=8` and :math:`J=1`. We select :math:`\beta=5` instead of :math:`\beta=10` to ease the problem slightly. For simplicity we fix the double-counting potential to :math:`\mu_{DC}=59` eV by::
DC_value = 59.0
SK.calc_dc(dm, U_interact=U, J_hund=J, orb=0, use_dc_value=DC_value)
For sensible results run this script in parallel on at least 20 cores. As a quick check of the results, we can compare the orbital occupation from the paper cited above (:math:`n_{eg} = 0.54` and :math:`n_{t2g}=1.0`) and those from the cthyb output (check lines `Orbital up_0 density:` for a t2g and `Orbital up_2 density:` for an eg orbital). They should coincide well.
Local lattice Green's function for all projected orbitals
---------------------------------------------------------
We calculate the local lattice Green's function - now also for the uncorrelated orbitals, i.e., the O p states, for what we use the script :ref:`NiO_local_lattice_GF.py`. The result is saved in the h5 file as `G_latt_orb_it<n_it>`, where `<n_it>` is the number of the last DMFT iteration.
Spectral function on real axis: MaxEnt
--------------------------------------
To compare to results from literature we make use of the `maxent triqs application <https://triqs.github.io/maxent/master/>`_ and calculate the spectral function on real axis. Use this script to perform a crude but quick calculation: :ref:`maxent.py` using a linear real axis and a line-fit analyzer to determine the optimal :math:`\alpha`. The output is saved in the h5 file in `DMFT_results/Iterations/G_latt_orb_w_o<n_o>_it<n_it>`, where `<n_o>` is the number of the orbital and `n_it` is again the number of the last iteration. The real axis information is stored in `DMFT_results/Iterations/w_it<n_it>`.
.. image:: images_scripts/nio_Aw.png
:width: 400
:align: center
Charge self-consistent DMFT
==================================================
In this part we will perform charge self-consistent DMFT calculations. To do so we have to adapt the VASP `INCAR` such that :program:`VASP` reads the updated charge density after each step. We add the lines::
ICHARG = 5
NELM = 1000
NELMIN = 1000
IMIX=0
which makes VASP wait after each step of its iterative diagonalization until the file vasp.lock is created. It then reads the update of the charge density in the file `GAMMA`. It is terminated by an external script after a desired amount of steps, such that we deactivate all automatic stoping criterion by setting the number of steps to a very high number.
We take the respective converged DFT and DMFT calculations from before as a starting point. I.e., we copy the `CHGCAR` and `nio.h5` together with the other :program:`VASP` input files and :file:`plo.cfg` in a new directory. We use a script called :program:`vasp_dmft` to invoke :program:`VASP` in the background and start the DMFT calculation together with :program:`plovasp` and the converter. This script assumes that the dmft sript contains a function `dmft_cycle()` and also the conversion from text files to the h5 file. Additionally we have to add a few lines to calculate the density correction and calculate the correlation energy. We adapt the script straightforardly (for a working example see :ref:`nio_csc.py`). The most important new lines are::
SK.chemical_potential = SK.calc_mu( precision = 0.000001 )
SK.calc_density_correction(dm_type='vasp')
correnerg = 0.5 * (S.G_iw * S.Sigma_iw).total_density()
where the chemical potential is determined to a greater precision than before, the correction to the dft density matrix is calculated and output to the file :file:`GAMMA`. The correlation energy is calculated via Migdal-Galitzki formula. We also slightly increase the tolerance for the detection of blocks since the DFT calculation now includes some QMC noise.
To help convergence, we keep the density (i.e., the GAMMA file) fixed for a few DFT iterations. We do so since VASP uses an iterative diagonalization.
We can start the whole machinery by excecuting::
vasp_dmft -n <n_procs> -i <n_iters> -j <n_iters_dft> -p <vasp_exec> nio_csc.py

23
doc/tutorials/svo_vasp/INCAR Normal file
View File

@ -0,0 +1,23 @@
SYSTEM = SrVO3
NCORE = 4
LMAXMIX=6
EDIFF = 1.E-10
# DOS energy window
NEDOS = 2001
! switch off symmetries
ISYM=-1
# Smearing procedure
ISMEAR = -5
# the energy window to optimize projector channels
EMIN = 3.9
EMAX = 7.1
# use the PAW channel optimization
LORBIT=14
# project to V d
LOCPROJ = 2 : d : Pr

6
doc/tutorials/svo_vasp/KPOINTS Normal file
View File

@ -0,0 +1,6 @@
Automatic Mesh
0
Gamma
9 9 9
0 0 0

13
doc/tutorials/svo_vasp/POSCAR Normal file
View File

@ -0,0 +1,13 @@
SrVO3
1.0
3.8420900000 0.0000000000 0.0000000000
0.0000000000 3.8420900000 0.0000000000
0.0000000000 0.0000000000 3.8420900000
Sr V O
1 1 3
Direct
0.000000000 0.000000000 0.000000000
0.500000000 0.500000000 0.500000000
0.500000000 0.500000000 0.000000000
0.000000000 0.500000000 0.500000000
0.500000000 0.000000000 0.500000000

5
doc/guide/images_scripts/srvo3.cfg → doc/tutorials/svo_vasp/plo.cfg
View File

@ -1,7 +1,10 @@
[General]
DOSMESH = -3.0 3.0 2001
[Shell 1]
LSHELL = 2
IONS = 2
EWINDOW = -1.45 1.8
EWINDOW = -1.4 2.0
TRANSFORM = 1.0 0.0 0.0 0.0 0.0
0.0 1.0 0.0 0.0 0.0

1248
doc/tutorials/svo_vasp/svo_dos.svg Normal file
View File

File diff suppressed because it is too large Load Diff
Side by Side

After

Width:  |  Height:  |  Size: 56 KiB

357
doc/tutorials/svo_vasp/svo_notebook.ipynb Normal file
View File

@ -0,0 +1,357 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"nbsphinx": "hidden"
},
"outputs": [],
"source": [
"import matplotlib\n",
"matplotlib.use(\"Pdf\")\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"\n",
"# set matplotlib parameters\n",
"params = {'backend': 'ps',\n",
" 'axes.labelsize': 13,\n",
" 'font.size': 13,\n",
" 'legend.fontsize': 13,\n",
" 'xtick.labelsize': 13,\n",
" 'ytick.labelsize': 13,\n",
" 'text.usetex': False,\n",
" 'text.latex.preamble': \"\\\\usepackage{mathpazo}, \\\\usepackage{amsmath}\",\n",
" 'font.family': 'sans-serif',\n",
" 'font.sans-serif': ['Computer Modern Sans serif']}\n",
"plt.rcParams.update(params)\n",
"\n",
"import warnings \n",
"warnings.filterwarnings(\"ignore\") #ignore some matplotlib warnings\n",
"\n",
"# numpy\n",
"import numpy as np"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this basic example we will perform a VASP calculation for SrVO$_3$, build PLOs for the Vanadium t$_{2g}$ orbitals, and load them as SumK object, which can be then used to perform a DMFT calculation.\n",
"\n",
"## VASP setup\n",
"\n",
"First we setup the VASP [INCAR link](INCAR) file by specifing the LOCPROJ, EMIN, EMAX and LORBIT flags:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```\n",
"SYSTEM = SrVO3\n",
"NCORE = 4\n",
"LMAXMIX=6\n",
"EDIFF = 1.E-10\n",
"\n",
"# DOS energy window\n",
"NEDOS = 2001\n",
"\n",
"! switch off symmetries\n",
"ISYM=-1\n",
"\n",
"# Smearing procedure\n",
"ISMEAR = -5\n",
"\n",
"# the energy window to optimize projector channels\n",
"EMIN = 3.9\n",
"EMAX = 7.1\n",
"\n",
"# use the PAW channel optimization\n",
"LORBIT=14\n",
"\n",
"# project to V d\n",
"LOCPROJ = 2 : d : Pr\n",
"```\n",
"Moreover we prepare a [KPOINTS link](KPOINTS), [POSCAR link](POSCAR), and a POTCAR file. For the POTCAR file please use the VASP provided PBE pseudopotentials: `Sr_sv`, `V`, and `O`. \n",
"\n",
"Now VASP is executed, which should converge in roughly 27 iterations. Afterwards you should find the files LOCPROJ and PROJCAR in you directory. \n",
"\n",
"## PLOVASP\n",
"\n",
"First import the PLOVASP module of DFTTools:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Starting run with 1 MPI rank(s) at : 2019-12-05 16:12:52.689539\n"
]
}
],
"source": [
"# import plovasp converter\n",
"import triqs_dft_tools.converters.plovasp.converter as plo_converter"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, create a configuration file for plovasp [plo.cfg link](plo.cfg):\n",
"\n",
"```\n",
"[General]\n",
"DOSMESH = -3.0 3.0 2001\n",
"\n",
"[Shell 1]\n",
"LSHELL = 2\n",
"IONS = 2\n",
"EWINDOW = -1.4 2.0\n",
"\n",
"TRANSFORM = 1.0 0.0 0.0 0.0 0.0\n",
" 0.0 1.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 1.0 0.0\n",
"\n",
"```\n",
"where the energy window of the t$_{2g}$ bands is specified by `EWINDOW` and the `TRANSFORM` flag picks the correct three orbitals out of the five Vanadium $d$ orbitals [see the guide for the ordering of orbitals](../../guide/conv_vasp.html). Before running PLOVASP, make sure that the Fermi energy is written in the first line of the `LOCPROJ` file, or copy it there (see the VASP [interface guide](../../guide/conv_vasp.html) for more information). The first line should look like\n",
"```\n",
"1 729 21 5 5.3834262 # of spin, # of k-points, # of bands, # of proj, Efermi\n",
"```\n",
"\n",
"Now run PLOVASP:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Read parameters:\n",
"0 -> {'m': 0, 'l': 2, 'isite': 2, 'label': 'dxy'}\n",
"1 -> {'m': 1, 'l': 2, 'isite': 2, 'label': 'dyz'}\n",
"2 -> {'m': 2, 'l': 2, 'isite': 2, 'label': 'dz2'}\n",
"3 -> {'m': 3, 'l': 2, 'isite': 2, 'label': 'dxz'}\n",
"4 -> {'m': 4, 'l': 2, 'isite': 2, 'label': 'dx2-y2'}\n",
" Found POSCAR, title line: SrVO3\n",
" Total number of ions: 5\n",
" Number of types: 3\n",
" Number of ions for each type: [1, 1, 3]\n",
"\n",
" Total number of k-points: 729\n",
" Total number of tetrahedra: 4374\n",
"eigvals from LOCPROJ\n",
"\n",
" Unorthonormalized density matrices and overlaps:\n",
" Spin: 1\n",
" Site: 2\n",
" Density matrix Overlap\n",
" 0.5875772 0.0015679 -0.0003707 0.0015674 0.0000000 0.9294791 -0.0000080 -0.0000078 -0.0000080 -0.0000001\n",
" 0.0015679 0.5876177 -0.0001854 -0.0016078 0.0003240 -0.0000080 0.9294790 -0.0000042 0.0000080 0.0000070\n",
" -0.0003707 -0.0001854 0.5815486 -0.0001854 -0.0000000 -0.0000078 -0.0000042 0.9715751 -0.0000038 0.0000003\n",
" 0.0015674 -0.0016078 -0.0001854 0.5876172 -0.0003240 -0.0000080 0.0000080 -0.0000038 0.9294791 -0.0000066\n",
" 0.0000000 0.0003240 -0.0000000 -0.0003240 0.5815487 -0.0000001 0.0000070 0.0000003 -0.0000066 0.9715748\n",
"\n",
" Generating 1 shell...\n",
"\n",
" Shell : 1\n",
" Orbital l : 2\n",
" Number of ions: 1\n",
" Dimension : 3\n",
" Correlated : True\n",
" Ion sort : [1]\n",
"Density matrix:\n",
" Shell 1\n",
" Site 1\n",
" 0.3332630 0.0021719 0.0021714\n",
" 0.0021719 0.3333128 -0.0022211\n",
" 0.0021714 -0.0022211 0.3333123\n",
" trace: 0.9998880790966638\n",
"\n",
" Impurity density: 0.9998880790966638\n",
"\n",
"Overlap:\n",
" Site 1\n",
"[[ 1. -0. 0.]\n",
" [-0. 1. 0.]\n",
" [ 0. 0. 1.]]\n",
"\n",
"Local Hamiltonian:\n",
" Shell 1\n",
" Site 1\n",
" 0.5633806 0.0007563 0.0007563\n",
" 0.0007563 0.5633801 -0.0007559\n",
" 0.0007563 -0.0007559 0.5633801\n",
"\n",
"Evaluating DOS...\n",
" Shell 1\n",
" Total number of states: [[[7.33737319 7.48285647 7.28002405]]]\n",
" Storing ctrl-file...\n",
" Storing PLO-group file 'vasp.pg1'...\n",
" Density within window: 0.9999741659673522\n"
]
}
],
"source": [
"# Generate and store PLOs\n",
"plo_converter.generate_and_output_as_text('plo.cfg', vasp_dir='./')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"PLOVASP created one shell with three orbitals, which are equally filled by 1/3, one electron in total. Additionally we calculated the density of states. Both in VASP and PLOVASP. The later stores the data in the file pdos_x.dat, which can be simply plotted with [matplotlib](https://matplotlib.org/). The result should look similar to:\n",
"\n",
"![](svo_dos.svg)\n",
"\n",
"Here the gray area highlights the energy window for the PLOs. The total DOS of VASP (blue) coincides with the PLO DOS in the window, as we re-orthonormalized the projector functions in the given window, picking up also Oxygen weight. This setting is closed to the result of maximally localized Wannier functions created with [wannier90](http://www.wannier.org/) without running the actual minimization of the spread. Note, for a proper comparison one can use the hydrogen projector in VASP by using the the line `LOCPROJ= 2 : d : Hy`, instead of `Pr`. \n",
"\n",
"\n",
"## Converting to hdf5 file\n",
"\n",
"Finally we can run the VASP converter to create a h5 file:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Reading input from vasp.ctrl...\n",
"{\n",
" \"ngroups\": 1,\n",
" \"nk\": 729,\n",
" \"ns\": 1,\n",
" \"nc_flag\": 0\n",
"}\n",
"\n",
" No. of inequivalent shells: 1\n"
]
}
],
"source": [
"# import VASPconverter\n",
"from triqs_dft_tools.converters.vasp_converter import *\n",
"\n",
"\n",
"# create Converter\n",
"Converter = VaspConverter(filename = 'vasp')\n",
"# run the converter\n",
"Converter.convert_dft_input()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The resulting h5 file `vasp.h5` can now be loaded as sumk object via:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# SumK\n",
"from triqs_dft_tools.sumk_dft_tools import SumkDFTTools\n",
"\n",
"SK = SumkDFTTools(hdf_file='vasp.h5', use_dft_blocks = False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here one should carefully determine the block structure manually. This is important to find degenerate orbitals and spin-channels:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"found 1 blocks of degenerate orbitals in shell 0\n",
"block 0 consists of orbitals:\n",
" up_2\n",
" up_0\n",
" up_1\n",
" down_2\n",
" down_1\n",
" down_0\n"
]
}
],
"source": [
"Sigma = SK.block_structure.create_gf(beta=40)\n",
"SK.put_Sigma([Sigma])\n",
"G = SK.extract_G_loc()\n",
"SK.analyse_block_structure_from_gf(G, threshold = 1e-3)\n",
"for i_sh in range(len(SK.deg_shells)):\n",
" num_block_deg_orbs = len(SK.deg_shells[i_sh])\n",
" mpi.report('found {0:d} blocks of degenerate orbitals in shell {1:d}'.format(num_block_deg_orbs, i_sh))\n",
" for iblock in range(num_block_deg_orbs):\n",
" mpi.report('block {0:d} consists of orbitals:'.format(iblock))\n",
" for keys in SK.deg_shells[i_sh][iblock].keys():\n",
" mpi.report(' '+keys)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This minimal example extracts the block structure by calculating once the local Green's functions and then finds degenerate orbitals with a certain threshold in `Gloc`. Afterwards the result is reported, where 1 block is found with size 6 (3x2 orbitals for spin), where a all 6 orbitals are marked as degenerate. This is indeed correct in cubic SrVO$_3$, as all 3 t$_{2g}$ orbitals are degenerate. Note: for a magnetic calculation one has to break the symmetry between up and down at this point manually. Moreover, one can reduce the threshold for example to `1e-5` and finds that then the degeneracy of the 3 t$_{2g}$ orbitals is not found anymore, resulting in only two degenerate blocks for spin up and down, each with size 3x3.\n",
"\n",
"Afterwards the exact same DMFT script as in the [Wien2k tutorial](../srvo3.html) can be used. For a more elaborate example including charge self-consistency take a look at the [VASP NiO example](../nio_csc.html)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.15+"
}
},
"nbformat": 4,
"nbformat_minor": 2
}

153
doc/vasp/Makefile
View File

@ -1,153 +0,0 @@
# Makefile for Sphinx documentation
#
# You can set these variables from the command line.
SPHINXOPTS =
SPHINXBUILD = sphinx-build
PAPER =
BUILDDIR = build
# Internal variables.
PAPEROPT_a4 = -D latex_paper_size=a4
PAPEROPT_letter = -D latex_paper_size=letter
ALLSPHINXOPTS = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) source
# the i18n builder cannot share the environment and doctrees with the others
I18NSPHINXOPTS = $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) source
.PHONY: help clean html dirhtml singlehtml pickle json htmlhelp qthelp devhelp epub latex latexpdf text man changes linkcheck doctest gettext
help:
@echo "Please use \`make <target>' where <target> is one of"
@echo " html to make standalone HTML files"
@echo " dirhtml to make HTML files named index.html in directories"
@echo " singlehtml to make a single large HTML file"
@echo " pickle to make pickle files"
@echo " json to make JSON files"
@echo " htmlhelp to make HTML files and a HTML help project"
@echo " qthelp to make HTML files and a qthelp project"
@echo " devhelp to make HTML files and a Devhelp project"
@echo " epub to make an epub"
@echo " latex to make LaTeX files, you can set PAPER=a4 or PAPER=letter"
@echo " latexpdf to make LaTeX files and run them through pdflatex"
@echo " text to make text files"
@echo " man to make manual pages"
@echo " texinfo to make Texinfo files"
@echo " info to make Texinfo files and run them through makeinfo"
@echo " gettext to make PO message catalogs"
@echo " changes to make an overview of all changed/added/deprecated items"
@echo " linkcheck to check all external links for integrity"
@echo " doctest to run all doctests embedded in the documentation (if enabled)"
clean:
-rm -rf $(BUILDDIR)/*
html:
$(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html
@echo
@echo "Build finished. The HTML pages are in $(BUILDDIR)/html."
dirhtml:
$(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) $(BUILDDIR)/dirhtml
@echo
@echo "Build finished. The HTML pages are in $(BUILDDIR)/dirhtml."
singlehtml:
$(SPHINXBUILD) -b singlehtml $(ALLSPHINXOPTS) $(BUILDDIR)/singlehtml
@echo
@echo "Build finished. The HTML page is in $(BUILDDIR)/singlehtml."
pickle:
$(SPHINXBUILD) -b pickle $(ALLSPHINXOPTS) $(BUILDDIR)/pickle
@echo
@echo "Build finished; now you can process the pickle files."
json:
$(SPHINXBUILD) -b json $(ALLSPHINXOPTS) $(BUILDDIR)/json
@echo
@echo "Build finished; now you can process the JSON files."
htmlhelp:
$(SPHINXBUILD) -b htmlhelp $(ALLSPHINXOPTS) $(BUILDDIR)/htmlhelp
@echo
@echo "Build finished; now you can run HTML Help Workshop with the" \
".hhp project file in $(BUILDDIR)/htmlhelp."
qthelp:
$(SPHINXBUILD) -b qthelp $(ALLSPHINXOPTS) $(BUILDDIR)/qthelp
@echo
@echo "Build finished; now you can run "qcollectiongenerator" with the" \
".qhcp project file in $(BUILDDIR)/qthelp, like this:"
@echo "# qcollectiongenerator $(BUILDDIR)/qthelp/plotools.qhcp"
@echo "To view the help file:"
@echo "# assistant -collectionFile $(BUILDDIR)/qthelp/plotools.qhc"
devhelp:
$(SPHINXBUILD) -b devhelp $(ALLSPHINXOPTS) $(BUILDDIR)/devhelp
@echo
@echo "Build finished."
@echo "To view the help file:"
@echo "# mkdir -p $$HOME/.local/share/devhelp/plotools"
@echo "# ln -s $(BUILDDIR)/devhelp $$HOME/.local/share/devhelp/plotools"
@echo "# devhelp"
epub:
$(SPHINXBUILD) -b epub $(ALLSPHINXOPTS) $(BUILDDIR)/epub
@echo
@echo "Build finished. The epub file is in $(BUILDDIR)/epub."
latex:
$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
@echo
@echo "Build finished; the LaTeX files are in $(BUILDDIR)/latex."
@echo "Run \`make' in that directory to run these through (pdf)latex" \
"(use \`make latexpdf' here to do that automatically)."
latexpdf:
$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
@echo "Running LaTeX files through pdflatex..."
$(MAKE) -C $(BUILDDIR)/latex all-pdf
@echo "pdflatex finished; the PDF files are in $(BUILDDIR)/latex."
text:
$(SPHINXBUILD) -b text $(ALLSPHINXOPTS) $(BUILDDIR)/text
@echo
@echo "Build finished. The text files are in $(BUILDDIR)/text."
man:
$(SPHINXBUILD) -b man $(ALLSPHINXOPTS) $(BUILDDIR)/man
@echo
@echo "Build finished. The manual pages are in $(BUILDDIR)/man."
texinfo:
$(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo
@echo
@echo "Build finished. The Texinfo files are in $(BUILDDIR)/texinfo."
@echo "Run \`make' in that directory to run these through makeinfo" \
"(use \`make info' here to do that automatically)."
info:
$(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo
@echo "Running Texinfo files through makeinfo..."
make -C $(BUILDDIR)/texinfo info
@echo "makeinfo finished; the Info files are in $(BUILDDIR)/texinfo."
gettext:
$(SPHINXBUILD) -b gettext $(I18NSPHINXOPTS) $(BUILDDIR)/locale
@echo
@echo "Build finished. The message catalogs are in $(BUILDDIR)/locale."
changes:
$(SPHINXBUILD) -b changes $(ALLSPHINXOPTS) $(BUILDDIR)/changes
@echo
@echo "The overview file is in $(BUILDDIR)/changes."
linkcheck:
$(SPHINXBUILD) -b linkcheck $(ALLSPHINXOPTS) $(BUILDDIR)/linkcheck
@echo
@echo "Link check complete; look for any errors in the above output " \
"or in $(BUILDDIR)/linkcheck/output.txt."
doctest:
$(SPHINXBUILD) -b doctest $(ALLSPHINXOPTS) $(BUILDDIR)/doctest
@echo "Testing of doctests in the sources finished, look at the " \
"results in $(BUILDDIR)/doctest/output.txt."

2
doc/vasp/build/.gitignore vendored
View File

@ -1,2 +0,0 @@
*
!/.gitignore

19
doc/vasp/source/adv_example.cfg
View File

@ -1,19 +0,0 @@
[General]
EFERMI = -0.6
[Group 1]
SHELLS = 1 2
EWINDOW = -7.6 2.7
[Shell 1]
# Ni shell
IONS = 5 6 7 8
LSHELL = 2
TRANSFORM =
0.0 0.0 0.0 0.0 1.0
0.0 0.0 1.0 0.0 0.0
[Shell 2]
# Oxygen shell
IONS = 9..20
LSHELL = 1

43
doc/vasp/source/charge_selfcons.rst
View File

@ -1,43 +0,0 @@
.. highlight:: python
#######################
Charge self-consistency
#######################
Introduction
************
Running a DFT+DMFT calculation in the charge self-consistent mode
requires additional process management whereby the control is given
to the VASP and to the TRIQS code in an alternating manner.
It is implemented by a lock system. A VASP process is launched in
the background and a self-consistency (SC) script in the foreground.
Once VASP reaches the point where the projectors are generated
it creates a lock file `vasp.lock` and waits until the lock file is
removed. The SC script, in turn, waits for the VASP process and once
the lock file is created it starts a DMFT iteration. The DMFT iteration
must finish by generating a Kohn-Sham (KS) density matrix (file `GAMMA`)
and removing the lock file. The VASP process then reads in `GAMMA`
and proceeds with the next iteration.
The DMFT iteration is performed using a user-defined script.
However, the control part is maintained by two universal scripts:
`sc_dmft.sh` and `sc_dmft.py`.
The first script, `sc_dmft.sh`, launches the VASP process
in the background mode, takes its pid, and launches `sc_dmft.py`
in the foreground mode.
Both processes are run within an MPI environment with an appropriate
number of MPI nodes.
The second script, `sc_dmft.py`, is responsible for the charge self-consitency
logic described in the first paragraph. It also combines the total
energy contributions coming from the DFT and DMFT parts.
The script imports a user-defined DMFT script `dmft_cycle.py` which
must produce KS density-matrix file `GAMMA` for the next DFT iteration
and also must return the correlation (including the double counting) energy
of the impurity model as well as a correlation correction to the
DFT band energy (resulting from the difference between the bare
and DMFT density matrices).

26
doc/vasp/source/code_struct.rst
View File

@ -1,26 +0,0 @@
Code Structure
##############
.. toctree::
vaspio
plotools
converter
charge_selfcons
The program consists of three main parts:
* :doc:`Import of data from VASP files </vaspio>`
* Processing of projectors according to an input config-file
* Conversion of the DFT data to TRIQS h5-file
Import of data from VASP files is implemented in `vaspio.py`. The data
is read from VASP files and then stored in objects in raw format
(i.e. practically no processing is done at this stage).
These objects are then combined into a dictionary which can be easily
passed to other routines.
The basic workflow is prescribed as follows:
* raw data is read from VASP files and passed to the main part
* in the main part the input config-file is read and the projectors are selected and process accordingly
* the processed data is stored into output text-files
* when the TRIQS-converter is requested the data is read from text-files and written into a h5-file in an appropriate format

244
doc/vasp/source/conf.py
View File

@ -1,244 +0,0 @@
# -*- coding: utf-8 -*-
#
# plotools documentation build configuration file, created by
# sphinx-quickstart on Tue Feb 3 20:25:36 2015.
#
# This file is execfile()d with the current directory set to its containing dir.
#
# Note that not all possible configuration values are present in this
# autogenerated file.
#
# All configuration values have a default; values that are commented out
# serve to show the default.
import sys, os
# If extensions (or modules to document with autodoc) are in another directory,
# add these directories to sys.path here. If the directory is relative to the
# documentation root, use os.path.abspath to make it absolute, like shown here.
sys.path.insert(0, os.path.abspath('../../../python/vasp'))
sys.path.insert(0, os.path.abspath('../../../c'))
# -- General configuration -----------------------------------------------------
# If your documentation needs a minimal Sphinx version, state it here.
#needs_sphinx = '1.0'
# Add any Sphinx extension module names here, as strings. They can be extensions
# coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
extensions = ['sphinx.ext.autodoc', 'sphinx.ext.doctest', 'sphinx.ext.pngmath', 'sphinx.ext.mathjax']
mathjax_path = 'http://cdn.mathjax.org/mathjax/latest/MathJax.js'
# Add any paths that contain templates here, relative to this directory.
templates_path = ['_templates']
# The suffix of source filenames.
source_suffix = '.rst'
# The encoding of source files.
#source_encoding = 'utf-8-sig'
# The master toctree document.
master_doc = 'index'
# General information about the project.
project = u'PLOVasp'
copyright = u'2015, Oleg E. Peil'
# The version info for the project you're documenting, acts as replacement for
# |version| and |release|, also used in various other places throughout the
# built documents.
#
# The short X.Y version.
version = '0.1'
# The full version, including alpha/beta/rc tags.
release = '0.1'
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.
#language = None
# There are two options for replacing |today|: either, you set today to some
# non-false value, then it is used:
#today = ''
# Else, today_fmt is used as the format for a strftime call.
#today_fmt = '%B %d, %Y'
# List of patterns, relative to source directory, that match files and
# directories to ignore when looking for source files.
exclude_patterns = []
# The reST default role (used for this markup: `text`) to use for all documents.
#default_role = None
# If true, '()' will be appended to :func: etc. cross-reference text.
#add_function_parentheses = True
# If true, the current module name will be prepended to all description
# unit titles (such as .. function::).
#add_module_names = True
# If true, sectionauthor and moduleauthor directives will be shown in the
# output. They are ignored by default.
#show_authors = False
# The name of the Pygments (syntax highlighting) style to use.
pygments_style = 'sphinx'
# A list of ignored prefixes for module index sorting.
#modindex_common_prefix = []
# -- Options for HTML output ---------------------------------------------------
# The theme to use for HTML and HTML Help pages. See the documentation for
# a list of builtin themes.
html_theme = 'default'
# Theme options are theme-specific and customize the look and feel of a theme
# further. For a list of options available for each theme, see the
# documentation.
#html_theme_options = {}
# Add any paths that contain custom themes here, relative to this directory.
#html_theme_path = []
# The name for this set of Sphinx documents. If None, it defaults to
# "<project> v<release> documentation".
#html_title = None
# A shorter title for the navigation bar. Default is the same as html_title.
#html_short_title = None
# The name of an image file (relative to this directory) to place at the top
# of the sidebar.
#html_logo = None
# The name of an image file (within the static path) to use as favicon of the
# docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32
# pixels large.
#html_favicon = None
# Add any paths that contain custom static files (such as style sheets) here,
# relative to this directory. They are copied after the builtin static files,
# so a file named "default.css" will overwrite the builtin "default.css".
html_static_path = ['_static']
# If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
# using the given strftime format.
#html_last_updated_fmt = '%b %d, %Y'
# If true, SmartyPants will be used to convert quotes and dashes to
# typographically correct entities.
#html_use_smartypants = True
# Custom sidebar templates, maps document names to template names.
#html_sidebars = {}
# Additional templates that should be rendered to pages, maps page names to
# template names.
#html_additional_pages = {}
# If false, no module index is generated.
#html_domain_indices = True
# If false, no index is generated.
#html_use_index = True
# If true, the index is split into individual pages for each letter.
#html_split_index = False
# If true, links to the reST sources are added to the pages.
#html_show_sourcelink = True
# If true, "Created using Sphinx" is shown in the HTML footer. Default is True.
#html_show_sphinx = True
# If true, "(C) Copyright ..." is shown in the HTML footer. Default is True.
#html_show_copyright = True
# If true, an OpenSearch description file will be output, and all pages will
# contain a <link> tag referring to it. The value of this option must be the
# base URL from which the finished HTML is served.
#html_use_opensearch = ''
# This is the file name suffix for HTML files (e.g. ".xhtml").
#html_file_suffix = None
# Output file base name for HTML help builder.
htmlhelp_basename = 'plotoolsdoc'
# -- Options for LaTeX output --------------------------------------------------
latex_elements = {
# The paper size ('letterpaper' or 'a4paper').
#'papersize': 'letterpaper',
# The font size ('10pt', '11pt' or '12pt').
#'pointsize': '10pt',
# Additional stuff for the LaTeX preamble.
#'preamble': '',
}
# Grouping the document tree into LaTeX files. List of tuples
# (source start file, target name, title, author, documentclass [howto/manual]).
latex_documents = [
('index', 'plotools.tex', u'plotools Documentation',
u'Oleg E. Peil', 'manual'),
]
# The name of an image file (relative to this directory) to place at the top of
# the title page.
#latex_logo = None
# For "manual" documents, if this is true, then toplevel headings are parts,
# not chapters.
#latex_use_parts = False
# If true, show page references after internal links.
#latex_show_pagerefs = False
# If true, show URL addresses after external links.
#latex_show_urls = False
# Documents to append as an appendix to all manuals.
#latex_appendices = []
# If false, no module index is generated.
#latex_domain_indices = True
# -- Options for manual page output --------------------------------------------
# One entry per manual page. List of tuples
# (source start file, name, description, authors, manual section).
man_pages = [
('index', 'plotools', u'plotools Documentation',
[u'Oleg E. Peil'], 1)
]
# If true, show URL addresses after external links.
#man_show_urls = False
# -- Options for Texinfo output ------------------------------------------------
# Grouping the document tree into Texinfo files. List of tuples
# (source start file, target name, title, author,
# dir menu entry, description, category)
texinfo_documents = [
('index', 'plotools', u'plotools Documentation',
u'Oleg E. Peil', 'plotools', 'One line description of project.',
'Miscellaneous'),
]
# Documents to append as an appendix to all manuals.
#texinfo_appendices = []
# If false, no module index is generated.
#texinfo_domain_indices = True
# How to display URL addresses: 'footnote', 'no', or 'inline'.
#texinfo_show_urls = 'footnote'

124
doc/vasp/source/config.rst
View File

@ -1,124 +0,0 @@
Input Config-file
=================
A config-file describes subsets of PLOs that are to be generated.
The PLOs are defined in terms of `shells` determined uniquely by an orbital
number `l` and a set of ions (nomrmally, of the same sort).
The shells are further combined into `groups` such that PLO in each group
are orthogonalized together. This, for instance, allows to create a group of several
mutually orthogonal subsets of PLOs.
A group is characterized by a single projection energy window.
A config-file contains three types of sections:
- **[General]** : providing information applicable to all projected shells
(e.g. Fermi level)
- **[Shell <Ns>]** : each section like this describes a PLO shell, with the index
`Ns` used for referencing
- **[Group <tag>]** : describes shell groups
The format requirements are relatively flexible. A config-file must contain
at least one `[Shell]` section. If there is only one shell defined, it is possible
to specify the energy window by providing parameters `EMIN`, `EMAX` (see below)
right in this section, in which case a group
will be created automatically and the `[Group]` section can be omitted.
If, nevertheless, a group referencing the single shell is explicitly given
the energy window parameters provided in the `[Group]` have higher priority
and in case of conflict with `[Shell]` section a warning is issued.
An example of a config-file:
.. literalinclude:: adv_example.cfg
Here two shells, one corresponding to `d` states on ions 5-8, another to `p`
states of ions 9-20, are created. They form a single group that, by default, will be
orthogonalized within a window `[-7.6, 2.7]` eV. Also Fermi level is explicitly
specified, which might be necessary sometimes, e.g., for non-self-consistent calcuation
of the band structure.
Below, the sections and their parameters are described.
All parameter names are case-insensitive.
Section [General]
-----------------
**Required parameters:**
In principle, there are unconditionally required parameters in this section.
However, if VASP data file do not contain a meaningful value of the Fermi level
it must be given here using parameter *EFERMI*. Note that if this parameter
is given it will override a value that might have been read from VASP files.
**Optional parameters:**
- *BASENAME* (string): a basename for output files. The filenames will be
of the form '<basename>.*'.
- *DOSMESH* (float, float, int): if this parameter is provided the projected
DOS for each projected shell will be generated, with the energy mesh parameters
given by the energy range (two floats) and the number of points (int). It is also
possible to omit the energy range, in which case it will be set to the energy window
of the corresponding projector group.
Section [Shell <Ns>]
--------------------
Defines a projected shell with an integer index `<Ns>`. Ideally, the indices should
run uniformly starting from 1. However, since they are used only to reference
shells in group sections, their values are not important. One should only
make sure that there are no sections with the same name, in which case one
of them will be ignored by config parser.
**Required parameters:**
- *IONS* ([int]): provides a list of ions. The list can be either given
by a explicit enumeration of ion indices or by a range `N1..N2` (ions `N1` and `N2`
will be included).
- *LSHELL* (int): orbital number `l` of the shell (`l = 0,1,2,3` for `s,p,d,f`, respectively).
**Optional parameters:**
- *TRANSFORM* (matrix of floats): transformation matrices
(real or complex) that are applied to projectors before orthogonalization.
The number of columns `Nc` must be consistent with the number of orbitals
(`Nc = 2l+1` for real matrices and `Nc = 2(2l+1)` for complex ones).
The dimension of the resulting orbital subspace is determined by the number of rows.
- *TRANSFILE* (str): file containing transformation matrices for individual ions.
The file must contain rows of floats. Empty lines and lines starting with '#' are ignored.
The data is interpreted as follows:
* The number of rows is divided by the number of ions to give the number
of rows per ion, `Nr`.
* The number of columns is divided by `Nm = 2 * l + 1` to give `nsize`.
There are, then, three options:
#. if `nc_flag = 0`, i.e. a calculation is collinear (no spin-orbit coupling),
and `nsize = 1` the matrices are considered to be real;
#. if `nc_flag = 0` and `nsize = 2` the matrices are considered to be complex;
#. if `nc_flag = 1`, i.e. a calculation is non-collinear, it is expected
that `nsize = 4`, and the matrices are considered to be complex and
spin-dependent.
* The subspace dimension is determined simply as `Ndim = Nr / nsize`.
In all cases when a division is performed the result must be integer. Otherwise,
the matrices are considered to be inconsistent.
- *EWINDOW* (float, float): energy window. Should be given only if no excplicit groups
is specified. Otherwise, the values are overriden by group parameters.
Section [Group <tag>]
---------------------
Defines a group of shells. Note that the group tag can be any string without whitespaces.
It will be used to tag intermediate output files.
**Required parameters:**
- *SHELLS* ([int]): indices refering to shells forming the group.
- *EWINDOW* (float, float): the bottom and top of the energy window with respect to the Fermi level.
**Optional parameters:**
- NORMALIZE (True/False): if True, orthogonalizetion is performed (default behavior).
- *NORMION* (True/False): if True, orthogonalization is performed on each site
separately; if False, all projectors of the group are orthogonalized together.

4
doc/vasp/source/converter.rst
View File

@ -1,4 +0,0 @@
TRIQS Converter
###############
The converter provides an interface between a DFT code and TRIQS solvers.

12
doc/vasp/source/example.cfg
View File

@ -1,12 +0,0 @@
[General]
[PLO Group 1]
IONS = 5 6 7 8
EMIN = -0.6
EMAX = 2.7
LSHELL = 2
RTRANSFORM =
0.0 0.0 0.0 0.0 1.0
0.0 0.0 1.0 0.0 0.0

27
doc/vasp/source/index.rst
View File

@ -1,27 +0,0 @@
.. plotools documentation master file, created by
sphinx-quickstart on Tue Feb 3 20:25:36 2015.
You can adapt this file completely to your liking, but it should at least
contain the root `toctree` directive.
Welcome to PLOVasp's documentation!
====================================
PLOVasp documentation.
Contents:
.. toctree::
:maxdepth: 1
users_guide
code_struct
Indices and tables
==================
* :ref:`genindex`
* :ref:`modindex`
* :ref:`search`

64
doc/vasp/source/notes.txt
View File

@ -1,64 +0,0 @@
Processing tools for PLOs generated by VASP.
The main purpose of this set of tools is to process raw data
generated by VASP and to convert data to the TRIQS format. The raw data
is read from files:
POSCAR (atomic positions)
EIGVAL (eignevalues)
PLOCAR (non-orthogonalized PLOs and Fermi-weights)
IBZKPT (k-points and tetrahedra)
An appropriate set of orthogonalized projectors is defined by
parameters defined in the config file (config-like syntax).
The config-file allows to define several groups of projectors.
Structure of PLOtools:
vaspio.py: reading of VASP-generated files
vaspio.py:
All VASP data are represented by objects which contain data read
from corresponding files. These objects will subsequently be used to
handle the data and convert it into a more functional internal representation.
Functions
read_lines(filename): generator yielding lines from a file <filename>
Classes:
Plocar: raw PLO data from PLOCAR file;
the data itself is read using an auxiliary C-routine 'read_plocar'
Poscar: structure data from POSCAR file
Kpoints: k-point data from IBZKPT file
note that k-point data is also contained in EIGENVAL file;
the two k-point sets will be checked for consistency.
Eigenval: Kohn-Sham eigenvalues as well as k-points with weights
Symmcar: symmetry operations stored by VASP into SYMMCAR file
ploortho.py (or projectors.py)
Set of routines for processing projectors. The functionality includes:
-- selecting a proper subset of non-orthogonalized projectors from the raw VASP input
-- transforming and orthogonalizing projectors
-- calculating DFT density matrix and local Hamiltonian
General design:
Input data: conf_pars (output of 'parse_input()'), VASP file descriptors (Eigenval, Plocar, etc.)
* A projector for a given k-point is described by class 'Projector'
PLOs project on a set or orbitals and a set of ions.
* Projector set is a container of 'Projector' objects.
Out of optimization purposes the projectors are stored in a multi-dimensional
array. A view in terms of Projector objects is, however, possible.

191
doc/vasp/source/plotools.rst
View File

@ -1,191 +0,0 @@
.. highlight:: python
#########
PLO tools
#########
Introduction
************
This set of tools is intended for processing of raw projectors read
from VASP. One of the main tasks is to generate an orthonormalized subset
of PLOs constructed according to the :doc:`config-file </config>`.
To produce the final output the following steps are undertaken:
* Parse input config-file
* Input raw VASP data
* Convert the raw VASP data into an internal representaion and check it for consistency.
* Generate a set of projector shells according to the parameters of the config-file
* Create a set of projector groups
* Perform necessary group operations (such as :ref:`orthogonalization<ortho>`) on the constituing shells
* Calculate and output some useful quantities (bare density matrix, DOS, etc.)
Initial Processing
******************
The raw data from VASP files is initially read in simple objects (containers).
Then these objects are combined in an another object containing all necessary
electronic structure data. At this stage simple consistency checks are performed:
* the basic dimensions, such as the number of bands, number of `k`-points, etc.,
are consistent for all data
* the `k`-point mesh are read both the IBZKPT and EIGENVAL and it is worth checking
that both sets are coinciding
* in case tetrahedron data is read from IBZKPT, the tetrahedron volume must be related
to the total volume of the unit cell as derived from POSCAR
All electronic structure from VASP is stored in a class ElectronicStructure:
.. autoclass:: triqs_dft_tools.converters.plovasp.elstruct.ElectronicStructure
:members:
Consistency with parameters
* parameters in the config-file should pass trivial checks such as that the ion
list does not contain non-existing ions (boundary check for ion indices)
.. function:: check_vasp_data_consistency(conf_pars, vasp_data)
**Parameters**:
- *conf_pars* (dict) : dictionary of input parameters from conf-file
- *vasp_data* (dict) : dictionary containing all VASP data
**Returns**:
*None*
**Raises**:
A meaningful exception indicating an inconsistency in the input data
Selecting projector subsets
===========================
The first step of PLO processing is to select subsets of projectors
corresponding to PLO groups. Each group contains a set of shells.
Each projector shell is represented by an object 'ProjectorShell'
that contains an array of projectors and information on the shell itself
(orbital number, ions, etc.). 'ProjectorShell's are contained in
both a list of shells (according to the original list as read
from config-file) and in a 'ProjectorGroup' object, the latter
also providing information about the energy window.
`[In fact, shell container can be a simple dictionary.]`
Order of operations:
- transform projectors (all bands) in each shell
- select transformed shell projectors for a given group within the window
- orthogonalize if necessary projectors within a group by performing
the following operations for each k-point:
* combine all projector shells into a single array
* orthogonalize the array
* distribute back the arrays assuming that the order is preserved
.. autoclass:: triqs_dft_tools.converters.plovasp.proj_shell.ProjectorShell
:members:
The class is using a helper function `select_bands()` for selecting a subset of bands.
.. function:: select_bands(eigvals, emin, emax)
**Parameters**:
- *eigvals* (numpy.array) : array of eigenvalues
- *emin*, *emax* (float) : energy window
**Returns**:
*ib_win*, *nb_min*, *nb_max* (numpy.array[int], int, int) :
lowest and highest indices of the selected bands
.. _ortho:
Orthogonalization
-----------------
At the second stage the selected projectors are orthogonalized (orthonormalized).
Orthogonalization can be performed in different ways if projection is made
on several ions or if several correlated shells per ion are considered.
In the case of several correlated ions per unit cell (and one correlated shell per ion)
at least two options can be considered:
#. Projectors are normalized for each ion separetely. In this case, corresponding
Wannier functions for different ions are generally not orthogonal.
#. Projectors are normalized for all ions in the unit cell simultaneously. This
ensures that the Wannier functions for different ions are mutually orthogonal.
The way the normalization of a PLO group is done is controlled by two group parameters:
- *NORMALIZE* (True/False) : indicates whether the PLO group is normalized (True by default)
- *NORMION* (True/False) : indicates whether the PLO group is normalized on a per-ion basis
(False by default)
Storing generated projectors
****************************
After the PLOs are generated they are stored to text files which can be subsequently
converted to TRIQS h5-files (using the converter). The general format of the file
is a JSON-header containing all necessary parameters followed by a set of arrays.
There is always one (control) file containing general information (`k`-kpoints, lattice vectors etc.)
and `at least` one file containing correlated groups (one file for each group).
Control file format
===================
Filename '<namebase>.ctrl'. Contains the data shared between all shells.
The JSON-header consists of the following elements:
* *nk*: number of `k`-points
* *ns*: number of spin channels
* *nc_flag*: collinear/noncollinear case (False/True)
* *ng*: number of projector groups
* Symmetry information (list of symmetry operations)
* *efermi*: Fermi level (optional)
* Lattice information
Projector-group file format
===========================
Projector-group files have names '<namebase>.plog<Ng>'.
They essentially contain serialized objects of class 'ProjectorGroup'.
The JSON-header has, thus, the following elements:
* *shells*: list of shells
* each shell is a dictionary:
- *lshell*: orbital number `l`
- *nion*: number of ions
- *ndim*: number of orbitals/ion
- *nbmax*: maxmimum number of bands (needed for array allocations)
* *emin*, *emax*: energy window

29
doc/vasp/source/users_guide.rst
View File

@ -1,29 +0,0 @@
User's Guide
############
.. toctree::
config
PLOVasp is a set of tools for performing DFT+DMFT calculations within
the fromalism of projected localized orbitals (PLO) using
VASP in combination with the TRIQS package.
It relies on projectors implemented in VASP >= 5.4.1 and generated with
the option LOCPROJ in INCAR. The main goal of PLOVasp is to use these raw
projectors generated by VASP to produce PLOs and to convert them to the standard
format of DFTTools, which allows a user to use the same DMFT scripts
as for Wien2TRIQS interface, with only the converter class being replaced.
Currently PLOVasp supports only one-shot DFT+DMFT calculations, with
the full charge self-consistency being in the process of implementation.
PLOVasp package consists of two main parts:
* A set of tools for processing VASP-generated projectors
* A converter class VaspConverter that prepares an HDF5 file for DFTTools
A typical one-shot DFT+DMFT calculation is generally performed as follows:
* Run VASP with LOCPROJ option to generate raw projectors
* Write an input file and run PLOVasp to produce data required by DFTTools
* Run a DMFT script based on DFTTools with VaspConverter class used as a converter
The format of the input file is described in section :doc:`Input Config-file <config>`.

39
doc/vasp/source/vaspio.rst
View File

@ -1,39 +0,0 @@
.. _vaspio:
VASP input-output
#################
The following VASP files are used by PLOtools:
* PROJCAR, LOCPROJ: raw projectors generated by VASP-PLO interface (VASP version >= 5.4.1)
* EIGENVAL: Kohn-Sham eigenvalues as well as `k`-points with weights and Fermi weights
* IBZKPT: `k`-point data (:math:`\Gamma`)
* POSCAR: crystal structure data
Starting from version 5.4.1 VASP now supports an official output of various types of
projectors that are requested in INCAR by specifying a set of sites, orbitals and types
of projectors. The calculated projectors are output into files **PROJCAR** and **LOCPROJ**.
The difference between these two files is that **LOCPROJ** contains raw matrices without
any reference to sites/orbitals, while **PROJCAR** is more detailed on that.
In particular, the information that can be obtained for each projector from **PROJCAR** is the following:
* site (and species) index
* for each `k`-point and band: a set of complex numbers for labeled orbitals
At the same time, **LOCPROJ** contains the total number of projectors (as well as the
number of `k`-points, bands, and spin channels) in the first line,
which can be used to allocate the arrays before parsing.
To enhance the performance of the parser, it is implemented in plain C. The idea is
that the python part of the parser first reads the first line of **LOCPROJ** and
then calls the C-routine with necessary parameters to parse **PROJCAR**.
The projectors are read in and stored in class `Plocar`. Two major data structures are stored:
* complex array `plo = nd.array((nproj, nspin, nk, nband))`
* list of projector descriptors `proj_params` containing the information on
the character of projectors
When a ProjectorShell is initialized it copies a subset of projectors corresponding
to selected sites/orbitals. This can be done by looping all shell sites/orbitals and
searching for the corresponding projector using the descriptor list `proj_params`.

11
python/converters/plovasp/converter.py
View File

@ -24,14 +24,16 @@
#
################################################################################
r"""
vasp.main
=========
plovasp.converter
=================
Main script of PLOVasp.
PLOVASP is a tool to transform raw, non-normalized
projectors generated by VASP into normalized projectors
corresponding to user-defined projected localized orbitals (PLO).
Runs routines in proper order to generate and store PLOs.
Usage: python main.py <conf-file> [<path-to-vasp-calcultaion>]
Usage: python converter.py <conf-file> [<path-to-vasp-calculation>]
"""
import sys
import vaspio
@ -85,4 +87,3 @@ def main():
if __name__ == '__main__':
main()

28
python/converters/plovasp/elstruct.py
View File

@ -24,8 +24,8 @@
#
################################################################################
r"""
vasp.elstruct
=============
plovasp.elstruct
================
Internal representation of VASP electronic structure data.
"""
@ -105,7 +105,7 @@ class ElectronicStructure:
print "eigvals from LOCPROJ"
self.eigvals = vasp_data.plocar.eigs
self.ferw = vasp_data.plocar.ferw.transpose((2, 0, 1))
self.efermi = vasp_data.plocar.efermi
self.efermi = vasp_data.doscar.efermi
# For later use it is more convenient to use a different order of indices
# [see ProjectorGroup.orthogonalization()]
@ -123,6 +123,21 @@ class ElectronicStructure:
# Concatenate coordinates grouped by type into one array
self.structure['qcoords'] = np.vstack(vasp_data.poscar.q_types)
self.structure['type_of_ion'] = vasp_data.poscar.type_of_ion
a = []
for ia in range(3):
a.append( vasp_data.poscar.a_brav[:,ia])
vol = np.dot(a[0],np.cross(a[1],a[2]))
b1 = 2.0*np.pi*np.cross(a[1],a[2])/vol
b2 = 2.0*np.pi*np.cross(a[2],a[0])/vol
b3 = 2.0*np.pi*np.cross(a[0],a[1])/vol
b = [b1,b2,b3]
self.kmesh['kpoints_cart'] = 0.0 * self.kmesh['kpoints']
for ik in range(self.nktot):
for ii in range(3):
self.kmesh['kpoints_cart'][ik] += self.kmesh['kpoints'][ik,ii]*b[ii]
# FIXME: This can be removed if ion coordinates are stored in a continuous array
## Construct a map to access coordinates by index
# self.structure['ion_index'] = []
@ -159,8 +174,9 @@ class ElectronicStructure:
# ov_min = np.minimum(ov, ov_min)
# Output only the site-diagonal parts of the matrices
print
print " Unorthonormalized density matrices and overlaps:"
for ispin in xrange(ns):
print
print " Spin:", ispin + 1
for io, ion in enumerate(ions):
print " Site:", ion
@ -173,11 +189,9 @@ class ElectronicStructure:
dm[iorb, iorb2] = den_mat[ispin, ind, ind2]
ov[iorb, iorb2] = overlap[ispin, ind, ind2]
print " Density matrix" + (12*norb - 12)*" " + "Overlap"
print " Density matrix" + (12*norb - 12 + 2)*" " + "Overlap"
for drow, dov in zip(dm, ov):
out = ''.join(map("{0:12.7f}".format, drow))
out += " "
out += ''.join(map("{0:12.7f}".format, dov))
print out

9
python/converters/plovasp/examples/ce/.gitignore vendored
View File

@ -1,9 +0,0 @@
*
!.gitignore
!INCAR
!POSCAR
!POTCAR
!KPOINTS
!plo.cfg
!*.sh
!*.py

49
python/converters/plovasp/examples/ce/INCAR
View File

@ -1,49 +0,0 @@
ADDGRID = TRUE
AMIX = 0.1
AMIN = 0.01
BMIX = 0.0001
#IMIX = 1
#MIXFIRST = True
ALGO = NORMAL
EDIFF = 0.00001
EDIFFG = -0.002
EMAX = 13.0
EMIN = -2.0
ENCUT = 600
ENAUG = 1100
IBRION = -1
ICHARG = 5
ISIF = 4
ISMEAR = 2
ISPIN = 2
ISYM = -1
LASPH = TRUE
LCHARG = TRUE
LDAU = True
LDAUJ = 0
LDAUL = 3
LDAUTYPE = 2
LDAUU = 0.0
LDAUPRINT = 2
LMAXMIX = 6
LORBIT = 13
LREAL = False
ROPT = 1e-3 1e-3 1e-3
LWAVE = TRUE
MAGMOM = 1.1
#NBANDS = 16
NEDOS = 501
NELM = 100
NELMDL = -7
NELMIN = 1
NPAR = 4
NSIM = 1
NSW = 0
POTIM = 0.4
PREC = Accurate
SIGMA = 0.1
SMASS = 0.5
SYMPREC = 1.0e-6
# Projectors
LOCPROJ = 1 : f : Pr 1

4
python/converters/plovasp/examples/ce/KPOINTS
View File

@ -1,4 +0,0 @@
Automatic kpoint scheme
0
Gamma
4 4 4

9
python/converters/plovasp/examples/ce/POSCAR
View File

@ -1,9 +0,0 @@
gamma-Ce
1.0
2.5790000000000000 2.5790000000000000 0.0000000000000000
0.0000000000000000 2.5790000000000000 2.5790000000000000
2.5790000000000000 0.0000000000000000 2.5790000000000000
Ce
1
Direct
0.0 0.0 0.0

3011
python/converters/plovasp/examples/ce/POTCAR
View File

File diff suppressed because it is too large Load Diff

252
python/converters/plovasp/examples/ce/U_matrix.py
View File

@ -1,252 +0,0 @@
from math import sqrt
from scipy.misc import factorial as fact
from itertools import product
import numpy as np
# The interaction matrix in desired basis
# U^{spherical}_{m1 m2 m3 m4} = \sum_{k=0}^{2l} F_k angular_matrix_element(l, k, m1, m2, m3, m4)
def U_matrix(l, radial_integrals=None, U_int=None, J_hund=None, basis="spherical", T=None):
"""Calculate U matrix being given either radial_integrals or U_int and J_hund.
l = angular momentum of shell being treated (l=2 for d shell, l=3 for f shell)
radial_integrals = [F0,F2,F4,..] (default None)
U_int, J_hund = values to use to compute radial_integrals (default None),
basis = "spherical", "cubic", or "other",
T = transformation matrix from spherical to desired basis, if basis='other' (default None)"""
# Check all necessary information is present and consistent
if radial_integrals is None and (U_int is None and J_hund is None):
raise ValueError("U_matrix: provide either the radial_integrals or U_int and J_hund.")
if radial_integrals is None and (U_int is not None and J_hund is not None):
radial_integrals = U_J_to_radial_integrals(l, U_int, J_hund)
if radial_integrals is not None and (U_int is not None and J_hund is not None):
if len(radial_integrals)-1 != l:
raise ValueError("U_matrix: inconsistency in l and number of radial_integrals provided.")
if (radial_integrals - U_J_to_radial_integrals(l, U_int, J_hund)).any() != 0.0:
print "Warning: U_matrix: radial_integrals provided do not match U_int and J_hund. Using radial_integrals to calculate spherical U_matrix."
# Full interaction matrix
# Basis of spherical harmonics Y_{-2}, Y_{-1}, Y_{0}, Y_{1}, Y_{2}
# U^{spherical}_{m1 m2 m3 m4} = \sum_{k=0}^{2l} F_k angular_matrix_element(l, k, m1, m2, m3, m4)
U_matrix = np.zeros((2*l+1,2*l+1,2*l+1,2*l+1),dtype=float)
m_range = range(-l,l+1)
for n, F in enumerate(radial_integrals):
k = 2*n
for m1, m2, m3, m4 in product(m_range,m_range,m_range,m_range):
U_matrix[m1+l,m2+l,m3+l,m4+l] += F * angular_matrix_element(l,k,m1,m2,m3,m4)
# Transform from spherical basis if needed
if basis == "cubic": T = spherical_to_cubic(l)
if basis == "other" and T is None:
raise ValueError("U_matrix: provide T for other bases.")
if T is not None: U_matrix = transform_U_matrix(U_matrix, T)
return U_matrix
# Convert full 4-index U matrix to 2-index density-density form
def reduce_4index_to_2index(U_4index):
"""Reduces the four-index matrix to two-index matrices."""
size = len(U_4index) # 2l+1
U = np.zeros((size,size),dtype=float) # matrix for same spin
Uprime = np.zeros((size,size),dtype=float) # matrix for opposite spin
m_range = range(size)
for m,mp in product(m_range,m_range):
U[m,mp] = U_4index[m,mp,m,mp].real - U_4index[m,mp,mp,m].real
Uprime[m,mp] = U_4index[m,mp,m,mp].real
return U, Uprime
# Construct the 2-index matrices for the density-density form
def U_matrix_kanamori(n_orb, U_int, J_hund):
"""Calculate the Kanamori U and Uprime matrices."""
U = np.zeros((n_orb,n_orb),dtype=float) # matrix for same spin
Uprime = np.zeros((n_orb,n_orb),dtype=float) # matrix for opposite spin
m_range = range(n_orb)
for m,mp in product(m_range,m_range):
if m == mp:
Uprime[m,mp] = U_int
else:
U[m,mp] = U_int - 3.0*J_hund
Uprime[m,mp] = U_int - 2.0*J_hund
return U, Uprime
# Get t2g or eg components
def t2g_submatrix(U, convention=''):
"""Return only the t2g part of the full d-manifold two- or four-index U matrix."""
if convention == 'wien2k':
return subarray(U, len(U.shape)*[(2,3,4)])
else:
return subarray(U, len(U.shape)*[(0,1,3)])
def eg_submatrix(U, convention=''):
"""Return only the eg part of the full d-manifold two- or four-index U matrix."""
if convention == 'wien2k':
return subarray(U, len(U.shape)*[(0,1)])
else:
return subarray(U, len(U.shape)*[(2,4)])
# Transform the interaction matrix into another basis
def transform_U_matrix(U_matrix, T):
"""Transform the interaction matrix into another basis by applying matrix T."""
return np.einsum("ij,kl,jlmo,mn,op",np.conj(T),np.conj(T),U_matrix,np.transpose(T),np.transpose(T))
# Rotation matrices: complex harmonics to cubic harmonics
# Complex harmonics basis: ..., Y_{-2}, Y_{-1}, Y_{0}, Y_{1}, Y_{2}, ...
def spherical_to_cubic(l, convention=''):
"""Returns the spherical harmonics to cubic harmonics rotation matrix."""
size = 2*l+1
T = np.zeros((size,size),dtype=complex)
if convention != 'wien2k' and l == 2:
raise ValueError("spherical_to_cubic: standard convention not implemented for l=1")
if l == 0:
cubic_names = ("s")
elif l == 1:
cubic_names = ("x","y","z")
T[0,0] = 1.0/sqrt(2); T[0,2] = -1.0/sqrt(2)
T[1,0] = 1j/sqrt(2); T[1,2] = 1j/sqrt(2)
T[2,1] = 1.0
elif l == 2:
if convention == 'wien2k':
cubic_names = ("z^2","x^2-y^2","xy","yz","xz")
T[0,2] = 1.0
T[1,0] = 1.0/sqrt(2); T[1,4] = 1.0/sqrt(2)
T[2,0] = -1j/sqrt(2); T[2,4] = 1j/sqrt(2)
T[3,1] = 1j/sqrt(2); T[3,3] = -1j/sqrt(2)
T[4,1] = 1.0/sqrt(2); T[4,3] = 1.0/sqrt(2)
else:
cubic_names = ("xy","yz","z^2","xz","x^2-y^2")
T[0,0] = 1j/sqrt(2); T[0,4] = -1j/sqrt(2)
T[1,1] = 1j/sqrt(2); T[1,3] = 1j/sqrt(2)
T[2,2] = 1.0
T[3,1] = 1.0/sqrt(2); T[3,3] = -1.0/sqrt(2)
T[4,0] = 1.0/sqrt(2); T[4,4] = 1.0/sqrt(2)
elif l == 3:
if convention == 'wien2k':
cubic_names = ("x(x^2-3y^2)","z(x^2-y^2)","xz^2","z^3","yz^2","xyz","y(3x^2-y^2)")
T[0,0] = 1.0/sqrt(2); T[0,6] = -1.0/sqrt(2)
T[1,1] = 1.0/sqrt(2); T[1,5] = 1.0/sqrt(2)
T[2,2] = 1.0/sqrt(2); T[2,4] = -1.0/sqrt(2)
T[3,5] = 1.0
T[4,2] = 1j/sqrt(2); T[4,4] = 1j/sqrt(2)
T[5,1] = 1j/sqrt(2); T[5,5] = -1j/sqrt(2)
T[6,0] = 1j/sqrt(2); T[6,6] = 1j/sqrt(2)
else:
cubic_names = ("y(3x^2-y^2)","xyz","yz^2","z^3","xz^2","z(x^2-y^2)","x(x^2-3y^2)")
# Verified with the matrix generated in VASP
T[0,0] = 1j/sqrt(2); T[0,6] = 1j/sqrt(2)
T[1,1] = 1j/sqrt(2); T[1,5] = -1j/sqrt(2)
T[2,2] = 1j/sqrt(2); T[2,4] = 1j/sqrt(2)
T[3,3] = 1.0
T[4,2] = 1.0/sqrt(2); T[4,4] = -1.0/sqrt(2)
T[5,1] = 1.0/sqrt(2); T[5,5] = 1.0/sqrt(2)
T[6,0] = 1.0/sqrt(2); T[6,6] = -1.0/sqrt(2)
else: raise ValueError("spherical_to_cubic: implemented only for l=0,1,2,3")
return T
# Names of cubic harmonics
def cubic_names(l):
if l == 0 or l == 's':
return ("s")
elif l == 1 or l == 'p':
return ("x","y","z")
elif l == 2 or l == 'd':
return ("xy","yz","z^2","xz","x^2-y^2")
elif l == 't2g':
return ("xy","yz","xz")
elif l == 'eg':
return ("z^2","x^2-y^2")
elif l == 3 or l == 'f':
return ("x(x^2-3y^2)","z(x^2-y^2)","xz^2","z^3","yz^2","xyz","y(3x^2-y^2)")
else: raise ValueError("cubic_names: implemented only for l=0,1,2,3")
# Convert U,J -> radial integrals F_k
def U_J_to_radial_integrals(l, U_int, J_hund):
"""Determines the radial integrals F_k from U_int and J_hund."""
F = np.zeros((l+1),dtype=float)
if l == 2:
F[0] = U_int
F[1] = J_hund * 14.0 / (1.0 + 0.63)
F[2] = 0.630 * F[1]
elif l == 3:
F[0] = U_int
F[1] = 6435.0 * J_hund / (286.0 + 195.0 * 451.0 / 675.0 + 250.0 * 1001.0 / 2025.0)
F[2] = 451.0 * F[1] / 675.0
F[3] = 1001.0 * F[1] / 2025.0
else: raise ValueError("U_J_to_radial_integrals: implemented only for l=2,3")
return F
# Convert radial integrals F_k -> U,J
def radial_integrals_to_U_J(l, F):
"""Determines U_int and J_hund from the radial integrals."""
if l == 2:
U_int = F[0]
J_hund = F[1] * (1.0 + 0.63) / 14.0
elif l == 3:
U_int = F[0]
J_hund = F[1] * (286.0 + 195.0 * 451.0 / 675.0 + 250.0 * 1001.0 / 2025.0) / 6435.0
else: raise ValueError("radial_integrals_to_U_J: implemented only for l=2,3")
return U_int,J_hund
# Angular matrix elements of particle-particle interaction
# (2l+1)^2 ((l 0) (k 0) (l 0))^2 \sum_{q=-k}^{k} (-1)^{m1+m2+q} ((l -m1) (k q) (l m3)) ((l -m2) (k -q) (l m4))
def angular_matrix_element(l, k, m1, m2, m3, m4):
result = 0
for q in range(-k,k+1):
result += three_j_symbol((l,-m1),(k,q),(l,m3))*three_j_symbol((l,-m2),(k,-q),(l,m4))*(-1.0 if (m1+q+m2) % 2 else 1.0)
result *= (2*l+1)**2 * (three_j_symbol((l,0),(k,0),(l,0))**2)
return result
# Wigner 3-j symbols
# ((j1 m1) (j2 m2) (j3 m3))
def three_j_symbol(jm1, jm2, jm3):
j1, m1 = jm1
j2, m2 = jm2
j3, m3 = jm3
if (m1+m2+m3 != 0 or
m1 < -j1 or m1 > j1 or
m2 < -j2 or m2 > j2 or
m3 < -j3 or m3 > j3 or
j3 > j1 + j2 or
j3 < abs(j1-j2)):
return .0
result = -1.0 if (j1-j2-m3) % 2 else 1.0
result *= sqrt(fact(j1+j2-j3)*fact(j1-j2+j3)*fact(-j1+j2+j3)/fact(j1+j2+j3+1))
result *= sqrt(fact(j1-m1)*fact(j1+m1)*fact(j2-m2)*fact(j2+m2)*fact(j3-m3)*fact(j3+m3))
t_min = max(j2-j3-m1,j1-j3+m2,0)
t_max = min(j1-m1,j2+m2,j1+j2-j3)
t_sum = 0
for t in range(t_min,t_max+1):
t_sum += (-1.0 if t % 2 else 1.0)/(fact(t)*fact(j3-j2+m1+t)*fact(j3-j1-m2+t)*fact(j1+j2-j3-t)*fact(j1-m1-t)*fact(j2+m2-t))
result *= t_sum
return result
# Clebsch-Gordan coefficients
# < j1 m1 j2 m2 | j3 m3 > = (-1)^{j1-j2+m3} \sqrt{2j3+1} ((j1 m1) (j2 m2) (j3 -m3))
def clebsch_gordan(jm1, jm2, jm3):
norm = sqrt(2*jm3[0]+1)*(-1 if jm1[0]-jm2[0]+jm3[1] % 2 else 1)
return norm*three_j_symbol(jm1,jm2,(jm3[0],-jm3[1]))
# Create subarray containing columns in idxlist
# e.g. idxlist = [(0),(2,3),(0,1,2,3)] gives
# column 0 for 1st dim,
# columns 2 and 3 for 2nd dim,
# columns 0,1,2 and 3 for 3rd dim.
def subarray(a,idxlist,n=None) :
if n is None: n = len(a.shape)-1
sa = a[tuple(slice(x) for x in a.shape[:n]) + (idxlist[n],)]
return subarray(sa,idxlist, n-1) if n > 0 else sa

335
python/converters/plovasp/examples/ce/hf_solver.py
View File

@ -1,335 +0,0 @@
################################################################################
#
# TRIQS: a Toolbox for Research in Interacting Quantum Systems
#
# Copyright (C) 2011 by M. Ferrero, O. Parcollet
#
# TRIQS is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.
#
# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# TRIQS. If not, see <http://www.gnu.org/licenses/>.
#
################################################################################
from types import *
#from triqs_dft_tools.U_matrix import *
from U_matrix import *
from pytriqs.gf import *
#from hubbard_I import gf_hi_fullu, sigma_atomic_fullu
import pytriqs.utility.mpi as mpi
from itertools import izip
import numpy
import copy
class Solver:
"""
Hartree-Fock Solver
"""
# initialisation:
def __init__(self, beta, l, n_msb=1025, use_spin_orbit=False, Nmoments=5, dudarev=False):
self.name = "Hartree-Fock"
self.beta = beta
self.l = l
self.Nmsb = n_msb
self.UseSpinOrbit = use_spin_orbit
self.Converged = False
self.Nspin = 2
self.Nmoments=Nmoments
self.dudarev = dudarev
assert use_spin_orbit == False, "Spin-orbit is not implemented"
self.Nlm = 2*l+1
if (use_spin_orbit):
# no blocks!
self.gf_struct = [ ('ud', range(2*self.Nlm)) ]
else:
# up/down blocks:
self.gf_struct = [ ('up', range(self.Nlm)), ('down', range(self.Nlm)) ]
# construct Greens functions:
self.a_list = [a for a,al in self.gf_struct]
glist = lambda : [ GfImFreq(indices = al, beta = self.beta, n_points = self.Nmsb) for a,al in self.gf_struct]
self.G = BlockGf(name_list = self.a_list, block_list = glist(),make_copies=False)
self.G_Old = self.G.copy()
self.G0 = self.G.copy()
self.Sigma = self.G.copy()
self.Sigma_Old = self.G.copy()
def solve(self, U_int, J_hund, T=None, verbosity=0, Iteration_Number=1, Test_Convergence=0.0001):
"""Calculation of the impurity Greens function using Hubbard-I"""
if self.Converged :
mpi.report("Solver %(name)s has already converged: SKIPPING"%self.__dict__)
return
if mpi.is_master_node():
self.verbosity = verbosity
else:
self.verbosity = 0
#self.Nmoments = 5
ur,ujmn,umn=self.__set_umatrix(U=U_int,J=J_hund,T=T)
M = [x for x in self.G.mesh]
self.zmsb = numpy.array([x for x in M],numpy.complex_)
# # for the tails:
# tailtempl={}
# for sig,g in self.G:
# tailtempl[sig] = copy.deepcopy(g.tail)
# for i in range(9): tailtempl[sig][i] *= 0.0
# self.__save_eal('eal.dat',Iteration_Number)
# mpi.report( "Starting Fortran solver %(name)s"%self.__dict__)
self.Sigma_Old <<= self.Sigma
self.G_Old <<= self.G
# # call the fortran solver:
# temp = 1.0/self.beta
# gf,tail,self.atocc,self.atmag = gf_hi_fullu(e0f=self.ealmat, ur=ur, umn=umn, ujmn=ujmn,
# zmsb=self.zmsb, nmom=self.Nmoments, ns=self.Nspin, temp=temp, verbosity = self.verbosity)
#self.sig = sigma_atomic_fullu(gf=self.gf,e0f=self.eal,zmsb=self.zmsb,ns=self.Nspin,nlm=self.Nlm)
def print_matrix(m):
for row in m:
print ''.join(map("{0:12.7f}".format, row))
# Hartree-Fock solver
self.Sigma.zero()
dm = self.G.density()
if mpi.is_master_node():
# print
# print " Reduced U-matrix:"
# print " U:"
# print_matrix(ujmn)
# print " Up:"
# print_matrix(umn)
#
## sig_test = {bl: numpy.zeros((self.Nlm, self.Nlm)) for bl in dm}
# sig_test = {}
# sig_test['up'] = numpy.dot(umn, dm['up'].real) + numpy.dot(ujmn, dm['down'].real)
# sig_test['down'] = numpy.dot(umn, dm['down'].real) + numpy.dot(ujmn, dm['up'].real)
# print " Sigma test:"
# print_matrix(sig_test['up'])
print
print " Density matrix (up):"
print_matrix(dm['up'])
print
print " Density matrix (down):"
print_matrix(dm['down'])
if self.dudarev:
Ueff = U_int - J_hund
corr_energy = 0.0
dft_dc = 0.0
for bl1 in dm:
# (U - J) * (1/2 - n)
self.Sigma[bl1] << Ueff * (0.5 * numpy.identity(self.Nlm) - dm[bl1])
# 1/2 (U - J) * \sum_{\sig} [\sum_{m} n_{m,m \sig} - \sum_{m1,m2} n_{m1,m2 \sig} n_{m2,m1 \sig}]
corr_energy += 0.5 * Ueff * (dm[bl1].trace() - (dm[bl1]*dm[bl1].conj()).sum()).real
# V[n] * n^{\dagger}
dft_dc += (self.Sigma[bl1](0) * dm[bl1].conj()).sum().real
else:
# !!!!!
# !!!!! Mind the order of indices in the 4-index matrix!
# !!!!!
for il1 in xrange(self.Nlm):
for il2 in xrange(self.Nlm):
for il3 in xrange(self.Nlm):
for il4 in xrange(self.Nlm):
for bl1 in dm:
for bl2 in dm:
self.Sigma[bl1][il1, il2] += ur[il1, il3, il2, il4] * dm[bl2][il3, il4]
if bl1 == bl2:
self.Sigma[bl1][il1, il2] -= ur[il1, il3, il4, il2] * dm[bl1][il3, il4]
if mpi.is_master_node() and self.verbosity > 0:
print
print " Sigma (up):"
print_matrix(self.Sigma['up'](0).real)
print
print " Sigma (down):"
print_matrix(self.Sigma['down'](0).real)
# if (self.verbosity==0):
# # No fortran output, so give basic results here
# mpi.report("Atomic occupancy in Hubbard I Solver : %s"%self.atocc)
# mpi.report("Atomic magn. mom. in Hubbard I Solver : %s"%self.atmag)
# transfer the data to the GF class:
if (self.UseSpinOrbit):
nlmtot = self.Nlm*2 # only one block in this case!
else:
nlmtot = self.Nlm
# M={}
# isp=-1
# for a,al in self.gf_struct:
# isp+=1
# M[a] = numpy.array(gf[isp*nlmtot:(isp+1)*nlmtot,isp*nlmtot:(isp+1)*nlmtot,:]).transpose(2,0,1).copy()
# for i in range(min(self.Nmoments,8)):
# tailtempl[a][i+1] = tail[i][isp*nlmtot:(isp+1)*nlmtot,isp*nlmtot:(isp+1)*nlmtot]
#
# #glist = lambda : [ GfImFreq(indices = al, beta = self.beta, n_points = self.Nmsb, data =M[a], tail =self.tailtempl[a])
# # for a,al in self.gf_struct]
# glist = lambda : [ GfImFreq(indices = al, beta = self.beta, n_points = self.Nmsb) for a,al in self.gf_struct]
# self.G = BlockGf(name_list = self.a_list, block_list = glist(),make_copies=False)
#
# self.__copy_Gf(self.G,M,tailtempl)
#
# # Self energy:
# self.G0 <<= iOmega_n
#
# M = [ self.ealmat[isp*nlmtot:(isp+1)*nlmtot,isp*nlmtot:(isp+1)*nlmtot] for isp in range((2*self.Nlm)/nlmtot) ]
# self.G0 -= M
# self.Sigma <<= self.G0 - inverse(self.G)
#
# # invert G0
# self.G0.invert()
def test_distance(G1,G2, dist) :
def f(G1,G2) :
#print abs(G1.data - G2.data)
dS = max(abs(G1.data - G2.data).flatten())
aS = max(abs(G1.data).flatten())
if mpi.is_master_node():
print " Distances:", dS, " vs ", aS * dist
return dS <= aS*dist
return reduce(lambda x,y : x and y, [f(g1,g2) for (i1,g1),(i2,g2) in izip(G1,G2)])
mpi.report("\nChecking Sigma for convergence...\nUsing tolerance %s"%Test_Convergence)
self.Converged = test_distance(self.Sigma,self.Sigma_Old,Test_Convergence)
if self.Converged :
mpi.report("Solver HAS CONVERGED")
else :
mpi.report("Solver has not yet converged")
return corr_energy, dft_dc
def GF_realomega(self, ommin, ommax, N_om, U_int, J_hund, T=None, verbosity=0, broadening=0.01):
"""Calculates the GF and spectral function on the real axis."""
delta_om = (ommax-ommin)/(1.0*(N_om-1))
omega = numpy.zeros([N_om],numpy.complex_)
ur,umn,ujmn=self.__set_umatrix(U=U_int,J=J_hund,T=T)
for i in range(N_om):
omega[i] = ommin + delta_om * i + 1j * broadening
tailtempl={}
for sig,g in self.G:
tailtempl[sig] = copy.deepcopy(g.tail)
for i in range(9): tailtempl[sig][i] *= 0.0
temp = 1.0/self.beta
gf,tail,self.atocc,self.atmag = gf_hi_fullu(e0f=self.ealmat, ur=ur, umn=umn, ujmn=ujmn,
zmsb=omega, nmom=self.Nmoments, ns=self.Nspin, temp=temp, verbosity = verbosity)
# transfer the data to the GF class:
if (self.UseSpinOrbit):
nlmtot = self.Nlm*2 # only one block in this case!
else:
nlmtot = self.Nlm
M={}
isp=-1
for a,al in self.gf_struct:
isp+=1
M[a] = numpy.array(gf[isp*nlmtot:(isp+1)*nlmtot,isp*nlmtot:(isp+1)*nlmtot,:]).transpose(2,0,1).copy()
for i in range(min(self.Nmoments,8)):
tailtempl[a][i+1] = tail[i][isp*nlmtot:(isp+1)*nlmtot,isp*nlmtot:(isp+1)*nlmtot]
#glist = lambda : [ GfReFreq(indices = al, window = (ommin, ommax), n_points = N_om, data = M[a], tail = self.tailtempl[a])
# for a,al in self.gf_struct] # Indices for the upfolded G
glist = lambda : [ GfReFreq(indices = al, window = (ommin, ommax), n_points = N_om) for a,al in self.gf_struct]
self.G = BlockGf(name_list = self.a_list, block_list = glist(),make_copies=False)
self.__copy_Gf(self.G,M,tailtempl)
# Self energy:
self.G0 = self.G.copy()
self.Sigma = self.G.copy()
self.G0 <<= Omega + 1j*broadening
M = [ self.ealmat[isp*nlmtot:(isp+1)*nlmtot,isp*nlmtot:(isp+1)*nlmtot] for isp in range((2*self.Nlm)/nlmtot) ]
self.G0 -= M
self.Sigma <<= self.G0 - inverse(self.G)
self.Sigma.note='ReFreq' # This is important for the put_Sigma routine!!!
#sigmamat = sigma_atomic_fullu(gf=gf,e0f=self.ealmat,zmsb=omega,nlm=self.Nlm,ns=self.Nspin)
#return omega,gf,sigmamat
def __save_eal(self,Filename,it):
f=open(Filename,'a')
f.write('\neff. atomic levels, Iteration %s\n'%it)
for i in range(self.Nlm*self.Nspin):
for j in range(self.Nlm*self.Nspin):
f.write("%10.6f %10.6f "%(self.ealmat[i,j].real,self.ealmat[i,j].imag))
f.write("\n")
f.close()
def __copy_Gf(self,G,data,tail):
""" Copies data and tail to Gf object GF """
for s,g in G:
g.data[:,:,:]=data[s][:,:,:]
for imom in range(1,min(self.Nmoments,8)):
g.tail.data[1+imom,:,:]=tail[s][imom]
def set_atomic_levels(self,eal):
""" Helps to set correctly the variables for the atomic levels from a dictionary."""
assert (type(eal)==DictType), "Give a dictionary to set_atomic_levels!"
cnt = 0
self.ealmat[:,:] *= 0.0
for ind in eal:
self.Eff_Atomic_Levels[ind] = copy.deepcopy(eal[ind])
if self.UseSpinOrbit:
for ii in range(self.Nlm*2):
for jj in range(self.Nlm*2):
self.ealmat[ii,jj] = self.Eff_Atomic_Levels[ind][ii,jj]
else:
for ii in range(self.Nlm):
for jj in range(self.Nlm):
self.ealmat[cnt*self.Nlm + ii,cnt*self.Nlm + jj] = self.Eff_Atomic_Levels[ind][ii,jj]
cnt += 1
def __set_umatrix(self,U,J,T=None):
# U matrix:
# l = (Nlm-1)/2
# If T is specified, it is used to transform the Basis set
Umat = U_matrix(l=self.l, U_int=U, J_hund=J, basis='cubic', T=T)
U, Up = reduce_4index_to_2index(Umat)
return Umat, U, Up

15
python/converters/plovasp/examples/ce/plo.cfg
View File

@ -1,15 +0,0 @@
[General]
#EFERMI =
#BASENAME = plo_vasp
#DOSMESH = -5.0 10.0 4001
[Group 1]
SHELLS = 1
NORMALIZE = False
EWINDOW = -41.45 41.8
[Shell 1]
LSHELL = 3
IONS = 1

2
python/converters/plovasp/examples/ce/run_build.sh
View File

@ -1,2 +0,0 @@
#!/bin/bash
PYTHONPATH=$HOME/Codes/triqs/triqs1.2/applications/dft_tools/dev_src/python:$HOME/Codes/triqs/triqs1.2/applications/dft_tools/dev_src/c:$PYTHONPATH pytriqs $@

2
python/converters/plovasp/examples/ce/run_plovasp.sh
View File

@ -1,2 +0,0 @@
PLOVASP_PATH=$HOME/Codes/triqs/triqs1.2/applications/dft_tools/dev_src
PYTHONPATH=$PLOVASP_PATH/python:$PLOVASP_PATH/c:$PYTHONPATH pytriqs $PLOVASP_PATH/python/converters/plovasp/main.py $@

164
python/converters/plovasp/examples/ce/sc_dmft.py
View File

@ -1,164 +0,0 @@
import os
import errno
import signal
import sys
import time
import pytriqs.utility.mpi as mpi
import converters.plovasp.main as plovasp
from test_ham_hf import dmft_cycle
debug = True
#
# Helper functions
#
def sigint_handler(signal, frame):
raise SystemExit(1)
def is_vasp_lock_present():
return os.path.isfile('./vasp.lock')
def is_vasp_running(vasp_pid):
"""
Tests if VASP initial process is still alive.
"""
pid_exists = False
if mpi.is_master_node():
try:
os.kill(vasp_pid, 0)
except OSError, e:
pid_exists = e.errno == errno.EPERM
else:
pid_exists = True
pid_exists = mpi.bcast(pid_exists)
return pid_exists
def get_dft_energy():
"""
Reads energy from the last line of OSZICAR.
"""
with open('OSZICAR', 'r') as f:
nextline = f.readline()
while nextline.strip():
line = nextline
nextline = f.readline()
# print "OSZICAR: ", line[:-1]
try:
dft_energy = float(line.split()[2])
except ValueError:
print "Cannot read energy from OSZICAR, setting it to zero"
dft_energy = 0.0
return dft_energy
class bcolors:
MAGENTA = '\033[95m'
BLUE = '\033[94m'
GREEN = '\033[92m'
YELLOW = '\033[93m'
RED = '\033[91m'
ENDC = '\033[0m'
#
# Main self-consistent cycle
#
def run_all(vasp_pid):
"""
"""
mpi.report(" Waiting for VASP lock to appear...")
while not is_vasp_lock_present():
time.sleep(1)
vasp_running = True
while vasp_running:
if debug: print bcolors.RED + "rank %s"%(mpi.rank) + bcolors.ENDC
mpi.report(" Waiting for VASP lock to disappear...")
mpi.barrier()
while is_vasp_lock_present():
time.sleep(1)
# if debug: print bcolors.YELLOW + " waiting: rank %s"%(mpi.rank) + bcolors.ENDC
if not is_vasp_running(vasp_pid):
mpi.report(" VASP stopped")
vasp_running = False
break
if debug: print bcolors.MAGENTA + "rank %s"%(mpi.rank) + bcolors.ENDC
err = 0
exc = None
try:
if debug: print bcolors.BLUE + "plovasp: rank %s"%(mpi.rank) + bcolors.ENDC
if mpi.is_master_node():
plovasp.generate_and_output_as_text('plo.cfg', vasp_dir='./')
# Read energy from OSZICAR
dft_energy = get_dft_energy()
except Exception, exc:
err = 1
err = mpi.bcast(err)
if err:
if mpi.is_master_node():
raise exc
else:
raise SystemExit(1)
mpi.barrier()
try:
if debug: print bcolors.GREEN + "rank %s"%(mpi.rank) + bcolors.ENDC
corr_energy, dft_dc = dmft_cycle()
except:
if mpi.is_master_node():
print " master forwarding the exception..."
raise
else:
print " rank %i exiting..."%(mpi.rank)
raise SystemExit(1)
mpi.barrier()
if mpi.is_master_node():
total_energy = dft_energy + corr_energy - dft_dc
print
print "="*80
print " Total energy: ", total_energy
print " DFT energy: ", dft_energy
print " Corr. energy: ", corr_energy
print " DFT DC: ", dft_dc
print "="*80
print
if mpi.is_master_node() and vasp_running:
open('./vasp.lock', 'a').close()
if mpi.is_master_node():
total_energy = dft_energy + corr_energy - dft_dc
with open('TOTENERGY', 'w') as f:
f.write(" Total energy: %s\n"%(total_energy))
f.write(" DFT energy: %s\n"%(dft_energy))
f.write(" Corr. energy: %s\n"%(corr_energy))
f.write(" DFT DC: %s\n"%(dft_dc))
f.write(" Energy correction: %s\n"%(corr_energy - dft_dc))
mpi.report("***Done")
if __name__ == '__main__':
try:
vasp_pid = int(sys.argv[1])
except (ValueError, KeyError):
if mpi.is_master_node():
print "VASP process pid must be provided as the first argument"
raise
# if len(sys.argv) > 1:
# vasp_path = sys.argv[1]
# else:
# try:
# vasp_path = os.environ['VASP_DIR']
# except KeyError:
# if mpi.is_master_node():
# print "Path to VASP must be specified either as an argument of in VASP_DIR"
# raise
signal.signal(signal.SIGINT, sigint_handler)
run_all(vasp_pid)

9
python/converters/plovasp/examples/ce/sc_dmft.sh
View File

@ -1,9 +0,0 @@
#/bin/bash
NPROC=8
rm -f vasp.lock
stdbuf -o 0 mpirun -np $NPROC ~/Codes/vasp/build_20151202/vasp.5.4.2/build/std/vasp &
mpirun -np $NPROC ./run_build.sh sc_dmft.py $(jobs -p) || kill %1

28
python/converters/plovasp/examples/ce/subm.sh
View File

@ -1,28 +0,0 @@
#!/bin/bash -l
#
#SBATCH --job-name="JOBNAME"
#SBATCH --time=36:00:00
##SBATCH --nodes=8
#SBATCH --ntasks=8
#SBATCH --partition=qmac
##SBATCH --mem-per-cpu=1024
#SBATCH --output=job.%j.o
#SBATCH --error=job.%j.e
#======START=====
# source $HOME/.profile
# export KMP_AFFINITY=none
# module load slurm
source /cm/shared/apps/intel/composer_xe/composer_xe_2015.0.090/mkl/bin/mklvars.sh intel64
# module switch PrgEnv-gnu PrgEnv-intel
# module load intel/14.0.2.144
#module load intel/mkl/64/11.2/2015.0.090
export OMP_NUM_THREADS=1
VASP_DIR=$HOME/Codes/vasp/build_20151202/vasp.5.4.2/build/std
#mpirun -np $SLURM_NTASKS $VASP_DIR/vasp
./sc_dmft.sh

437
python/converters/plovasp/examples/ce/test_ham_hf.py
View File

@ -1,437 +0,0 @@
#from triqs_dft_tools.sumk_dft import *
from sumk_dft import *
#from triqs_dft_tools.converters.wien2k_converter import *
from converters.vasp_converter import *
#from pytriqs.applications.impurity_solvers.hubbard_I.hubbard_solver import Solver
from hf_solver import Solver
import shutil
class TestSumkDFT(SumkDFT):
# calculate and save occupancy matrix in the Bloch basis for VASP charge denity recalculation
def calc_density_correction(self, filename='GAMMA', dm_type='wien2k'):
r"""
Calculates the charge density correction and stores it into a file.
The charge density correction is needed for charge-self-consistent DFT+DMFT calculations.
It represents a density matrix of the interacting system defined in Bloch basis
and it is calculated from the sum over Matsubara frequecies of the full GF,
..math:: N_{\nu\nu'}(k) = \sum_{i\omega_{n}} G_{\nu\nu'}(k, i\omega_{n})
The density matrix for every `k`-point is stored into a file.
Parameters
----------
filename : string
Name of the file to store the charge density correction.
Returns
-------
(deltaN, dens) : tuple
Returns a tuple containing the density matrix `deltaN` and
the corresponing total charge `dens`.
"""
assert type(filename) == StringType, "calc_density_correction: filename has to be a string!"
assert dm_type in ('vasp', 'wien2k'), "'type' must be either 'vasp' or 'wienk'"
ntoi = self.spin_names_to_ind[self.SO]
spn = self.spin_block_names[self.SO]
dens = {sp: 0.0 for sp in spn}
# Fetch Fermi weights and energy window band indices
if dm_type == 'vasp':
fermi_weights = 0
band_window = 0
if mpi.is_master_node():
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)
# Convert Fermi weights to a density matrix
dens_mat_dft = {}
for sp in spn:
dens_mat_dft[sp] = [fermi_weights[ik, ntoi[sp], :].astype(numpy.complex_) for ik in xrange(self.n_k)]
# Set up deltaN:
deltaN = {}
for sp in spn:
deltaN[sp] = [numpy.zeros([self.n_orbitals[ik,ntoi[sp]],self.n_orbitals[ik,ntoi[sp]]], numpy.complex_) for ik in range(self.n_k)]
ikarray = numpy.array(range(self.n_k))
for ik in mpi.slice_array(ikarray):
G_latt_iw = self.lattice_gf(ik = ik, mu = self.chemical_potential, iw_or_w = "iw")
for bname,gf in G_latt_iw:
deltaN[bname][ik][:, :] = G_latt_iw[bname].density()
dens[bname] += self.bz_weights[ik] * G_latt_iw[bname].total_density()
if dm_type == 'vasp':
# In 'vasp'-mode subtract the DFT density matrix
diag_inds = numpy.diag_indices(self.n_orbitals[ik, ntoi[bname]])
deltaN[bname][ik][diag_inds] -= dens_mat_dft[bname][ik]
dens[bname] -= self.bz_weights[ik] * dens_mat_dft[bname][ik].sum().real
# mpi reduce:
for bname in deltaN:
for ik in range(self.n_k):
deltaN[bname][ik] = mpi.all_reduce(mpi.world, deltaN[bname][ik], lambda x,y : x+y)
dens[bname] = mpi.all_reduce(mpi.world, dens[bname], lambda x,y : x+y)
mpi.barrier()
# now save to file:
if dm_type == 'wien2k':
if mpi.is_master_node():
if self.SP == 0:
f = open(filename,'w')
else:
f = open(filename+'up','w')
f1 = open(filename+'dn','w')
# write chemical potential (in Rydberg):
f.write("%.14f\n"%(self.chemical_potential/self.energy_unit))
if self.SP != 0: f1.write("%.14f\n"%(self.chemical_potential/self.energy_unit))
# write beta in rydberg-1
f.write("%.14f\n"%(G_latt_iw.mesh.beta*self.energy_unit))
if self.SP != 0: f1.write("%.14f\n"%(G_latt_iw.mesh.beta*self.energy_unit))
if self.SP == 0: # no spin-polarization
for ik in range(self.n_k):
f.write("%s\n"%self.n_orbitals[ik,0])
for inu in range(self.n_orbitals[ik,0]):
for imu in range(self.n_orbitals[ik,0]):
valre = (deltaN['up'][ik][inu,imu].real + deltaN['down'][ik][inu,imu].real) / 2.0
valim = (deltaN['up'][ik][inu,imu].imag + deltaN['down'][ik][inu,imu].imag) / 2.0
f.write("%.14f %.14f "%(valre,valim))
f.write("\n")
f.write("\n")
f.close()
elif self.SP == 1: # with spin-polarization
# dict of filename: (spin index, block_name)
if self.SO == 0: to_write = {f: (0, 'up'), f1: (1, 'down')}
if self.SO == 1: to_write = {f: (0, 'ud'), f1: (0, 'ud')}
for fout in to_write.iterkeys():
isp, sp = to_write[fout]
for ik in range(self.n_k):
fout.write("%s\n"%self.n_orbitals[ik,isp])
for inu in range(self.n_orbitals[ik,isp]):
for imu in range(self.n_orbitals[ik,isp]):
fout.write("%.14f %.14f "%(deltaN[sp][ik][inu,imu].real,deltaN[sp][ik][inu,imu].imag))
fout.write("\n")
fout.write("\n")
fout.close()
elif dm_type == 'vasp':
# assert self.SP == 0, "Spin-polarized density matrix is not implemented"
if mpi.is_master_node():
with open(filename, 'w') as f:
f.write(" %i -1 ! Number of k-points, default number of bands\n"%(self.n_k))
for sp in spn:
for ik in xrange(self.n_k):
ib1 = band_window[0][ik, 0]
ib2 = band_window[0][ik, 1]
f.write(" %i %i %i\n"%(ik + 1, ib1, ib2))
for inu in xrange(self.n_orbitals[ik, 0]):
for imu in xrange(self.n_orbitals[ik, 0]):
if self.SP == 0:
valre = (deltaN['up'][ik][inu, imu].real + deltaN['down'][ik][inu, imu].real) / 2.0
valim = (deltaN['up'][ik][inu, imu].imag + deltaN['down'][ik][inu, imu].imag) / 2.0
else:
valre = deltaN[sp][ik][inu, imu].real
valim = deltaN[sp][ik][inu, imu].imag
f.write(" %.14f %.14f"%(valre, valim))
f.write("\n")
else:
raise NotImplementedError("Unknown density matrix type: '%s'"%(dm_type))
return deltaN, dens
def calc_hamiltonian_correction(self, filename='GAMMA'):
r"""
Calculates the charge density correction and stores it into a file.
The charge density correction is needed for charge-self-consistent DFT+DMFT calculations.
It represents a density matrix of the interacting system defined in Bloch basis
and it is calculated from the sum over Matsubara frequecies of the full GF,
..math:: N_{\nu\nu'}(k) = \sum_{i\omega_{n}} G_{\nu\nu'}(k, i\omega_{n})
The density matrix for every `k`-point is stored into a file.
Parameters
----------
filename : string
Name of the file to store the charge density correction.
Returns
-------
(deltaN, dens) : tuple
Returns a tuple containing the density matrix `deltaN` and
the corresponing total charge `dens`.
"""
assert type(filename) == StringType, "calc_density_correction: filename has to be a string!"
ntoi = self.spin_names_to_ind[self.SO]
spn = self.spin_block_names[self.SO]
dens = {sp: 0.0 for sp in spn}
# Fetch Fermi weights and energy window band indices
fermi_weights = 0
band_window = 0
if mpi.is_master_node():
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)
# Set up deltaH:
deltaH = {}
for sp in spn:
deltaH[sp] = [numpy.zeros([self.n_orbitals[ik,ntoi[sp]],self.n_orbitals[ik,ntoi[sp]]], numpy.complex_) for ik in range(self.n_k)]
ikarray = numpy.array(range(self.n_k))
for ik in mpi.slice_array(ikarray):
sigma_minus_dc = [s.copy() for s in self.Sigma_imp_iw]
sigma_minus_dc = self.add_dc('iw')
beta = self.Sigma_imp_iw[0].mesh.beta # override beta if Sigma_iw is present
n_iw = len(self.Sigma_imp_iw[0].mesh)
block_structure = [ range(self.n_orbitals[ik,ntoi[sp]]) for sp in spn ]
gf_struct = [ (spn[isp], block_structure[isp]) for isp in range(self.n_spin_blocks[self.SO]) ]
block_ind_list = [block for block,inner in gf_struct]
glist = lambda : [ GfImFreq(indices=inner,beta=beta,n_points=n_iw) for block,inner in gf_struct]
G_latt = BlockGf(name_list = block_ind_list, block_list = glist(), make_copies = False)
G_latt.zero()
for icrsh in range(self.n_corr_shells):
for bname, gf in G_latt:
gf += self.upfold(ik,icrsh,bname,sigma_minus_dc[icrsh][bname],gf)
for sp in spn:
deltaH[sp][ik][:, :] = G_latt[sp](0) # Any Matsubara frequency will do
# G_latt_iw = self.lattice_gf(ik = ik, mu = self.chemical_potential, iw_or_w = "iw")
# for bname,gf in G_latt_iw:
# deltaN[bname][ik][:, :] = G_latt_iw[bname].density()
# dens[bname] += self.bz_weights[ik] * G_latt_iw[bname].total_density()
# if dm_type == 'vasp':
## In 'vasp'-mode subtract the DFT density matrix
# diag_inds = numpy.diag_indices(self.n_orbitals[ik, ntoi[bname]])
# deltaN[bname][ik][diag_inds] -= dens_mat_dft[bname][ik]
# dens[bname] -= self.bz_weights[ik] * dens_mat_dft[bname][ik].sum().real
# mpi reduce:
for bname in deltaH:
for ik in range(self.n_k):
deltaH[bname][ik] = mpi.all_reduce(mpi.world, deltaH[bname][ik], lambda x,y : x+y)
mpi.barrier()
# now save to file:
if mpi.is_master_node():
with open(filename, 'w') as f:
f.write("H %i -1 ! Number of k-points, default number of bands\n"%(self.n_k))
for sp in spn:
for ik in xrange(self.n_k):
ib1 = band_window[0][ik, 0]
ib2 = band_window[0][ik, 1]
f.write(" %i %i %i\n"%(ik + 1, ib1, ib2))
for inu in xrange(self.n_orbitals[ik, 0]):
for imu in xrange(self.n_orbitals[ik, 0]):
if self.SP == 0:
valre = (deltaH['up'][ik][inu, imu].real + deltaH['down'][ik][inu, imu].real) / 2.0
valim = (deltaH['up'][ik][inu, imu].imag + deltaH['down'][ik][inu, imu].imag) / 2.0
else:
valre = deltaH[sp][ik][inu, imu].real
valim = deltaH[sp][ik][inu, imu].imag
f.write(" %.14f %.14f"%(valre, valim))
f.write("\n")
return deltaH
def dmft_cycle():
lda_filename = 'vasp'
beta = 400
U_int = 4.00
J_hund = 0.70
Loops = 1 # Number of DMFT sc-loops
Mix = 1.0 # Mixing factor in QMC
DC_type = 0 # 0...FLL, 1...Held, 2... AMF, 3...Lichtenstein
DC_Mix = 1.0 # 1.0 ... all from imp; 0.0 ... all from Gloc
useBlocs = False # use bloc structure from LDA input
useMatrix = True # use the U matrix calculated from Slater coefficients instead of (U+2J, U, U-J)
chemical_potential_init=-0.0 # initial chemical potential
use_dudarev = True
HDFfilename = lda_filename+'.h5'
# Convert DMFT input:
# Can be commented after the first run
Converter = VaspConverter(filename=lda_filename)
Converter.convert_dft_input()
#check if there are previous runs:
previous_runs = 0
previous_present = False
if mpi.is_master_node():
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)
previous_present = mpi.bcast(previous_present)
# Init the SumK class
#SK=SumkDFT(hdf_file=lda_filename+'.h5',use_dft_blocks=False)
SK=TestSumkDFT(hdf_file=lda_filename+'.h5',use_dft_blocks=False)
Norb = SK.corr_shells[0]['dim']
l = SK.corr_shells[0]['l']
# Init the Hubbard-I solver:
S = Solver(beta = beta, l = l, dudarev=use_dudarev)
# DEBUG
#SK.put_Sigma(Sigma_imp=[S.Sigma])
#dH = SK.calc_hamiltonian_correction(filename='GAMMA')
# END DEBUG
chemical_potential=chemical_potential_init
# load previous data: old self-energy, chemical potential, DC correction
if (previous_present):
mpi.report("Using stored data for initialisation")
if (mpi.is_master_node()):
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]
SK.dc_imp=old_data[1]
S.Sigma = mpi.bcast(S.Sigma)
chemical_potential=mpi.bcast(chemical_potential)
SK.dc_imp=mpi.bcast(SK.dc_imp)
# DMFT loop:
for Iteration_Number in range(1,Loops+1):
itn = Iteration_Number + previous_runs
# put Sigma into the SumK class:
S.Sigma.zero() # !!!!
SK.put_Sigma(Sigma_imp = [ S.Sigma ])
# Compute the SumK, possibly fixing mu by dichotomy
if SK.density_required and (Iteration_Number > 1):
chemical_potential = SK.calc_mu( precision = 0.001 )
else:
mpi.report("No adjustment of chemical potential\nTotal density = %.3f"%SK.total_density(mu=chemical_potential))
# Density:
S.G <<= SK.extract_G_loc()[0]
mpi.report("Total charge of Gloc : %.6f"%S.G.total_density())
# calculated DC at the first run to have reasonable initial non-interacting atomic level positions
if ((Iteration_Number==1)and(previous_present==False)):
if use_dudarev:
dc_value_init = 0.0
else:
dc_value_init=U_int/2.0
dm=S.G.density()
SK.calc_dc( dm, U_interact = U_int, J_hund = J_hund, orb = 0, use_dc_formula = DC_type, use_dc_value=dc_value_init)
# calculate non-interacting atomic level positions:
# eal = SK.eff_atomic_levels()[0]
# S.set_atomic_levels( eal = eal )
# solve it:
corr_energy, dft_dc = S.solve(U_int = U_int, J_hund = J_hund, verbosity = 1)
# Now mix Sigma and G:
if ((itn>1)or(previous_present)):
if (mpi.is_master_node()and (Mix<1.0)):
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)
# after the Solver has finished, set new double counting:
# dm = S.G.density()
# if not use_dudarev:
# SK.calc_dc( dm, U_interact = U_int, J_hund = J_hund, orb = 0, use_dc_formula = DC_type )
# correlation energy calculations:
correnerg = 0.5 * (S.G * S.Sigma).total_density()
mpi.report("Corr. energy = %s"%correnerg)
# store the impurity self-energy, GF as well as correlation energy in h5
if (mpi.is_master_node()):
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']
SK.save(things_to_save)
# print out occupancy matrix of Ce 4f
# mpi.report("Orbital densities of impurity Green function:")
# for s in dm:
# mpi.report("Block %s: "%s)
# for ii in range(len(dm[s])):
# str = ''
# for jj in range(len(dm[s])):
# if (dm[s][ii,jj].real>0):
# str += " %.4f"%(dm[s][ii,jj].real)
# else:
# str += " %.4f"%(dm[s][ii,jj].real)
# mpi.report(str)
mpi.report("Total charge of impurity problem : %.6f"%S.G.total_density())
# find exact chemical potential
#if (SK.density_required):
# SK.chemical_potential = SK.calc_mu( precision = 0.000001 )
# SK.chemical_potential = SK.calc_mu( precision = 0.01 )
#dN, d = SK.calc_density_correction(filename='GAMMA', dm_type='vasp')
#mpi.report("Trace of Density Matrix: %s"%d)
mpi.report("Storing Hamiltonian correction GAMMA...")
SK.put_Sigma(Sigma_imp=[S.Sigma])
dH = SK.calc_hamiltonian_correction(filename='GAMMA')
# shutil.copyfile('GAMMA', 'it%i.GAMMA'%(itn))
# store correlation energy contribution to be read by Wien2ki and then included to DFT+DMFT total energy
if (mpi.is_master_node()):
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)
f.close()
return corr_energy, dft_dc
if __name__ == '__main__':
dmft_cycle()

11
python/converters/plovasp/examples/lunio3/.gitignore vendored
View File

@ -1,11 +0,0 @@
*
!.gitignore
!INCAR
!POSCAR
!POTCAR
!KPOINTS
!plo.cfg
!run_plovasp.sh
!conv_example.py
!rotations
!rot_dz2_dx2

39
python/converters/plovasp/examples/lunio3/INCAR
View File

@ -1,39 +0,0 @@
ADDGRID = TRUE
ALGO = NORMAL
EDIFF = 0.00001
EDIFFG = -0.01
EMAX = 13.0
EMIN = -3.0
ENCUT = 500
IBRION = -1
ICHARG = 11
ISIF = 2
ISMEAR = -5
ISPIN = 1
ISYM = -1
LASPH = FALSE
LCHARG = TRUE
LDAU = FALSE
LDAUJ = 0 0 0 0 0
LDAUL = 0 2 0 0 0
LDAUTYPE = 2
LDAUU = 0 6.2 0 0 0
LMAXMIX = 6
LORBIT = 0
LREAL = Auto
ROPT = 1e-3 1e-3 1e-3
LWAVE = FALSE
NEDOS = 1001
NELM = 100
NELMDL = -9
NELMIN = 5
NPAR = 4
NSIM = 1
NSW = 0
POTIM = 0.1
PREC = Accurate
SIGMA = 0.05
SMASS = 0.5
SYMPREC = 1.0e-12
LOCPROJ = 5 6 7 8 : d : Hy

4
python/converters/plovasp/examples/lunio3/KPOINTS
View File

@ -1,4 +0,0 @@
Automatic kpoint scheme
0
Gamma
5 5 3

28
python/converters/plovasp/examples/lunio3/POSCAR
View File

@ -1,28 +0,0 @@
LuNiO3, exp. low-T
1.0
5.1234998703 0.0000000000 0.0000000000
0.0000000000 5.5089001656 0.0000000000
-0.0166880521 0.0000000000 7.3551808822
Lu Ni O
4 4 12
Direct
0.977199972 0.077000000 0.252999991
0.022800028 0.922999978 0.746999979
0.522800028 0.577000022 0.247000009
0.477199972 0.423000008 0.753000021
0.500000000 0.000000000 0.000000000
0.000000000 0.500000000 0.500000000
0.500000000 0.000000000 0.500000000
0.000000000 0.500000000 0.000000000
0.110100001 0.462700009 0.244100004
0.889899969 0.537299991 0.755900025
0.389899999 0.962700009 0.255899996
0.610100031 0.037299991 0.744099975
0.693300009 0.313699991 0.053900000
0.306699991 0.686300039 0.946099997
0.806699991 0.813699961 0.446099997
0.193300009 0.186300009 0.553900003
0.185100004 0.201600000 0.943799973
0.814899981 0.798399985 0.056200027
0.314899981 0.701600015 0.556200027
0.685100019 0.298399985 0.443799973

7683
python/converters/plovasp/examples/lunio3/POTCAR
View File

File diff suppressed because it is too large Load Diff

12
python/converters/plovasp/examples/lunio3/plo.cfg
View File

@ -1,12 +0,0 @@
[General]
DOSMESH = 401
[Shell 1]
EWINDOW = -8.6 2.7
#EWINDOW = -0.6 2.7
LSHELL = 2
IONS = 5..8
NORMION = True
TRANSFILE = rotations
#TRANSFILE = rot_dz2_dx2
NORMALIZE = True

18
python/converters/plovasp/examples/lunio3/rot_dz2_dx2
View File

@ -1,18 +0,0 @@
# Rotations to the local frames
# Atom 1, Euler angles: 44.57639968 16.85783381 -11.27701880
0.02793493813 -0.09400347991 0.87385177595 -0.47142545676 0.06726141074
-0.87830274777 -0.28356545426 0.00107688253 0.06072810673 0.38011294872
# Atom 2, Euler angles: 47.67919913 16.85783423 11.27702026
-0.02793494291 0.09400349403 0.87385176987 -0.47142546468 0.06726141256
-0.84565748799 -0.27808790765 -0.00680142052 -0.08187952499 -0.44808482755
# Atom 3, Euler angles:-133.92682764 17.32502138 165.73386489
0.03668336697 0.12133769713 0.86698007031 0.47720500874 0.06747170631
-0.85519672598 0.29102955535 -0.00287628125 -0.06301365234 0.42421853379
# Atom 4, Euler angles:-135.60182669 17.32502172 -165.73386378
-0.03668337098 -0.12133770845 0.86698006522 0.47720501465 0.06747170748
-0.84809293709 0.29001625337 0.00161325481 0.06758002057 -0.43824568573

30
python/converters/plovasp/examples/lunio3/rotations
View File

@ -1,30 +0,0 @@
# Rotations to the local frames
# Atom 1, Euler angles: 44.57639968 16.85783381 -11.27701880
0.38045444863 0.05847209742 0.07282374934 0.27130935338 0.87916060116
-0.26548421431 0.78270891327 0.33738843326 0.43922773606 -0.10066245585
0.02793493813 -0.09400347991 0.87385177595 -0.47142545676 0.06726141074
0.11214484811 -0.54286439359 0.34241445738 0.71241867354 -0.26064104853
-0.87830274777 -0.28356545426 0.00107688253 0.06072810673 0.38011294872
# Atom 2, Euler angles: 47.67919913 16.85783423 11.27702026
-0.44835779602 -0.08059449489 0.07251344244 0.26569233153 0.84653954413
-0.10160916307 0.51163832768 0.35542791179 0.72915358208 -0.26440093208
-0.02793494291 0.09400349403 0.87385176987 -0.47142546468 0.06726141256
0.26969226294 -0.80346823694 0.32365046757 0.41084840677 -0.09032628101
-0.84565748799 -0.27808790765 -0.00680142052 -0.08187952499 -0.44808482755
# Atom 3, Euler angles:-133.92682764 17.32502138 165.73386489
0.42474364412 -0.05919623743 0.07674521204 -0.27806951500 0.85608186496
0.27929719729 0.78784591591 -0.35463192502 0.41101628718 0.08120158747
0.03668336697 0.12133769713 0.86698007031 0.47720500874 0.06747170631
-0.09422739956 -0.52571244938 -0.34158989281 0.72252862625 0.27571062205
-0.85519672598 0.29102955535 -0.00287628125 -0.06301365234 0.42421853379
# Atom 4, Euler angles:-135.60182669 17.32502172 -165.73386378
-0.43875616246 0.06397634469 0.07678214913 -0.27699516110 0.84898659561
0.09192658044 0.51921347838 -0.34449569877 0.72588538389 0.27636930122
-0.03668337098 -0.12133770845 0.86698006522 0.47720501465 0.06747170748
-0.28006291089 -0.79214400488 -0.35180986204 0.40505852046 0.07893071297
-0.84809293709 0.29001625337 0.00161325481 0.06758002057 -0.43824568573

2
python/converters/plovasp/examples/lunio3/run_plovasp.sh
View File

@ -1,2 +0,0 @@
PLOVASP_PATH=/path/to/triqs1.2/applications/dft_tools/src
PYTHONPATH=$PLOVASP_PATH/python:$PLOVASP_PATH/c:$PYTHONPATH pytriqs $PLOVASP_PATH/python/converters/plovasp/main.py $@

7
python/converters/plovasp/examples/simple/.gitignore vendored
View File

@ -1,7 +0,0 @@
*
!example.cfg
!INCAR
!KPOINTS
!POSCAR
!POTCAR
!conv_example.py

42
python/converters/plovasp/examples/simple/INCAR
View File

@ -1,42 +0,0 @@
ADDGRID = TRUE
ALGO = NORMAL
EDIFF = 0.00001
EDIFFG = -0.002
EMAX = 13.0
EMIN = -2.0
ENCUT = 550
ENAUG = 1100
IBRION = -1
ICHARG = 11
ISIF = 4
ISMEAR = -5
ISPIN = 1
ISYM = -1
LASPH = TRUE
LCHARG = TRUE
LDAU = False
LDAUJ = 0 1.0 0 0 0
LDAUL = -1 2 -1
LDAUTYPE = 1
LDAUU = 0 5.0 0 0 0
LDAUPRINT = 1
LMAXMIX = 6
LORBIT = 0
LREAL = False
ROPT = 1e-3 1e-3 1e-3
LWAVE = FALSE
NBANDS = 16
NEDOS = 501
NELM = 100
NELMDL = -7
NELMIN = 7
NPAR = 4
NSIM = 8
NSW = 0
POTIM = 0.4
PREC = Accurate
SIGMA = 0.1
SMASS = 0.5
SYMPREC = 1.0e-6
LOCPROJ = 1 : d : Hy

4
python/converters/plovasp/examples/simple/KPOINTS
View File

@ -1,4 +0,0 @@
Automatic kpoint scheme
0
Gamma
2 2 2

10
python/converters/plovasp/examples/simple/POSCAR
View File

@ -1,10 +0,0 @@
V SF test
2.9878
-0.5 0.5 0.5
0.5 -0.5 0.5
0.5 0.5 -0.5
V
1
Direct
0.0 0.0 0.0

2615
python/converters/plovasp/examples/simple/POTCAR
View File

File diff suppressed because it is too large Load Diff

7
python/converters/plovasp/examples/simple/conv_example.py
View File

@ -1,7 +0,0 @@
from vasp_converter import VaspConverter
if __name__ == '__main__':
conv = VaspConverter('vasp')
conv.convert_dft_input()

7
python/converters/plovasp/examples/simple/example.cfg
View File

@ -1,7 +0,0 @@
[Shell 1]
LSHELL = 2
IONS = 1
EWINDOW = -15.0 5.0

9
python/converters/plovasp/examples/srvo3/.gitignore vendored
View File

@ -1,9 +0,0 @@
*
!.gitignore
!INCAR
!POSCAR
!POTCAR
!KPOINTS
!plo.cfg
!run_plovasp.sh
!conv_example.py

43
python/converters/plovasp/examples/srvo3/INCAR
View File

@ -1,43 +0,0 @@
ADDGRID = TRUE
ALGO = NORMAL
EDIFF = 0.00001
EDIFFG = -0.002
EMAX = 13.0
EMIN = -2.0
ENCUT = 550
ENAUG = 1100
IBRION = -1
ICHARG = 11
ISIF = 4
ISMEAR = -5
ISPIN = 1
ISYM = -1
LASPH = TRUE
LCHARG = TRUE
LDAU = False
LDAUJ = 0 1.0 0 0 0
LDAUL = -1 2 -1
LDAUTYPE = 1
LDAUU = 0 5.0 0 0 0
LDAUPRINT = 1
LMAXMIX = 6
LORBIT = 0
LREAL = False
ROPT = 1e-3 1e-3 1e-3
LWAVE = FALSE
NBANDS = 32
NEDOS = 501
NELM = 100
NELMDL = -7
NELMIN = 7
NPAR = 1
NSIM = 8
NSW = 0
POTIM = 0.4
PREC = Accurate
SIGMA = 0.1
SMASS = 0.5
SYMPREC = 1.0e-6
# Projectors
LOCPROJ = 2 : d : Hy

4
python/converters/plovasp/examples/srvo3/KPOINTS
View File

@ -1,4 +0,0 @@
Automatic kpoint scheme
0
Gamma
3 3 3

13
python/converters/plovasp/examples/srvo3/POSCAR
View File

@ -1,13 +0,0 @@
SrVO3
3.841
1.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 1.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 1.0000000000000000
Sr V O
1 1 3
Direct
0.5 0.5 0.5
0.0 0.0 0.0
0.5 0.0 0.0
0.0 0.5 0.0
0.0 0.0 0.5

7739
python/converters/plovasp/examples/srvo3/POTCAR
View File

File diff suppressed because it is too large Load Diff

51
python/converters/plovasp/examples/srvo3/conv_example.py
View File

@ -1,51 +0,0 @@
import numpy as np
from pytriqs.gf import *
#from sumk_dft import SumkDFT
from sumk_dft_tools import SumkDFTTools
from converters.vasp_converter import VaspConverter
np.set_printoptions(suppress=True)
def density_matrix_and_overlap(sk):
glist = [GfImFreq(indices=inds,beta=1000.0) for bl, inds in sk.gf_struct_sumk[0]]
sigma = BlockGf(name_list=[bl for bl, inds in sk.gf_struct_sumk[0]], block_list=glist)
sigma.zero()
sk.put_Sigma(Sigma_imp=[sigma])
print "Overlap matrix:"
dm_blocks = sk.check_projectors()
print dm_blocks[0].real
sk.calc_mu(precision=0.001)
print
print "Denisty matrix:"
dm_blocks = sk.density_matrix(method='using_gf', beta=1000.0)
ntot = 0.0
for bl in dm_blocks[0]:
print
print bl
print dm_blocks[0][bl].real
ntot += dm_blocks[0][bl].real.trace()
print
print " Impurity density:", ntot
def density_of_states(sk):
glist = [GfReFreq(indices=inds, window=(-10.0, 10.0), n_points=2000) for bl, inds in sk.gf_struct_sumk[0]]
sigma = BlockGf(name_list=[bl for bl, inds in sk.gf_struct_sumk[0]], block_list=glist)
sigma.zero()
sk.put_Sigma(Sigma_imp=[sigma])
print " Evaluating DOS..."
sk.dos_wannier_basis(broadening=0.03, with_dc=False)
if __name__ == '__main__':
conv = VaspConverter('vasp')
conv.convert_dft_input()
sk = SumkDFTTools(hdf_file='vasp.h5')
density_matrix_and_overlap(sk)
# density_of_states(sk)

16
python/converters/plovasp/examples/srvo3/plo.cfg
View File

@ -1,16 +0,0 @@
[General]
DOSMESH = -3.0 5.0 4001
[Group 1]
SHELLS = 1
NORMALIZE = True
EWINDOW = -1.45 1.8
[Shell 1]
LSHELL = 2
IONS = 2
TRANSFORM = 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

2
python/converters/plovasp/examples/srvo3/run_plovasp.sh
View File

@ -1,2 +0,0 @@
PLOVASP_PATH=/path/to/triqs1.2/applications/dft_tools/src
PYTHONPATH=$PLOVASP_PATH/python:$PLOVASP_PATH/c:$PYTHONPATH pytriqs $PLOVASP_PATH/python/converters/plovasp/main.py $@

42
python/converters/plovasp/examples/two_site/INCAR
View File

@ -1,42 +0,0 @@
ADDGRID = TRUE
ALGO = NORMAL
EDIFF = 0.00001
EDIFFG = -0.002
EMAX = 13.0
EMIN = -2.0
ENCUT = 550
ENAUG = 1100
IBRION = -1
ICHARG = 11
ISIF = 4
ISMEAR = -5
ISPIN = 1
ISYM = -1
LASPH = TRUE
LCHARG = TRUE
LDAU = False
LDAUJ = 0 1.0 0 0 0
LDAUL = -1 2 -1
LDAUTYPE = 1
LDAUU = 0 5.0 0 0 0
LDAUPRINT = 1
LMAXMIX = 6
LORBIT = 0
LREAL = False
ROPT = 1e-3 1e-3 1e-3
LWAVE = FALSE
NBANDS = 32
NEDOS = 501
NELM = 100
NELMDL = -7
NELMIN = 7
NPAR = 4
NSIM = 8
NSW = 0
POTIM = 0.4
PREC = Accurate
SIGMA = 0.1
SMASS = 0.5
SYMPREC = 1.0e-6
LOCPROJ = 1 2 : d : Hy

4
python/converters/plovasp/examples/two_site/KPOINTS
View File

@ -1,4 +0,0 @@
Automatic kpoint scheme
0
Gamma
2 2 2

11
python/converters/plovasp/examples/two_site/POSCAR
View File

@ -1,11 +0,0 @@
V SF test
2.9878
1.0 0.0 0.0
0.0 1.0 0.0
0.0 0.0 1.0
V
2
Direct
0.0 0.0 0.0
0.5 0.5 0.5

2615
python/converters/plovasp/examples/two_site/POTCAR
View File

File diff suppressed because it is too large Load Diff

7
python/converters/plovasp/examples/two_site/example.cfg
View File

@ -1,7 +0,0 @@
[Shell 1]
LSHELL = 2
IONS = 1
EWINDOW = -15.0 5.0

55
python/converters/plovasp/inpconf.py
View File

@ -24,8 +24,8 @@
#
################################################################################
r"""
vasp.inpconfig
==============
plovasp.inpconfig
=================
Module for parsing and checking an input config-file.
"""
@ -84,7 +84,9 @@ class ConfigParameters:
self.sh_optional = {
'transform': ('tmatrix', lambda s: self.parse_string_tmatrix(s, real=True)),
'transfile': ('tmatrices', self.parse_file_tmatrix)}
'transfile': ('tmatrices', self.parse_file_tmatrix),
'sort': ('ion_sort', self.parse_string_int,None),
'corr': ('corr', self.parse_string_logical, True)}
self.gr_required = {
'shells': ('shells', lambda s: map(int, s.split())),
@ -92,12 +94,16 @@ class ConfigParameters:
self.gr_optional = {
'normalize' : ('normalize', self.parse_string_logical, True),
'normion' : ('normion', self.parse_string_logical, True)}
'normion' : ('normion', self.parse_string_logical, True),
'complement' : ('complement', self.parse_string_logical, False),
'bands': ('bands', self.parse_band_window)}
self.gen_optional = {
'basename' : ('basename', str, 'vasp'),
'efermi' : ('efermi', float),
'dosmesh': ('dosmesh', self.parse_string_dosmesh)}
'dosmesh': ('dosmesh', self.parse_string_dosmesh),
'hk': ('hk', self.parse_string_logical, False)}
#
# Special parsers
@ -190,6 +196,18 @@ class ConfigParameters:
assert first_char in 'tf', "Logical parameters should be given by either 'True' or 'False'"
return first_char == 't'
################################################################################
#
# parse_string_int()
#
################################################################################
def parse_string_int(self, par_str):
"""
int parameters
"""
return int(par_str)
################################################################################
#
# parse_energy_window()
@ -205,6 +223,21 @@ class ConfigParameters:
assert ftmp[0] < ftmp[1], "The first float in EWINDOW must be smaller than the second one"
return tuple(ftmp)
################################################################################
#
# parse_band_window()
#
################################################################################
def parse_band_window(self, par_str):
"""
Band window is given by two ints, with the first one being smaller
than the second one.
"""
ftmp = map(int, par_str.split())
assert len(ftmp) == 2, "BANDS must be specified by exactly two ints"
assert ftmp[0] < ftmp[1], "The first int in BANDS must be smaller than the second one"
return tuple(ftmp)
################################################################################
#
# parse_string_tmatrix()
@ -461,12 +494,11 @@ class ConfigParameters:
################################################################################
def groups_shells_consistency(self):
"""
Ensures consistency between groups and shells.
In particular:
- if no groups are explicitly defined and only shell is defined create
a group automatically
- check the existance of all shells referenced in the groups
- check that all shells are referenced in the groups
Ensures consistency between groups and shells. In particular:
- if no groups are explicitly defined and only shell is defined create a group automatically
- check the existance of all shells referenced in the groups
- check that all shells are referenced in the groups
"""
# Special case: no groups is defined
if self.ngroups == 0:
@ -627,4 +659,3 @@ if __name__ == '__main__':
doscar = vaspio.Doscar()
doscar.from_file(vasp_dir)
# pars = parse_input(filename)

265
python/converters/plovasp/plotools.py
View File

@ -24,8 +24,8 @@
#
################################################################################
r"""
vasp.plotools
=============
plovasp.plotools
================
Set of routines for processing and outputting PLOs.
@ -33,11 +33,31 @@ r"""
the consistency of the input data, generation of projected localized
orbitals (PLOs) out of raw VASP projectors, and outputting data
required by DFTTools.
The first step of PLO processing is to select subsets of projectors
corresponding to PLO groups. Each group contains a set of shells. Each
projector shell is represented by an object 'ProjectorShell' that contains
an array of projectors and information on the shell itself (orbital number,
ions, etc.). 'ProjectorShell's are contained in both a list of shells
(according to the original list as read from config-file) and in a
'ProjectorGroup' object, the latter also providing information about the
energy window.
Order of operations:
- transform projectors (all bands) in each shell
- select transformed shell projectors for a given group within the window
- orthogonalize if necessary projectors within a group by performing
the following operations for each k-point:
* combine all projector shells into a single array
* orthogonalize the array
* distribute back the arrays assuming that the order is preserved
"""
import itertools as it
import numpy as np
from proj_group import ProjectorGroup
from proj_shell import ProjectorShell
from proj_shell import ComplementShell
np.set_printoptions(suppress=True)
@ -81,6 +101,7 @@ def check_data_consistency(pars, el_struct):
errmsg = "Projector for isite = %s, l = %s does not match PROJCAR"%(ion + 1, lshell)
raise Exception(errmsg)
################################################################################
#
# generate_plo()
@ -105,7 +126,8 @@ def generate_plo(conf_pars, el_struct):
# eigvals(nktot, nband, ispin) are defined with respect to the Fermi level
eigvals = el_struct.eigvals - efermi
# check if at least one shell is correlated
assert np.any([shell['corr'] for shell in conf_pars.shells]), 'at least one shell has be CORR = True'
nshell = len(conf_pars.shells)
print
print " Generating %i shell%s..."%(nshell, '' if nshell == 1 else 's')
@ -117,39 +139,59 @@ def generate_plo(conf_pars, el_struct):
print " Orbital l : %i"%(pshell.lorb)
print " Number of ions: %i"%(pshell.nion)
print " Dimension : %i"%(pshell.ndim)
print " Correlated : %r"%(pshell.corr)
print " Ion sort : %r"%(pshell.ion_sort)
pshells.append(pshell)
pgroups = []
for gr_par in conf_pars.groups:
pgroup = ProjectorGroup(gr_par, pshells, eigvals)
pgroup.orthogonalize()
if pgroup.complement:
pgroup.calc_complement(eigvals)
if conf_pars.general['hk']:
pgroup.calc_hk(eigvals)
#testout = 'hk.out.h5'
#from pytriqs.archive import HDFArchive
#with HDFArchive(testout, 'w') as h5test:
# h5test['hk'] = pgroup.hk
# DEBUG output
print "Density matrix:"
dm_all, ov_all = pshells[pgroup.ishells[0]].density_matrix(el_struct)
nimp = 0.0
spin_fac = 2 if dm_all.shape[0] == 1 else 1
for io in xrange(dm_all.shape[1]):
print " Site %i"%(io + 1)
dm = spin_fac * dm_all[:, io, : ,:].sum(0)
for row in dm:
print ''.join(map("{0:12.7f}".format, row))
ndm = dm.trace()
nimp += ndm
print " trace: ", ndm
ov_all = []
for ish in pgroup.ishells:
if not isinstance(pshells[pgroup.ishells[ish]],ComplementShell):
print " Shell %i"%(ish + 1)
dm_all, ov_all_ = pshells[ish].density_matrix(el_struct)
ov_all.append(ov_all_[0])
spin_fac = 2 if dm_all.shape[0] == 1 else 1
for io in xrange(dm_all.shape[1]):
print " Site %i"%(io + 1)
dm = spin_fac * dm_all[:, io, : ,:].sum(0)
for row in dm:
print ''.join(map("{0:14.7f}".format, row))
ndm = dm.trace()
if pshells[ish].corr:
nimp += ndm
print " trace: ", ndm
print
print " Impurity density:", nimp
print
print "Overlap:"
for io, ov in enumerate(ov_all[0]):
for io, ov in enumerate(ov_all):
print " Site %i"%(io + 1)
print ov
print ov[0,...]
print
print "Local Hamiltonian:"
loc_ham = pshells[pgroup.ishells[0]].local_hamiltonian(el_struct)
for io in xrange(loc_ham.shape[1]):
print " Site %i"%(io + 1)
for row in loc_ham[:, io, :, :].sum(0):
print ''.join(map("{0:12.7f}".format, row))
for ish in pgroup.ishells:
if not isinstance(pshells[pgroup.ishells[ish]],ComplementShell):
print " Shell %i"%(ish + 1)
loc_ham = pshells[pgroup.ishells[ish]].local_hamiltonian(el_struct)
for io in xrange(loc_ham.shape[1]):
print " Site %i"%(io + 1)
for row in loc_ham[:, io, :, :].sum(0):
print ''.join(map("{0:14.7f}".format, row))
# END DEBUG output
if 'dosmesh' in conf_pars.general:
print
@ -164,12 +206,15 @@ def generate_plo(conf_pars, el_struct):
n_points = mesh_pars['n_points']
emesh = np.linspace(dos_emin, dos_emax, n_points)
dos = pshells[pgroup.ishells[0]].density_of_states(el_struct, emesh)
de = emesh[1] - emesh[0]
ntot = (dos[1:,...] + dos[:-1,...]).sum(0) / 2 * de
print " Total number of states:", ntot
for io in xrange(dos.shape[2]):
np.savetxt('pdos_%i.dat'%(io), np.vstack((emesh.T, dos[:, 0, io, :].T)).T)
for ish in pgroup.ishells:
if not isinstance(pshells[pgroup.ishells[ish]],ComplementShell) or True:
print " Shell %i"%(ish + 1)
dos = pshells[pgroup.ishells[ish]].density_of_states(el_struct, emesh)
de = emesh[1] - emesh[0]
ntot = (dos[1:,...] + dos[:-1,...]).sum(0) / 2 * de
print " Total number of states:", ntot
for io in xrange(dos.shape[2]):
np.savetxt('pdos_%i_%i.dat'%(ish,io), np.vstack((emesh.T, dos[:, 0, io, :].T)).T)
pgroups.append(pgroup)
@ -183,9 +228,12 @@ def generate_plo(conf_pars, el_struct):
def output_as_text(pars, el_struct, pshells, pgroups):
"""
Output all information necessary for the converter as text files.
"""
ctrl_output(pars, el_struct, len(pgroups))
plo_output(pars, el_struct, pshells, pgroups)
if pars.general['hk']:
hk_output(pars, el_struct, pgroups)
# TODO: k-points with weights should be stored once and for all
@ -232,6 +280,27 @@ def kpoints_output(basename, el_struct):
def ctrl_output(conf_pars, el_struct, ng):
"""
Outputs a ctrl-file.
Control file format
""""""""""""""""""""""""""""""
Filename '<namebase>.ctrl'. Contains the data shared between all shells.
The JSON-header consists of the following elements:
* *nk*: number of `k`-points
* *ns*: number of spin channels
* *nc_flag*: collinear/noncollinear case (False/True)
* *ng*: number of projector groups
* Symmetry information (list of symmetry operations)
* *efermi*: Fermi level (optional)
* Lattice information
"""
ctrl_fname = conf_pars.general['basename'] + '.ctrl'
head_dict = {}
@ -259,6 +328,13 @@ def ctrl_output(conf_pars, el_struct, ng):
tmp1 = "".join(map("{0:15.10f}".format, kp))
out = tmp1 + "{0:16.10f}".format(el_struct.kmesh['kweights'][ik])
f.write(out + "\n")
f.write("# k-points and weights cartesian\n")
labels = ['kx', 'ky', 'kz']
out = "".join(map(lambda s: s.center(15), labels))
f.write("#" + out + "\n")
for ik, kp in enumerate(el_struct.kmesh['kpoints_cart']):
out = "".join(map("{0:15.10f}".format, kp))
f.write(out + "\n")
################################################################################
@ -276,31 +352,31 @@ def plo_output(conf_pars, el_struct, pshells, pgroups):
Each group is stored in a '<basename>.plog<Ng>' file. The format is the
following:
# Energy window: emin, emax
ib_min, ib_max
nelect
# Eigenvalues
isp, ik1, kx, ky, kz, kweight
ib1, ib2
eig1
eig2
...
eigN
ik2, kx, ky, kz, kweight
...
| # Energy window: emin, emax
| ib_min, ib_max
| nelect
| # Eigenvalues
| isp, ik1, kx, ky, kz, kweight
| ib1, ib2
| eig1
| eig2
| ...
| eigN
| ik2, kx, ky, kz, kweight
| ...
# Projected shells
Nshells
# Shells: <shell indices>
# Shell <1>
Shell 1
ndim
# complex arrays: plo(ns, nion, ndim, nb)
...
# Shells: <shell indices>
# Shell <2>
Shell 2
...
| # Projected shells
| Nshells
| # Shells: <shell indices>
| # Shell <1>
| Shell 1
| ndim
| # complex arrays: plo(ns, nion, ndim, nb)
| ...
| # Shells: <shell indices>
| # Shell <2>
| Shell 2
| ...
"""
for ig, pgroup in enumerate(pgroups):
@ -308,9 +384,14 @@ def plo_output(conf_pars, el_struct, pshells, pgroups):
print " Storing PLO-group file '%s'..."%(plo_fname)
head_dict = {}
head_dict['ewindow'] = (pgroup.emin, pgroup.emax)
head_dict['nb_max'] = pgroup.nb_max
if 'bands' in conf_pars.groups[ig]:
head_dict['bandwindow'] = (pgroup.ib_min, pgroup.ib_max)
else:
head_dict['ewindow'] = (pgroup.emin, pgroup.emax)
# Number of electrons within the window
head_dict['nelect'] = pgroup.nelect_window(el_struct)
print " Density within window:", head_dict['nelect']
@ -322,6 +403,7 @@ def plo_output(conf_pars, el_struct, pshells, pgroups):
sh_dict['shell_index'] = ish
sh_dict['lorb'] = shell.lorb
sh_dict['ndim'] = shell.ndim
sh_dict['corr'] = shell.corr
# Convert ion indices from the internal representation (starting from 0)
# to conventional VASP representation (starting from 1)
ion_output = [io + 1 for io in shell.ion_list]
@ -342,7 +424,11 @@ def plo_output(conf_pars, el_struct, pshells, pgroups):
f.write("#END OF HEADER\n")
# Eigenvalues within the window
f.write("# Eigenvalues within the energy window: %s, %s\n"%(pgroup.emin, pgroup.emax))
if 'bands' in conf_pars.groups[ig]:
f.write("# Eigenvalues within the band window: %s, %s\n"%(pgroup.ib_min+1, pgroup.ib_max+1))
else:
f.write("# Eigenvalues within the energy window: %s, %s\n"%(pgroup.emin, pgroup.emax))
nk, nband, ns_band = el_struct.eigvals.shape
for isp in xrange(ns_band):
f.write("# is = %i\n"%(isp + 1))
@ -376,4 +462,79 @@ def plo_output(conf_pars, el_struct, pshells, pgroups):
f.write("{0:16.10f}{1:16.10f}\n".format(p.real, p.imag))
f.write("\n")
################################################################################
#
# plo_output
#
################################################################################
def hk_output(conf_pars, el_struct, pgroups):
"""
Outputs HK into text file.
Filename is defined by <basename> that is passed from config-file.
The Hk for each groups is stored in a '<basename>.hk<Ng>' file. The format is
similar as defined in the Hk dft_tools format, but does not store info
about correlated shells and irreps
nk # number of k-points
n_el # electron density
n_sh # number of total atomic shells
at sort l dim # atom, sort, l, dim
at sort l dim # atom, sort, l, dim
After these header lines, the file has to contain the Hamiltonian matrix
in orbital space. The standard convention is that you give for each k-point
first the matrix of the real part, then the matrix of the imaginary part,
and then move on to the next k-point.
"""
for ig, pgroup in enumerate(pgroups):
hk_fname = conf_pars.general['basename'] + '.hk%i'%(ig + 1)
print " Storing HK-group file '%s'..."%(hk_fname)
head_shells = []
for ish in pgroup.ishells:
shell = pgroup.shells[ish]
ion_output = [io + 1 for io in shell.ion_list]
for iion in ion_output:
sh_dict = {}
sh_dict['shell_index'] = ish
sh_dict['lorb'] = shell.lorb
sh_dict['ndim'] = shell.ndim
# Convert ion indices from the internal representation (starting from 0)
# to conventional VASP representation (starting from 1)
# Derive sorts from equivalence classes
sh_dict['ion_list'] = ion_output
sh_dict['ion_sort'] = shell.ion_sort
head_shells.append(sh_dict)
with open(hk_fname, 'wt') as f:
# Eigenvalues within the window
nk, nband, ns_band = el_struct.eigvals.shape
f.write('%i # number of kpoints\n'%nk)
f.write('{0:0.4f} # electron density\n'.format(pgroup.nelect_window(el_struct)))
f.write('%i # number of shells\n'%len(head_shells))
for head in head_shells:
f.write('%i %i %i %i # atom sort l dim\n'%(head['ion_list'][0],head['ion_sort'][0],head['lorb'],head['ndim']))
norbs = pgroup.hk.shape[2]
for isp in xrange(ns_band):
for ik in xrange(nk):
for io in xrange(norbs):
for iop in xrange(norbs):
f.write(" {0:14.10f}".format(pgroup.hk[isp,ik,io,iop].real))
f.write("\n")
for io in xrange(norbs):
for iop in xrange(norbs):
f.write(" {0:14.10f}".format(pgroup.hk[isp,ik,io,iop].imag))
f.write("\n")

166
python/converters/plovasp/proj_group.py
View File

@ -24,13 +24,13 @@
#
################################################################################
r"""
vasp.proj_group
===============
plovasp.proj_group
==================
Storage and manipulation of projector groups.
"""
import numpy as np
from proj_shell import ComplementShell
np.set_printoptions(suppress=True)
################################################################################
@ -48,10 +48,9 @@ class ProjectorGroup:
the parameters from the config-file (passed in `pars`).
Parameters:
- gr_pars (dict) : group parameters from the config-file
- shells ([ProjectorShell]) : array of ProjectorShell objects
- eigvals (numpy.array) : array of KS eigenvalues
- gr_pars (dict) : group parameters from the config-file
- shells ([ProjectorShell]) : array of ProjectorShell objects
- eigvals (numpy.array) : array of KS eigenvalues
"""
def __init__(self, gr_pars, shells, eigvals):
@ -62,16 +61,38 @@ class ProjectorGroup:
self.ishells = gr_pars['shells']
self.ortho = gr_pars['normalize']
self.normion = gr_pars['normion']
self.complement = gr_pars['complement']
self.shells = shells
# Determine the minimum and maximum band numbers
ib_win, ib_min, ib_max = self.select_bands(eigvals)
if 'bands' in gr_pars:
nk, nband, ns_band = eigvals.shape
ib_win = np.zeros((nk, ns_band, 2), dtype=np.int32)
ib_win[:,:,0] = gr_pars['bands'][0]-1
ib_win[:,:,1] = gr_pars['bands'][1]-1
ib_min = gr_pars['bands'][0] - 1
ib_max = gr_pars['bands'][1] - 1
else:
ib_win, ib_min, ib_max = self.select_bands(eigvals)
self.ib_win = ib_win
self.ib_min = ib_min
self.ib_max = ib_max
self.nb_max = ib_max - ib_min + 1
if self.complement:
n_bands = self.ib_win[:,:,1] - self.ib_win[:,:,0]+1
n_orbs = sum([x.ndim for x in self.shells])
assert np.all( n_bands == n_bands[0,0] ), "At each band the same number of bands has to be selected for calculating the complement (to end up with an equal number of orbitals at each k-point)."
if n_orbs == n_bands[0,0]:
self.complement = False
print "\nWARNING: The total number of orbitals in this group is "
print "equal to the number of bands. Setting COMPLEMENT to FALSE!\n"
# Select projectors within the energy window
for ish in self.ishells:
shell = self.shells[ish]
@ -156,6 +177,130 @@ class ProjectorGroup:
shell = self.shells[ish]
shell.proj_win[ion, isp, ik, :nlm, :nb] = p_orth[i1:i2, :nb]
################################################################################
#
# calc_hk
#
################################################################################
def calc_hk(self, eigvals):
"""
Calculate H(k) for a group by applying the projectors P
to the eigenvalues eps.
H_ij(k) = sum_l P*_il eps_l P_lj
"""
# here we abuse the get_block_matrix_map(), however, it only works
# if self.normion is false
temp = self.normion
self.normion = False
block_maps, ndim = self.get_block_matrix_map()
self.normion = temp
_, ns, nk, _, _ = self.shells[0].proj_win.shape
self.hk = np.zeros((ns,nk,ndim,ndim), dtype=np.complex128)
# Note that 'ns' and 'nk' are the same for all shells
for isp in xrange(ns):
for ik in xrange(nk):
bmin = self.ib_win[ik, isp, 0]
bmax = self.ib_win[ik, isp, 1]+1
nb = bmax - bmin
p_mat = np.zeros((ndim, nb), dtype=np.complex128)
#print(bmin,bmax,nb)
# Combine all projectors of the group to one block projector
for bl_map in block_maps:
p_mat[:, :] = 0.0j # !!! Clean-up from the last k-point and block!
for ibl, block in enumerate(bl_map):
i1, i2 = block['bmat_range']
ish, ion = block['shell_ion']
nlm = i2 - i1 + 1
shell = self.shells[ish]
p_mat[i1:i2, :nb] = shell.proj_win[ion, isp, ik, :nlm, :nb]
self.hk[isp,ik,:,:] = np.dot(p_mat*eigvals[ik,bmin:bmax,isp],
p_mat.transpose().conjugate())
################################################################################
#
# complement
#
################################################################################
def calc_complement(self,eigvals):
"""
Calculate the complement for a group of projectors.
This leads to quadtratic projectors P = <l|n> by using a Gram-Schmidt.
The projector on the orthogonal complement of the existing projectors
|l> is P^u = 1 - sum_l |l><l|
We get candidates for complement projectors by applying P^u to a Bloch
state |n>: |l*> = P^u |n>. For numerical stability we select that Bloch
state which leads to the |l*> with the largest norm (that corresponds to
that Bloch state with the smallest overlap with the space spanned by |l>)
We normalize |l*> and add it to |l>. We do so untill we have as many
|l> states as we have |n> states.
"""
print '\nCalculating complement\n'
block_maps, ndim = self.get_block_matrix_map()
_, ns, nk, _, _ = self.shells[0].proj_win.shape
p_mat = np.zeros((ndim, self.nb_max), dtype=np.complex128)
p_full = np.zeros((1,ns,nk,self.nb_max, self.nb_max), dtype=np.complex128)
# Note that 'ns' and 'nk' are the same for all shells
for isp in xrange(ns):
for ik in xrange(nk):
bmin = self.ib_win[ik, isp, 0]
bmax = self.ib_win[ik, isp, 1]+1
nb = bmax - bmin
# Combine all projectors of the group to one block projector
for bl_map in block_maps:
p_mat[:, :] = 0.0j # !!! Clean-up from the last k-point and block!
for ibl, block in enumerate(bl_map):
i1, i2 = block['bmat_range']
ish, ion = block['shell_ion']
nlm = i2 - i1 + 1
shell = self.shells[ish]
p_mat[i1:i2, :nb] = shell.proj_win[ion, isp, ik, :nlm, :nb]
orbs_done = 1*ndim
p_full[0,isp,ik,:ndim,:] = p_mat
while orbs_done < self.nb_max:
#We calculate the overlap of all bloch states: sum_l <n|l><l|m>
overlap = np.dot(p_full[0,isp,ik,:orbs_done,:].transpose().conjugate(),p_full[0,isp,ik,:orbs_done,:])
# work is the projector onto the orthogonal complment <n| ( 1 - sum_l |l><l| ) |m>
work = np.eye(self.nb_max) - overlap
# calculate the norm of the projected bloch function
norm = np.sqrt(np.sum(work*work.transpose(),axis=1))
# select the bloch function leading to the largest norm
max_ind = np.argmax(norm)
# normalize and put it to the projectors
p_full[0,isp,ik,orbs_done,:] = work[:,max_ind].conjugate()/norm[max_ind]
orbs_done += 1
sh_pars = {}
sh_pars['lshell'] = -1
sh_pars['ions'] = {'nion':1,'ion_list':[[1]]}
sh_pars['user_index'] = 'complement'
sh_pars['corr'] = False
sh_pars['ib_min'] = bmin
sh_pars['ib_max'] = bmax
sh_pars['ib_win'] = self.ib_win
self.shells.append(ComplementShell(sh_pars,p_full[:,:,:,ndim:,:],False))
self.ishells.append(self.ishells[-1]+1)
################################################################################
#
# gen_block_matrix_map
@ -299,7 +444,8 @@ class ProjectorGroup:
Returns
-------
ib_win, nb_min, nb_max :
ib_win, nb_min, nb_max : lowest and highest indices of the selected bands
"""
# Sanity check
if self.emin > eigvals.max() or self.emax < eigvals.min():
@ -335,5 +481,3 @@ class ProjectorGroup:
ib_max = max(ib_max, ib2)
return ib_win, ib_min, ib_max

75
python/converters/plovasp/proj_shell.py
View File

@ -24,8 +24,8 @@
#
################################################################################
r"""
vasp.proj_shell
===============
plovasp.proj_shell
==================
Storage and manipulation on projector shells.
"""
@ -72,6 +72,8 @@ class ProjectorShell:
self.lorb = sh_pars['lshell']
self.ions = sh_pars['ions']
self.user_index = sh_pars['user_index']
self.corr = sh_pars['corr']
self.ion_sort = [sh_pars['ion_sort']]
self.nc_flag = nc_flag
# try:
# self.tmatrix = sh_pars['tmatrix']
@ -84,12 +86,14 @@ class ProjectorShell:
self.nion = self.ions['nion']
# Extract ion list and equivalence classes (ion sorts)
self.ion_list = sorted(it.chain(*self.ions['ion_list']))
self.ion_sort = []
for ion in self.ion_list:
for icl, eq_cl in enumerate(self.ions['ion_list']):
if ion in eq_cl:
self.ion_sort.append(icl + 1) # Enumerate classes starting from 1
break
if self.ion_sort[0] is None:
self.ion_sort = []
for ion in self.ion_list:
for icl, eq_cl in enumerate(self.ions['ion_list']):
if ion in eq_cl:
self.ion_sort.append(icl + 1) # Enumerate classes starting from 1
break
self.ndim = self.extract_tmatrices(sh_pars)
@ -449,5 +453,60 @@ class ProjectorShell:
return dos
################################################################################
################################################################################
#
# class ProjectorShell
#
################################################################################
################################################################################
class ComplementShell(ProjectorShell):
"""
Container of projectors related to a complement shell.
Parameters:
- sh_pars (dict) : shell parameters from the config-file
- proj_compl (numpy.array) : array of complement projectors
"""
def __init__(self, sh_pars, proj_compl, nc_flag):
self.lorb = sh_pars['lshell']
self.ions = sh_pars['ions']
self.user_index = sh_pars['user_index']
self.corr = sh_pars['corr']
self.nc_flag = nc_flag
self.ib_min = sh_pars['ib_min']
self.ib_max = sh_pars['ib_max']
self.ib_win = sh_pars['ib_win']
#self.lm1 = self.lorb**2
#self.lm2 = (self.lorb+1)**2
self.nion = self.ions['nion']
# Extract ion list and equivalence classes (ion sorts)
self.ion_list = sorted(it.chain(*self.ions['ion_list']))
self.ion_sort = []
for ion in self.ion_list:
for icl, eq_cl in enumerate(self.ions['ion_list']):
if ion in eq_cl:
self.ion_sort.append(icl + 1) # Enumerate classes starting from 1
break
self.ndim = proj_compl.shape[3]
self.proj_win = proj_compl
def extract_tmatrices(self, sh_pars):
raise Exception('not implemented')
def local_hamiltonian(self, el_struct, site_diag=True, spin_diag=True):
raise Exception('not implemented')
def density_matrix(self, el_struct, site_diag=True, spin_diag=True):
raise Exception('not implemented')
#def density_of_states(self, el_struct, emesh):
# raise Exception('not implemented')

72
python/converters/plovasp/sc_dmft.py
View File

@ -32,6 +32,7 @@ import signal
import sys
import pytriqs.utility.mpi as mpi
import converter
from shutil import copyfile
xch = sys.excepthook
def excepthook(typ, value, traceback):
@ -40,7 +41,7 @@ def excepthook(typ, value, traceback):
mpi.MPI.COMM_WORLD.Abort(1)
sys.excepthook = excepthook
debug = True
debug = False
#
# Helper functions
#
@ -99,7 +100,7 @@ class bcolors:
#
# Main self-consistent cycle
#
def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter):
def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter, n_iter_dft, vasp_version):
"""
"""
mpi.report(" Waiting for VASP lock to appear...")
@ -120,15 +121,9 @@ def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter):
mpi.report(" VASP stopped")
vasp_running = False
break
# Tell VASP to stop if the maximum number of iterations is reached
iter += 1
if iter == n_iter:
if mpi.is_master_node():
print "\n Maximum number of iterations reached."
print " Aborting VASP iterations...\n"
f_stop = open('STOPCAR', 'wt')
f_stop.write("LABORT = .TRUE.\n")
f_stop.close()
if debug: print bcolors.MAGENTA + "rank %s"%(mpi.rank) + bcolors.ENDC
err = 0
@ -155,9 +150,37 @@ def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter):
print "="*80
print
if mpi.is_master_node() and vasp_running:
open('./vasp.lock', 'a').close()
# check if we should do additional VASP calculations
# in the standard VASP version, VASP writes out GAMMA itself
# so that if we want to keep GAMMA fixed we have to copy it to
# GAMMA_recent and copy it back after VASP has completed an iteration
# if we are using a hacked Version of VASP where the write out is
# disabled we can skip this step.
# the hack consists of removing the call of LPRJ_LDApU in VASP src file
# electron.F around line 644
iter_dft = 0
if vasp_version == 'standard':
copyfile(src='GAMMA',dst='GAMMA_recent')
while iter_dft < n_iter_dft:
if mpi.is_master_node():
open('./vasp.lock', 'a').close()
while is_vasp_lock_present():
time.sleep(1)
if not is_vasp_running(vasp_pid):
mpi.report(" VASP stopped")
vasp_running = False
break
iter_dft += 1
if vasp_version == 'standard':
copyfile(src='GAMMA_recent',dst='GAMMA')
iter += 1
if iter == n_iter:
print "\n Maximum number of iterations reached."
print " Aborting VASP iterations...\n"
f_stop = open('STOPCAR', 'wt')
f_stop.write("LABORT = .TRUE.\n")
f_stop.close()
if mpi.is_master_node():
total_energy = dft_energy + corr_energy - dft_dc
with open('TOTENERGY', 'w') as f:
@ -170,7 +193,9 @@ def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter):
mpi.report("***Done")
def main():
import importlib
try:
vasp_pid = int(sys.argv[1])
except (ValueError, KeyError):
@ -186,18 +211,33 @@ def main():
raise
try:
dmft_script = re.sub("\.py$", "", sys.argv[3])
n_iter_dft = int(sys.argv[3])
except (ValueError, KeyError):
if mpi.is_master_node():
print "ERROR: Number of VASP iterations with fixed charge density must be provided as the third argument"
raise
try:
dmft_script = re.sub("\.py$", "", sys.argv[4])
except:
if mpi.is_master_node():
print "ERROR: User-defined DMFT script must be provided as the third argument"
print "ERROR: User-defined DMFT script must be provided as the fourth argument"
raise
# Optional parameter: config-file name
try:
cfg_file = sys.argv[4]
cfg_file = sys.argv[5]
except KeyError:
cfg_file = 'plo.cfg'
try:
vasp_version = sys.argv[6]
except KeyError:
vasp_version = 'standard'
if vasp_version != 'standard' and vasp_version != 'no_gamma_write':
raise Exception('vasp_version has to be standard or no_gamma_write')
# if len(sys.argv) > 1:
# vasp_path = sys.argv[1]
# else:
@ -211,7 +251,7 @@ def main():
dmft_mod = importlib.import_module(dmft_script)
run_all(vasp_pid, dmft_mod.dmft_cycle, cfg_file, n_iter)
run_all(vasp_pid, dmft_mod.dmft_cycle, cfg_file, n_iter, n_iter_dft, vasp_version)
if __name__ == '__main__':
main()

73
python/converters/plovasp/vaspio.py
View File

@ -24,8 +24,8 @@
#
################################################################################
r"""
vasp.vaspio
===========
plovasp.vaspio
==============
Input of required VASP data.
@ -88,8 +88,12 @@ class VaspData:
self.doscar.from_file(vasp_dir)
except (IOError, StopIteration):
if efermi_required:
# raise Exception("Efermi cannot be read from DOSCAR")
pass
print "!!! WARNING !!!: Error reading from Efermi from DOSCAR, trying LOCPROJ"
try:
self.plocar.efermi
self.doscar.efermi = self.plocar.efermi
except NameError:
raise Exception("Efermi cannot be read from DOSCAR or LOCPROJ")
else:
# TODO: This a hack. Find out a way to determine ncdij without DOSCAR
print "!!! WARNING !!!: Error reading from DOSCAR, taking Efermi from config"
@ -103,16 +107,18 @@ class VaspData:
################################################################################
################################################################################
class Plocar:
r"""
"""
Class containing raw PLO data from VASP.
Properties
----------
Properties:
- *plo* (numpy.array((nion, ns, nk, nb, nlmmax))) : raw projectors
- *params* (dict) : parameters read from PLOCAR
- *ferw* (array(nion, ns, nk, nb)) : Fermi weights from VASP
- *plo* (numpy.array((nion, ns, nk, nb, nlmmax))) : raw projectors
- *params* (dict) : parameters read from PLOCAR
- *ferw* (array(nion, ns, nk, nb)) : Fermi weights from VASP
"""
def __init__(self):
self.plo = None
self.proj_params = None
def from_file(self, vasp_dir='./', plocar_filename='PLOCAR'):
r"""
@ -240,7 +246,10 @@ class Plocar:
self.nspin = 1 if self.ncdij == 1 else 2
self.nspin_band = 2 if self.ncdij == 2 else 1
self.efermi = float(sline[4])
try:
self.efermi = float(sline[4])
except:
print "!!! WARNING !!!: Error reading E-Fermi from LOCPROJ, trying DOSCAR"
plo = np.zeros((nproj, self.nspin, nk, self.nband), dtype=np.complex128)
proj_params = [{} for i in xrange(nproj)]
@ -323,14 +332,12 @@ class Poscar:
"""
Class containing POSCAR data from VASP.
Properties
----------
nq (int) : total number of ions
ntypes ([int]) : number of ion types
nions (int) : a list of number of ions of each type
a_brav (numpy.array((3, 3), dtype=float)) : lattice vectors
q_types ([numpy.array((nions, 3), dtype=float)]) : a list of
Properties:
- nq (int) : total number of ions
- ntypes ([int]) : number of ion types
- nions (int) : a list of number of ions of each type
- a_brav (numpy.array((3, 3), dtype=float)) : lattice vectors
- q_types ([numpy.array((nions, 3), dtype=float)]) : a list of
arrays each containing fractional coordinates of ions of a given type
"""
def __init__(self):
@ -350,7 +357,7 @@ class Poscar:
"""
# Convenince local function
def readline_remove_comments():
return f.next().split('!')[0].strip()
return f.next().split('!')[0].split('#')[0].strip()
# Add a slash to the path name if necessary
if vasp_dir[-1] != '/':
@ -443,16 +450,18 @@ class Kpoints:
"""
Class describing k-points and optionally tetrahedra.
Properties
----------
- nktot (int) : total number of k-points in the IBZ
- kpts (numpy.array((nktot, 3), dtype=float)) : k-point vectors (fractional coordinates)
- ntet (int) : total number of k-point tetrahedra
- itet (numpy.array((ntet, 5), dtype=float) : array of tetrahedra
- volt (float) : volume of a tetrahedron (the k-grid is assumed to
be uniform)
Properties:
- nktot (int) : total number of k-points in the IBZ
- kpts (numpy.array((nktot, 3), dtype=float)) : k-point vectors (fractional coordinates)
- ntet (int) : total number of k-point tetrahedra
- itet (numpy.array((ntet, 5), dtype=float) : array of tetrahedra
- volt (float) : volume of a tetrahedron (the k-grid is assumed to
be uniform)
"""
def __init__(self):
self.kpts = None
self.nktot = None
self.kwghts = None
#
# Reads IBZKPT file
#
@ -610,6 +619,10 @@ class Doscar:
"""
Class containing some data from DOSCAR
"""
def __init__(self):
self.ncdij = None
self.efermi = None
def from_file(self, vasp_dir='./', dos_filename='DOSCAR'):
"""
Reads only E_Fermi from DOSCAR.
@ -705,5 +718,3 @@ def read_symmcar(vasp_dir, symm_filename='SYMMCAR'):
data.update({ 'nrot': nrot, 'ntrans': ntrans,
'lmax': lmax, 'nion': nion,
'sym_rots': rot_mats, 'perm_map': rot_map })

162
python/converters/vasp_converter.py
View File

@ -43,9 +43,36 @@ class VaspConverter(ConverterTools):
dft_subgrp = 'dft_input', symmcorr_subgrp = 'dft_symmcorr_input',
parproj_subgrp='dft_parproj_input', symmpar_subgrp='dft_symmpar_input',
bands_subgrp = 'dft_bands_input', misc_subgrp = 'dft_misc_input',
transp_subgrp = 'dft_transp_input', repacking = False):
transp_subgrp = 'dft_transp_input', repacking = False,
proj_or_hk='proj'):
"""
Init of the class. Variable filename gives the root of all filenames, e.g. case.ctqmcout, case.h5, and so on.
Parameters
----------
filename : string
Base name of DFT files.
hdf_filename : string, optional
Name of hdf5 archive to be created.
dft_subgrp : string, optional
Name of subgroup storing necessary DFT data.
symmcorr_subgrp : string, optional
Name of subgroup storing correlated-shell symmetry data.
parproj_subgrp : string, optional
Name of subgroup storing partial projector data.
symmpar_subgrp : string, optional
Name of subgroup storing partial-projector symmetry data.
bands_subgrp : string, optional
Name of subgroup storing band data.
misc_subgrp : string, optional
Name of subgroup storing miscellaneous DFT data.
transp_subgrp : string, optional
Name of subgroup storing transport data.
repacking : boolean, optional
Does the hdf5 archive need to be repacked to save space?
proj_or_hk : string, optional
Select scheme to convert between KS bands and localized orbitals.
"""
assert type(filename)==StringType, "Please provide the DFT files' base name as a string."
@ -61,15 +88,22 @@ class VaspConverter(ConverterTools):
self.bands_subgrp = bands_subgrp
self.misc_subgrp = misc_subgrp
self.transp_subgrp = transp_subgrp
assert (proj_or_hk == 'proj') or (proj_or_hk == 'hk'), "proj_or_hk has to be 'proj' of 'hk'"
self.proj_or_hk = proj_or_hk
# Checks if h5 file is there and repacks it if wanted:
if (os.path.exists(self.hdf_file) and repacking):
ConverterTools.repack(self)
# this is to test pull request
def read_data(self, fh):
"""
Generator for reading plain data.
Parameters
----------
fh : file object
file object which is read in.
"""
for line in fh:
line_ = line.strip()
@ -82,6 +116,11 @@ class VaspConverter(ConverterTools):
def read_header_and_data(self, filename):
"""
Opens a file and returns a JSON-header and the generator for the plain data.
Parameters
----------
filename : string
file name of the file to read.
"""
fh = open(filename, 'rt')
header = ""
@ -121,12 +160,16 @@ class VaspConverter(ConverterTools):
SO = ctrl_head['nc_flag']
kpts = numpy.zeros((n_k, 3))
kpts_cart = numpy.zeros((n_k, 3))
bz_weights = numpy.zeros(n_k)
try:
for ik in xrange(n_k):
kx, ky, kz = rf.next(), rf.next(), rf.next()
kpts[ik, :] = kx, ky, kz
bz_weights[ik] = rf.next()
for ik in xrange(n_k):
kx, ky, kz = rf.next(), rf.next(), rf.next()
kpts_cart[ik, :] = kx, ky, kz
except StopIteration:
raise "VaspConverter: error reading %s"%self.ctrl_file
@ -148,17 +191,23 @@ class VaspConverter(ConverterTools):
jheader, rf = self.read_header_and_data(gr_file)
gr_head = json.loads(jheader)
e_win = gr_head['ewindow']
nb_max = gr_head['nb_max']
p_shells = gr_head['shells']
density_required = gr_head['nelect']
charge_below = 0.0 # This is not defined in VASP interface
# Note that in the DftTools convention each site gives a separate correlated shell!
n_shells = sum([len(sh['ion_list']) for sh in p_shells])
n_corr_shells = sum([len(sh['ion_list']) for sh in p_shells])
shells = []
corr_shells = []
shion_to_corr_shell = [[] for ish in xrange(len(p_shells))]
shion_to_shell = [[] for ish in xrange(len(p_shells))]
cr_shion_to_shell = [[] for ish in xrange(len(p_shells))]
shorbs_to_globalorbs = [[] for ish in xrange(len(p_shells))]
last_dimension = 0
crshorbs_to_globalorbs = []
icsh = 0
for ish, sh in enumerate(p_shells):
ion_list = sh['ion_list']
@ -168,16 +217,25 @@ class VaspConverter(ConverterTools):
# We set all sites inequivalent
pars['sort'] = sh['ion_sort'][i]
pars['l'] = sh['lorb']
#pars['corr'] = sh['corr']
pars['dim'] = sh['ndim']
#pars['ion_list'] = sh['ion_list']
pars['SO'] = SO
# TODO: check what 'irep' entry does (it seems to be very specific to dmftproj)
pars['irep'] = 0
corr_shells.append(pars)
shion_to_corr_shell[ish].append(i)
shells.append(pars)
shion_to_shell[ish].append(ish)
shorbs_to_globalorbs[ish].append([last_dimension,
last_dimension+sh['ndim']])
last_dimension = last_dimension+sh['ndim']
if sh['corr']:
corr_shells.append(pars)
# TODO: generalize this to the case of multiple shell groups
n_shells = n_corr_shells # No non-correlated shells at the moment
shells = corr_shells
n_corr_shells = len(corr_shells)
n_orbs = sum([sh['dim'] for sh in shells])
# FIXME: atomic sorts in Wien2K are not the same as in VASP.
# A symmetry analysis from OUTCAR or symmetry file should be used
@ -217,6 +275,7 @@ class VaspConverter(ConverterTools):
band_window = [numpy.zeros((n_k, 2), dtype=int) for isp in xrange(n_spin_blocs)]
n_orbitals = numpy.zeros([n_k, n_spin_blocs], numpy.int)
for isp in xrange(n_spin_blocs):
for ik in xrange(n_k):
ib1, ib2 = int(rf.next()), int(rf.next())
@ -227,10 +286,33 @@ class VaspConverter(ConverterTools):
hopping[ik, isp, ib, ib] = rf.next()
f_weights[ik, isp, ib] = rf.next()
if self.proj_or_hk == 'hk':
hopping = numpy.zeros([n_k, n_spin_blocs, n_orbs, n_orbs], numpy.complex_)
# skip header lines
hk_file = self.basename + '.hk%i'%(ig + 1)
f_hk = open(hk_file, 'rt')
# skip the header (1 line for n_kpoints, n_electrons, n_shells)
# and one line per shell
count = 0
while count < 3 + n_shells:
f_hk.readline()
count += 1
rf_hk = self.read_data(f_hk)
for isp in xrange(n_spin_blocs):
for ik in xrange(n_k):
n_orbitals[ik, isp] = n_orbs
for ib in xrange(n_orbs):
for jb in xrange(n_orbs):
hopping[ik, isp, ib, jb] = rf_hk.next()
for ib in xrange(n_orbs):
for jb in xrange(n_orbs):
hopping[ik, isp, ib, jb] += 1j*rf_hk.next()
rf_hk.close()
# Projectors
# print n_orbitals
# print [crsh['dim'] for crsh in corr_shells]
proj_mat = numpy.zeros([n_k, n_spin_blocs, n_corr_shells, max([crsh['dim'] for crsh in corr_shells]), numpy.max(n_orbitals)], numpy.complex_)
proj_mat_csc = numpy.zeros([n_k, n_spin_blocs, sum([sh['dim'] for sh in shells]), numpy.max(n_orbitals)], numpy.complex_)
# TODO: implement reading from more than one projector group
# In 'dmftproj' each ion represents a separate correlated shell.
@ -249,14 +331,37 @@ class VaspConverter(ConverterTools):
for isp in xrange(n_spin_blocs):
for ik in xrange(n_k):
for ion in xrange(len(sh['ion_list'])):
icsh = shion_to_corr_shell[ish][ion]
for ilm in xrange(sh['ndim']):
for ilm in xrange(shorbs_to_globalorbs[ish][ion][0],shorbs_to_globalorbs[ish][ion][1]):
for ib in xrange(n_orbitals[ik, isp]):
# This is to avoid confusion with the order of arguments
pr = rf.next()
pi = rf.next()
proj_mat[ik, isp, icsh, ilm, ib] = complex(pr, pi)
proj_mat_csc[ik, isp, ilm, ib] = complex(pr, pi)
# now save only projectors with flag 'corr' to proj_mat
proj_mat = numpy.zeros([n_k, n_spin_blocs, n_corr_shells, max([crsh['dim'] for crsh in corr_shells]), numpy.max(n_orbitals)], numpy.complex_)
if self.proj_or_hk == 'proj':
for ish, sh in enumerate(p_shells):
if sh['corr']:
for isp in xrange(n_spin_blocs):
for ik in xrange(n_k):
for ion in xrange(len(sh['ion_list'])):
icsh = shion_to_shell[ish][ion]
for iclm,ilm in enumerate(xrange(shorbs_to_globalorbs[ish][ion][0],shorbs_to_globalorbs[ish][ion][1])):
for ib in xrange(n_orbitals[ik, isp]):
proj_mat[ik,isp,icsh,iclm,ib] = proj_mat_csc[ik,isp,ilm,ib]
elif self.proj_or_hk == 'hk':
for ish, sh in enumerate(p_shells):
if sh['corr']:
for ion in xrange(len(sh['ion_list'])):
icsh = shion_to_shell[ish][ion]
for isp in xrange(n_spin_blocs):
for ik in xrange(n_k):
for iclm,ilm in enumerate(xrange(shorbs_to_globalorbs[ish][ion][0],shorbs_to_globalorbs[ish][ion][1])):
proj_mat[ik,isp,icsh,iclm,ilm] = 1.0
#corr_shell.pop('ion_list')
things_to_set = ['n_shells','shells','n_corr_shells','corr_shells','n_spin_blocs','n_orbitals','n_k','SO','SP','energy_unit']
for it in things_to_set:
# print "%s:"%(it), locals()[it]
@ -268,6 +373,8 @@ class VaspConverter(ConverterTools):
rf.close()
proj_or_hk = self.proj_or_hk
# Save it to the HDF:
with HDFArchive(self.hdf_file,'a') as ar:
if not (self.dft_subgrp in ar): ar.create_group(self.dft_subgrp)
@ -275,7 +382,9 @@ class VaspConverter(ConverterTools):
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']
'n_inequiv_shells', 'corr_to_inequiv', 'inequiv_to_corr','proj_or_hk','kpts','kpts_cart']
if self.proj_or_hk == 'hk' or True:
things_to_save.append('proj_mat_csc')
for it in things_to_save: ar[self.dft_subgrp][it] = locals()[it]
# Store Fermi weights to 'dft_misc_input'
@ -296,6 +405,24 @@ class VaspConverter(ConverterTools):
Reads input for the band window from bandwin_file, which is case.oubwin,
structure from struct_file, which is case.struct,
symmetries from outputs_file, which is case.outputs.
Parameters
----------
bandwin_file : string
filename of .oubwin/up/dn file.
struct_file : string
filename of .struct file.
outputs_file : string
filename of .outputs file.
misc_subgrp : string
name of the subgroup in which to save
SO : boolean
spin-orbit switch
SP : int
spin
n_k : int
number of k-points
"""
if not (mpi.is_master_node()): return
@ -390,6 +517,15 @@ class VaspConverter(ConverterTools):
def convert_symmetry_input(self, ctrl_head, orbits, symm_subgrp):
"""
Reads input for the symmetrisations from symm_file, which is case.sympar or case.symqmc.
Parameters
----------
ctrl_head : dict
dictionary of header of .ctrl file
orbits : list of shells
contains all shells
symm_subgrp : name of symmetry group in h5 archive
"""
# In VASP interface the symmetries are read directly from *.ctrl file

Some files were not shown because too many files have changed in this diff Show More