RESPACK: An ab initio tool for derivation of effective low-energy model of material

RESPACK is a first-principles calculation software for evaluating the interaction parameters of materials and is able to calculate maximally localized Wannier functions, response functions based on the random phase approximation and related optical properties, and frequency-dependent electronic interaction parameters. RESPACK receives its input data from a band-calculation code using norm-conserving pseudopotentials with plane-wave basis sets. Automatic generation scripts that convert the band-structure results to the RESPACK inputs are prepared for xTAPP and Quantum ESPRESSO. An input file for specifying the RESPACK calculation conditions is designed pursuing simplicity and is given in the Fortran namelist format. RESPACK supports hybrid parallelization using OpenMP and MPI and can treat large systems including a few hundred atoms in the calculation cell.


Introduction
First-principles calculations based on density functional theory (DFT) [1,2] have currently been established and are widely used not only by theoretical researchers but also by experimental researchers. The users only input the crystal structure, and can easily evaluate the electronic and structural properties of materials. Density-functional calculations based on the local density approximation (LDA) are attractive due to its reasonable accuracy and low computational cost. Thanks to the recent advances in computer powers, it has successfully been applied to large-scale systems [3][4][5][6][7]. There are also many efforts to predict new materials by performing a huge number of DFT calculations [8,9]. On the other hand, it is well known that the LDA often fails to describe strongly-correlated electron systems [10]. Low-energy excitation and quantum fluctuations due to local interactions (on the order of several eV) dominate the low-energy properties of strongly correlated electron systems; the LDA cannot describe such quantum fluctuations. Therefore, the development of ab initio methods for strongly correlated materials has been an active research subject .
A well-known procedure for this attempt is the first-principle effective-model approach [10,32]. In this method, an effective model of a strongly correlated electron system is derived from first principles, and then the resulting effective model is numerically analyzed. Applying the highly accurate solvers which are able to evaluate electronic correlation accurately (exact diagonalization [33], dynamical mean field theory [34,35], manyvariable variational Monte Carlo method [36][37][38], configuration interaction method [39], etc.) to the ab initio derived effective model enables quantitative understanding of real strongly correlated materials. It is important to establish a reliable ab initio derivation method for the effective models, and many studies have been performed for this purpose . So far, development of the program was often closed at the laboratory level, and therefore, in order to spread the first-principle effectivemodel approach, software development and release are necessary.
In this paper, we introduce a software RESPACK [40] which contains a program group (maximally localized Wannier function and ab initio many body perturbation calculation) for deriving an effective low-energy model from first principles. The present paper focuses on a feature of the derivation tool for the effective model, but RESPACK has currently been extended to include the GW calculation, the spin-orbit interaction, and the electron-lattice coupling evaluation, which will be reported in the future. The present paper is organized as follows: In Section 2, we describe the methodological background. Some program details are given in Section 3. We give in Section 4 calculation procedure of RESPACK, and demonstrate in Section 5 how to install and compile of the source codes. In Section 6, we show a quantitative check for derived effective-model parameters of a cubic perovskite oxide SrVO 3 . A summary is given in Section 7. We give in Appendix A input details for the RESPACK calculation. In Appendix B, we introduce a utility tool for transfer analysis after the RESPACK calculations.

Effective low-energy model
We consider the derivation of the following extended Hubbard model within the two-center integrals written as H = σ RR i j t iR jR a σ † iR a σ jR + 1 2 σρ RR i j U iR jR a σ † iR a ρ † jR a ρ jR a σ iR + J iR jR a σ † iR a ρ † jR a ρ iR a σ jR + a σ † iR a ρ † iR a ρ jR a σ jR , where a σ † iR and a σ iR are creation and annihilation operators, respectively, of an electron with spin σ in the ith Wannier orbital in the lattice R. In this expression, the Wannier orbital is taken to be real. t iR jR is a transfer integral defined as t iR jR = φ iR |H 0 |φ jR = V drφ * iR (r)H 0 (r)φ jR (r ) (2) with |φ iR = a † iR |0 and the diagonal term (i = j, R = R ) being onsite energy. H 0 in Eq. (2) is the one-body part of H and is often taken to be the Kohn-Sham (KS) Hamiltonian H KS . The integral in Eq. (2) is taken over the crystal volume V.
The static limit of the screened direct-Coulomb U iR jR (ω) and exchange J iR jR (ω) integrals gives interaction parameters in the Hamiltonian H, which are given by and respectively. t iR jR , U iR jR (ω), and J iR jR (ω) have the lattice translational symmetry of t i0 jR −R = t iR jR , U i0 jR −R (ω) = U iR jR (ω), where we used φ iR (r) = φ i0 (r − R). In this paper, we focus on an ab initio derivation of these parameters.

Wannier function
The calculation of the Wannier function follows the algorithm for the maximally-localized Wannier function [41,42]. The ith Wannier function of the lattice R is defined as where k is a wave vector in the first Brillouin zone and N k is the total number of the Monkhorst-Pack k mesh. The Wannier function is constructed from the N k b KS bands [from (N k s + 1)-th to (N k s + N k b )-th bands]. N k s and N k b are determined from the energy-window information. U k αi is a matrix that transforms the αth Bloch wave function into the ith Wannier function. The αth Bloch wave function is defined as where G is a reciprocal lattice vector, and N ψ G is the total number of the plane waves used for the expansion of the wave function, which is determined by the cutoff energy E ψ cut from the inequality 1 2 |k + G| 2 ≤ E ψ cut . Ω is the volume of the unit cell. Both the Wannier function and the Bloch function are normalized for the crystal volume V = N k Ω. C Gα (k) is the expansion coefficient of the plane wave e i(k+G)·r / √ Ω. By inserting Eq. (11) into Eq. (10), the Wannier function at the home cell (R = 0) is written as withC Here,C Gi (k) is the expansion coefficient of the plane wave for the Wannier function. The center of the ith Wannier orbital at the lattice R is defined as Similarly, the spread of the Wannier orbital is defined as with r 2 iR = φ iR |r 2 |φ iR .

Response function with random phase approximation
In the RPA and constrained RPA, the polarization function in the plane-wave basis is written as with M G αβ (k + q, k) = ψ αk+q |e i(q+G)·r |ψ βk (17) and Here, q is a wave vector in the first Brillouin zone, ω is frequency, and indices α and β specify the unoccupied and occupied bands, respectively. The interstate matrix M G αβ (k + q, k) is evaluated using the fast Fourier transformation technique. E αk and δ in Eq. (18) are the energy of the Bloch state and the broadening factor, respectively. The quantity X αk+q,βk (ω) is calculated with the generalized tetrahedron technique [43,44] as where δ a and δ r in Eq.
where N w is the total number of the Wannier orbitals and U k αi is defined in Eq. (10). The quantity T αk is introduced to calculate constrained polarization [45], and, in the usual RPA, the T αk are set to zero.
The symmetric dielectric function [46] can be written with using the polarization function as When |q + G| and |q + G | are large, the contribution from the second term in the right hand side is negligible and GG (q, ω) ≈ δ GG . Therefore, the matrix is the number of the plane waves used to expand the polarization function, which is determined by the cutoff energy E cut as an inequality 1 2 In the calculation of GG (q, ω), we note on our special treatment of the head component which corresponds to the G = G = 0 component in the q → 0 limit. In the usual RPA case, we calculate the following [47] lim q→0 Here, the last term results from the intraband transition, and δ is the broadening factor introduced in Eq. (18). ω pl is the bare plasma frequency calculated via the Fermi-surface integral as where E F is the Fermi energy determined in the DFT band calculation, and p ν ααk is the diagonal element of the transitionmoment matrix with respect to the bands as with x µ being the Cartesian coordinate. On the above evaluation, we ignore the contribution from the non-local part of the pseudopotential, V NL . We note that this neglect is not so serious in the evaluation of the effective interaction, because the effective interaction is written in terms of the sum over the q points [see Eq. (30)]; in this case, the contribution from the q = 0 to the effective interaction becomes small relatively. On the other hand, in the optical response, the non-local pseudopotential contribution may manifest itself as a significant effect, especially for the transition metals [48]. This is because the optical properties are completely the q = 0 quantity [see Eqs. (26), which is obtained by dropping the last term in the right hand side of Eq. (22). The Wannier functions are constructed to include the low-energy bands near the Fermi level, and thus, T αk = 1 is expected and the bare plasma frequency ω µν pl in Eq. (23) becomes zero. This is why we dropped the last term in the right hand side of Eq. (22). The second derivative of the polarization function with respect to the wavenumber q in Eqs. (22) or (25) can be calculated analytically, as has been done for insulators [46].
The optical properties such as the macroscopic dielectric function M (ω), the electronic energy loss spectrum (EELS) L(ω), the real part of the optical conductivity σ(ω), and the reflectance spectrum R(ω) are also calculated from the inverse of the matrix GG (q, ω) in Eqs.(21) and (22) or (25) as and respectively.

Direct-Coulomb and exchange integrals
The evaluations of the interaction integrals U i0 jR (ω) in Eq. (8) and J i0 jR (ω) in Eq. (9) proceed as follows: First, the screened Coulomb interaction is written in the reciprocal space by using the Fourier transform as Note that the q-grid is the same as the k-grid; thus, N q = N k . W GG (q, ω) in the right-hand side is written with using the inverse dielectric matrix as follows: Note that the inverse dielectric matrix has off-diagonal elements in the first N G × N G block and becomes diagonal matrix with unity in the area beyond this block. By inserting Eqs. (30) and (31) into Eq. (3) and noting the lattice-translational symmetry [Eqs. (8) and (9)], we obtain the form of with and |φ ik = N R R |φ iR e ik·R with N R being the total number of the lattices in the system. The divergence of Eq. (33) in q → 0 with G=G'=0 is removed by following the prescription of Ref. [49].
The dependence of the static (ω = 0) direct-Coulomb integral on the distance between the two Wannier functions is evaluated via where r i jR is the distance between the two Wannier centers, Matrix elements of the bare (or unscreened) Coulomb interaction, U bare i0 jR = φ i0 φ i0 |v|φ jR φ jR , are calculated with replacing −1 GG (q, ω) of Eq. (32) by δ GG as The dependence of U bare i0 jR on the distance between two Wannier functions is obtained from The parallel argument can be applied to the derivation of the screened exchange integrals in Eq. (9). The result is with The bare exchange integral J bare i0 jR = φ i0 φ j0 |v|φ jR φ iR is given as 3. Some technical aspects 3.1. Symmetry RESPACK makes use of the space group symmetries of an input crystal structure and requires irreducible data obtained from a band calculation. The Bloch function at a reducible k point, k, is calculated from its irreducible part via where k * is the corresponding irreducible k point.R andT are the operators of a rotation and a fractional translation, respectively, which are represented by the 3×3 matrix R and threedimensional vector T. Equations (11) and (41) leads to with G * being reciprocal lattice vector for expansion of the wave function at the irreducible k point. In the above expression, we remove the global phase which does not depend on G * . By comparing Eq. (42) with Eq. (11), we find the following relations: Here, ∆ rw is a rewind vector which is introduced to pull back the rotated k * vector to the first Brillouin zone. From Eq. (44), G * and G have the following relationship and with Eqs. (45) and (46), we obtain More specifically, in the code, we treat S = (R −1 ) t instead of R, so the following expression is practically implemented Similarly, the inverse dielectric matrix at a reducible q point is generated from the irreducible one as follows:

Frequency grid
The frequency grid of the polarization function is generated as a logarithmic grid: Here, ω i is the ith frequency, N Ω is the total number of the frequency grids, and N max is the total number of the frequency grids in the frequency range 0 ≤ ω i ≤ E max . By default, N Ω and N max are set to 70 and (9N Ω )/10. The parameter s in Eq. (50) is determined by solving the equation s N max −1 + E max ∆ω (s−1) = 1. The resulting grids satisfy the following boundary conditions: (i) ω 1 = 0, (ii) ω 2 = ∆ω with ∆ω being 0.05 eV by default, (iii) ω N max = E max , and (iv) ω N Ω = 3E max . An example of the generated grid is shown in Fig. 1.

Interpolation treatment
Using the resulting transfer data t i0 jR in Sec. 2.2, the onebody Hamiltonian matrix at an arbitrary k-point k is calculated as Here, k is the k point employed in the band dispersion or the k point used in the Monkhorst-Pack mesh for the density of state calculation: Also, w R in Eq. 51 is a weight factor at the lattice R, which is introduced to avoid the double counting of the transfer at the boundary edge of the system with the periodic boundary condition. Note that w R satisfies the following sum rule R w R = N k . By diagonalizing the matrix H i j (k ), we obtain eigenvectors {C jα (k )} and eigenvalues αk . With αk , we can calculate the density of state as with being the partial density of state associated with the Wannier orbital φ i0 . The factor of 2 comes from the sum over spin degrees of freedom. The Brillouin Zone integral is performed with the generalized tetrahedron technique [43,44]. As a similar quantity, a density matrix is calculated as follows: This is convenient to monitor occupancy of each Wannier orbital or bond order between the Wannier orbitals. The Fermi surface is also calculated as constant energy surface F(k ) The output can be visualized by software Fermisurfer [50].

Parallel calculation
The polarization function can be calculated in parallel. There are two parallelization levels; one over the irreducible q points and the other over band pairs. First, let us consider the parallel calculation over the band pairs. To see this treatment, we rewrite the polarization function [Eq. (16)] as follows: The band sums are divided, and each χ αβ GG (q, ω) can be calculated by an independent MPI process. Now, we write this process as follows: Here, n and m specify an index of an MPI process. The divided virtual-state data are interchanged among the MPI processes to compute the partial-sum contribution to the polarization function. MPI SENDRECV routine is used for this data interchange. Figure 2 is a practical procedure for the case of N MPI = 2. The occupied-state and virtual-stat data are divided into two, and each data are stored in each MPI process (Step 1). After performing the polarization calculations in each MPI process (Step 2), only the virtual-state data {α 1 } and {α 2 } are interchanged between the two MPI processes (Step 3). Then, the polarization calculation is performed again (Step 4). Finally, the data stored in each MPI process is collected in the master process (Step 5).
A parallel calculation over the q points is more trivial. Consider the case where the total number of the MPI processes is 64, and the number of irreducible q points is 4. We first divide all the 64 MPI processes into 4 communities, and thus each Step 0: initial situation Step 1: data decomposition Step 2: calculate χ MPI#1 MPI#2 Occupied bands Virtual bands Step 3: data rotation Step 4: calculate MPI#1 MPI#2 χ α1β1 +χ α2β1 χ α2β2 +χ α1β2

MPI#1 MPI#2
Step 5: reduce on MPI#1 MPI#1 MPI#2 χ α1β1 +χ α2β1 χ α2β2 +χ α1β2 χ α2β2 +χ α1β2 χ + Figure 2: Practical procedure for an MPI calculation of polarization function. On the right side, work images are displayed, and the work at each step is highlighted in red. In this example, the total number of the MPI processes is 2. The band data are divided into two and accommodated in each MPI process (step 1). After performing polarization calculations in each MPI process (step 2), only virtual data are interchanged among the MPI processes (step 3). The polarization calculations are performed again (step 4). Finally, the partial-sum data of the polarization functions stored in each MPI process are collected in the master node (step 5).
6 community consists of 16 MPI processes. One community performs the polarization calculation of one q point. Figure. 3 is a schematic figure showing this procedure. MPI COMM SPLIT routine is used for splitting to the communities. The 16 MPI processes in each community are assigned to perform the parallel calculation over the band pairs mentioned above.

Calculation flow
We next describe the practical procedure of a RESPACK calculation [40]. Figure 4 shows an overall flow diagram of calculation processes; first, we perform band-structure calculations with xTAPP [51] or Quantum Espresso [52,53]. Next, with using the interface script, we convert the band-structure results to the inputs of RESPACK. Then, with the obtained bandcalculation data and an input file that specifies the RESPACKcalculation condition, we perform the Wannier-function calculation. We call this calculation wannier. Then, we calculate the polarization and dielectric functions, and this calculation is called chiqw. In the constrained RPA, the polarization process is restricted by using the information of the Wannier function [see Eqs. (16) and (20)], so one has to perform the wannier calculation before the chiqw calculation. Lastly, with the wannier and chiqw outputs, we evaluate the matrix elements of the screened interaction. This calculation is called calc int. In the following subsections, we describe details.

Preparation for RESPACK
We show in Fig. 5(a) a preparation process from bandstructure calculation with xTAPP or Quantum Espresso to RESPACK. In RESPACK, interface scripts that convert outputs of the band calculation to inputs for RESPACK are prepared for these two codes. In the case of xTAPP, xtapp2respack.sh generates a directory dir-wfn, in which the following 9 files are created.  Figure 5(b) is a flow diagram of wannier. With the data in dir-wfn and an input input.in that describes conditions of the wannier calculation, the calculation is performed with an executable file calc wannier. Details of input.in are described in Appendix A. After the calculation, two directories dir-wan and dir-model are generated, in which the calculation results are saved. For details of the generated output files, see the descriptions in Fig. 5

Chiqw calculation
We next show in Fig. 5(c) a flow diagram of chiqw. This code calculates the polarization and dielectric functions. With the data in dir-wfn and dir-wan and the input file input.in, the calculation is performed with an executable file calc chiqw. After the calculation, a directory dir-eps is generated and, under this directory, subdirectories q001, q002, ... qNirr are generated, where 001, 002, and Nirr are the numbers of the irreducible q points. The calculation results of every q points are saved in each subdirectory.

Figures 5(d) and (e) show flow diagrams for the direct-
Coulomb-integral and exchange-integral calculations, respectively. A common namelist &param calc int described in input.in can be used for the two calculation programs (see Appendix A.15). Executable files are calc w3d for the direct-Coulomb integral and calc j3d for the exchange integral. The calculation results are saved in the directories dir-intW, dir-intJ, and dir-model.

Download source files and compile
The source code of RESPACK can be obtained from the official website https://sites.google.com/view/kazuma7k6r. A gzipped tar file RESPACK.tar.gz contains everything necessary for installation. When moving to the directory src, one finds three source directories: calc int, chiqw, wannier. Makefile is prepared in each source directory, and one executes the make command to compile these source codes. In the case of chiqw, the work so far is as follows: > tar -zxvf RESPACK.tar.gz > cd RESPACK/src/chiqw/ > make Here, > is a prompt character. After make, an executable file calc chiqw is generated. This procedure is the same for the other programs.

Compile using cmake
RESPACK can also be compiled using CMake. In CMake, one needs to make a temporary directory for compilation, and executes the cmake and make commands from that directory as follows: If make is successful, executable files are generated in each directory under RESPACK/build/src.
By executing make install, the executable files will be installed in the bin directory under the directory specified by the -DCMAKE INSTALL PREFIX option.
If the -DCMAKE INSTALL PREFIX option is omitted, it will be installed under /usr/local/bin. The -DCONFIG option is used for reading the CMake configuration files stored in the RESPACK/config directory. For -DCONFIG, the following options are available: • intel: Intel compiler • gcc: GNU compiler If one wants to execute the cmake command again, it is recommended that one deletes the temporary directory build and restart from scratch, because the previous settings may remain.

Typical outputs
In this section, we show benchmark results for a t 2g -model derivation of perovskite oxide SrVO 3 with a simple cubic structure having a lattice constant of 3.8425 Å. Density functional calculations with plane-wave basis sets were performed using the xTAPP code [51], where the norm-conserving pseudopotential [54,55] and the generalized gradient approximation to the exchange correlation energy were employed [56]. The calculation condition is set to 8 × 8 × 8 k-point sampling, 100-Ry wavefunction cutoff, and 400-Ry charge-density cutoff. The 50 bands are considered for the polarization function, which corresponds to considering the excitation from the Fermi level to 35 eV. The numbers of the doubly-occupied, partially-occupied, and unoccupied bands are 12, 3, and 35, respectively. Table 1 is an input file input.in for RESPACK calculations. Details of input.in are described in Appendix A. We construct the t 2g -type Wannier functions from the low-energy bands near the Fermi level. The calculations are performed with both of the constrained RPA and usual RPA to show the difference between these two. The cutoff for the polarization function is set to 10 Ry, and the broadening factor δ of the generalized tetrahedron calculation is set to 0.1 eV. Figure 6 is a comparison between the original KS band (redsolid curves) and the Wannier-interpolated band (green-dashed curves). A region between the two blue-dashed horizontal lines indicates the energy window used to construct the t 2g -type Wannier functions. We also show in Fig. 7 the calculated d xy Wannier function in realspace.
We show in Table 2 important transfers for the t 2g band, where the definition of t 1 , t 2 , t 3 , and t 4 are illustrated in Fig. 8. In this figure, we depict the d xy Wannier function as an example. Since the lattice of the system is simple cubic, there exist equivalent transfers for the d yz and d zx orbitals. RESPACK provides a utility code that searches the equivalent transfers, which is described in Appendix B in more detail. We note that the original band structure in Fig. 6 are well reproduced by these four transfers. Figure 9 shows calculated macroscopic dielectric functions [Eq. (26)] with the RESPACK-chiqw code. Panels (a) and (b) describe the constrained RPA and usual RPA results, respectively. Red-solid and green-dashed curves describe the real and imaginary parts, respectively, and circles represent calculation values. The difference between the constrained RPA and usual RPA spectra is appreciable in the low-energy excitation region less than 2-3 eV. In the constrained RPA, a metallic charge excitation is excluded by the polarization constraint described in Sec. 2.3, and then the real part of M (ω) converges to the finite  dir-eps     (2), where we show main 4 transfer integrals. Definition for t 1 , t 2 , t 3 , and t 4 are given in Fig. 8. These four transfers successfully reproduce the original t 2g -band structure in Fig. 6. The unit of transfer integral is eV. value in the ω → 0 limit, while, in the usual RPA, the real part of M (ω) diverges negatively due to the metallic charge excitation [58], thus leading to the Drude behavior of the imaginary part of M (ω). It should be noted here that the present spectra neglect the transition moment contributed from the commutation relation between the non-local pseudopotential and electronic position, [V NL , r], in Eq. (24). For transition metals, this contribution is known to affect the spectral property in the low-energy excitation region [48]. In the present SrVO 3 , we checked that this contribution is not significant. Support for this contribution is a future issue in the RESPACK project.
We next show in Fig. 10 other optical properties calculated with the usual RPA and constrained RPA. Panels (a), (b), and (c) display EELS L(ω) [Eq. (27)], the real part of the optical con-ductivity σ(ω) [Eq. (28)], and the reflectance spectrum R(ω) [Eq. (29)], respectively. Red-solid and green-dashed curves represent the results based on the usual RPA and constrained RPA, respectively. Circles denote the calculation values. There is a difference between the constrained RPA and usual RPA in the EELS around the low-energy excitation region less than 2 eV; in the usual RPA case, an additional peak appears in L(ω), which is due to the low-energy plasmon excitation in the t 2g band [14,59]. This plasma excitation is also observed in the reflectance spectrum R(ω); we see a sharp drop from around 1 to 0 in the RPA reflectance spectrum. On the optical conductivity, the spectral trend is basically the same as the macroscopic dielectric function M (ω); the usual RPA spectrum exhibits the Drude behavior characteristic of a metallic system in the low-excitation region less than 2 eV, while the constrained RPA spectrum has no intensity in this frequency region, which is a characteristic aspect of the insulating system.  Figure 12 shows a frequency dependence of the onsite direct-Coulomb integral [U 1010 (ω) in Eq. (32) with orbital index 1 denoting d xy -type orbital]. Panels (a) and (b) represent the constrained RPA and usual RPA results, respectively. Red-solid and green-dashed curves describe the real and imaginary parts, respectively. As well as the optical data, the difference between the constrained RPA and usual RPA occurs in the low-energy excitation region below 2 eV due to the low-energy plasmon excitation considered in the usual RPA calculation. Figure 13 is a frequency dependence of the onsite exchange integral. [J 1020 (ω) in Eq. (38) with orbital indices 1 and 2 denoting d xy -type and d yz -type orbitals, respectively]. In exchange integral, differences between the constrained RPA (a) and usual RPA (b) can also be observed around ω ∼ 2 eV.

Convergence check to calculation conditions
We next show a convergence behavior of the calculated interaction parameters of the t 2g model of SrVO 3 with respect to the various computational conditions. Table 3 shows a convergence of the static interaction parameters with increasing the sampling k-point density, where we list the static onsite-intraorbital interaction U, onsite-interorbital U , onsite-exchage J, and nearestneighbor V interactions averaged over orbitals. We see that the convergence is achieved around 8 × 8 × 8. Table 4 is a dependence of the static interaction parameters of SrVO 3 on the total number of bands. We see that the 50 bands is enough to obtain the converged results. Table 5 gives a convergence behavior of the static interaction parameters with respect to the cutoff energy E cut of the polarization function. This parameter is important for the convergence and it is desirable to take large enough. The convergence within       Table 6 is a dependence of the static onsite interaction parameter on the broadening factor δ introduced in the polarizationfunction calculation of Eqs. (16), (18), and (19). We see that the δ does not affect the static cRPA and RPA results. We note that the δ may affect dynamical properties, especially in the low-energy collective excitation. Figure 14 is the δ dependence of EELS function L(ω) in Eq. (27). As the δ value increases [0.001 eV (black curves), 0.01 eV (purple curves), 0.05 eV (blue curves), 0.1 eV (green curves), and 0.2 eV (red curves)], the intensity of the plasmon peak around 1-2 eV decreases and eventually the peak position shifts to the lower energy. The plasma excitation is also sensitive to the k-point density [47]. Therefore, one should be careful about computational conditions for the quantitative discussion of the dynamical properties.

Pseudopotential dependence
Here, we discuss the pseudopotential dependence of the interaction parameters in details. The pseudopotential depends  (19)]. The calculation is based on the usual RPA. The spectra drawn by black, purple, blue, green, and red indicates the spectra with δ = 0.001 eV, 0.01 eV, 0.05 eV, 0.1 eV, and 0.2 eV, respectively. A spectral behavior around 1-2 eV results from the plasmon excitation within the t 2g band. As the δ increases, the peak intensity decreases and the peak position shifts to the lower-energy side eventually. mainly on the cutoff radius r loc for local pseudopotential. The pseudopotential with a small r loc parameter makes deeper potential, and the resulting pseudo wavefunction tends to be more localized near the ion core. We constructed the five pseudopotentials with the different r loc values (0.8, 1.0, 1.5, 1.8 and 2.1 bohr) for a vanadium atom. The pseudopotential with r loc = 0.8, 1.0, and 1.5 bohr were constructed for an ionic semicore configuration of (3s) 2 (3p) 6 (3d) 3 . The pseudopotentials with r loc = 1.8 and 2.1 bohr were constructed with an ionic valence configuration of (3d) 3 (4s) 0 (4p) 0 . The Troullier-Martins (TM) type was adopted as a function form of the pseudo wavefunction [55]. 8 × 8 × 8 k-point sampling and wavefunction cutoff E ψ cut of 196 Ry are employed. The total number of bands N band is 50, and the broadening factor δ was set to be 0.1 eV. A polarizationfunction cutoff E cut is important for the effective interaction parameters, so the convergence behavior are discussed for this parameter. We note that the t 2g -band structures obtained with the above 5 pseudopotentials are in almost perfect agreements. Figure 15 compares the atomic pseudo wavefunctions obtained with the different r loc parameters with the atomic allelectron wavefunction (black curve). As the r loc parameter is reduced from 2.1 bohr (light blue) → 1.8 bohr (purple) → 1.5 bohr (blue) → 1.0 bohr (green) → 0.8 bohr (red), the maximum amplitude position of the pseudo wavefunction is shifted to the ion-core side. The pseudo wavefunctions with the r loc = 0.8 and 1.0 bohr are almost the same as the all-electron wave function.
This trend can affect the localization of the Wannier function; the Wannier function generated with a pseudopotential with a small r loc cutoff tends to be more localized. As a result, it can give a large bare (unscreened) direct-Coulomb integral. Table 7 shows the r loc dependence of the bare interaction parameters for the t 2g Wannier function of SrVO 3 . We see from the table that the r loc and the Wannier spread S i0 in Eq. (15) has clear positive correlation; a smaller r loc leads to a smaller Wannier spread. As a result, the r loc and bare onsite direct-Coulomb integrals U bare and U bare correlate negatively (i.e., the small r loc brings about the large U bare and U bare ). There are no discernible effects on the bare exchange J bare and the bare nearest-neighbor direct integrals V bare . Tables 8 is the pseudopotential dependence of constrained-RPA and usual-RPA interaction parameters. An interesting trend can be seen in the table; details of the pseudopotential hardly affect the interaction values in contrast to the bare (unscreened) interaction parameters. In general, however, sufficiently large E cut would be desirable for a safer quantitative discussion about effective interaction parameters, especially for the usual RPA case.
Finally, we mention the effects of the pseudopotential type. As the famous pseudopotential types, besides the TM type, there are the ONCV (Optimized Norm-Conserving Vanderbilt) type [60] and the RRKJ (Rappe-Rabe-Kaxiras-Joannopoulos) type [61]. Even if the same cutoff r loc is employed, the results can be quantitatively different due to the difference in the functional form of the pseudopotential, so the user should be Table 8: Dependence of static (ω=0) interaction parameters of SrVO 3 on pseudopotential, where we change the r loc parameter which is a cutoff radius for the local pseudopotential. Note that the pseudopotentials with r loc = 0.8, 1.0, and 1.5 bohr are constructed for the semicore configuration and those with r loc = 1.8 and 2.1 bohr are constructed for the valence configuration (see the text). Under each pseudopotential, the E cut cutoff dependence on the interaction parameters is investigated. Other calculation condition is as follows: 8 × 8 × 8 k-point sampling, N band = 50, E ψ cut = 196 Ry, and δ = 0.1 eV. The unit of the interaction parameter is eV. careful about this point. Table 9 compares the calculated interaction parameters based on the TM-type, ONCV-type [62] and RRKJ-type [63] pseudopotentials. The calculations with the TM-type pseudopotential were performed with xTAPP, which is referred to as TM-xTAPP. The calculations with the ONCVtype and RRKJ-type pseudopotentials were performed with Quantum Espresso, which are referred to as ONCV-QE and RRKJ-QE, respectively. For the TM-xTAPP, we give two results TM(v)-xTAPP and TM(s)-xTAPP, for which the former is the results based on the pseudopotential constructed with the valence-electron configuration, and the latter is the results based on the pseudopotential with the semicore configuration for V and Sr. The ONCV and RRKJ pseudopotentials are also constructed for the semicore configurations for V and Sr. We found that, for the screened direct-Coulomb interaction U and U , RRKJ-QE gives significantly smaller values than others.

Conclusion
In conclusion, we present a new software RESPACK for deriving effective low-energy models of materials from first principles. The software contains programs for computing maximally-localized Wannier functions, RPA response functions, and matrix elements for screened interaction with respect to the Wannier functions. RESPACK is freely available under the GNU General Public Licence [64]. RESPACK is written in Fortran90 and supports plane-wave DFT codes xTAPP [65] and Quantum Espresso [52,53] , for which an interface script that converts the band-calculation results to the RESPACK inputs is provided. As an important notice, RESPACK currently supports ab initio codes with the norm-conserving pseudopotential. The present paper focuses on the derivation of the ef-fective model, but it has currently been extended to the GW calculation, effective-model derivation with spin-orbit interaction, and electron-lattice coupling evaluation. We will report on these additional features in near future.

Acknowledgments
We thank Yoshiro Nohara for providing a module for generalized tetrahedron calculation. We also acknowledge Maxime Charlebois and Jean-Baptiste Morée for useful discussions about the code development. We thank Masatoshi Imada, Ryotaro Arita, and Takashi

Appendix A. Input files
In this appendix, we describe details of the input file input.in for specifying the calculation condition of RESPACK. The input.in is given in the Fortran namelist format. Table A.10 lists six namelists available in input.in. The namelist &param wannier contains variables for the Wannier-function calculation. The namelists &param interpolation and &param visualization describe variables for the Wannier-interpolation band and the visualization of the realspace Wannier functions, respectively. These three namelists are read by a common executable file calc wannier. This executable file is generated by compilation of source code in the directory RESPACK/src/wannier. The namelist &param chiqw contains parameters for calculations of polarization and dielectric functions and are read by an executable file calc chiqw which is created by compilation of source code in RESPACK/src/chiqw. The namelist &param calc int describes variables for direct-Coulomb and exchange-integral calculation and are read by executable files calc w3d and calc j3d. These are made by compilation in the source code in RESPACK/src/calc int.  Table A.11 shows main four variables in the namelist &param wannier. N wannier is the total number of the Wannier function to be calculated. Lower energy window and Upper energy window specify the energy window and the Bloch bands in the energy window are used for the Wannierfunction calculation. The user has to set proper values by seeing the band dispersion obtained with DFT band calculations. N initial guess is the number of initial Gaussian for the Wannier-function calculation. The user must enter initial-guess information just below this namelist (see Table 1). vec ini(1:N initial guess) in Table A.11 are variables to specify the initial guess and is treated as type in Fortran. vec ini(i)%orb specifies the orbital type of the ith initial guess, where %orb can take a character of s, px, py, pz, dxy, dyz, dzx, dx2, dz2. The initial guess is treated with the Gaussian function, exp(−α i |r − R i | 2 ), where α i =%a is orbital exponent of the ith Gaussian. The Gaussian position R i = xa 1 + ya 2 + za 3 with x=%x, y=%y, and z=%z are given in the fractional coordinates, where a 1 , a 2 , and a 3 are the lattice vectors in the calculation cell. Note that one has to enter N initial guess pieces of the vec init information. Specific descriptions are given in Table 1. Table A.11: Main variables in namelist &param wannier to construct maximally-localized Wannier functions. The variable type and its brief explanation are shown in the second and third columns, respectively. The lower five variables in vec ini are information about initial guesses. Note that vec ini is treated as type of Fortran. The initial guess is specified by the orbital type "s, px, py, pz, dxy, dyz, dzx, dx2, dz2" and parameters in the Gaussian function, exp(−α i |r − R i | 2 ). Here, α i is an orbital exponent of the ith Gaussian, and R i = xa 1 +ya 2 +za 3 is the Gaussian position in the fractional coordinates, where a 1 , a 2 , and a 3 are the lattice vectors in the calculation cell.  Tables A.12 summarize main variables in the namelist &param interpolation to draw the Wannier-interpolated band. The user specifies the total number of symmetric k points in the dispersion line as N sym points, and, just below this namelist, describes the coordinates {(s 1 , s 2 , s 3 )} of the symmetric k vectors (see the Table 1 for concrete descriptions). Here, k i = s 1 b 1 + s 2 b 2 + s 3 b 3 with b 1 , b 2 , and b 3 being the basic reciprocal lattice vectors. These variables are read as array SK sym pts(1:3,1:N sym points). Usually, the same k points as those employed in the DFT band calculation are adopted.  We show in Table A.13 parameters in the namelist &param visualization. By default, RESPACK skips the visualization calculation of the realspace Wannier function. So, the user must set Flg vis wannier to 1 when calculating visualization data. Structural data of the same spatial range as the Wannier-function data is also output as cif (Crystallographic Information File) format. All the results can be drawn by software VESTA [57]. Next, we list in Table A.14 main parameters of the namelist &param chiqw. In this namelist, the user specifies the calculation condition of the polarization and dielectric functions. Ecut for eps is the cutoff energy for the polarization function. Default value is set to 1/10 of wave-function cutoff E ψ cut . Num freq grid is the total number of the frequency grid. By default, a log grid of 70 points is generated. The maximum excitation energy E max is estimated from the total bands considered in the polarization function, and the range of the frequency grid is set to 3E max to consider the Lorentz tail of the polarization function. See Section 3.2 for details. Green func delt is smearing value used in the generalized tetrahedron calculation (0.1 eV by default). In RESPACK, the polarization calculation is performed in parallel. To this end, two variables MPI num qcomm and MPI num proc per qcomm are prepared. MPI num qcomm is a variable for parallel computation on the q points and defines the number of qcommunity. MPI num proc per qcomm specifies the total number of MPI processes in each q-community. Note that when mod(N MPI,MPI num qcomm)/=0, the chiqw program stops, where N MPI is the total number of the MPI processes and this value is set in the job script. By default, MPI num qcomm=1 and MPI num proc per qcomm=N MPI. Details of the parallel calculation can be found in Section. 3.4. Lastly, Flg cRPA is a variable to specify the constrained RPA calculation. In the default setting, Flg cRPA=0 and the program performs a usual RPA calculation. transfers, and k-grid size (--kgrd) considered in the banddispersion and density-of-states calculations. Fermi-surface data generation for Fermisurfer --diff D threshold for detecting equivalent transfers (0.01 eV) --ecut E cut transfers below E (0.00 eV) --rcut R cut transfers over R (100 Å) --elec N set electron numbers in unit cell to N (no setting) --delta δ set broadening for DOS calculation to δ (0.01 eV) --kgrd='n m l set k-grid to n × m × l (the same as the band calculation)  Table 2), and therefore the band dispersion and density of states are appreciably changed from others. As another example, we show in Fig. B.18 the k-grid density dependence of the Fermi surface. This is an Al result. Accurate description for a small Fermi surface colored by red needs dense k-grid density. The calculation is performed as > python tr.py --frm --kgrd='n m l' with n, m, and l being numbers of k-grid along the basic reciprocal lattice vectors b 1 , b 2 , and b 3 , respectively.