Branch data Line data Source code
1 : : !!****m* ABINIT/m_dyson_solver
2 : : !! NAME
3 : : !! m_dyson_solver
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains procedures to solve the Dyson equation to find QP energies.
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2008-2026 ABINIT group (MG)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : : module m_dyson_solver
23 : :
24 : : use defs_basis
25 : : use m_xmpi
26 : : use m_errors
27 : : use m_abicore
28 : : use m_dtfil
29 : :
30 : : use m_time, only : timab
31 : : use m_gwdefs, only : sigparams_t
32 : : use m_numeric_tools, only : linfit, pade, dpade, newrap_step
33 : : use m_io_tools, only : open_file
34 : : use m_fstrings, only : int2char10
35 : : use m_hide_lapack, only : xheev
36 : : use m_dtset, only : dataset_type
37 : : use m_bz_mesh, only : kmesh_t
38 : : use m_sigma, only : sigma_t
39 : : use m_melemts, only : melements_t
40 : :
41 : : implicit none
42 : :
43 : : private
44 : : !!***
45 : :
46 : : public :: solve_dyson ! Solve the Dyson equation for the QP energies.
47 : :
48 : : !----------------------------------------------------------------------
49 : :
50 : : !!****t* m_dyson_solver/sigma_pade_t
51 : : !! NAME
52 : : !! sigma_pade_t
53 : : !!
54 : : !! FUNCTION
55 : : !! Object to perform the analytic continuation with Pade' and
56 : : !! find the QP solution with Newton-Rapson method.
57 : : !!
58 : : !! SOURCE
59 : :
60 : : type, public :: sigma_pade_t
61 : :
62 : : integer :: npts
63 : : ! Number of points
64 : : integer :: nsig_ab
65 : : ! Number of spinor components
66 : :
67 : : real(dp) :: betar_pm(2), zcut_pm(2)
68 : : complex(dp) :: alphac_pm(2)
69 : : logical :: do_sigma_fit
70 : :
71 : : complex(dp),allocatable :: zmesh(:)
72 : : ! input mesh
73 : :
74 : : complex(dp),allocatable :: sigc_cvals(:,:)
75 : : ! values on mesh
76 : :
77 : : contains
78 : :
79 : : procedure :: init_default => sigma_pade_init
80 : : procedure :: init_spinor => sigma_pade_init_spinor
81 : : generic :: init => init_default, init_spinor
82 : : ! Init object
83 : :
84 : : procedure :: free => sigma_pade_free
85 : : ! Free memory
86 : :
87 : : procedure :: eval => sigma_pade_eval
88 : : ! Evaluate self-energy and derivative along the real axis.
89 : :
90 : : procedure :: qp_solve => sigma_pade_qp_solve
91 : : ! Find the QP solution with Newton-Rapson method
92 : :
93 : : end type sigma_pade_t
94 : : !!***
95 : :
96 : : integer,private,parameter :: NR_MAX_NITER = 1000
97 : : ! Max no of iterations in the Newton-Raphson method.
98 : :
99 : : real(dp),private,parameter :: NR_ABS_ROOT_ERR = 0.0001/Ha_eV
100 : : ! Tolerance on the absolute error on the Newton-Raphson root.
101 : :
102 : : CONTAINS !====================================================================
103 : : !!***
104 : :
105 : : !!****f* m_dyson_solver/solve_dyson
106 : : !! NAME
107 : : !! solve_dyson
108 : : !!
109 : : !! FUNCTION
110 : : !! Solve the Dyson equation for the QP energies. Two different methods are coded:
111 : : !! The first one is based on the standard perturbative approach in which the self-energy
112 : : !! is linearly expanded around the previous single-particle energy (KS energy if one-shot)
113 : : !! and the derivative is evaluated by finite differences.
114 : : !! In the second method (AC), the values of the self-energy operator on the real axis are obtained
115 : : !! by means of an analytic continuation based on the Pade extrapolation.
116 : : !!
117 : : !! INPUTS
118 : : !! ikcalc=Index of the considered k-point in the Sigp%kptgw2bz array.
119 : : !! nomega_sigc=Number of frequencies used to evaluate the correlation part of Sigma.
120 : : !! Sigp<sigparams_t>=Structure gathering parameters on the calculation of Sigma.
121 : : !! Sr<sigma_t>=Structure containing the matrix elements of the self-energy INOUT
122 : : !! Kmesh<kmesh_t>=Info on the K-mesh for the wavefunctions.
123 : : !! sigcme=(nomega_sigc,ib1:ib2,ib1:ib2,nsppol)=Matrix elements of Sigma_c.
124 : : !! qp_ene(nbnds,nkibz,nsppol)= KS or QP energies, only used in case of calculation with scissor operator.
125 : : !! comm=MPI communicator.
126 : : !!
127 : : !! OUTPUT
128 : : !! Sr<sigma_t>=Structure containing the matrix elements of the self-energy:
129 : : !! %sigxme(ib1:ib2,jkibz,nsspol)=Diagonal elements of Sigma_x
130 : : !! %sigcmee0(ib1:ib2,jkibz,nsppol)=Matrix elements of Sigma_c at the initial energy E0.
131 : : !! %dsigmee0(jb,ib1:ib2,nsppol)=Derivate of sigma at the energy E0.
132 : : !! %ze0(ib1:ib2,jkibz,is)=Renormalization factor at the energy E0.
133 : : !! %degw(ib1:ib2,jkibz,is)= QP correction i.e DeltaE_GW=E-E0
134 : : !! %egw(ib1:ib2,jkibz,is)=QP energy
135 : : !! %sigmee(ib1:ib2,jkibz,is)=Self-energy evaluated at the QP energy.
136 : : !! %sigcme (ib1:ib2,jkibz,io,is)= Sigma_c as a function of frequency.
137 : : !! %sigxcme(ib1:ib2,jkibz,io,is)= Sigma_xc as a function of frequency.
138 : : !! %sigcme4sd (ib1:ib2,jkibz,io,is)= Diagonal matrix elements of \Sigma_c at frequencies around the KS eigenvalue
139 : : !! %sigxcme4sd(ib1:ib2,jkibz,io,is)= Diagonal matrix elements of \Sigma_xc at frequencies around the KS eigenvalue
140 : : !! where ib1 and ib2 are the band indices included in the GW calculation for this k-point.
141 : : !!
142 : : !! SOURCE
143 : :
144 : 614 : subroutine solve_dyson(ikcalc, minbnd, maxbnd, nomega_sigc, dtset, Sigp, Kmesh, sigcme, qp_ene, Sr, ks_me, Dtfil, comm)
145 : :
146 : : !Arguments ------------------------------------
147 : : !scalars
148 : : integer,intent(in) :: ikcalc,nomega_sigc,minbnd,maxbnd,comm
149 : : type(dataset_type),intent(in) :: dtset
150 : : type(sigparams_t),intent(in) :: Sigp
151 : : type(kmesh_t),intent(in) :: Kmesh
152 : : type(Datafiles_type),intent(in) :: Dtfil
153 : : type(sigma_t),intent(inout) :: Sr
154 : : type(melements_t),intent(in) :: ks_me
155 : : !arrays
156 : : real(dp),intent(in) :: qp_ene(Sr%nbnds,Sr%nkibz,Sr%nsppol)
157 : : complex(dp),intent(in) :: sigcme(nomega_sigc,minbnd:maxbnd,minbnd:maxbnd,Sigp%nsppol*Sigp%nsig_ab)
158 : :
159 : : !Local variables-------------------------------
160 : : !scalars
161 : : integer,parameter :: master=0
162 : : integer :: iab,ib1,ib2,ikbz_gw,io,spin,is_idx,isym,itim,jb, ie0, ierr
163 : : integer :: ik_ibz,kb,ld_matrix,mod10,nsploop,my_rank, units(2)
164 : : real(dp) :: alpha, beta, smrt, vxc_val, vu, v_meanf, sigx
165 : : complex(dp) :: dsigc, sigc, sigc_zsc, zz, zsc, phase
166 : : logical :: ltest
167 : : character(len=500) :: msg
168 : 614 : type(sigma_pade_t) :: spade
169 : : !arrays
170 : : real(dp) :: kbz_gw(3),tsec(2), betar_pm(2), zcut_pm(2)
171 : 614 : real(dp),allocatable :: e0pde(:),eig(:),scme(:)
172 : : complex(dp) :: alphac_pm(2)
173 : 614 : complex(dp),allocatable :: hdp(:,:),tmpcdp(:),hhartree(:,:,:),htotal(:,:,:),h_tmp1(:,:),h_tmp2(:,:)
174 : : ! *************************************************************************
175 : :
176 : : DBG_ENTER("COLL")
177 : :
178 : 614 : units = [std_out, ab_out]
179 : :
180 : 614 : call timab(490,1,tsec) ! csigme(Dyson)
181 : :
182 : 614 : my_rank = xmpi_comm_rank(comm)
183 : :
184 : 614 : mod10 = MOD(Sigp%gwcalctyp,10)
185 : :
186 : 614 : ltest=(nomega_sigc==Sr%nomega_r+Sr%nomega4sd)
187 [ + + ]: 614 : if (mod10==1) ltest=(nomega_sigc==Sr%nomega_i)
188 [ - + ]: 614 : ABI_CHECK(ltest,'Wrong number of frequencies')
189 : :
190 : : ! Index of the KS or QP energy.
191 : : !ioe0j=Sr%nomega4sd/2+1
192 : :
193 : : ! min and Max band index for GW corrections (for this k-point).
194 [ + + ]: 1236 : ib1 = MINVAL(Sigp%minbnd(ikcalc,:))
195 [ + + ]: 1236 : ib2 = MAXVAL(Sigp%maxbnd(ikcalc,:))
196 : :
197 : : ! Find the index of the k-point for sigma in the IBZ array.
198 : 614 : ikbz_gw = Sigp%kptgw2bz(ikcalc)
199 : 614 : call kmesh%get_BZ_item(ikbz_gw, kbz_gw, ik_ibz, isym, itim, phase)
200 : :
201 : 614 : sigc = czero; dsigc = czero
202 : :
203 : : ! ===========================================================
204 : : ! ==== Solve the Dyson Equation and store results in Sr% ====
205 : : ! ===========================================================
206 : :
207 [ + + ]: 614 : if (mod10 /= 1) then
208 : : ! ===============================
209 : : ! ==== Perturbative approach ====
210 : : ! ===============================
211 : :
212 : : ! Index of the KS or QP energy in sigme_tmp
213 : 607 : ie0 = sr%nomega_r + Sr%nomega4sd/2+1
214 : :
215 [ + + ]: 1222 : do spin=1,Sr%nsppol
216 [ + + ]: 6414 : do jb=ib1,ib2
217 : : ! Get matrix elements of Sigma_c at energy E0.
218 : : ! SigC(w) is linearly interpolated and the slope alpha is assumed as dSigC/dE
219 [ + + ]: 10816 : do iab=1,Sr%nsig_ab
220 [ + + ]: 5624 : is_idx = spin; if (Sr%nsig_ab>1) is_idx=iab
221 : :
222 : 5624 : Sr%sigcmee0(jb,ik_ibz,is_idx) = sigcme(ie0,jb,jb,is_idx)
223 : :
224 [ - + + - : 16872 : ABI_MALLOC(scme, (Sr%nomega4sd))
+ - + - -
+ - + ]
225 [ - + + - : 11248 : ABI_MALLOC(e0pde, (Sr%nomega4sd))
- + - + ]
226 [ + + ]: 25440 : e0pde(:) = Sr%omega4sd(jb,ik_ibz,:,spin)
227 [ + + ]: 25440 : scme(:) = REAL(sigcme(Sr%nomega_r+1:Sr%nomega_r+Sr%nomega4sd,jb,jb,is_idx))
228 : :
229 [ + + ]: 5624 : if (Sr%nomega4sd==1) then
230 : 3846 : smrt = zero; alpha = zero
231 : : else
232 : 1778 : smrt = linfit(Sr%nomega4sd,e0pde(:),scme(:),alpha,beta)
233 : : end if
234 : :
235 [ + + ]: 5624 : if (smrt > 0.1/Ha_eV) then
236 : : write(msg,'(3a,i0,a,i0,2a,2(f22.15,2a))')&
237 : 65 : 'WARNING: Values of Re Sig_c(omega) are not linear ',ch10,&
238 : 65 : 'band index: ',jb,' spin|component: ',is_idx,ch10,&
239 : 65 : 'root mean square: ',smrt,ch10,&
240 : 65 : 'estimated slope: ',alpha,ch10,&
241 : 130 : 'Omega [eV] SigC [eV]'
242 : 65 : ABI_WARNING(msg)
243 [ + + ]: 650 : do io=1,Sr%nomega4sd
244 : 585 : write(msg, '(2f8.4)')e0pde(io)*Ha_eV,scme(io)*Ha_eV
245 : 650 : call wrtout(std_out, msg)
246 : : end do
247 : : end if
248 : :
249 [ - + ]: 5624 : ABI_FREE(scme)
250 [ - + ]: 5624 : ABI_FREE(e0pde)
251 : : !
252 : : ! === Evaluate renormalization factor and QP correction ===
253 : : ! * Z=(1-dSigma/domega(E0))^-1
254 : : ! * DeltaE_GW=E-E0= (Sigma(E0)-V_xc)/(1-dSigma/domega)
255 : : ! * If nspinor==2, this part is done at the end.
256 : : !
257 : 5624 : Sr%dsigmee0(jb,ik_ibz,is_idx)=CMPLX(alpha,zero)
258 : :
259 [ + + ]: 5624 : if (Sr%nsig_ab==1) then
260 : 5048 : Sr%ze0(jb,ik_ibz,spin)= one / (one-Sr%dsigmee0(jb,ik_ibz,spin))
261 : :
262 [ + + ]: 5048 : if (ABS(Sigp%mbpt_sciss) < tol6) then
263 : : Sr%degw(jb,ik_ibz,spin) = Sr%ze0(jb,ik_ibz,spin) * &
264 : : (Sr%sigxme(jb,ik_ibz,spin) + Sr%sigcmee0(jb,ik_ibz,spin) - Sr%e0(jb,ik_ibz,spin) + &
265 : 5044 : Sr%hhartree(jb,jb,ik_ibz,spin))
266 : :
267 : 5044 : Sr%egw(jb,ik_ibz,spin) = Sr%e0(jb,ik_ibz,spin) + Sr%degw(jb,ik_ibz,spin)
268 : :
269 : : ! Estimate Sigma at the QP-energy: Sigma(E_qp)=Sigma(E0)+(E_qp-E0)*dSigma/dE
270 : : Sr%sigmee(jb,ik_ibz,spin) = &
271 : 5044 : Sr%sigxme(jb,ik_ibz,spin)+Sr%sigcmee0(jb,ik_ibz,spin)+Sr%degw(jb,ik_ibz,spin)*Sr%dsigmee0(jb,ik_ibz,spin)
272 : :
273 : : else
274 : : ! If GW+scissor: e0 is replaced by qp_ene which contains the updated energy eigenvalue
275 : : Sr%degw(jb,ik_ibz,spin)= Sr%ze0(jb,ik_ibz,spin) * &
276 : : (Sr%sigxme(jb,ik_ibz,spin) + Sr%sigcmee0(jb,ik_ibz,spin) - qp_ene(jb,ik_ibz,spin) + &
277 : 4 : Sr%hhartree(jb,jb,ik_ibz,spin))
278 : :
279 : 4 : Sr%egw(jb,ik_ibz,spin) = qp_ene(jb,ik_ibz,spin) + Sr%degw(jb,ik_ibz,spin)
280 : :
281 : : ! Estimate Sigma at the QP-energy: Sigma(E_qp)=Sigma(E0)+(E_qp-E0)*dSigma/dE
282 : : Sr%sigmee(jb,ik_ibz,spin)= &
283 : : Sr%sigxme(jb,ik_ibz,spin) + Sr%sigcmee0(jb,ik_ibz,spin) + &
284 : 4 : Sr%degw(jb,ik_ibz,spin) * Sr%dsigmee0(jb,ik_ibz,spin)
285 : :
286 : : ! RS: In the output, the gw corr with respect to e0 without mbpt_sciss is reported.
287 : 4 : Sr%degw(jb,ik_ibz,spin) = Sr%egw(jb,ik_ibz,spin) - Sr%e0(jb,ik_ibz,spin)
288 : : end if
289 : : end if !Sigp%nsig_ab==1
290 : :
291 : : ! Spectrum of Sigma
292 [ + + ]: 8614 : do io=1,Sr%nomega_r
293 : 2990 : Sr%sigcme (jb,ik_ibz,io,is_idx)= sigcme(io,jb,jb,is_idx)
294 : 8614 : Sr%sigxcme(jb,ik_ibz,io,is_idx)= Sr%sigxme(jb,ik_ibz,is_idx)+Sr%sigcme(jb,ik_ibz,io,is_idx)
295 : : end do
296 [ + + ]: 30632 : do io=1,Sr%nomega4sd
297 : 19816 : Sr%sigcme4sd (jb,ik_ibz,io,is_idx)= sigcme(Sr%nomega_r+io,jb,jb,is_idx)
298 : 25440 : Sr%sigxcme4sd(jb,ik_ibz,io,is_idx)= Sr%sigxme(jb,ik_ibz,is_idx)+Sr%sigcme4sd(jb,ik_ibz,io,is_idx)
299 : : end do
300 : : end do !iab
301 : :
302 [ + + ]: 5807 : if (Sr%nsig_ab > 1) then
303 [ - + ]: 144 : ABI_CHECK(ABS(Sigp%mbpt_sciss)<0.1d-4,'Scissor with spinor not coded')
304 : : !TODO this should be allocated with nsppol, recheck this part
305 : :
306 : : ! Evaluate renormalization factor and QP correction.
307 : : ! Z=(1-dSigma/domega(E0))^-1
308 : : ! DeltaE_GW=E-E0= (Sigma(E0)-V_xc)/(1-dSigma/domega)
309 : : !write(std_out,'(a,i2,10f8.3)')' Correlation',jb,Sr%sigcmee0(jb,ik_ibz,:)*Ha_eV,SUM(Sr%sigcmee0(jb,ik_ibz,:))*Ha_eV
310 : :
311 [ + + ]: 720 : Sr%ze0 (jb,ik_ibz,1) = one/(one-SUM(Sr%dsigmee0(jb,ik_ibz,:)))
312 : :
313 : : Sr%degw(jb,ik_ibz,1) = Sr%ze0(jb,ik_ibz,1) * &
314 [ + + ]: 720 : (SUM(Sr%sigxme(jb,ik_ibz,:)+Sr%sigcmee0(jb,ik_ibz,:)+Sr%hhartree(jb,jb,ik_ibz,:))-Sr%e0(jb,ik_ibz,1))
315 : :
316 : 144 : Sr%egw(jb,ik_ibz,1)=Sr%e0(jb,ik_ibz,1)+Sr%degw(jb,ik_ibz,1)
317 : :
318 : : ! Estimate Sigma at the QP-energy.
319 [ + + ]: 720 : do iab=1,Sr%nsig_ab
320 : : Sr%sigmee(jb,ik_ibz,iab)= &
321 : 720 : Sr%sigxme(jb,ik_ibz,iab)+Sr%sigcmee0(jb,ik_ibz,iab)+Sr%degw(jb,ik_ibz,1)*Sr%dsigmee0(jb,ik_ibz,iab)
322 : : end do
323 : : end if
324 : :
325 : : end do ! jb
326 : : end do ! spin
327 : :
328 : : else
329 : : ! =============================
330 : : ! === Analytic Continuation ===
331 : : ! =============================
332 : : ! ABI_CHECK(Sr%nsig_ab == 1, "AC with spinor not implemented")
333 : :
334 : : ! Index of the KS or QP energy in sigme_tmp
335 : : !ie0 = sr%nomega_r + Sr%nomega4sd/2+1
336 : :
337 [ + + ]: 14 : do spin=1,Sr%nsppol
338 [ + + ]: 67 : do jb=ib1,ib2
339 [ - + + - : 159 : ABI_MALLOC(tmpcdp,(Sr%nomega_i))
+ - + - -
+ - + ]
340 : : ! Calculate Sigc(E0), dSigc(E0)
341 : 53 : zz = CMPLX(Sr%e0(jb,ik_ibz,spin), zero)
342 : :
343 [ - + ]: 53 : if (Sigp%mbpt_sciss > 0.1d-4) then
344 : : ! e0 is replaced by qp_ene which contains the updated energy eigenvalue
345 : 0 : zz = CMPLX(qp_ene(jb,ik_ibz,spin), zero)
346 : : end if
347 : :
348 : : ! Diagonal elements of sigcme
349 [ + + ]: 106 : do iab=1,Sr%nsig_ab
350 [ - + ]: 53 : is_idx=spin; if (Sr%nsig_ab>1) is_idx=iab
351 : 53 : alphac_pm = zero; betar_pm = zero; zcut_pm = zero
352 : 53 : call spade%init(sr%nomega_i, sr%omega_i, sigcme(:,jb,jb,is_idx), alphac_pm, betar_pm, zcut_pm)
353 : 53 : call spade%eval(zz, Sr%sigcmee0(jb,ik_ibz,is_idx), dvdz=Sr%dsigmee0(jb,ik_ibz,is_idx))
354 : 106 : call spade%free()
355 : : end do !iab
356 : :
357 : : ! Z = (1 - dSigma / domega(E0))^{-1}
358 [ + - ]: 53 : if (Sr%nsig_ab == 1) then
359 : 53 : Sr%ze0(jb,ik_ibz,spin) = one / (one - Sr%dsigmee0(jb,ik_ibz,spin))
360 : : else
361 [ # # ]: 0 : Sr%ze0(jb,ik_ibz,1) = one / (one - SUM(Sr%dsigmee0(jb,ik_ibz,:)))
362 : : end if
363 : :
364 : : !#define _DEV_PERTURBATIVE
365 : : !#ifdef _DEV_PERTURBATIVE
366 : : #if 0
367 : : call wrtout(units, "COMMENT: Using perturbative approach with Z.")
368 : :
369 : : ! Note vxc[n_val] instead of vxc[n_val + n_nlcc] with the model core charge.
370 : : ! With gw_sigxcore=1 (PAW), the core contribution is treated with Fock exchange: use vxc - sxcore,
371 : : ! consistent with Sr%hhartree used by the other solvers.
372 : : if (Sigp%use_sigxcore == 1) then
373 : : vxc_val = ks_me%vxc(jb, jb, ik_ibz, spin) - ks_me%sxcore(jb, jb, ik_ibz, spin)
374 : : if (Sr%nsig_ab > 1) vxc_val = SUM(ks_me%vxc(jb, jb, ik_ibz, :) - ks_me%sxcore(jb, jb, ik_ibz, :))
375 : : else
376 : : vxc_val = ks_me%vxcval(jb, jb, ik_ibz, spin)
377 : : if (Sr%nsig_ab > 1) vxc_val = SUM(ks_me%vxcval(jb, jb, ik_ibz, :))
378 : : end if
379 : : vu = zero; if (dtset%usepawu /= 0) vu = ks_me%vu(jb, jb, ik_ibz, spin)
380 : : v_meanf = vxc_val + vu
381 : :
382 : : ! qp_ene = e0 + z_e0 * (sigc_e0__ + sigx - v_meanf)
383 : : if (Sr%nsig_ab == 1) then
384 : : Sr%egw(jb,ik_ibz,spin) = Sr%e0(jb,ik_ibz,spin) + Sr%ze0(jb,ik_ibz,spin) * &
385 : : (Sr%sigcmee0(jb,ik_ibz,spin) + Sr%sigxme(jb,ik_ibz,spin) - v_meanf)
386 : :
387 : : Sr%degw(jb,ik_ibz,spin) = Sr%egw(jb,ik_ibz,spin) - Sr%e0(jb,ik_ibz,spin)
388 : :
389 : : ! Estimate Sigma at the QP-energy: Sigma(E_qp)=Sigma(E0)+(E_qp-E0)*dSigma/dE
390 : : Sr%sigmee(jb,ik_ibz,spin) = &
391 : : Sr%sigxme(jb,ik_ibz,spin)+Sr%sigcmee0(jb,ik_ibz,spin)+Sr%degw(jb,ik_ibz,spin)*Sr%dsigmee0(jb,ik_ibz,spin)
392 : : else
393 : : Sr%egw(jb,ik_ibz,1) = Sr%e0(jb,ik_ibz,1) + Sr%ze0(jb,ik_ibz,1) * &
394 : : (SUM(Sr%sigcmee0(jb,ik_ibz,:)+Sr%sigxme(jb,ik_ibz,:))-v_meanf)
395 : :
396 : : Sr%degw(jb,ik_ibz,1) = Sr%egw(jb,ik_ibz,1) - Sr%e0(jb,ik_ibz,1)
397 : :
398 : : Sr%sigmee(jb,ik_ibz,1) = &
399 : : SUM(Sr%sigxme(jb,ik_ibz,:)+Sr%sigcmee0(jb,ik_ibz,:))+Sr%degw(jb,ik_ibz,1)*SUM(Sr%dsigmee0(jb,ik_ibz,:))
400 : : end if
401 : :
402 : : #else
403 : : ! MG FIXME: Here we are solving the non-linear QP equation using the Pade' continuation + root finding
404 : : ! but this is very misleading because in the output file we are still reporting the Z factor
405 : : ! and there's no mention that the QP energies have been obtained from the non-linear equation!!
406 : : ! One should change the format used to print the results or at least warn the user!
407 : :
408 : 53 : zz = CMPLX(Sr%e0(jb,ik_ibz,spin), zero)
409 : :
410 [ - + ]: 53 : if (Sigp%mbpt_sciss>0.1d-4) then
411 : : ! e0 is replaced by qp_ene which contains the updated energy eigenvalue.
412 : 0 : zz = CMPLX(qp_ene(jb,ik_ibz,spin), zero)
413 : : end if
414 : :
415 : : ! Solve the QP equation with Newton-Rapson starting from e0
416 : : ! Find root of E^0-V_xc-V_U+Sig_x+Sig_c(z)-z, i.e E^qp.
417 : 53 : alphac_pm = zero; betar_pm = zero; zcut_pm = zero
418 [ - + ]: 53 : if (Sr%nsig_ab > 1) then
419 [ # # # # : 0 : call spade%init(sr%nomega_i, sr%omega_i, Sr%nsig_ab, sigcme(:,jb,jb,:), alphac_pm, betar_pm, zcut_pm)
# # # # #
# ]
420 : : else
421 : 53 : call spade%init(sr%nomega_i, sr%omega_i, sigcme(:,jb,jb,spin), alphac_pm, betar_pm, zcut_pm)
422 : : end if
423 : :
424 : : ! Note vxc[n_val] instead of vxc[n_val + n_nlcc] with the model core charge.
425 : : ! With gw_sigxcore=1 (PAW), the core contribution is treated with Fock exchange: use vxc - sxcore,
426 : : ! consistent with Sr%hhartree used by the other solvers.
427 [ - + ]: 53 : if (Sigp%use_sigxcore == 1) then
428 : 0 : vxc_val = ks_me%vxc(jb, jb, ik_ibz, spin) - ks_me%sxcore(jb, jb, ik_ibz, spin)
429 [ # # # # ]: 0 : if (Sr%nsig_ab > 1) vxc_val = SUM(ks_me%vxc(jb, jb, ik_ibz, :) - ks_me%sxcore(jb, jb, ik_ibz, :))
430 : : else
431 : 53 : vxc_val = ks_me%vxcval(jb, jb, ik_ibz, spin)
432 [ - + - - ]: 53 : if (Sr%nsig_ab > 1) vxc_val = SUM(ks_me%vxcval(jb, jb, ik_ibz, :))
433 : : end if
434 [ - + ]: 53 : vu = zero; if (dtset%usepawu /= 0) vu = ks_me%vu(jb, jb, ik_ibz, spin)
435 : 53 : v_meanf = vxc_val + vu
436 [ + - ]: 53 : if (Sr%nsig_ab == 1) then
437 : 53 : sigx = Sr%sigxme(jb,ik_ibz,spin)
438 : : else
439 [ # # ]: 0 : sigx = SUM(Sr%sigxme(jb,ik_ibz,:))
440 : : end if
441 : :
442 : 53 : call spade%qp_solve(sr%e0(jb,ik_ibz,spin), v_meanf, sigx, zz, zsc, sigc_zsc, msg, ierr)
443 : 53 : call spade%free()
444 : : !qpe_pade_kcalc(ibc, ikcalc, spin) = zsc
445 : : !qp_solver_ierr(ibc, ikcalc, spin) = ierr
446 [ - + ]: 53 : if (ierr /= 0) then
447 : 0 : ABI_WARNING(msg)
448 : : end if
449 : :
450 : : ! Store the final result (self-consistent result for zz and Sigma_c(zz_scf)
451 : 53 : Sr%egw(jb,ik_ibz,spin) = zsc
452 : 53 : Sr%degw(jb,ik_ibz,spin) = Sr%egw(jb,ik_ibz,spin) - Sr%e0(jb,ik_ibz,spin)
453 : 53 : Sr%sigmee(jb,ik_ibz,spin) = Sr%sigxme(jb,ik_ibz,spin) + sigc_zsc
454 : : #endif
455 : :
456 : : ! Spectra of Sigma, remember that Sr%nomega_r does not contain the frequencies
457 : : ! used to evaluate the derivative each frequency is obtained using the pade_expression
458 : : ! In sigma indeed we have:
459 : : ! nomega_sigc=Sr%nomega_r+Sr%nomega4sd; if (mod10==SIG_GW_AC) nomega_sigc=Sr%nomega_i
460 [ + + ]: 803 : do io=1,Sr%nomega_r
461 : 750 : zz=Sr%omega_r(io)
462 [ + + ]: 750 : if (REAL(zz) > zero) then
463 [ + + + + ]: 7875 : tmpcdp(:) = SUM(sigcme(:,jb,jb,:), DIM=2)
464 : 375 : Sr%sigcme(jb,ik_ibz,io,spin) = pade(Sr%nomega_i, Sr%omega_i, tmpcdp, zz)
465 : : else
466 [ + + + + ]: 7875 : tmpcdp(:) = CONJG(SUM(sigcme(:,jb,jb,:), DIM=2))
467 [ + - + + : 4125 : Sr%sigcme(jb,ik_ibz,io,spin) = pade(Sr%nomega_i, CONJG(Sr%omega_i), tmpcdp, zz)
- + ]
468 : : end if
469 : 803 : Sr%sigxcme(jb,ik_ibz,io,spin) = Sr%sigxme(jb,ik_ibz,spin) + Sr%sigcme(jb,ik_ibz,io,spin)
470 : : end do
471 : :
472 : : ! Save sigma values along the imaginary axis
473 [ + + ]: 106 : do iab=1,Sr%nsig_ab
474 [ - + ]: 53 : is_idx=spin; if (Sr%nsig_ab > 1) is_idx = iab
475 [ + + ]: 716 : do io=1,Sr%nomega_i
476 : 610 : Sr%sigcmesi (jb,ik_ibz,io,is_idx) = sigcme(io,jb,jb,is_idx)
477 : 663 : Sr%sigxcmesi(jb,ik_ibz,io,is_idx) = Sr%sigxme(jb,ik_ibz,is_idx) + Sr%sigcmesi(jb,ik_ibz,io,is_idx)
478 : : end do
479 : : end do
480 : :
481 [ - + ]: 60 : ABI_FREE(tmpcdp)
482 : : end do !jb
483 : : end do !is
484 : :
485 : : end if ! Analytic continuation.
486 : :
487 : : ! === Diagonalize the QP Hamiltonian (forced to be Hermitian) ===
488 : : ! Calculate Sr%en_qp_diago and Sr%eigvec_qp to be written in the QPS file.
489 : : ! TODO in case of AC results are wrong.
490 : :
491 [ + + ]: 614 : if (mod10 /= 1) then
492 [ - + + - : 3035 : ABI_MALLOC(hhartree, (ib1:ib2,ib1:ib2,Sr%nsppol*Sr%nsig_ab))
- + + - -
+ + - + -
+ - - + -
+ ]
493 [ + - - + : 66114 : hhartree = Sr%hhartree(ib1:ib2,ib1:ib2,ik_ibz,:)
- - + + +
+ + + ]
494 : :
495 : : ! If non self-consistent erase all off-diagonal elements
496 [ + + ]: 607 : if (Sigp%gwcalctyp<20) then
497 [ + + ]: 1778 : do jb=ib1,ib2
498 [ + + ]: 13612 : do kb=ib1,ib2
499 [ + + ]: 11834 : if (jb == kb) CYCLE
500 [ + + ]: 27416 : hhartree(jb,kb,:) = czero
501 : : end do
502 : : end do
503 : : end if
504 : :
505 [ - + + - : 2428 : ABI_MALLOC(htotal, (ib1:ib2,ib1:ib2,Sr%nsppol*Sr%nsig_ab))
- + + - -
+ + - - +
+ - ]
506 [ + + ]: 1255 : do spin=1,Sr%nsppol*Sr%nsig_ab
507 [ + + ]: 6879 : do jb=ib1,ib2
508 [ + + ]: 64900 : do kb=ib1,ib2
509 : 64252 : htotal(kb,jb,spin) = hhartree(kb,jb,spin) + Sr%x_mat(kb,jb,ik_ibz,spin) + sigcme(ie0,kb,jb,spin)
510 : : end do
511 : : end do
512 : : end do
513 : :
514 : : ! Get the Hermitian part of htotal
515 : : ! In the noncollinear case A_{12}^{ab} = A_{21}^{ba}^* if A is Hermitian.
516 [ - + + - : 2428 : ABI_MALLOC(h_tmp1, (ib1:ib2,ib1:ib2))
- + + - +
- + - - +
- + ]
517 [ - + + - : 1821 : ABI_MALLOC(h_tmp2, (ib1:ib2,ib1:ib2))
- + + - -
+ - + ]
518 : :
519 [ + + ]: 607 : nsploop=Sr%nsppol; if (Sr%nsig_ab/=1) nsploop=2
520 [ + + ]: 1233 : do spin=1,nsploop
521 [ + - - + : 61424 : h_tmp1 = CONJG(htotal(:,:,spin))
- - + + +
+ ]
522 [ + - - + : 61424 : h_tmp2 = TRANSPOSE(h_tmp1)
- - + + +
+ ]
523 [ + - - + : 61424 : h_tmp1 = htotal(:,:,spin)
- - + + +
+ ]
524 [ + + + + ]: 61405 : htotal(:,:,spin)= half * (h_tmp1 + h_tmp2)
525 : : end do
526 : :
527 : : ! Print the different matrix elements of sigma if QPSC and prtvol>9
528 [ + + - + : 607 : if (Sigp%gwcalctyp >=20 .and. mod10 /= 1 .and. dtset%prtvol>9 .and. my_rank==master) then
- - ]
529 : : call print_sigma_melems(ikcalc,ib1,ib2,Sr%nsppol*Sr%nsig_ab,htotal,hhartree,&
530 [ # # # # : 0 : Sr%x_mat(ib1:ib2,ib1:ib2,ik_ibz,:),sigcme(ie0,:,:,:),Dtfil%filnam_ds(4))
# # # # #
# # # # #
# # # # #
# # # #
# ]
531 : : end if
532 : :
533 [ + + ]: 607 : if (Sr%nsig_ab==4) then
534 [ + - - + : 2062 : h_tmp1 = CONJG(htotal(:,:,4))
- - + + +
+ ]
535 [ + - - + : 2062 : h_tmp2 = TRANSPOSE(h_tmp1)
- - + + +
+ ]
536 [ + - - + : 2062 : h_tmp1 = htotal(:,:,3)
- - + + +
+ ]
537 [ + + + + ]: 2051 : htotal(:,:,3)= half * (h_tmp1 + h_tmp2)
538 : :
539 [ + - - + : 2062 : h_tmp1 = CONJG(htotal(:,:,3))
- - + + +
+ ]
540 [ + - - + : 2062 : h_tmp2 = TRANSPOSE(h_tmp1)
- - + + +
+ ]
541 [ + + + + ]: 2051 : htotal(:,:,4) = h_tmp2
542 : : end if
543 : :
544 : : ! Solve Herm(htotal)*U = E*U
545 : 607 : ld_matrix = ib2 - ib1 + 1
546 [ - + + - : 2428 : ABI_MALLOC(hdp, (ld_matrix, ld_matrix))
- + + - +
- + - - +
- + ]
547 [ - + + - : 1821 : ABI_MALLOC(eig, (ld_matrix))
+ - + - -
+ - + ]
548 : :
549 [ + + ]: 1222 : do spin=1,Sr%nsppol
550 [ + + ]: 615 : if (Sr%nsig_ab==1) then
551 [ + - - + : 57300 : hdp=htotal(ib1:ib2,ib1:ib2,spin)
- - + + +
+ ]
552 : : else
553 [ + - - + : 9646 : hdp = SUM(htotal(ib1:ib2,ib1:ib2,:), DIM=3)
- - + + +
+ + + ]
554 : : end if
555 : 615 : call xheev("Vectors","Upper", ld_matrix, hdp, eig)
556 : :
557 [ + + ]: 615 : if (Sr%needs_eigvec_qp) then
558 [ + + + + ]: 46183 : Sr%eigvec_qp(ib1:ib2,ib1:ib2,ik_ibz,spin)=hdp(:,:)
559 : : end if
560 [ + + ]: 6414 : Sr%en_qp_diago(ib1:ib2,ik_ibz,spin)=eig(:)
561 : : end do
562 : :
563 [ - + ]: 607 : ABI_FREE(hdp)
564 [ - + ]: 607 : ABI_FREE(eig)
565 [ - + ]: 607 : ABI_FREE(htotal)
566 [ - + ]: 607 : ABI_FREE(hhartree)
567 [ - + ]: 607 : ABI_FREE(h_tmp1)
568 [ - + ]: 607 : ABI_FREE(h_tmp2)
569 : : end if ! (mod10 /= 1)
570 : :
571 : 614 : call timab(490,2,tsec)
572 : :
573 : : DBG_EXIT("COLL")
574 : :
575 [ - + - + ]: 614 : end subroutine solve_dyson
576 : : !!***
577 : :
578 : : !----------------------------------------------------------------------
579 : :
580 : : !!****f* m_dyson_solver/print_sigma_melems
581 : : !! NAME
582 : : !! print_sigma_melems
583 : : !!
584 : : !! FUNCTION
585 : : !! This routine prints the Hermitian and the non-hermitian part of the matrix
586 : : !! elements of Sigma, as well as the individual contributions.
587 : : !! The first 14x14 are printed to screen, and the full matrices are printed
588 : : !! to files: sigma_melems_, sigma_nonH_melems_, sigma_Hart_melems_,
589 : : !! sigma_x_melems, and sigma_c_melems
590 : : !!
591 : : !! INPUTS
592 : : !! ikcalc : index of k-point
593 : : !! ib1,ib2 : starting and ending band indices
594 : : !! nsp : no. of spin elements
595 : : !! htotal : Hermitianised matrix elements of Sigma
596 : : !! hhartree : Hartree contribution to matrix elements
597 : : !! sigxme : Sigma_x contribution to matrix elements
598 : : !! sigcme : Sigma_c contribution to matrix elements
599 : : !! prefix : prefix for output files.
600 : : !!
601 : : !! OUTPUT
602 : : !!
603 : : !! SOURCE
604 : :
605 : 0 : subroutine print_sigma_melems(ikcalc, ib1, ib2, nsp, htotal, hhartree, sigxme, sigcme, prefix)
606 : :
607 : : ! Arguments ------------------------------------
608 : : !scalars
609 : : integer,intent(in) :: ikcalc,ib1,ib2,nsp
610 : : character(len=*),intent(in) :: prefix
611 : : !arrays
612 : : complex(dp),intent(in) :: htotal(ib1:ib2,ib1:ib2,nsp),hhartree(ib1:ib2,ib1:ib2,nsp)
613 : : complex(dp),intent(in) :: sigxme(ib1:ib2,ib1:ib2,nsp),sigcme(ib1:ib2,ib1:ib2,nsp)
614 : :
615 : : ! Local variables ------------------------------
616 : : integer,parameter :: MAX_NCOLS = 14
617 : : integer :: isp,mc,mr,jj,ii,temp_unit,ount
618 : : character(len=10) :: sidx
619 : : character(len=500) :: msg
620 : : character(len=100) :: fmth,fmt1,fmt2,fmthh,kpt_index,fmtfile
621 : : character(len=fnlen) :: filename
622 : : ! *************************************************************************
623 : :
624 [ # # ]: 0 : if (nsp==3.or.nsp>4) then
625 : 0 : ABI_ERROR('nsp has wrong value in print_sigma_melems')
626 : : end if
627 : :
628 : 0 : ount = std_out
629 : :
630 [ # # ]: 0 : mc = ib2-ib1+1; if (mc>MAX_NCOLS) mc = MAX_NCOLS
631 : 0 : mr = mc
632 : :
633 : 0 : write(fmthh,*)'(2(a),2(I2,a))'
634 : 0 : write(fmth,*)'(7x,',mc,'(i2,8x))'
635 : 0 : write(fmt1,*)'(3x,i2,',mc,'f10.5)'
636 : 0 : write(fmt2,*)'(5x ,',mc,'f10.5,a)'
637 : :
638 : : ! First print to screen
639 [ # # ]: 0 : do isp=1,nsp
640 : 0 : write(msg,'(a)') ''
641 : 0 : call wrtout(ount,msg)
642 : 0 : write(msg,fmthh) ch10,' Hermitianised matrix elements of Sigma (spin ',isp,' of ',nsp,'):'
643 : 0 : call wrtout(ount,msg)
644 [ # # # # ]: 0 : write(msg,fmth)(jj,jj=1,mc)
645 : 0 : call wrtout(ount,msg) !header
646 [ # # ]: 0 : do ii=ib1,ib1+mr-1
647 [ # # ]: 0 : write(msg,fmt1)ii-ib1+1,DBLE(htotal(ii,ib1:(ib1+mc-1),isp))
648 : 0 : call wrtout(ount,msg) !real part
649 [ # # ]: 0 : write(msg,fmt2) AIMAG(htotal(ii,ib1:(ib1+mc-1),isp)),ch10
650 : 0 : call wrtout(ount,msg) !imag part
651 : : end do
652 : : end do !nsp
653 : :
654 : 0 : write(msg,'(a,i2,a)')" Max. ",MAX_NCOLS," elements printed. Full matrix output in _HTOTAL files"
655 : 0 : call wrtout(ount,msg)
656 : :
657 [ # # ]: 0 : do isp=1,nsp
658 : 0 : write(msg,fmthh) ch10,' H_Hartree matrix elements (spin ',isp,' of ',nsp,'):'
659 : 0 : call wrtout(ount,msg)
660 [ # # # # ]: 0 : write(msg,fmth)(jj,jj=1,mc)
661 : 0 : call wrtout(ount,msg) !header
662 [ # # ]: 0 : do ii=ib1,ib1+mr-1
663 [ # # ]: 0 : write(msg,fmt1)ii-ib1+1,DBLE(hhartree(ii,ib1:(ib1+mc-1),isp))
664 : 0 : call wrtout(ount,msg) !real part
665 [ # # ]: 0 : write(msg,fmt2) AIMAG(hhartree(ii,ib1:(ib1+mc-1),isp)),ch10
666 : 0 : call wrtout(ount,msg) !imag part
667 : : end do
668 : : end do !nsp
669 : :
670 : 0 : write(msg,'(a,i2,a)')" Max. ",MAX_NCOLS," elements printed. Full matrix output in _HHARTREE files"
671 : 0 : call wrtout(ount,msg)
672 : :
673 [ # # ]: 0 : do isp=1,nsp
674 : 0 : write(msg,fmthh) ch10,' Sigma_x matrix elements (spin ',isp,' of ',nsp,'):'
675 : 0 : call wrtout(ount,msg)
676 [ # # # # ]: 0 : write(msg,fmth)(jj,jj=1,mc)
677 : 0 : call wrtout(ount,msg) !header
678 [ # # ]: 0 : do ii=ib1,ib1+mr-1
679 [ # # ]: 0 : write(msg,fmt1)ii-ib1+1,DBLE(sigxme(ii,ib1:(ib1+mc-1),isp))
680 : 0 : call wrtout(ount,msg) !real part
681 [ # # ]: 0 : write(msg,fmt2) AIMAG(sigxme(ii,ib1:(ib1+mc-1),isp)),ch10
682 : 0 : call wrtout(ount,msg) !imag part
683 : : end do
684 : : end do !nsp
685 : :
686 : 0 : write(msg,'(a,i2,a)')" Max. ",MAX_NCOLS," elements printed. Full matrix output _SIGX files"
687 : 0 : call wrtout(ount,msg)
688 : :
689 [ # # ]: 0 : do isp=1,nsp
690 : 0 : write(msg,fmthh) ch10,' Sigma_c matrix elements (spin ',isp,' of ',nsp,'):'
691 : 0 : call wrtout(ount,msg)
692 [ # # # # ]: 0 : write(msg,fmth)(jj,jj=1,mc)
693 : 0 : call wrtout(ount,msg) !header
694 [ # # ]: 0 : do ii=ib1,ib1+mr-1
695 [ # # ]: 0 : write(msg,fmt1)ii-ib1+1,DBLE(sigcme(ii,ib1:(ib1+mc-1),isp))
696 : 0 : call wrtout(ount,msg) !real part
697 [ # # ]: 0 : write(msg,fmt2) AIMAG(sigcme(ii,ib1:(ib1+mc-1),isp)),ch10
698 : 0 : call wrtout(ount,msg) !imag part
699 : : end do
700 : : end do !nsp
701 : :
702 : 0 : write(msg,'(a,i2,a)')" Max ",MAX_NCOLS," elements printed. Full matrix output _SIGC files"
703 : 0 : call wrtout(ount,msg)
704 : :
705 : : ! Then print to file
706 : : ! Format is: row, column, value; with a blank space for each full
707 : : ! set of columns for easy plotting with the gnuplot splot command
708 : 0 : write(fmtfile,*)'(3X,I6,2X,I6,',nsp,'(2(ES28.16E3,3x)))'
709 : :
710 : 0 : call int2char10(ikcalc,sidx)
711 : 0 : kpt_index = "_KPT"//TRIM(sidx)
712 : :
713 [ # # ]: 0 : filename = TRIM(prefix)//'_HTOTAL'//TRIM(kpt_index)
714 : :
715 [ # # ]: 0 : if (open_file(filename,msg,newunit=temp_unit,form="formatted",status="replace",action="write") /= 0) then
716 : 0 : ABI_ERROR(msg)
717 : : end if
718 : :
719 : 0 : msg = '# row col. Re(htotal(r,c)) Im(htotal(r,c)) for spin11 ... spin22 ... spin12 ... spin13'
720 : 0 : call wrtout(temp_unit,msg)
721 [ # # ]: 0 : do ii=ib1,ib2
722 [ # # ]: 0 : do jj=ib1,ib2
723 : 0 : write(msg,fmtfile) ii,jj,(htotal(jj,ii,isp),isp=1,nsp)
724 : 0 : call wrtout(temp_unit,msg)
725 : : end do
726 : 0 : call wrtout(temp_unit,"")
727 : : end do
728 : 0 : close(temp_unit)
729 : :
730 [ # # ]: 0 : filename = TRIM(prefix)//'_HHARTREE'//TRIM(kpt_index)
731 [ # # ]: 0 : if (open_file(filename,msg,newunit=temp_unit,form="formatted",status="replace",action="write") /= 0) then
732 : 0 : ABI_ERROR(msg)
733 : : end if
734 : :
735 : 0 : msg = '# row col. Re(hhartree(r,c)) Im(hhartree(r,c) for spin11 ... spin22 ... spin12 ... spin13'
736 : 0 : call wrtout(temp_unit,msg)
737 [ # # ]: 0 : do ii=ib1,ib2
738 [ # # ]: 0 : do jj=ib1,ib2
739 : 0 : write(msg,fmtfile) ii,jj,(hhartree(jj,ii,isp),isp=1,nsp)
740 : 0 : call wrtout(temp_unit,msg)
741 : : end do
742 : 0 : call wrtout(temp_unit,"")
743 : : end do
744 : 0 : close(temp_unit)
745 : :
746 [ # # ]: 0 : filename = TRIM(prefix)//'_SIGX'//TRIM(kpt_index)
747 [ # # ]: 0 : if (open_file(filename,msg,newunit=temp_unit,form="formatted",status="replace",action="write") /= 0) then
748 : 0 : ABI_ERROR(msg)
749 : : end if
750 : :
751 : 0 : write(msg,'(a)')'# row col. Re(Sigx(r,c)) Im(Sigx(r,c) for spin11 ... spin22 ... spin12 ... spin13'
752 : 0 : call wrtout(temp_unit,msg)
753 [ # # ]: 0 : do ii=ib1,ib2
754 [ # # ]: 0 : do jj=ib1,ib2
755 : 0 : write(msg,fmtfile) ii,jj,(sigxme(jj,ii,isp),isp=1,nsp)
756 : 0 : call wrtout(temp_unit,msg)
757 : : end do
758 : 0 : call wrtout(temp_unit,"")
759 : : end do
760 : 0 : close(temp_unit)
761 : :
762 [ # # ]: 0 : filename = TRIM(prefix)//'_SIGC'//TRIM(kpt_index)
763 [ # # ]: 0 : if (open_file(filename,msg,newunit=temp_unit,form="formatted",status="replace",action="write") /= 0) then
764 : 0 : ABI_ERROR(msg)
765 : : end if
766 : :
767 : 0 : write(msg,'(a)')'# row col. Re(Sigc(r,c)) Im(Sigc(r,c) for spin11 ... spin22 ... spin12 ... spin21'
768 : 0 : call wrtout(temp_unit,msg)
769 [ # # ]: 0 : do ii=ib1,ib2
770 [ # # ]: 0 : do jj=ib1,ib2
771 : 0 : write(msg,fmtfile) ii,jj,(sigcme(jj,ii,isp),isp=1,nsp)
772 : 0 : call wrtout(temp_unit,msg)
773 : : end do
774 : 0 : call wrtout(temp_unit,"")
775 : : end do
776 : :
777 : 0 : close(temp_unit)
778 : :
779 : 0 : end subroutine print_sigma_melems
780 : : !!***
781 : :
782 : : !----------------------------------------------------------------------
783 : :
784 : : !!****f* m_dyson_solver/sigma_pade_init
785 : : !! NAME
786 : : !! sigma_pade_init
787 : : !!
788 : : !! FUNCTION
789 : : !! Initialize the Pade' from the `npts` values of Sigma_c(iw) given on the mesh `zmesh`.
790 : : !!
791 : : !! SOURCE
792 : :
793 [ + - ]: 106 : subroutine sigma_pade_init(self, npts, zmesh, sigc_cvals, alphac_pm, betar_pm, zcut_pm)
794 : :
795 : : !Arguments ------------------------------------
796 : : class(sigma_pade_t),intent(out) :: self
797 : : integer,intent(in) :: npts
798 : : complex(dp),target,intent(in) :: zmesh(npts), sigc_cvals(npts), alphac_pm(2)
799 : : real(dp),intent(in) :: betar_pm(2), zcut_pm(2)
800 : : ! *************************************************************************
801 : :
802 : 106 : self%npts = npts
803 [ - + + - : 318 : ABI_MALLOC(self%zmesh, (npts))
+ - + - -
+ - + -
+ ]
804 [ - + + - : 212 : ABI_MALLOC(self%sigc_cvals, (npts,1))
- + - + -
+ ]
805 [ - + + + ]: 1432 : self%zmesh = zmesh
806 [ + + ]: 1326 : self%sigc_cvals(:,1) = sigc_cvals
807 : 106 : self%nsig_ab = 1
808 : :
809 [ + + ]: 318 : self%alphac_pm = alphac_pm
810 [ + + ]: 318 : self%betar_pm = betar_pm
811 [ + + ]: 318 : self%zcut_pm = zcut_pm
812 : 106 : self%do_sigma_fit = .False.
813 : :
814 : 106 : end subroutine sigma_pade_init
815 : : !!***
816 : : !----------------------------------------------------------------------
817 : :
818 : : !!****f* m_dyson_solver/sigma_pade_init_spinor
819 : : !! NAME
820 : : !! sigma_pade_init_spinor
821 : : !!
822 : : !! FUNCTION
823 : : !! Initialize the Pade' from the `npts` values of Sigma_c(iw) given on the mesh `zmesh`.
824 : : !!
825 : : !! SOURCE
826 : :
827 [ # # ]: 0 : subroutine sigma_pade_init_spinor(self, npts, zmesh, nsig_ab, sigc_cvals, alphac_pm, betar_pm, zcut_pm)
828 : :
829 : : !Arguments ------------------------------------
830 : : class(sigma_pade_t),intent(out) :: self
831 : : integer,intent(in) :: npts, nsig_ab
832 : : complex(dp),target,intent(in) :: zmesh(npts), sigc_cvals(npts,nsig_ab), alphac_pm(2)
833 : : real(dp),intent(in) :: betar_pm(2), zcut_pm(2)
834 : : ! *************************************************************************
835 : :
836 : 0 : self%npts = npts
837 [ # # # # : 0 : ABI_MALLOC(self%zmesh, (npts))
# # # # #
# # # #
# ]
838 [ # # # # : 0 : ABI_MALLOC(self%sigc_cvals, (npts,nsig_ab))
# # # # #
# # # # #
# # # # ]
839 [ # # # # ]: 0 : self%zmesh = zmesh
840 [ # # # # : 0 : self%sigc_cvals = sigc_cvals
# # ]
841 : 0 : self%nsig_ab = nsig_ab
842 : :
843 [ # # ]: 0 : self%alphac_pm = alphac_pm
844 [ # # ]: 0 : self%betar_pm = betar_pm
845 [ # # ]: 0 : self%zcut_pm = zcut_pm
846 : 0 : self%do_sigma_fit = .False.
847 : :
848 : 0 : end subroutine sigma_pade_init_spinor
849 : : !!***
850 : :
851 : :
852 : 106 : subroutine sigma_pade_free(self)
853 : : !Arguments ------------------------------------
854 : : class(sigma_pade_t),intent(inout) :: self
855 : : ! *************************************************************************
856 : :
857 [ + - ]: 106 : ABI_SFREE(self%zmesh)
858 [ + - ]: 106 : ABI_SFREE(self%sigc_cvals)
859 : :
860 : 106 : end subroutine sigma_pade_free
861 : : !!***
862 : :
863 : : !----------------------------------------------------------------------
864 : :
865 : : !!****f* m_dyson_solver/sigma_pade_eval
866 : : !! NAME
867 : : !! sigma_pade_eval
868 : : !!
869 : : !! FUNCTION
870 : : !! Evaluate the Pade' at the complex point `zz`.
871 : : !! Return result in `val` and, optionally, the derivative at zz in `dvdz`
872 : : !!
873 : : !! SOURCE
874 : :
875 : 215 : subroutine sigma_pade_eval(self, zz, val, &
876 : : dvdz) ! optional
877 : :
878 : : !Arguments ------------------------------------
879 : : class(sigma_pade_t),intent(in) :: self
880 : : complex(dp),intent(in) :: zz
881 : : complex(dp),intent(out) :: val
882 : : complex(dp),optional,intent(out) :: dvdz
883 : : integer :: iab
884 : : ! *************************************************************************
885 : :
886 : : ! if zz in 2 or 3 quadrant, avoid branch cut in the complex plane using Sigma(-iw) = Sigma(iw)*.
887 : 215 : val = czero
888 [ + - ]: 215 : if (present(dvdz)) dvdz = czero
889 : :
890 [ + + ]: 430 : do iab = 1, self%nsig_ab
891 [ + + ]: 430 : if (real(zz) > zero) then
892 : 104 : val = val + pade(self%npts, self%zmesh, self%sigc_cvals(:,iab), zz)
893 [ + - ]: 104 : if (present(dvdz)) then
894 : 104 : dvdz = dvdz + dpade(self%npts, self%zmesh, self%sigc_cvals(:,iab), zz)
895 : : end if
896 : : else
897 [ + - + + : 2651 : val = val + pade(self%npts, -self%zmesh, conjg(self%sigc_cvals(:,iab)), zz)
+ - + + -
+ ]
898 [ + - ]: 111 : if (present(dvdz)) then
899 [ + - + + : 2651 : dvdz = dvdz + dpade(self%npts, -self%zmesh, conjg(self%sigc_cvals(:,iab)), zz)
+ - + + -
+ ]
900 : : end if
901 : : end if
902 : : end do
903 : :
904 [ - + ]: 215 : if (self%do_sigma_fit) then
905 : : ! Add the analytic transform of the model removed in imaginary time:
906 : : !
907 : : ! Sigma(+tau) = alpha_1 e^{-beta_1 tau} --> alpha_1 / (beta_1 - z)
908 : : ! Sigma(-tau) = alpha_2 e^{-beta_2 tau} --> alpha_2 / (beta_2 + z)
909 : : !
910 : : ! with tau > 0, consistent with the cos/sin transforms used for the residual.
911 : : val = val + self%alphac_pm(1) / (self%betar_pm(1) - zz) &
912 : 0 : + self%alphac_pm(2) / (self%betar_pm(2) + zz)
913 : :
914 [ # # ]: 0 : if (present(dvdz)) then
915 : : dvdz = dvdz + self%alphac_pm(1) / ((self%betar_pm(1) - zz) ** 2) &
916 : 0 : - self%alphac_pm(2) / ((self%betar_pm(2) + zz) ** 2)
917 : : end if
918 : : end if
919 : :
920 : 215 : end subroutine sigma_pade_eval
921 : : !!***
922 : :
923 : : !----------------------------------------------------------------------
924 : :
925 : : !!****f* m_dyson_solver/sigma_pade_qp_solve
926 : : !! NAME
927 : : !! sigma_pade_qp_solve
928 : : !!
929 : : !! FUNCTION
930 : : !! Use the Pade' approximant and Newton-Rapson method to solve the QP equation
931 : : !! in the complex plane starting from the initial guess `z_guess`.
932 : : !!
933 : : !! INPUTS
934 : : !! e0: KS energy
935 : : !! v_meanf: matrix element of the mean-field Hamiltonian
936 : : !! sigx: matrix element of the exchange self-energy.
937 : : !! z_guess: Initial guess for the QP energy
938 : : !!
939 : : !! OUTPUT
940 : : !! zsc: root.
941 : : !! sigc: Sigma_c(zsc)
942 : : !! msg: Error message if ierr /= 0.
943 : : !! ierr: Exit status.
944 : : !!
945 : : !! SOURCE
946 : :
947 : 53 : subroutine sigma_pade_qp_solve(self, e0, v_meanf, sigx, z_guess, zsc, sigc_zsc, msg, ierr)
948 : :
949 : : !Arguments ------------------------------------
950 : : class(sigma_pade_t),intent(in) :: self
951 : : real(dp),intent(in) :: e0, v_meanf, sigx
952 : : complex(dp),intent(in) :: z_guess
953 : : complex(dp),intent(out) :: zsc, sigc_zsc
954 : : integer,intent(out) :: ierr
955 : :
956 : : !Local variables-------------------------------
957 : : !scalars
958 : : integer :: iter
959 : : logical :: converged
960 : : complex(dp) :: ctdpc, dct, dsigc
961 : : character(len=500) :: msg
962 : : ! *************************************************************************
963 : :
964 : : ! Use Newton-Rapson to find the root of:
965 : : ! f(z) = e0 - zz + Sigma_xc(z) - v_meanf
966 : : ! f'(z) = -1 + Sigma_c'(z)
967 : :
968 : 53 : iter = 0; converged = .FALSE.; ctdpc = cone
969 : 53 : zsc = z_guess
970 [ - + - - ]: 162 : do while (abs(ctdpc) > NR_ABS_ROOT_ERR .or. iter < NR_MAX_NITER)
971 : 162 : iter = iter + 1
972 : :
973 : 162 : call self%eval(zsc, sigc_zsc, dvdz=dsigc)
974 : 162 : ctdpc = e0 - v_meanf + sigx + sigc_zsc - zsc
975 : :
976 [ + + ]: 162 : if (abs(ctdpc) < NR_ABS_ROOT_ERR) then
977 : : converged=.TRUE.; EXIT
978 : : end if
979 : 109 : dct = dsigc - one
980 : 109 : zsc = newrap_step(zsc, ctdpc, dct)
981 : : end do
982 : :
983 : 53 : ierr = 0; msg = ""
984 [ - + ]: 53 : if (.not. converged) then
985 : : write(msg,'(a,i0,3a,f8.4,a,f8.4)')&
986 : 0 : 'Newton-Raphson method did not converge after: ', NR_MAX_NITER,' iterations.',ch10,&
987 : 0 : 'Absolute error: ', abs(ctdpc), ' > ', NR_ABS_ROOT_ERR
988 : 0 : ierr = 1
989 : : end if
990 : :
991 : 53 : end subroutine sigma_pade_qp_solve
992 : : !!***
993 : :
994 : : !----------------------------------------------------------------------
995 : :
996 [ - + - - : 212 : end module m_dyson_solver
+ + - + +
- - + -
+ ]
997 : : !!***
|
