diff --git a/doc/analytic_continuation.rst b/doc/analytic_continuation.rst new file mode 100644 index 00000000..b438f62d --- /dev/null +++ b/doc/analytic_continuation.rst @@ -0,0 +1,19 @@ +Analytic continuation +===================== + +The analytic continuation is a method to obtain the real-frequency Green's function from the imaginary-frequency Green's function. +The program ``dcore_anacont`` performs the analytic continuation with the Pade approximation and the sparse-modeling method. + +.. toctree:: + :maxdepth: 1 + + analytic_continuation/pade + analytic_continuation/spm + + +Users can also use external codes for the analytic continuation. + +.. toctree:: + :maxdepth: 1 + + analytic_continuation/external_code diff --git a/doc/analytic_continuation/external_code.rst b/doc/analytic_continuation/external_code.rst new file mode 100644 index 00000000..7fefda0d --- /dev/null +++ b/doc/analytic_continuation/external_code.rst @@ -0,0 +1,55 @@ +How to use an external AC code +=============================== + +You can use an external code to perform AC of the self-energy from the Matsubara frequency to the real frequency. +The only thing the code needs to do is to write the self-energy in the real frequency as a NumPy binary file, ``post/sigma_w.npz`` +from the Matsubara frequency self-energy stored in ``seedname_sigma_iw.npz``. + +Input file: ``seedname_sigma_iw.npz`` +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +``seedname_sigma_iw.npz`` is a NumPy binary file, so you can load it with ``numpy.load`` function. + +.. code-block:: python3 + + import numpy as np + npz = np.load("./seedname_sigma_iw.npz") + +The returned object, ``npz``, is a dictionary. The keys are as follows: + +.. csv-table:: + :header: Key, Type, Description + :widths: 5, 10, 20 + + "beta", float, "Inverse temperature" + "iwn", array of complex, "Matsubara frequency" + "data#", array of complex, "Self-energy of #-th inequivalent shell" + "hartree\_fock#", array of complex, "Hartree-Fock term of #-th inequivalent shell" + +Here, "#" is 0, 1, 2, ... . + +"data#" is a :math:`N_{i\omega} \times N_\text{orb} \times N_\text{orb}` array, where :math:`N_{i\omega}` is the number of Matsubara frequencies, and :math:`N_\text{orb}` is the number of orbitals. + +"hartree\_fock#" is a Hartree-Fock term of "#"-th inequivalent shell, :math:`H^\text{f}`. + +.. math:: + + H^\text{f}_{ik} = \sum_{jl} U_{ijkl} \left\langle c^\dagger_j c_l \right\rangle + +The data format is a :math:`N_\text{orb} \times N_\text{orb}` array, where :math:`N_\text{orb}` is the number of orbitals. + +Output file: ``post/sigma_w.npz`` +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The output file, ``post/sigma_w.npz``, is also a NumPy binary file storing one dictionary with the following keys: + +.. csv-table:: + :header: Key, Type, Description + :widths: 5, 10, 20 + + "omega", array of real, "frequency" + "data#", array of complex, "Self-energy of #-th inequivalent shell" + +Here, "#" is 0, 1, 2, ... . + +"data#" is a :math:`N_{\omega} \times N_\text{orb} \times N_\text{orb}` array, where :math:`N_{\omega}` is the number of frequencies, and :math:`N_\text{orb}` is the number of orbitals. \ No newline at end of file diff --git a/doc/analytic_continuation/pade.ini b/doc/analytic_continuation/pade.ini new file mode 100644 index 00000000..d13ca9bb --- /dev/null +++ b/doc/analytic_continuation/pade.ini @@ -0,0 +1,39 @@ +[model] +seedname = square +lattice = square +norb = 1 +nelec = 1.0 +t = -1.0 +kanamori = [(4.0, 0.0, 0.0)] +nk0 = 8 +nk1 = 8 +nk2 = 1 + +[system] +T = 0.1 +n_iw = 1000 +fix_mu = True +mu = 2.0 + +[impurity_solver] +name = pomerol +exec_path{str} = pomerol2dcore +n_bath{int} = 3 +fit_gtol{float} = 1e-6 + +[control] +max_step = 100 +sigma_mix = 0.5 +time_reversal = true +converge_tol = 1e-5 + +[post.anacont] +solver = pade +show_result = false +save_result = true +omega_max = 6.0 +omega_min = -6.0 +Nomega = 401 + +[post.anacont.pade] +eta = 0.1 diff --git a/doc/analytic_continuation/pade.png b/doc/analytic_continuation/pade.png new file mode 100644 index 00000000..17789ed7 Binary files /dev/null and b/doc/analytic_continuation/pade.png differ diff --git a/doc/analytic_continuation/pade.rst b/doc/analytic_continuation/pade.rst new file mode 100644 index 00000000..7f3ea6f1 --- /dev/null +++ b/doc/analytic_continuation/pade.rst @@ -0,0 +1,63 @@ +Analytic continuation with the Pade approximation +==================================================== + +The Pade approximation is a simple and widely used method for the analytic continuation. +The Pade approximation is based on the assumption that the Green's function can be approximated by a rational function of the frequency. +The Pade approximation is implemented in the ``dcore_anacont`` program. + +The common parameters of ``dcore_anacont`` is specified in the ``[post.anacont]`` block as follows: + +.. code-block:: + + [post.anacont] + solver = pade + omega_min = -6.0 + omega_max = 6.0 + Nomega = 101 + show_result = false + save_result = true + +``solver`` specifies the method for the analytic continuation. +``omega_min``, ``omega_max``, and ``Nomega`` specify the frequency range and the number of frequency points for the output; +the frequency points are linearly spaced between ``omega_min`` and ``omega_max``. +``show_result`` specifies whether the obtained self-energies :math:`\mathrm{Re}\Sigma_{ij}(\omega)` and :math:`\mathrm{Im}\Sigma_{ij}(\omega)` are displayed on the screen by using Matplotlib. +``save_result`` specifies whether the plots of the obtained self-energies are saved in the ``work/post`` directory as ``sigma_w_{ish}_{iorb}_{jorb}.png``, +where ``ish`` is the shell index, and ``iorb`` and ``jorb`` are the orbital indices. +The output directory ``work/post`` can be changed by the parameter ``dir_work`` in the ``[post]`` block. + + +The parameters for the Pade approximation are specified in the ``[post.anacont.pade]`` block as follows: + +.. code-block:: + + [post.anacont.pade] + iomega_max = 10000 + n_min = 0 + n_max = 100 + eta = 0.01 + +``iomega_max``, ``n_min``, and ``n_max`` specify how many Matsubara frequencies are used. +First, ``iomega_max`` specifies the cutoff of Matsubara frequency. +When ``iomega_max > 0``, solve :math:`\omega_N = \pi/\beta (2N+1) \le \text{iomega_max} \lt \omega_{N+1}` and obtain :math:`N`. +If :math:`N < \text{n_min}` or :math:`N > \text{n_max}`, :math:`N` is replaced by :math:`\text{n_min}` or :math:`\text{n_max}`. +Otherwise (``iomega_max == 0``, default value), :math:`N = \text{n_max}`. +Then, the Matsubara frequencies :math:`\omega_n` for :math:`0 \le n \le N` are used. +When the number of self-energies calculated by the DMFT loop is less than :math:`N`, all the data are used. + +``eta`` is the imaginary frequency shift for the analytic continuation, that is, the analytic continuation is performed as :math:`i\omega_n \to \omega + i\eta`. + +Example +-------- + +The following input file is used to perform a DMFT calculation of a Hubbard model on a square lattice with the Pomerol solver and an analytic continuation of the self-energy with the Pade approximation: + +:download:`pade.ini ` + +.. literalinclude:: pade.ini + :language: ini + +The figure below shows the self-energy in real frequency obtained by analytic continuation (``work/post/sigma_w_0_0_0.png``). + +.. image:: pade.png + :width: 700 + :align: center diff --git a/doc/analytic_continuation/spm.ini b/doc/analytic_continuation/spm.ini new file mode 100644 index 00000000..8c969d1e --- /dev/null +++ b/doc/analytic_continuation/spm.ini @@ -0,0 +1,47 @@ +[model] +seedname = square +lattice = square +norb = 1 +nelec = 1.0 +t = -1.0 +kanamori = [(4.0, 0.0, 0.0)] +nk0 = 8 +nk1 = 8 +nk2 = 1 + +[system] +T = 0.1 +n_iw = 1000 +fix_mu = True +mu = 2.0 + +[impurity_solver] +name = pomerol +exec_path{str} = pomerol2dcore +n_bath{int} = 3 +fit_gtol{float} = 1e-6 + +[control] +max_step = 100 +sigma_mix = 0.5 +time_reversal = true +converge_tol = 1e-5 + +[post.anacont] +solver = spm +show_result = false +save_result = true +omega_max = 6.0 +omega_min = -6.0 +Nomega = 401 + +[post.anacont.spm] +n_matsubara = 1000 +n_tau = 101 +n_tail = 5 +n_sv = 30 +lambda = 1e-5 +solver = ECOS + +[post.anacont.spm.solver] +max_iters{int} = 100 diff --git a/doc/analytic_continuation/spm.png b/doc/analytic_continuation/spm.png new file mode 100644 index 00000000..5cb15fe4 Binary files /dev/null and b/doc/analytic_continuation/spm.png differ diff --git a/doc/analytic_continuation/spm.rst b/doc/analytic_continuation/spm.rst new file mode 100644 index 00000000..0691127c --- /dev/null +++ b/doc/analytic_continuation/spm.rst @@ -0,0 +1,59 @@ +Analytic continuation with the sparse-modeling method +======================================================== + +The sparse-modeling (SpM) method is a method for the analytic continuation of the self-energy by solving the integral equation + +.. math:: + + \Sigma(\tau) = \int_{-\infty}^{\infty} d\omega \frac{e^{-\tau \omega}}{1+e^{-\beta\omega}} \Sigma(\omega) + +where :math:`\Sigma(\tau)` is the self-energy in imaginary time and :math:`\Sigma(\omega)` is the self-energy in real frequency. + +In the SpM method, the kernel matrix of the integral equation, :math:`e^{-\tau\omega}/(1+e^{-\beta\omega})` is decomposed by the singular value decomposition (SVD), +and the self-energies :math:`\Sigma(\tau)` and :math:`\Sigma(\omega)` are transformed by the left and right singular vectors, respectively. + +The SpM method is implemented in the ``dcore_anacont`` program. +While the SpM method in the original paper uses the ADMM for optimization, the SpM method in the ``dcore_anacont`` program can use more solvers via `the CVXPY package `_. + +The common parameters of ``dcore_anacont`` is described in :doc:`the section of the Pade approximation `. +The SpM method is specified by ``solver = spm`` in the ``[post.anacont]`` block. +The parameters for the SpM method are specified in the ``[post.anacont.spm]`` block as follows: + +.. code-block:: + + [post.anacont.spm] + solver = ECOS + n_matsubara = 1000 + n_tau = 101 + n_tail = 5 + n_sv = 30 + lambda = 1e-5 + +``n_matsubara`` is the number of Matsubara frequencies used. When it is larger than the number of the data obtained by the DMFT calculation, the number of the data is used. +``n_tau`` is the number of imaginary time points. +``n_tail`` is the number of the last Matsubara frequencies used for the tail fitting, :math:`\Sigma(i\omega_n) \sim a/i\omega_n`. +``n_sv`` is the number of singular values used after truncation. +``lambda`` is the coefficient of the L1 regularization term in the optimization. +``solver`` specifies the solver for the optimization. + +The parameter of the solver can be specified in the ``[post.anacont.spm.solver]`` block. +In this block, the format of ``name{type} = value`` should be used, for example, ``max_iters{int} = 100``. +Available types are ``int``, ``float``, and ``str``. +Available parameters depend on the solver used. +For details of solvers, see `the CVXPY documentation `_. + +Example +-------- + +The following input file is used to perform a DMFT calculation of a Hubbard model on a square lattice with the Pomerol solver and an analytic continuation of the self-energy with the SpM method: + +:download:`spm.ini ` + +.. literalinclude:: spm.ini + :language: ini + +The figure below shows the self-energy in real frequency obtained by analytic continuation (``work/post/sigma_w_0_0_0.png``). + +.. image:: spm.png + :width: 700 + :align: center diff --git a/doc/basicnotions/structure.rst b/doc/basicnotions/structure.rst index 78010126..b12d6134 100644 --- a/doc/basicnotions/structure.rst +++ b/doc/basicnotions/structure.rst @@ -15,7 +15,7 @@ The structure of programs and data flow for the DMFT calculation is summarized b The DMFT calculation includes two **DCore** programs: (i) ``dcore_pre`` and (ii) ``dcore`` as described later. -After the DMFT loop is finished, one can compute dynamical physical quantities such as the density of states and the momentum-resolved spectrum functions using the post-processing tool. +After the DMFT loop (``dcore``) is finished, one can compute dynamical physical quantities such as the density of states and the momentum-resolved spectrum functions using the post-processing tool. The structure of the post-processing is shown below. .. image:: images/structure_post.png @@ -24,51 +24,13 @@ The structure of the post-processing is shown below. The post-processing tool consists of two **DCore** programs: (iii) ``dcore_anacont`` and (iv) ``dcore_spectrum``. -.. - **DCore** is a set of DMFT (Dynamical Mean Field Theory) programs which works together with other first-principles calculation packages. - **DCore** supports input hopping parameters in the wannier90 format. - Simple preset models such as a tight-binding model on the Bethe lattice are also available. - After the DMFT loop is finished, one can compute physical quantities such as the density of states and the momentum-resolved spectrum functions using the post-processing tool. - -.. - **DCore** consists of three layers: (i) interface layer, (ii) DMFT loop, and (iii) post-processing. - Those are respectively performed by the executables ``dcore_pre``, ``dcore``, and ``dcore_anacont`` and ``dcore_spectrum``. - Input parameters are provided by a single text file, which is read by all the three programs. - Data generated by ``dcore_pre`` and ``dcore`` are severally stored in a file with HDF5 format and passed to the next process. - - (i) The interface layer ``dcore_pre`` ------------------------------------- -.. - The pre-processing tool, ``dcore_pre`` can generate models from the wannier orbitals - as well as intrinsic model-generator (Standard interface). - ``dcore_pre`` generates a HDF5 file necessary for the DMFT loop. Users specify parameters defining a model such as hopping parameters on a certain lattice, and interactions. The hopping parameters are given either for **preset models** (e.g., square lattice, Bethe lattice) or using **Wannier90 format** -.. - There are two types of interfaces: **standard interface** for models and **Wannier90 interface** for materials. - - - Standard interface - - The preset models include - - * Bethe lattice - * Chain lattice - * Square lattice - * Cubic lattice - - .. - For more details, please see :ref:`inputformat`. - - - Wannier90 interface - - Results of ab-initio calculations can be imported into **DCore**. - ``dcore_pre`` reads the hopping parameters in the Wannier90 format and transform to HDF5. - - (ii) DMFT loop ``dcore`` ------------------------ diff --git a/doc/impuritysolvers/triqs_cthyb/cthyb.rst b/doc/impuritysolvers/triqs_cthyb/cthyb.rst index f1879475..2888098d 100644 --- a/doc/impuritysolvers/triqs_cthyb/cthyb.rst +++ b/doc/impuritysolvers/triqs_cthyb/cthyb.rst @@ -350,7 +350,7 @@ Finally, we choose the following setting: [control] max_step = 1 - [tool] + [post.check] omega_check = 30.0 and obtain diff --git a/doc/reference.rst b/doc/reference.rst index a07f130f..1d43f4e2 100644 --- a/doc/reference.rst +++ b/doc/reference.rst @@ -8,4 +8,5 @@ Reference Manual reference/input reference/output impuritysolvers + analytic_continuation tools diff --git a/doc/reference/input.rst b/doc/reference/input.rst index 278405d0..f91fe90b 100644 --- a/doc/reference/input.rst +++ b/doc/reference/input.rst @@ -1,19 +1,13 @@ Input-file format ================= -The input file consists of six parameter blocks named [model], [system], [impurity\_solver], [control], [tool], and [mpi]. - -.. - The example of the input file is shown as follows: - - .. literalinclude:: ../tutorial/nis/nis.ini - :language: ini +The input file consists of several parameter blocks. The following table shows which blocks are used by each program. .. csv-table:: - :header: "", ``dcore_pre``, ``dcore``, ``dcore_check``, ``dcore_anacont``, ``dcore_spectrum`` - :widths: 5, 5, 5, 5, 5, 5 + :header: "Block name", ``dcore_pre``, ``dcore``, ``dcore_check``, ``dcore_anacont``, ``dcore_spectrum`` + :widths: 10, 5, 5, 5, 5, 5 [model] , Yes, Yes, Yes, Yes, Yes [pre] , Yes, , , , diff --git a/doc/reference/output.rst b/doc/reference/output.rst index 49727408..cd4943eb 100644 --- a/doc/reference/output.rst +++ b/doc/reference/output.rst @@ -9,6 +9,7 @@ The output files generated by each program are summarized in the table below: ``dcore_pre``, "*seedname*.h5" ``dcore``, *seedname*.out.h5 + , *seedname*\_sigma\_iw.npz, :doc:`../analytic_continuation/external_code` , work/imp_shell#\_iter#/*solver_dependent_output* ``dcore_check``, check/iter\_mu.dat , check/iter\_mu.png, :doc:`programs` @@ -20,10 +21,12 @@ The output files generated by each program are summarized in the table below: , check/iter\_spin-ish0.png, :doc:`programs` , check/sigma.dat , check/sigma\_ave.png, :doc:`programs` - ``dcore_post``, post/*seedname*\_dos.dat - , post/*seedname*\_akw.dat - , post/*seedname*\_akw.gp, :doc:`programs` - , post/*seedname*\_momdist.dat + ``dcore_anacont``, post/sigma\_w.npz, :doc:`../analytic_continuation/external_code` + , post/sigma\_w.npz + ``dcore_spectrum``, post/dos.dat + , post/akw.dat + , post/akw.gp, :doc:`programs` + , post/momdist.dat Files that have empty Reference in this table are explained in the following. @@ -129,10 +132,10 @@ Files that have empty Reference in this table are explained in the following. -156.686934 0.994751 0.006371 0.994751 0.006371 : -``dcore_post`` -~~~~~~~~~~~~~~ +``dcore_spectrum`` +~~~~~~~~~~~~~~~~~~~~ -- **post/**\ *seedname*\ **_dos.dat** +- **post/dos.dat** The k-integrated single-particle excitation spectrum :math:`A(\omega)` (density of states). @@ -150,7 +153,7 @@ Files that have empty Reference in this table are explained in the following. -4.889724 0.011126 0.011126 0.011126 0.011126 : -- **post/**\ *seedname*\ **_akw.dat** +- **post/akw.dat** The single-particle excitation spectrum :math:`A(\boldsymbol{k}, \omega)` on the given k-path. See :ref:`program_dcore_post` for how to plot this data. @@ -169,7 +172,7 @@ Files that have empty Reference in this table are explained in the following. 0.000000 -4.751880 0.145834 : -- **post/**\ *seedname*\ **_momdist.dat** +- **post/momdist.dat** The momentum distribution function. diff --git a/doc/reference/programs.rst b/doc/reference/programs.rst index 4855c5db..062b57ef 100644 --- a/doc/reference/programs.rst +++ b/doc/reference/programs.rst @@ -1,25 +1,11 @@ Programs ======== -**DCore** consists of four main programs, ``dcore_pre``, ``dcore``, ``dcore_check``, and ``dcore_post``. +**DCore** consists of five main programs, ``dcore_pre``, ``dcore``, ``dcore_check``, ``dcore_anacont``, and ``dcore_spectrum``. All programs read a single input file (say *input.ini*) which contains parameters classified into blocks. It is also possible to pass multiple input files to the program. Parameters from all input files are merged. Duplicate paramegers are overwritten by the later one. See :doc:`input`, for the list of parameters and detailed descriptions. -.. - All programs can read input files of the same type and get the information by using blocks. - For details of input parameters defined in each block, see the next section. - -.. - ================= ======================================================= ==================== - Program Blocks to read from the input file Output HDF files - ================= ======================================================= ==================== - ``dcore_pre`` [model], [system] *seedname*.h5 - ``dcore`` [model], [system], [impurity-solver], [control], [mpi] *seedname*.out.h5 - ``dcore_check`` [model], [tool] --- - ``dcore_post`` [model], [system], [impurity-solver], [tool], [mpi] --- - ================= ================================================== ==================== - In the following, brief explanations are given for each program. Pre-processing : ``dcore_pre`` @@ -140,33 +126,46 @@ Three kinds of figures will be included: /\left[\sum_i^{\rm shell} N_{\rm orb}^{i}\right], The maximum frequency of this plot is specified with the parameter ``omega_check`` - in the ``[tool]`` block. + in the ``[post.check]`` block. .. Here, the average is taken over the shell index *i* and the orbital indices *a*, *b*. .. .. image:: ../tutorial/square/convergence.png -.. _program_dcore_post: +.. _program_dcore_anacont: -Post-processing : ``dcore_post`` -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Analytic continuation : ``dcore_anacont`` +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -This program computes the total DOS (*seedname*\_dos.dat) and momenum-resolved spectral function (*seedname*\_akw.dat) reading the DMFT results in *seedname*.out.h5. +This program performs the analytic continuation of the self-energy of Matsubara frequencies to the real frequency. +The self-energy of Matsubara frequencies is stored in the NumPy binary file *seedname*\_sigma\_iw.npz, which is one of the output files of the main program ``dcore``. -.. - This program reads the parameters defined in the ``[model]``, ``[system]``, ``[impurity-solver]`` and ``[tool]`` blocks. +.. code-block:: bash + + $ dcore_anacont input.ini + +The obtained self-energy on the real axis is stored in the NumPy binary file sigma\_w.npz in the ``post`` directory; the name of the directory can be changed by the parameter ``dir_post`` in the ``[post]`` block. + +.. _program_dcore_spectrum: + +Spectral functions : ``dcore_spectrum`` +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This program computes the total DOS (post/dos.dat) and momenum-resolved spectral function (post/akw.dat) reading the self-energy on the real axis stored in post/sigma\_w.npz. .. code-block:: bash - $ dcore_post input.ini --np 4 + $ dcore_spectrum input.ini --np 4 Here, please specify the number of MPI processes. +The output files are saved in the directory ``post``; the name of the directory can be changed by the parameter ``dir_post`` in the ``[post]`` block. The computed spectral function can be drawn by .. code-block:: bash - $ gnuplot [seedname]_akw.gp + $ cd post + $ gnuplot akw.gp Using this gnuplot script, you can also see the original (DFT) band structure as follows if either *seedname*\_band.dat (Wannier90 output) or dir-wan/dat.iband (RESPACK output) exists. diff --git a/doc/tutorial/afm/afm.rst b/doc/tutorial/afm/afm.rst index ef2788fa..93cc575d 100644 --- a/doc/tutorial/afm/afm.rst +++ b/doc/tutorial/afm/afm.rst @@ -87,7 +87,8 @@ Now, DMFT calculations can be done as usual. export MPIRUN="mpirun" dcore_pre cubic.ini > output-pre dcore cubic.ini --np 48 > output - dcore_post cubic.ini --np 48 > output-post + dcore_anacont cubic.ini > output-anacont + dcore_spectrum cubic.ini --np 48 > output-spectrum In the standard output of ``dcore``, you will see that the magnetic moments converge to :math:`\simeq 0.28` (57 % of the saturated moment).