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 609 : 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 609 : type(sigma_pade_t) :: spade
169 : !arrays
170 : real(dp) :: kbz_gw(3),tsec(2), betar_pm(2), zcut_pm(2)
171 609 : real(dp),allocatable :: e0pde(:),eig(:),scme(:)
172 : complex(dp) :: alphac_pm(2)
173 609 : complex(dp),allocatable :: hdp(:,:),tmpcdp(:),hhartree(:,:,:),htotal(:,:,:),h_tmp1(:,:),h_tmp2(:,:)
174 : ! *************************************************************************
175 :
176 : DBG_ENTER("COLL")
177 :
178 609 : units = [std_out, ab_out]
179 :
180 609 : call timab(490,1,tsec) ! csigme(Dyson)
181 :
182 609 : my_rank = xmpi_comm_rank(comm)
183 :
184 609 : mod10 = MOD(Sigp%gwcalctyp,10)
185 :
186 609 : ltest=(nomega_sigc==Sr%nomega_r+Sr%nomega4sd)
187 609 : if (mod10==1) ltest=(nomega_sigc==Sr%nomega_i)
188 609 : 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 1226 : ib1 = MINVAL(Sigp%minbnd(ikcalc,:))
195 1226 : ib2 = MAXVAL(Sigp%maxbnd(ikcalc,:))
196 :
197 : ! Find the index of the k-point for sigma in the IBZ array.
198 609 : ikbz_gw = Sigp%kptgw2bz(ikcalc)
199 609 : call kmesh%get_BZ_item(ikbz_gw, kbz_gw, ik_ibz, isym, itim, phase)
200 :
201 609 : sigc = czero; dsigc = czero
202 :
203 : ! ===========================================================
204 : ! ==== Solve the Dyson Equation and store results in Sr% ====
205 : ! ===========================================================
206 :
207 609 : if (mod10 /= 1) then
208 : ! ===============================
209 : ! ==== Perturbative approach ====
210 : ! ===============================
211 :
212 : ! Index of the KS or QP energy in sigme_tmp
213 602 : ie0 = sr%nomega_r + Sr%nomega4sd/2+1
214 :
215 1212 : do spin=1,Sr%nsppol
216 6390 : 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 10788 : do iab=1,Sr%nsig_ab
220 5610 : is_idx = spin; if (Sr%nsig_ab>1) is_idx=iab
221 :
222 5610 : Sr%sigcmee0(jb,ik_ibz,is_idx) = sigcme(ie0,jb,jb,is_idx)
223 :
224 16830 : ABI_MALLOC(scme, (Sr%nomega4sd))
225 11220 : ABI_MALLOC(e0pde, (Sr%nomega4sd))
226 25300 : e0pde(:) = Sr%omega4sd(jb,ik_ibz,:,spin)
227 25300 : scme(:) = REAL(sigcme(Sr%nomega_r+1:Sr%nomega_r+Sr%nomega4sd,jb,jb,is_idx))
228 :
229 5610 : if (Sr%nomega4sd==1) then
230 3846 : smrt = zero; alpha = zero
231 : else
232 1764 : smrt = linfit(Sr%nomega4sd,e0pde(:),scme(:),alpha,beta)
233 : end if
234 :
235 5610 : 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 5610 : ABI_FREE(scme)
250 5610 : 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 5610 : Sr%dsigmee0(jb,ik_ibz,is_idx)=CMPLX(alpha,zero)
258 :
259 5610 : if (Sr%nsig_ab==1) then
260 5034 : Sr%ze0(jb,ik_ibz,spin)= one / (one-Sr%dsigmee0(jb,ik_ibz,spin))
261 :
262 5034 : 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 5030 : Sr%hhartree(jb,jb,ik_ibz,spin))
266 :
267 5030 : 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 5030 : 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 8600 : do io=1,Sr%nomega_r
293 2990 : Sr%sigcme (jb,ik_ibz,io,is_idx)= sigcme(io,jb,jb,is_idx)
294 8600 : 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 30478 : do io=1,Sr%nomega4sd
297 19690 : Sr%sigcme4sd (jb,ik_ibz,io,is_idx)= sigcme(Sr%nomega_r+io,jb,jb,is_idx)
298 25300 : 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 5788 : 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 : vxc_val = ks_me%vxcval(jb, jb, ik_ibz, spin)
371 : if (Sr%nsig_ab > 1) vxc_val = SUM(ks_me%vxcval(jb, jb, ik_ibz, :))
372 : vu = zero; if (dtset%usepawu /= 0) vu = ks_me%vu(jb, jb, ik_ibz, spin)
373 : v_meanf = vxc_val + vu
374 :
375 : ! qp_ene = e0 + z_e0 * (sigc_e0__ + sigx - v_meanf)
376 : if (Sr%nsig_ab == 1) then
377 : Sr%egw(jb,ik_ibz,spin) = Sr%e0(jb,ik_ibz,spin) + Sr%ze0(jb,ik_ibz,spin) * &
378 : (Sr%sigcmee0(jb,ik_ibz,spin) + Sr%sigxme(jb,ik_ibz,spin) - v_meanf)
379 :
380 : Sr%degw(jb,ik_ibz,spin) = Sr%egw(jb,ik_ibz,spin) - Sr%e0(jb,ik_ibz,spin)
381 :
382 : ! Estimate Sigma at the QP-energy: Sigma(E_qp)=Sigma(E0)+(E_qp-E0)*dSigma/dE
383 : Sr%sigmee(jb,ik_ibz,spin) = &
384 : Sr%sigxme(jb,ik_ibz,spin)+Sr%sigcmee0(jb,ik_ibz,spin)+Sr%degw(jb,ik_ibz,spin)*Sr%dsigmee0(jb,ik_ibz,spin)
385 : else
386 : Sr%egw(jb,ik_ibz,1) = Sr%e0(jb,ik_ibz,1) + Sr%ze0(jb,ik_ibz,1) * &
387 : (SUM(Sr%sigcmee0(jb,ik_ibz,:)+Sr%sigxme(jb,ik_ibz,:))-v_meanf)
388 :
389 : Sr%degw(jb,ik_ibz,1) = Sr%egw(jb,ik_ibz,1) - Sr%e0(jb,ik_ibz,1)
390 :
391 : Sr%sigmee(jb,ik_ibz,1) = &
392 : SUM(Sr%sigxme(jb,ik_ibz,:)+Sr%sigcmee0(jb,ik_ibz,:))+Sr%degw(jb,ik_ibz,1)*SUM(Sr%dsigmee0(jb,ik_ibz,:))
393 : end if
394 :
395 : #else
396 : ! MG FIXME: Here we are solving the non-linear QP equation using the Pade' continuation + root finding
397 : ! but this is very misleading because in the output file we are still reporting the Z factor
398 : ! and there's no mention that the QP energies have been obtained from the non-linear equation!!
399 : ! One should change the format used to print the results or at least warn the user!
400 :
401 53 : zz = CMPLX(Sr%e0(jb,ik_ibz,spin), zero)
402 :
403 53 : if (Sigp%mbpt_sciss>0.1d-4) then
404 : ! e0 is replaced by qp_ene which contains the updated energy eigenvalue.
405 0 : zz = CMPLX(qp_ene(jb,ik_ibz,spin), zero)
406 : end if
407 :
408 : ! Solve the QP equation with Newton-Rapson starting from e0
409 : ! Find root of E^0-V_xc-V_U+Sig_x+Sig_c(z)-z, i.e E^qp.
410 53 : alphac_pm = zero; betar_pm = zero; zcut_pm = zero
411 53 : if (Sr%nsig_ab > 1) then
412 0 : call spade%init(sr%nomega_i, sr%omega_i, Sr%nsig_ab, sigcme(:,jb,jb,:), alphac_pm, betar_pm, zcut_pm)
413 : else
414 53 : call spade%init(sr%nomega_i, sr%omega_i, sigcme(:,jb,jb,spin), alphac_pm, betar_pm, zcut_pm)
415 : end if
416 :
417 : ! Note vxc[n_val] instead of vxc[n_val + n_nlcc] with the model core charge.
418 53 : vxc_val = ks_me%vxcval(jb, jb, ik_ibz, spin)
419 53 : if (Sr%nsig_ab > 1) vxc_val = SUM(ks_me%vxcval(jb, jb, ik_ibz, :))
420 53 : vu = zero; if (dtset%usepawu /= 0) vu = ks_me%vu(jb, jb, ik_ibz, spin)
421 53 : v_meanf = vxc_val + vu
422 53 : if (Sr%nsig_ab == 1) then
423 53 : sigx = Sr%sigxme(jb,ik_ibz,spin)
424 : else
425 0 : sigx = SUM(Sr%sigxme(jb,ik_ibz,:))
426 : end if
427 :
428 53 : call spade%qp_solve(sr%e0(jb,ik_ibz,spin), v_meanf, sigx, zz, zsc, sigc_zsc, msg, ierr)
429 53 : call spade%free()
430 : !qpe_pade_kcalc(ibc, ikcalc, spin) = zsc
431 : !qp_solver_ierr(ibc, ikcalc, spin) = ierr
432 53 : if (ierr /= 0) then
433 0 : ABI_WARNING(msg)
434 : end if
435 :
436 : ! Store the final result (self-consistent result for zz and Sigma_c(zz_scf)
437 53 : Sr%egw(jb,ik_ibz,spin) = zsc
438 53 : Sr%degw(jb,ik_ibz,spin) = Sr%egw(jb,ik_ibz,spin) - Sr%e0(jb,ik_ibz,spin)
439 53 : Sr%sigmee(jb,ik_ibz,spin) = Sr%sigxme(jb,ik_ibz,spin) + sigc_zsc
440 : #endif
441 :
442 : ! Spectra of Sigma, remember that Sr%nomega_r does not contain the frequencies
443 : ! used to evaluate the derivative each frequency is obtained using the pade_expression
444 : ! In sigma indeed we have:
445 : ! nomega_sigc=Sr%nomega_r+Sr%nomega4sd; if (mod10==SIG_GW_AC) nomega_sigc=Sr%nomega_i
446 803 : do io=1,Sr%nomega_r
447 750 : zz=Sr%omega_r(io)
448 750 : if (REAL(zz) > zero) then
449 7875 : tmpcdp(:) = SUM(sigcme(:,jb,jb,:), DIM=2)
450 375 : Sr%sigcme(jb,ik_ibz,io,spin) = pade(Sr%nomega_i, Sr%omega_i, tmpcdp, zz)
451 : else
452 7875 : tmpcdp(:) = CONJG(SUM(sigcme(:,jb,jb,:), DIM=2))
453 4125 : Sr%sigcme(jb,ik_ibz,io,spin) = pade(Sr%nomega_i, CONJG(Sr%omega_i), tmpcdp, zz)
454 : end if
455 803 : Sr%sigxcme(jb,ik_ibz,io,spin) = Sr%sigxme(jb,ik_ibz,spin) + Sr%sigcme(jb,ik_ibz,io,spin)
456 : end do
457 :
458 : ! Save sigma values along the imaginary axis
459 106 : do iab=1,Sr%nsig_ab
460 53 : is_idx=spin; if (Sr%nsig_ab > 1) is_idx = iab
461 716 : do io=1,Sr%nomega_i
462 610 : Sr%sigcmesi (jb,ik_ibz,io,is_idx) = sigcme(io,jb,jb,is_idx)
463 663 : Sr%sigxcmesi(jb,ik_ibz,io,is_idx) = Sr%sigxme(jb,ik_ibz,is_idx) + Sr%sigcmesi(jb,ik_ibz,io,is_idx)
464 : end do
465 : end do
466 :
467 60 : ABI_FREE(tmpcdp)
468 : end do !jb
469 : end do !is
470 :
471 : end if ! Analytic continuation.
472 :
473 : ! === Diagonalize the QP Hamiltonian (forced to be Hermitian) ===
474 : ! Calculate Sr%en_qp_diago and Sr%eigvec_qp to be written in the QPS file.
475 : ! TODO in case of AC results are wrong.
476 :
477 609 : if (mod10 /= 1) then
478 3010 : ABI_MALLOC(hhartree, (ib1:ib2,ib1:ib2,Sr%nsppol*Sr%nsig_ab))
479 66041 : hhartree = Sr%hhartree(ib1:ib2,ib1:ib2,ik_ibz,:)
480 :
481 : ! If non self-consistent erase all off-diagonal elements
482 602 : if (Sigp%gwcalctyp<20) then
483 1759 : do jb=ib1,ib2
484 13549 : do kb=ib1,ib2
485 11790 : if (jb == kb) CYCLE
486 27342 : hhartree(jb,kb,:) = czero
487 : end do
488 : end do
489 : end if
490 :
491 2408 : ABI_MALLOC(htotal, (ib1:ib2,ib1:ib2,Sr%nsppol*Sr%nsig_ab))
492 1245 : do spin=1,Sr%nsppol*Sr%nsig_ab
493 6855 : do jb=ib1,ib2
494 64837 : do kb=ib1,ib2
495 64194 : htotal(kb,jb,spin) = hhartree(kb,jb,spin) + Sr%x_mat(kb,jb,ik_ibz,spin) + sigcme(ie0,kb,jb,spin)
496 : end do
497 : end do
498 : end do
499 :
500 : ! Get the Hermitian part of htotal
501 : ! In the noncollinear case A_{12}^{ab} = A_{21}^{ba}^* if A is Hermitian.
502 2408 : ABI_MALLOC(h_tmp1, (ib1:ib2,ib1:ib2))
503 1806 : ABI_MALLOC(h_tmp2, (ib1:ib2,ib1:ib2))
504 :
505 602 : nsploop=Sr%nsppol; if (Sr%nsig_ab/=1) nsploop=2
506 1223 : do spin=1,nsploop
507 61356 : h_tmp1 = CONJG(htotal(:,:,spin))
508 61356 : h_tmp2 = TRANSPOSE(h_tmp1)
509 61356 : h_tmp1 = htotal(:,:,spin)
510 61337 : htotal(:,:,spin)= half * (h_tmp1 + h_tmp2)
511 : end do
512 :
513 : ! Print the different matrix elements of sigma if QPSC and prtvol>9
514 602 : if (Sigp%gwcalctyp >=20 .and. mod10 /= 1 .and. dtset%prtvol>9 .and. my_rank==master) then
515 : call print_sigma_melems(ikcalc,ib1,ib2,Sr%nsppol*Sr%nsig_ab,htotal,hhartree,&
516 0 : Sr%x_mat(ib1:ib2,ib1:ib2,ik_ibz,:),sigcme(ie0,:,:,:),Dtfil%filnam_ds(4))
517 : end if
518 :
519 602 : if (Sr%nsig_ab==4) then
520 2062 : h_tmp1 = CONJG(htotal(:,:,4))
521 2062 : h_tmp2 = TRANSPOSE(h_tmp1)
522 2062 : h_tmp1 = htotal(:,:,3)
523 2051 : htotal(:,:,3)= half * (h_tmp1 + h_tmp2)
524 :
525 2062 : h_tmp1 = CONJG(htotal(:,:,3))
526 2062 : h_tmp2 = TRANSPOSE(h_tmp1)
527 2051 : htotal(:,:,4) = h_tmp2
528 : end if
529 :
530 : ! Solve Herm(htotal)*U = E*U
531 602 : ld_matrix = ib2 - ib1 + 1
532 2408 : ABI_MALLOC(hdp, (ld_matrix, ld_matrix))
533 1806 : ABI_MALLOC(eig, (ld_matrix))
534 :
535 1212 : do spin=1,Sr%nsppol
536 610 : if (Sr%nsig_ab==1) then
537 57232 : hdp=htotal(ib1:ib2,ib1:ib2,spin)
538 : else
539 9646 : hdp = SUM(htotal(ib1:ib2,ib1:ib2,:), DIM=3)
540 : end if
541 610 : call xheev("Vectors","Upper", ld_matrix, hdp, eig)
542 :
543 610 : if (Sr%needs_eigvec_qp) then
544 46183 : Sr%eigvec_qp(ib1:ib2,ib1:ib2,ik_ibz,spin)=hdp(:,:)
545 : end if
546 6390 : Sr%en_qp_diago(ib1:ib2,ik_ibz,spin)=eig(:)
547 : end do
548 :
549 602 : ABI_FREE(hdp)
550 602 : ABI_FREE(eig)
551 602 : ABI_FREE(htotal)
552 602 : ABI_FREE(hhartree)
553 602 : ABI_FREE(h_tmp1)
554 602 : ABI_FREE(h_tmp2)
555 : end if ! (mod10 /= 1)
556 :
557 609 : call timab(490,2,tsec)
558 :
559 : DBG_EXIT("COLL")
560 :
561 609 : end subroutine solve_dyson
562 : !!***
563 :
564 : !----------------------------------------------------------------------
565 :
566 : !!****f* m_dyson_solver/print_sigma_melems
567 : !! NAME
568 : !! print_sigma_melems
569 : !!
570 : !! FUNCTION
571 : !! This routine prints the Hermitian and the non-hermitian part of the matrix
572 : !! elements of Sigma, as well as the individual contributions.
573 : !! The first 14x14 are printed to screen, and the full matrices are printed
574 : !! to files: sigma_melems_, sigma_nonH_melems_, sigma_Hart_melems_,
575 : !! sigma_x_melems, and sigma_c_melems
576 : !!
577 : !! INPUTS
578 : !! ikcalc : index of k-point
579 : !! ib1,ib2 : starting and ending band indices
580 : !! nsp : no. of spin elements
581 : !! htotal : Hermitianised matrix elements of Sigma
582 : !! hhartree : Hartree contribution to matrix elements
583 : !! sigxme : Sigma_x contribution to matrix elements
584 : !! sigcme : Sigma_c contribution to matrix elements
585 : !! prefix : prefix for output files.
586 : !!
587 : !! OUTPUT
588 : !!
589 : !! SOURCE
590 :
591 0 : subroutine print_sigma_melems(ikcalc, ib1, ib2, nsp, htotal, hhartree, sigxme, sigcme, prefix)
592 :
593 : ! Arguments ------------------------------------
594 : !scalars
595 : integer,intent(in) :: ikcalc,ib1,ib2,nsp
596 : character(len=*),intent(in) :: prefix
597 : !arrays
598 : complex(dp),intent(in) :: htotal(ib1:ib2,ib1:ib2,nsp),hhartree(ib1:ib2,ib1:ib2,nsp)
599 : complex(dp),intent(in) :: sigxme(ib1:ib2,ib1:ib2,nsp),sigcme(ib1:ib2,ib1:ib2,nsp)
600 :
601 : ! Local variables ------------------------------
602 : integer,parameter :: MAX_NCOLS = 14
603 : integer :: isp,mc,mr,jj,ii,temp_unit,ount
604 : character(len=10) :: sidx
605 : character(len=500) :: msg
606 : character(len=100) :: fmth,fmt1,fmt2,fmthh,kpt_index,fmtfile
607 : character(len=fnlen) :: filename
608 : ! *************************************************************************
609 :
610 0 : if (nsp==3.or.nsp>4) then
611 0 : ABI_ERROR('nsp has wrong value in print_sigma_melems')
612 : end if
613 :
614 0 : ount = std_out
615 :
616 0 : mc = ib2-ib1+1; if (mc>MAX_NCOLS) mc = MAX_NCOLS
617 0 : mr = mc
618 :
619 0 : write(fmthh,*)'(2(a),2(I2,a))'
620 0 : write(fmth,*)'(7x,',mc,'(i2,8x))'
621 0 : write(fmt1,*)'(3x,i2,',mc,'f10.5)'
622 0 : write(fmt2,*)'(5x ,',mc,'f10.5,a)'
623 :
624 : ! First print to screen
625 0 : do isp=1,nsp
626 0 : write(msg,'(a)') ''
627 0 : call wrtout(ount,msg)
628 0 : write(msg,fmthh) ch10,' Hermitianised matrix elements of Sigma (spin ',isp,' of ',nsp,'):'
629 0 : call wrtout(ount,msg)
630 0 : write(msg,fmth)(jj,jj=1,mc)
631 0 : call wrtout(ount,msg) !header
632 0 : do ii=ib1,ib1+mr-1
633 0 : write(msg,fmt1)ii-ib1+1,DBLE(htotal(ii,ib1:(ib1+mc-1),isp))
634 0 : call wrtout(ount,msg) !real part
635 0 : write(msg,fmt2) AIMAG(htotal(ii,ib1:(ib1+mc-1),isp)),ch10
636 0 : call wrtout(ount,msg) !imag part
637 : end do
638 : end do !nsp
639 :
640 0 : write(msg,'(a,i2,a)')" Max. ",MAX_NCOLS," elements printed. Full matrix output in _HTOTAL files"
641 0 : call wrtout(ount,msg)
642 :
643 0 : do isp=1,nsp
644 0 : write(msg,fmthh) ch10,' H_Hartree matrix elements (spin ',isp,' of ',nsp,'):'
645 0 : call wrtout(ount,msg)
646 0 : write(msg,fmth)(jj,jj=1,mc)
647 0 : call wrtout(ount,msg) !header
648 0 : do ii=ib1,ib1+mr-1
649 0 : write(msg,fmt1)ii-ib1+1,DBLE(hhartree(ii,ib1:(ib1+mc-1),isp))
650 0 : call wrtout(ount,msg) !real part
651 0 : write(msg,fmt2) AIMAG(hhartree(ii,ib1:(ib1+mc-1),isp)),ch10
652 0 : call wrtout(ount,msg) !imag part
653 : end do
654 : end do !nsp
655 :
656 0 : write(msg,'(a,i2,a)')" Max. ",MAX_NCOLS," elements printed. Full matrix output in _HHARTREE files"
657 0 : call wrtout(ount,msg)
658 :
659 0 : do isp=1,nsp
660 0 : write(msg,fmthh) ch10,' Sigma_x matrix elements (spin ',isp,' of ',nsp,'):'
661 0 : call wrtout(ount,msg)
662 0 : write(msg,fmth)(jj,jj=1,mc)
663 0 : call wrtout(ount,msg) !header
664 0 : do ii=ib1,ib1+mr-1
665 0 : write(msg,fmt1)ii-ib1+1,DBLE(sigxme(ii,ib1:(ib1+mc-1),isp))
666 0 : call wrtout(ount,msg) !real part
667 0 : write(msg,fmt2) AIMAG(sigxme(ii,ib1:(ib1+mc-1),isp)),ch10
668 0 : call wrtout(ount,msg) !imag part
669 : end do
670 : end do !nsp
671 :
672 0 : write(msg,'(a,i2,a)')" Max. ",MAX_NCOLS," elements printed. Full matrix output _SIGX files"
673 0 : call wrtout(ount,msg)
674 :
675 0 : do isp=1,nsp
676 0 : write(msg,fmthh) ch10,' Sigma_c matrix elements (spin ',isp,' of ',nsp,'):'
677 0 : call wrtout(ount,msg)
678 0 : write(msg,fmth)(jj,jj=1,mc)
679 0 : call wrtout(ount,msg) !header
680 0 : do ii=ib1,ib1+mr-1
681 0 : write(msg,fmt1)ii-ib1+1,DBLE(sigcme(ii,ib1:(ib1+mc-1),isp))
682 0 : call wrtout(ount,msg) !real part
683 0 : write(msg,fmt2) AIMAG(sigcme(ii,ib1:(ib1+mc-1),isp)),ch10
684 0 : call wrtout(ount,msg) !imag part
685 : end do
686 : end do !nsp
687 :
688 0 : write(msg,'(a,i2,a)')" Max ",MAX_NCOLS," elements printed. Full matrix output _SIGC files"
689 0 : call wrtout(ount,msg)
690 :
691 : ! Then print to file
692 : ! Format is: row, column, value; with a blank space for each full
693 : ! set of columns for easy plotting with the gnuplot splot command
694 0 : write(fmtfile,*)'(3X,I6,2X,I6,',nsp,'(2(ES28.16E3,3x)))'
695 :
696 0 : call int2char10(ikcalc,sidx)
697 0 : kpt_index = "_KPT"//TRIM(sidx)
698 :
699 0 : filename = TRIM(prefix)//'_HTOTAL'//TRIM(kpt_index)
700 :
701 0 : if (open_file(filename,msg,newunit=temp_unit,form="formatted",status="replace",action="write") /= 0) then
702 0 : ABI_ERROR(msg)
703 : end if
704 :
705 0 : msg = '# row col. Re(htotal(r,c)) Im(htotal(r,c)) for spin11 ... spin22 ... spin12 ... spin13'
706 0 : call wrtout(temp_unit,msg)
707 0 : do ii=ib1,ib2
708 0 : do jj=ib1,ib2
709 0 : write(msg,fmtfile) ii,jj,(htotal(jj,ii,isp),isp=1,nsp)
710 0 : call wrtout(temp_unit,msg)
711 : end do
712 0 : call wrtout(temp_unit,"")
713 : end do
714 0 : close(temp_unit)
715 :
716 0 : filename = TRIM(prefix)//'_HHARTREE'//TRIM(kpt_index)
717 0 : if (open_file(filename,msg,newunit=temp_unit,form="formatted",status="replace",action="write") /= 0) then
718 0 : ABI_ERROR(msg)
719 : end if
720 :
721 0 : msg = '# row col. Re(hhartree(r,c)) Im(hhartree(r,c) for spin11 ... spin22 ... spin12 ... spin13'
722 0 : call wrtout(temp_unit,msg)
723 0 : do ii=ib1,ib2
724 0 : do jj=ib1,ib2
725 0 : write(msg,fmtfile) ii,jj,(hhartree(jj,ii,isp),isp=1,nsp)
726 0 : call wrtout(temp_unit,msg)
727 : end do
728 0 : call wrtout(temp_unit,"")
729 : end do
730 0 : close(temp_unit)
731 :
732 0 : filename = TRIM(prefix)//'_SIGX'//TRIM(kpt_index)
733 0 : if (open_file(filename,msg,newunit=temp_unit,form="formatted",status="replace",action="write") /= 0) then
734 0 : ABI_ERROR(msg)
735 : end if
736 :
737 0 : write(msg,'(a)')'# row col. Re(Sigx(r,c)) Im(Sigx(r,c) for spin11 ... spin22 ... spin12 ... spin13'
738 0 : call wrtout(temp_unit,msg)
739 0 : do ii=ib1,ib2
740 0 : do jj=ib1,ib2
741 0 : write(msg,fmtfile) ii,jj,(sigxme(jj,ii,isp),isp=1,nsp)
742 0 : call wrtout(temp_unit,msg)
743 : end do
744 0 : call wrtout(temp_unit,"")
745 : end do
746 0 : close(temp_unit)
747 :
748 0 : filename = TRIM(prefix)//'_SIGC'//TRIM(kpt_index)
749 0 : if (open_file(filename,msg,newunit=temp_unit,form="formatted",status="replace",action="write") /= 0) then
750 0 : ABI_ERROR(msg)
751 : end if
752 :
753 0 : write(msg,'(a)')'# row col. Re(Sigc(r,c)) Im(Sigc(r,c) for spin11 ... spin22 ... spin12 ... spin21'
754 0 : call wrtout(temp_unit,msg)
755 0 : do ii=ib1,ib2
756 0 : do jj=ib1,ib2
757 0 : write(msg,fmtfile) ii,jj,(sigcme(jj,ii,isp),isp=1,nsp)
758 0 : call wrtout(temp_unit,msg)
759 : end do
760 0 : call wrtout(temp_unit,"")
761 : end do
762 :
763 0 : close(temp_unit)
764 :
765 0 : end subroutine print_sigma_melems
766 : !!***
767 :
768 : !----------------------------------------------------------------------
769 :
770 : !!****f* m_dyson_solver/sigma_pade_init
771 : !! NAME
772 : !! sigma_pade_init
773 : !!
774 : !! FUNCTION
775 : !! Initialize the Pade' from the `npts` values of Sigma_c(iw) given on the mesh `zmesh`.
776 : !!
777 : !! SOURCE
778 :
779 106 : subroutine sigma_pade_init(self, npts, zmesh, sigc_cvals, alphac_pm, betar_pm, zcut_pm)
780 :
781 : !Arguments ------------------------------------
782 : class(sigma_pade_t),intent(out) :: self
783 : integer,intent(in) :: npts
784 : complex(dp),target,intent(in) :: zmesh(npts), sigc_cvals(npts), alphac_pm(2)
785 : real(dp),intent(in) :: betar_pm(2), zcut_pm(2)
786 : ! *************************************************************************
787 :
788 106 : self%npts = npts
789 318 : ABI_MALLOC(self%zmesh, (npts))
790 212 : ABI_MALLOC(self%sigc_cvals, (npts,1))
791 1432 : self%zmesh = zmesh
792 1326 : self%sigc_cvals(:,1) = sigc_cvals
793 106 : self%nsig_ab = 1
794 :
795 318 : self%alphac_pm = alphac_pm
796 318 : self%betar_pm = betar_pm
797 318 : self%zcut_pm = zcut_pm
798 106 : self%do_sigma_fit = .False.
799 :
800 106 : end subroutine sigma_pade_init
801 : !!***
802 : !----------------------------------------------------------------------
803 :
804 : !!****f* m_dyson_solver/sigma_pade_init_spinor
805 : !! NAME
806 : !! sigma_pade_init_spinor
807 : !!
808 : !! FUNCTION
809 : !! Initialize the Pade' from the `npts` values of Sigma_c(iw) given on the mesh `zmesh`.
810 : !!
811 : !! SOURCE
812 :
813 0 : subroutine sigma_pade_init_spinor(self, npts, zmesh, nsig_ab, sigc_cvals, alphac_pm, betar_pm, zcut_pm)
814 :
815 : !Arguments ------------------------------------
816 : class(sigma_pade_t),intent(out) :: self
817 : integer,intent(in) :: npts, nsig_ab
818 : complex(dp),target,intent(in) :: zmesh(npts), sigc_cvals(npts,nsig_ab), alphac_pm(2)
819 : real(dp),intent(in) :: betar_pm(2), zcut_pm(2)
820 : ! *************************************************************************
821 :
822 0 : self%npts = npts
823 0 : ABI_MALLOC(self%zmesh, (npts))
824 0 : ABI_MALLOC(self%sigc_cvals, (npts,nsig_ab))
825 0 : self%zmesh = zmesh
826 0 : self%sigc_cvals = sigc_cvals
827 0 : self%nsig_ab = nsig_ab
828 :
829 0 : self%alphac_pm = alphac_pm
830 0 : self%betar_pm = betar_pm
831 0 : self%zcut_pm = zcut_pm
832 0 : self%do_sigma_fit = .False.
833 :
834 0 : end subroutine sigma_pade_init_spinor
835 : !!***
836 :
837 :
838 106 : subroutine sigma_pade_free(self)
839 : !Arguments ------------------------------------
840 : class(sigma_pade_t),intent(inout) :: self
841 : ! *************************************************************************
842 :
843 106 : ABI_SFREE(self%zmesh)
844 106 : ABI_SFREE(self%sigc_cvals)
845 :
846 106 : end subroutine sigma_pade_free
847 : !!***
848 :
849 : !----------------------------------------------------------------------
850 :
851 : !!****f* m_dyson_solver/sigma_pade_eval
852 : !! NAME
853 : !! sigma_pade_eval
854 : !!
855 : !! FUNCTION
856 : !! Evaluate the Pade' at the complex point `zz`.
857 : !! Return result in `val` and, optionally, the derivative at zz in `dvdz`
858 : !!
859 : !! SOURCE
860 :
861 215 : subroutine sigma_pade_eval(self, zz, val, &
862 : dvdz) ! optional
863 :
864 : !Arguments ------------------------------------
865 : class(sigma_pade_t),intent(in) :: self
866 : complex(dp),intent(in) :: zz
867 : complex(dp),intent(out) :: val
868 : complex(dp),optional,intent(out) :: dvdz
869 : integer :: iab
870 : ! *************************************************************************
871 :
872 : ! if zz in 2 or 3 quadrant, avoid branch cut in the complex plane using Sigma(-iw) = Sigma(iw)*.
873 215 : val = czero
874 215 : if (present(dvdz)) dvdz = czero
875 :
876 430 : do iab = 1, self%nsig_ab
877 430 : if (real(zz) > zero) then
878 104 : val = val + pade(self%npts, self%zmesh, self%sigc_cvals(:,iab), zz)
879 104 : if (present(dvdz)) then
880 104 : dvdz = dvdz + dpade(self%npts, self%zmesh, self%sigc_cvals(:,iab), zz)
881 : end if
882 : else
883 2651 : val = val + pade(self%npts, -self%zmesh, conjg(self%sigc_cvals(:,iab)), zz)
884 111 : if (present(dvdz)) then
885 2651 : dvdz = dvdz + dpade(self%npts, -self%zmesh, conjg(self%sigc_cvals(:,iab)), zz)
886 : end if
887 : end if
888 : end do
889 :
890 215 : if (self%do_sigma_fit) then
891 : ! Add analytic expression.
892 : val = val + self%alphac_pm(1) / (self%betar_pm(1) + zz) &
893 0 : + self%alphac_pm(2) / (self%betar_pm(2) - zz)
894 :
895 0 : if (present(dvdz)) then
896 : ! Add analytic expression.
897 : dvdz = dvdz - self%alphac_pm(1) / ((self%betar_pm(1) + zz) ** 2) &
898 0 : - self%alphac_pm(2) / ((self%betar_pm(2) - zz) ** 2)
899 : end if
900 : end if
901 :
902 215 : end subroutine sigma_pade_eval
903 : !!***
904 :
905 : !----------------------------------------------------------------------
906 :
907 : !!****f* m_dyson_solver/sigma_pade_qp_solve
908 : !! NAME
909 : !! sigma_pade_qp_solve
910 : !!
911 : !! FUNCTION
912 : !! Use the Pade' approximant and Newton-Rapson method to solve the QP equation
913 : !! in the complex plane starting from the initial guess `z_guess`.
914 : !!
915 : !! INPUTS
916 : !! e0: KS energy
917 : !! v_meanf: matrix element of the mean-field Hamiltonian
918 : !! sigx: matrix element of the exchange self-energy.
919 : !! z_guess: Initial guess for the QP energy
920 : !!
921 : !! OUTPUT
922 : !! zsc: root.
923 : !! sigc: Sigma_c(zsc)
924 : !! msg: Error message if ierr /= 0.
925 : !! ierr: Exit status.
926 : !!
927 : !! SOURCE
928 :
929 53 : subroutine sigma_pade_qp_solve(self, e0, v_meanf, sigx, z_guess, zsc, sigc_zsc, msg, ierr)
930 :
931 : !Arguments ------------------------------------
932 : class(sigma_pade_t),intent(in) :: self
933 : real(dp),intent(in) :: e0, v_meanf, sigx
934 : complex(dp),intent(in) :: z_guess
935 : complex(dp),intent(out) :: zsc, sigc_zsc
936 : integer,intent(out) :: ierr
937 :
938 : !Local variables-------------------------------
939 : !scalars
940 : integer :: iter
941 : logical :: converged
942 : complex(dp) :: ctdpc, dct, dsigc
943 : character(len=500) :: msg
944 : ! *************************************************************************
945 :
946 : ! Use Newton-Rapson to find the root of:
947 : ! f(z) = e0 - zz + Sigma_xc(z) - v_meanf
948 : ! f'(z) = -1 + Sigma_c'(z)
949 :
950 53 : iter = 0; converged = .FALSE.; ctdpc = cone
951 53 : zsc = z_guess
952 162 : do while (abs(ctdpc) > NR_ABS_ROOT_ERR .or. iter < NR_MAX_NITER)
953 162 : iter = iter + 1
954 :
955 162 : call self%eval(zsc, sigc_zsc, dvdz=dsigc)
956 162 : ctdpc = e0 - v_meanf + sigx + sigc_zsc - zsc
957 :
958 162 : if (abs(ctdpc) < NR_ABS_ROOT_ERR) then
959 : converged=.TRUE.; EXIT
960 : end if
961 109 : dct = dsigc - one
962 109 : zsc = newrap_step(zsc, ctdpc, dct)
963 : end do
964 :
965 53 : ierr = 0; msg = ""
966 53 : if (.not. converged) then
967 : write(msg,'(a,i0,3a,f8.4,a,f8.4)')&
968 0 : 'Newton-Raphson method did not converge after: ', NR_MAX_NITER,' iterations.',ch10,&
969 0 : 'Absolute error: ', abs(ctdpc), ' > ', NR_ABS_ROOT_ERR
970 0 : ierr = 1
971 : end if
972 :
973 53 : end subroutine sigma_pade_qp_solve
974 : !!***
975 :
976 : !----------------------------------------------------------------------
977 :
978 212 : end module m_dyson_solver
979 : !!***
|
