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 : 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 the analytic transform of the model removed in imaginary time:
892 : : !
893 : : ! Sigma(+tau) = alpha_1 e^{-beta_1 tau} --> alpha_1 / (beta_1 - z)
894 : : ! Sigma(-tau) = alpha_2 e^{-beta_2 tau} --> alpha_2 / (beta_2 + z)
895 : : !
896 : : ! with tau > 0, consistent with the cos/sin transforms used for the residual.
897 : : val = val + self%alphac_pm(1) / (self%betar_pm(1) - zz) &
898 : 0 : + self%alphac_pm(2) / (self%betar_pm(2) + zz)
899 : :
900 [ # # ]: 0 : if (present(dvdz)) then
901 : : dvdz = dvdz + self%alphac_pm(1) / ((self%betar_pm(1) - zz) ** 2) &
902 : 0 : - self%alphac_pm(2) / ((self%betar_pm(2) + zz) ** 2)
903 : : end if
904 : : end if
905 : :
906 : 215 : end subroutine sigma_pade_eval
907 : : !!***
908 : :
909 : : !----------------------------------------------------------------------
910 : :
911 : : !!****f* m_dyson_solver/sigma_pade_qp_solve
912 : : !! NAME
913 : : !! sigma_pade_qp_solve
914 : : !!
915 : : !! FUNCTION
916 : : !! Use the Pade' approximant and Newton-Rapson method to solve the QP equation
917 : : !! in the complex plane starting from the initial guess `z_guess`.
918 : : !!
919 : : !! INPUTS
920 : : !! e0: KS energy
921 : : !! v_meanf: matrix element of the mean-field Hamiltonian
922 : : !! sigx: matrix element of the exchange self-energy.
923 : : !! z_guess: Initial guess for the QP energy
924 : : !!
925 : : !! OUTPUT
926 : : !! zsc: root.
927 : : !! sigc: Sigma_c(zsc)
928 : : !! msg: Error message if ierr /= 0.
929 : : !! ierr: Exit status.
930 : : !!
931 : : !! SOURCE
932 : :
933 : 53 : subroutine sigma_pade_qp_solve(self, e0, v_meanf, sigx, z_guess, zsc, sigc_zsc, msg, ierr)
934 : :
935 : : !Arguments ------------------------------------
936 : : class(sigma_pade_t),intent(in) :: self
937 : : real(dp),intent(in) :: e0, v_meanf, sigx
938 : : complex(dp),intent(in) :: z_guess
939 : : complex(dp),intent(out) :: zsc, sigc_zsc
940 : : integer,intent(out) :: ierr
941 : :
942 : : !Local variables-------------------------------
943 : : !scalars
944 : : integer :: iter
945 : : logical :: converged
946 : : complex(dp) :: ctdpc, dct, dsigc
947 : : character(len=500) :: msg
948 : : ! *************************************************************************
949 : :
950 : : ! Use Newton-Rapson to find the root of:
951 : : ! f(z) = e0 - zz + Sigma_xc(z) - v_meanf
952 : : ! f'(z) = -1 + Sigma_c'(z)
953 : :
954 : 53 : iter = 0; converged = .FALSE.; ctdpc = cone
955 : 53 : zsc = z_guess
956 [ - + - - ]: 162 : do while (abs(ctdpc) > NR_ABS_ROOT_ERR .or. iter < NR_MAX_NITER)
957 : 162 : iter = iter + 1
958 : :
959 : 162 : call self%eval(zsc, sigc_zsc, dvdz=dsigc)
960 : 162 : ctdpc = e0 - v_meanf + sigx + sigc_zsc - zsc
961 : :
962 [ + + ]: 162 : if (abs(ctdpc) < NR_ABS_ROOT_ERR) then
963 : : converged=.TRUE.; EXIT
964 : : end if
965 : 109 : dct = dsigc - one
966 : 109 : zsc = newrap_step(zsc, ctdpc, dct)
967 : : end do
968 : :
969 : 53 : ierr = 0; msg = ""
970 [ - + ]: 53 : if (.not. converged) then
971 : : write(msg,'(a,i0,3a,f8.4,a,f8.4)')&
972 : 0 : 'Newton-Raphson method did not converge after: ', NR_MAX_NITER,' iterations.',ch10,&
973 : 0 : 'Absolute error: ', abs(ctdpc), ' > ', NR_ABS_ROOT_ERR
974 : 0 : ierr = 1
975 : : end if
976 : :
977 : 53 : end subroutine sigma_pade_qp_solve
978 : : !!***
979 : :
980 : : !----------------------------------------------------------------------
981 : :
982 [ - + - - : 212 : end module m_dyson_solver
+ + - + +
- - + -
+ ]
983 : : !!***
|
