Branch data Line data Source code
1 : : !!****m* ABINIT/m_haydock
2 : : !! NAME
3 : : !! m_haydock
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 2008-2026 ABINIT group (M.Giantomassi, Y. Gillet, L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida)
9 : : !! This file is distributed under the terms of the
10 : : !! GNU General Public License, see ~abinit/COPYING
11 : : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : : !!
13 : : !! SOURCE
14 : :
15 : : #if defined HAVE_CONFIG_H
16 : : #include "config.h"
17 : : #endif
18 : :
19 : : #include "abi_common.h"
20 : :
21 : : MODULE m_haydock
22 : :
23 : : use defs_basis
24 : : use m_abicore
25 : : use m_bs_defs
26 : : use m_xmpi
27 : : use m_errors
28 : : use m_haydock_io
29 : : use m_linalg_interfaces
30 : : use netcdf
31 : : use m_nctk
32 : :
33 : : use m_time, only : timab
34 : : use m_fstrings, only : strcat, sjoin, itoa, int2char4
35 : : use m_io_tools, only : file_exists, open_file
36 : : use defs_datatypes, only : pseudopotential_type
37 : : use m_geometry, only : normv
38 : : use m_hide_blas, only : xdotc, xgemv
39 : : use m_hide_lapack, only : matrginv
40 : : use m_numeric_tools, only : print_arr, symmetrize, hermitianize, continued_fract, wrap2_pmhalf, iseven
41 : : use m_kpts, only : listkk
42 : : use m_crystal, only : crystal_t
43 : : use m_ebands, only : ebands_t
44 : : use m_hdr, only : hdr_type
45 : : use m_bz_mesh, only : kmesh_t, findqg0
46 : : use m_double_grid, only : double_grid_t, get_kpt_from_indices_coarse, compute_corresp
47 : : use m_paw_hr, only : pawhur_t
48 : : use m_wfd, only : wfdgw_t
49 : : use m_bse_io, only : exc_write_optme
50 : : use m_pawtab, only : pawtab_type
51 : : use m_vcoul, only : vcoul_t
52 : : use m_hexc, only : hexc_t, hexc_interp_t
53 : : use m_exc_spectra, only : exc_write_data, exc_eps_rpa, exc_write_tensor, mdfs_ncwrite
54 : : use m_eprenorms, only : eprenorms_t, renorm_bst
55 : : use m_wfd_optic, only : calc_optical_mels
56 : :
57 : : implicit none
58 : :
59 : : private
60 : : !!***
61 : :
62 : : public :: exc_haydock_driver ! Driver for the Haydock method (main entry point for client code).
63 : :
64 : : CONTAINS !=======================================================================
65 : : !!***
66 : :
67 : : !!****f* m_haydock/exc_haydock_driver
68 : : !! NAME
69 : : !! exc_haydock_driver
70 : : !!
71 : : !! FUNCTION
72 : : !! Calculate the imaginary part of the macroscopic dielectric function with the Haydock recursive method.
73 : : !!
74 : : !! INPUTS
75 : : !! BSp<type(excparam)=The parameter for the Bethe-Salpeter run.
76 : : !! BS_files<excparam>=Files associated to the bethe_salpeter code.
77 : : !! Cryst<crystal_t>=Info on the crystalline structure.
78 : : !! Kmesh<type(kmesh_t)>=The list of k-points in the BZ, IBZ and symmetry tables.
79 : : !! Cryst<type(crystal_t)>=Info on the crystalline structure.
80 : : !! Hdr_bse
81 : : !! KS_BSt=The KS energies.
82 : : !! QP_BSt=The QP energies.
83 : : !! Wfd<wfdgw_t>=Wavefunction descriptor (input k-mesh)
84 : : !! Psps <type(pseudopotential_type)>=variables related to pseudopotentials.
85 : : !! Pawtab(Cryst%ntypat*usepaw)<pawtab_type>=PAW tabulated starting data.
86 : : !! Hur(Cryst%natom*usepaw)<type(pawhur_t)>=Only for PAW and DFT+U, quantities used to evaluate the commutator [H_u,r].
87 : : !!
88 : : !! OUTPUT
89 : : !! The imaginary part of the macroscopic dielectric function is written on the external file _EXC_MDF
90 : : !!
91 : : !! SOURCE
92 : :
93 : 23 : subroutine exc_haydock_driver(BSp,BS_files,Cryst,Kmesh,Hdr_bse,KS_BSt,QP_Bst,Wfd,Psps,Pawtab,Hur,Epren, &
94 : : Kmesh_dense, KS_BSt_dense, QP_BSt_dense, Wfd_dense, Vcp_dense, grid) ! Optional args
95 : :
96 : : !Arguments ------------------------------------
97 : : !scalars
98 : : type(excparam),intent(in) :: BSp
99 : : type(excfiles),intent(in) :: BS_files
100 : : type(kmesh_t),intent(in) :: Kmesh
101 : : type(crystal_t),intent(in) :: Cryst
102 : : type(Hdr_type),intent(in) :: Hdr_bse
103 : : type(wfdgw_t),intent(inout) :: Wfd
104 : : type(pseudopotential_type),intent(in) :: Psps
105 : : type(ebands_t),intent(in) :: KS_BSt,QP_Bst
106 : : type(double_grid_t),intent(in),optional :: grid
107 : : type(kmesh_t),intent(in),optional :: Kmesh_dense
108 : : type(wfdgw_t),intent(inout),optional :: Wfd_dense
109 : : type(ebands_t),intent(in),optional :: KS_BSt_dense, QP_Bst_dense
110 : : type(vcoul_t),intent(in),optional :: Vcp_dense
111 : : type(eprenorms_t),intent(in) :: Epren
112 : : !arrays
113 : : type(pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd%usepaw)
114 : : type(pawhur_t),intent(in) :: Hur(Cryst%natom*Wfd%usepaw)
115 : :
116 : : !Local variables ------------------------------
117 : : !scalars
118 : : integer,parameter :: master=0
119 : : integer :: io,my_rank,iq,itt,ierr
120 : : integer :: hsize,comm,my_t1,my_t2,nsppol,nkets,nproc,ncid
121 : : integer :: spin,spad,ik_bz,iv,ic,trans_idx,lomo_min,max_band
122 : : real(dp) :: omegaev,rand_phi !,norm
123 : : complex(dp) :: ks_avg,gw_avg,exc_avg
124 : : logical :: use_mpio,prtdos
125 : : character(len=500) :: msg
126 : 23 : type(hexc_t) :: hexc
127 : 23 : type(hexc_interp_t) :: hexc_i
128 : : !arrays
129 : : real(dp) :: tsec(2)
130 : 23 : real(dp),allocatable :: dos(:),dos_gw(:),dos_ks(:)
131 : 23 : complex(dp),allocatable :: green(:,:)
132 : 23 : complex(dp),allocatable :: opt_cvk(:,:,:,:,:),kets(:,:)
133 : 23 : complex(dp),allocatable :: eps_rpanlf(:,:),eps_gwnlf(:,:)
134 : 23 : complex(dp),allocatable :: tensor_cart(:,:),tensor_cart_rpanlf(:,:),tensor_cart_gwnlf(:,:)
135 : 23 : complex(dp),allocatable :: tensor_red(:,:),tensor_red_rpanlf(:,:),tensor_red_gwnlf(:,:)
136 : :
137 : : !Temperature
138 : : integer :: ep_ik, ik, ireh, isppol, itemp, ntemp
139 : : real(dp) :: dksqmax, en
140 : 23 : integer,allocatable :: bs2eph(:,:)
141 : : integer :: sppoldbl, timrev
142 : : logical :: do_ep_renorm, do_ep_lifetime
143 : : character(len=4) :: ts
144 : : character(len=fnlen) :: prefix, path
145 : 23 : complex(dp),allocatable :: ep_renorms(:)
146 : 23 : type(ebands_t) :: EPBSt, EP_QPBSt
147 : : !************************************************************************
148 : :
149 : 23 : call timab(690,1,tsec) ! exc_haydock_driver
150 : 23 : call timab(691,1,tsec) ! exc_haydock_driver(read)
151 : :
152 [ - + ]: 23 : if (BSp%have_complex_ene) then
153 : 0 : ABI_ERROR("Complex energies are not supported yet")
154 : : end if
155 : :
156 : 23 : my_rank = Wfd%my_rank
157 : 23 : comm = Wfd%comm
158 : 23 : nsppol = Wfd%nsppol
159 : 23 : nproc = Wfd%nproc
160 : :
161 : : use_mpio=.FALSE.
162 : : #ifdef HAVE_MPI_IO
163 : : use_mpio = (nproc > 1)
164 : : !use_mpio = .TRUE.
165 : : #endif
166 : 23 : use_mpio=.FALSE.
167 : : !use_mpio = .TRUE.
168 : :
169 : : ! Hsize refers to the size of the individual blocks (resonant and coupling).
170 : : ! Thanks to the symmetry property of the starting vector, the Haydock method
171 : : ! can be reformulated in terms of matrix-vector multiplication involving the
172 : : ! blocks thus avoiding to allocation of the full matrix ( R C )
173 : : ! -C* -R*)
174 [ + + ]: 46 : hsize=SUM(BSp%nreh)
175 : :
176 : : !YG2014
177 : 23 : call hexc%init(BSp, BS_files, Cryst, Kmesh, Wfd, KS_BSt, QP_BSt, comm)
178 : :
179 : : !YG2014
180 [ + + ]: 23 : if(BSp%use_interp) then
181 : : call hexc_i%init(hexc, BSp%interp_m3_width, BSp%interp_method,&
182 : : Kmesh_dense, Vcp_dense, grid, Wfd_dense, &
183 : 4 : KS_BSt_dense, QP_BSt_dense, Psps, Pawtab)
184 : : end if
185 : :
186 : 23 : call timab(691,2,tsec) ! exc_haydock_driver(read)
187 : 23 : call timab(692,1,tsec) ! exc_haydock_driver(prep)
188 : :
189 : : ! Prepare the starting vectors for the Lanczos chain.
190 : 23 : nkets=Bsp%nq
191 : :
192 : 23 : prtdos=.FALSE. !prtdos=.TRUE.
193 : : if (prtdos) then
194 : : nkets=nkets+1
195 : : if (Bsp%use_coupling>0) then
196 : : ABI_ERROR("DOS with coupling not coded")
197 : : nkets=nkets+1
198 : : end if
199 : : end if
200 : :
201 : : !YG2014
202 [ - + + - : 92 : ABI_MALLOC_OR_DIE(kets,(hexc%hsize,nkets), ierr)
- + + - +
- + - + -
- + - + ]
203 [ + + + + ]: 64289 : kets=czero
204 : : !
205 : : ! Prepare the kets for the macroscopic dielectric function.
206 : 23 : lomo_min=Bsp%lomo_min; max_band=Bsp%nbnds
207 : :
208 : : !YG2014
209 [ - + + - : 161 : ABI_MALLOC_OR_DIE(opt_cvk,(lomo_min:max_band,lomo_min:max_band,hexc%nbz,Wfd%nsppol,BSp%nq), ierr)
- + + - -
+ + - - +
+ - - + +
- + - + -
+ - - + -
+ ]
210 : :
211 [ + + ]: 161 : do iq=1,Bsp%nq
212 : : ! Note KS_BSt is used here to calculate the commutator.
213 : : call calc_optical_mels(hexc%Wfd,hexc%Kmesh,hexc%KS_BSt,Cryst,Psps,Pawtab,Hur, &
214 : 138 : & BSp%inclvkb,BSp%lomo_spin,lomo_min,max_band,hexc%nbz,BSp%q(:,iq),opt_cvk(:,:,:,:,iq))
215 : :
216 : : ! Fill ket0 using the same ordering for the indices as the one used for the excitonic Hamiltonian.
217 : : ! Note that only the resonant part is used here.
218 [ + + ]: 299 : do spin=1,nsppol
219 : :
220 [ + + ]: 138 : if(BSp%use_interp) then
221 : 24 : spad=(spin-1)*BSp%nreh_interp(spin)
222 : : else
223 : 114 : spad=(spin-1)*BSp%nreh(spin)
224 : : end if
225 : :
226 [ + + ]: 7524 : do ik_bz=1,hexc%nbz
227 [ + + ]: 28794 : do iv=BSp%lomo_spin(spin),BSp%homo_spin(spin)
228 [ + + ]: 92784 : do ic=BSp%lumo_spin(spin),BSp%nbnds
229 : :
230 [ + + ]: 64128 : if(BSp%use_interp) then
231 : 18432 : trans_idx = BSp%vcks2t_interp(iv,ic,ik_bz,spin)
232 : : else
233 : 45696 : trans_idx = BSp%vcks2t(iv,ic,ik_bz,spin)
234 : : end if
235 : :
236 [ + - ]: 85536 : if (trans_idx>0) kets(trans_idx+spad,iq)=opt_cvk(ic,iv,ik_bz,spin,iq)
237 : : end do
238 : : end do
239 : : end do
240 : : end do
241 : :
242 : : end do
243 : :
244 : : ! ========================================================
245 : : ! === Write the Optical Matrix Elements to NetCDF file ===
246 : : ! ========================================================
247 : :
248 : : !if (.false.) then
249 : : ! ome_fname='test_OME.nc'
250 : : ! call exc_write_optme(ome_fname,minb,maxb,BSp%nkbz,Wfd%nsppol,BSp%nq,opt_cvk,ierr)
251 : : !end if
252 : :
253 : : ! Free WFD descriptor, we don't need ur and ug anymore !
254 : : ! We make space for interpolated hamiltonian
255 : 23 : call wfd%wave_free("All")
256 [ + + ]: 23 : if(BSp%use_interp) call wfd_dense%wave_free("All")
257 : :
258 : : ! Build interpolated hamiltonian
259 [ + + ]: 23 : if(BSp%use_interp) then
260 [ + + ]: 4 : if (any(BSp%interp_mode == [2,3,4])) call hexc%build_hinterp(hexc_i)
261 : : end if
262 : :
263 : 23 : call timab(692,2,tsec) ! exc_haydock_driver(prep)
264 : 23 : call timab(693,1,tsec) ! exc_haydock_driver(wo lf) - that is, without local field
265 : :
266 : 23 : do_ep_renorm = .FALSE.
267 : 23 : ntemp = 1
268 : 23 : do_ep_lifetime = .FALSE.
269 : :
270 [ + + ]: 23 : if(BSp%do_ep_renorm) then
271 [ - + ]: 1 : if (BSp%nsppol == 2) then
272 : 0 : ABI_ERROR('Elphon renorm with nsppol == 2 not yet coded !')
273 : : end if
274 : 1 : do_ep_renorm = .TRUE.
275 : 1 : ntemp = Epren%ntemp
276 : : if(BSp%do_lifetime) do_ep_lifetime = .TRUE.
277 : :
278 : : ! Force elphon linewidth
279 : 1 : do_ep_lifetime = .TRUE.
280 : :
281 : : ! Map points from BSE to elphon kpoints
282 : 1 : sppoldbl = 1 !; if (any(Cryst%symafm == -1) .and. Epren%nsppol == 1) nsppoldbl=2
283 [ - + + - : 3 : ABI_MALLOC(bs2eph, (Kmesh%nbz*sppoldbl, 6))
+ - + - -
+ - + ]
284 [ - + + - : 3 : ABI_MALLOC(ep_renorms, (hsize))
+ - + - -
+ - + ]
285 : 1 : timrev = 1
286 : : call listkk(dksqmax, Cryst%gmet, bs2eph, Epren%kpts, Kmesh%bz, Epren%nkpt, Kmesh%nbz, Cryst%nsym, &
287 : 1 : sppoldbl, Cryst%symafm, Cryst%symrel, timrev, comm, use_symrec=.False.)
288 : : end if
289 : :
290 : 23 : call timab(693,2,tsec) ! exc_haydock_driver(wo lf - that is, without local field
291 : 23 : call timab(694,1,tsec) ! exc_haydock_driver(apply
292 : :
293 : 23 : prefix = ""
294 [ + + ]: 55 : do itemp = 1, ntemp
295 : :
296 : 32 : call hexc%KS_BSt%copy(EPBSt)
297 : 32 : call hexc%QP_BSt%copy(EP_QPBSt)
298 : :
299 : : ! =================================================
300 : : ! == Calculate elphon vector in transition space ==
301 : : ! =================================================
302 [ + + ]: 32 : if (do_ep_renorm) then
303 : :
304 : : ! Will perform elphon renormalization for itemp
305 : :
306 : 10 : call int2char4(itemp,ts)
307 : 10 : prefix = TRIM("_T") // ts
308 : :
309 : : ! No scissor with KS Bands
310 : 10 : call renorm_bst(Epren, EPBSt, Cryst, itemp, do_lifetime=.TRUE.,do_check=.TRUE.)
311 : 10 : call renorm_bst(Epren, EP_QPBSt, Cryst, itemp, do_lifetime=.TRUE.,do_check=.FALSE.)
312 : :
313 [ + + ]: 20 : do isppol = 1, BSp%nsppol
314 [ + + ]: 980 : do ireh = 1, BSp%nreh(isppol)
315 : 960 : ic = BSp%Trans(ireh,isppol)%c
316 : 960 : iv = BSp%Trans(ireh,isppol)%v
317 : 960 : ik = BSp%Trans(ireh,isppol)%k ! In the full bz
318 : 960 : en = BSp%Trans(ireh,isppol)%en
319 : :
320 : 960 : ep_ik = bs2eph(ik,1)
321 : :
322 : : !TODO support multiple spins !
323 [ - + ]: 960 : if(ABS(en - (Epren%eigens(ic,ep_ik,isppol)-Epren%eigens(iv,ep_ik,isppol)+BSp%mbpt_sciss)) > tol3) then
324 : 0 : ABI_ERROR("Eigen from the transition does not correspond to the EP file !")
325 : : end if
326 : 960 : ep_renorms(ireh) = (Epren%renorms(1,ic,ik,isppol,itemp) - Epren%renorms(1,iv,ik,isppol,itemp))
327 : :
328 : : ! Add linewith
329 : 10 : if(do_ep_lifetime) then
330 : : ep_renorms(ireh) = ep_renorms(ireh) - j_dpc*(Epren%linewidth(1,ic,ik,isppol,itemp) +&
331 : 960 : Epren%linewidth(1,iv,ik,isppol,itemp))
332 : : end if
333 : :
334 : : end do
335 : : end do
336 : : end if
337 : :
338 : : ! =======================================================
339 : : ! === Make EPS RPA and GW without local-field effects ===
340 : : ! =======================================================
341 [ - + + - : 128 : ABI_MALLOC(eps_rpanlf,(BSp%nomega,BSp%nq))
- + + - +
- + - - +
- + ]
342 [ - + + - : 96 : ABI_MALLOC(dos_ks,(BSp%nomega))
+ - + - -
+ - + ]
343 [ - + + - : 96 : ABI_MALLOC(eps_gwnlf ,(BSp%nomega,BSp%nq))
- + + - -
+ - + ]
344 [ - + + - : 64 : ABI_MALLOC(dos_gw,(BSp%nomega))
- + - + ]
345 : :
346 : 32 : call wrtout(std_out," Calculating RPA NLF and QP NLF epsilon")
347 : :
348 : : call exc_eps_rpa(BSp%nbnds,BSp%lomo_spin,BSp%lomo_min,BSp%homo_spin,hexc%Kmesh,EPBSt,BSp%nq,nsppol,&
349 : 32 : opt_cvk,Cryst%ucvol,BSp%broad,BSp%nomega,BSp%omega,eps_rpanlf,dos_ks)
350 : :
351 : : call exc_eps_rpa(BSp%nbnds,BSp%lomo_spin,BSp%lomo_min,BSp%homo_spin,hexc%Kmesh,EP_QPBSt,BSp%nq,nsppol,&
352 : 32 : opt_cvk,Cryst%ucvol,Bsp%broad,BSp%nomega,BSp%omega,eps_gwnlf,dos_gw)
353 : :
354 [ + + ]: 32 : if (my_rank==master) then ! Only master works.
355 : : !
356 : : ! Master node writes final results on file.
357 : 28 : call exc_write_data(BSp,BS_files,"RPA_NLF_MDF",eps_rpanlf,prefix=prefix,dos=dos_ks)
358 : 28 : call exc_write_data(BSp,BS_files,"GW_NLF_MDF",eps_gwnlf,prefix=prefix,dos=dos_gw)
359 : :
360 : : ! Computing and writing tensor in files
361 : :
362 : : ! RPA_NLF
363 [ - + + - : 84 : ABI_MALLOC(tensor_cart_rpanlf,(BSp%nomega,6))
+ - + - -
+ - + ]
364 [ - + + - : 56 : ABI_MALLOC(tensor_red_rpanlf,(BSp%nomega,6))
- + - + ]
365 : :
366 : 28 : call wrtout(std_out," Calculating RPA NLF dielectric tensor")
367 : 28 : call haydock_mdf_to_tensor(BSp,Cryst,eps_rpanlf,tensor_cart_rpanlf, tensor_red_rpanlf, ierr)
368 : :
369 [ + - ]: 28 : if (ierr == 0) then
370 : : ! Writing tensor
371 : 28 : call exc_write_tensor(BSp,BS_files,"RPA_NLF_TSR_CART",tensor_cart_rpanlf)
372 : 28 : call exc_write_tensor(BSp,BS_files,"RPA_NLF_TSR_RED",tensor_red_rpanlf)
373 : : else
374 : : write(msg,'(3a)')&
375 : 0 : 'The RPA_NLF dielectric complex tensor cannot be computed',ch10,&
376 : 0 : 'There must be 6 different q-points in long wavelength limit (see gw_nqlwl)'
377 : 0 : ABI_COMMENT(msg)
378 : : end if
379 : :
380 [ - + ]: 28 : ABI_FREE(tensor_cart_rpanlf)
381 [ - + ]: 28 : ABI_FREE(tensor_red_rpanlf)
382 : :
383 : : ! GW_NLF
384 [ - + + - : 56 : ABI_MALLOC(tensor_cart_gwnlf,(BSp%nomega,6))
- + - + ]
385 [ - + + - : 56 : ABI_MALLOC(tensor_red_gwnlf,(BSp%nomega,6))
- + - + ]
386 : :
387 : 28 : call wrtout(std_out," Calculating GW NLF dielectric tensor")
388 : :
389 : 28 : call haydock_mdf_to_tensor(BSp,Cryst,eps_gwnlf,tensor_cart_gwnlf, tensor_red_gwnlf, ierr)
390 : :
391 [ + - ]: 28 : if(ierr == 0) then
392 : : ! Writing tensor
393 : 28 : call exc_write_tensor(BSp,BS_files,"GW_NLF_TSR_CART",tensor_cart_gwnlf)
394 : 28 : call exc_write_tensor(BSp,BS_files,"GW_NLF_TSR_RED",tensor_red_gwnlf)
395 : : else
396 : : write(msg,'(3a)')&
397 : 0 : 'The GW_NLF dielectric complex tensor cannot be computed',ch10,&
398 : 0 : 'There must be 6 different q-points in long wavelength limit (see gw_nqlwl)'
399 : 0 : ABI_COMMENT(msg)
400 : : end if
401 : :
402 [ - + ]: 28 : ABI_FREE(tensor_cart_gwnlf)
403 [ - + ]: 84 : ABI_FREE(tensor_red_gwnlf)
404 : :
405 : : !call wrtout(std_out," Checking Kramers Kronig on Excitonic Macroscopic Epsilon")
406 : : !call check_kramerskronig(BSp%nomega,REAL(BSp%omega),eps_exc(:,1))
407 : :
408 : : !call wrtout(std_out," Checking Kramers Kronig on RPA NLF Macroscopic Epsilon")
409 : : !call check_kramerskronig(BSp%nomega,REAL(BSp%omega),eps_rpanlf(:,1))
410 : :
411 : : !call wrtout(std_out," Checking Kramers Kronig on GW NLF Macroscopic Epsilon")
412 : : !call check_kramerskronig(BSp%nomega,REAL(BSp%omega),eps_gwnlf(:,1))
413 : :
414 : : !call wrtout(std_out," Checking f-sum rule on Excitonic Macroscopic Epsilon")
415 : :
416 : : !if (BSp%exchange_term>0) then
417 : : ! ABI_COMMENT(' f-sum rule should be checked without LF')
418 : : !end if
419 : : !call check_fsumrule(BSp%nomega,REAL(BSp%omega),AIMAG(eps_exc(:,1)),drude_plsmf)
420 : :
421 : : !call wrtout(std_out," Checking f-sum rule on RPA NLF Macroscopic Epsilon")
422 : : !call check_fsumrule(BSp%nomega,REAL(BSp%omega),AIMAG(eps_rpanlf(:,1)),drude_plsmf)
423 : :
424 : : !call wrtout(std_out," Checking f-sum rule on GW NLF Macroscopic Epsilon")
425 : : !call check_fsumrule(BSp%nomega,REAL(BSp%omega),AIMAG(eps_gwnlf(:,1)),drude_plsmf)
426 : : end if ! my_rank==master
427 : :
428 : : !call xmpi_barrier(comm)
429 : : !
430 : : ! The ket for the approximated DOS.
431 : : if (prtdos) then
432 : : ABI_WARNING("Calculating DOS with Haydock method")
433 : : ABI_CHECK(BSp%use_coupling==0,"DOS with coupling not coded")
434 : : iq = BSp%nq + 1
435 : : if (my_rank==master) then
436 : : !call random_seed()
437 : : do itt=1,SUM(Bsp%nreh)
438 : : call RANDOM_NUMBER(rand_phi)
439 : : rand_phi = two_pi*rand_phi
440 : : kets(itt,iq) = CMPLX( COS(rand_phi), SIN(rand_phi) )
441 : : end do
442 : : ! Normalize the vector.
443 : : !norm = SQRT( DOT_PRODUCT(kets(:,iq), kets(:,iq)) )
444 : : !kets(:,iq) = kets(:,iq)/norm
445 : : end if
446 : : call xmpi_bcast(kets(:,iq),master,comm,ierr)
447 : : end if
448 : :
449 [ - + + - : 128 : ABI_MALLOC(green,(BSp%nomega,nkets))
- + + - +
- + - - +
- + ]
450 : :
451 [ + - ]: 32 : if (BSp%use_coupling==0) then
452 [ + + ]: 32 : if(do_ep_renorm) then
453 : 10 : call haydock_bilanczos(BSp,BS_files,Cryst,Hdr_bse,hexc,hexc_i,hsize,hexc%my_t1,hexc%my_t2,nkets,kets,ep_renorms,green,comm)
454 : : else
455 : : !YG2014
456 : 22 : call haydock_herm(BSp,BS_files,Cryst,Hdr_bse,hexc%my_t1,hexc%my_t2,nkets,kets,green,hexc,hexc_i,comm)
457 : : end if
458 : : else
459 [ # # ]: 0 : if (BSp%use_interp) then
460 : 0 : ABI_ERROR("BSE Interpolation with coupling is not supported")
461 : : else
462 : 0 : call haydock_psherm(BSp,BS_files,Cryst,Hdr_bse,hexc,hexc_i,hsize,my_t1,my_t2,nkets,kets,green,comm)
463 : : end if
464 : : end if
465 : :
466 : : ! Add 1 to have the real part right.
467 [ + + + + ]: 87374 : green = one + green
468 : :
469 [ + + ]: 32 : if (my_rank==master) then ! Master writes the final results.
470 : :
471 : : if (prtdos) then
472 : : ABI_MALLOC(dos,(BSp%nomega))
473 : : dos = -AIMAG(green(:,BSp%nq+1))
474 : : call exc_write_data(BSp,BS_files,"EXC_MDF",green,prefix=prefix,dos=dos)
475 : : ABI_FREE(dos)
476 : : else
477 : 28 : call exc_write_data(BSp,BS_files,"EXC_MDF",green,prefix=prefix)
478 : : end if
479 : :
480 : : ! =========================
481 : : ! === Write out Epsilon ===
482 : : ! =========================
483 : :
484 [ - + + - : 84 : ABI_MALLOC(tensor_cart,(BSp%nomega,6))
+ - + - -
+ - + ]
485 [ - + + - : 56 : ABI_MALLOC(tensor_red,(BSp%nomega,6))
- + - + ]
486 : :
487 : 28 : call wrtout(std_out," Calculating EXC dielectric tensor")
488 : 28 : call haydock_mdf_to_tensor(BSp,Cryst,green,tensor_cart,tensor_red,ierr)
489 : :
490 [ + - ]: 28 : if (ierr == 0) then
491 : : ! Writing tensor
492 : 28 : call exc_write_tensor(BSp,BS_files,"EXC_TSR_CART",tensor_cart)
493 : 28 : call exc_write_tensor(BSp,BS_files,"EXC_TSR_RED",tensor_red)
494 : : else
495 : : write(msg,'(3a)')&
496 : 0 : 'The EXC dielectric complex tensor cannot be computed',ch10,&
497 : 0 : 'There must be 6 different q-points in long wavelength limit (see gw_nqlwl)'
498 : 0 : ABI_COMMENT(msg)
499 : : end if
500 : :
501 [ - + ]: 28 : ABI_FREE(tensor_cart)
502 [ - + ]: 28 : ABI_FREE(tensor_red)
503 : : !
504 : : ! This part will be removed when fldiff will be able to compare two mdf files.
505 : 28 : write(ab_out,*)" "
506 : 28 : write(ab_out,*)"Macroscopic dielectric function:"
507 : 28 : write(ab_out,*)"omega [eV] <KS_RPA_nlf> <GW_RPA_nlf> <BSE> "
508 [ + + ]: 308 : do io=1,MIN(BSp%nomega,10)
509 : 280 : omegaev = REAL(BSp%omega(io))*Ha_eV
510 [ + + ]: 1960 : ks_avg = SUM( eps_rpanlf(io,:)) / Bsp%nq
511 [ + + ]: 1960 : gw_avg = SUM( eps_gwnlf (io,:)) / Bsp%nq
512 [ + + ]: 1960 : exc_avg = SUM( green (io,:)) / BSp%nq
513 : 308 : write(ab_out,'(7f9.4)')omegaev,ks_avg,gw_avg,exc_avg
514 : : end do
515 : 28 : write(ab_out,*)" "
516 : :
517 : : ! Write MDF file with the final results.
518 : : ! FIXME: It won't work if prtdos == True
519 [ + - ]: 28 : path = strcat(BS_files%out_basename,strcat(prefix,"_MDF.nc"))
520 [ - + ]: 28 : NCF_CHECK(nctk_open_create(ncid, path, xmpi_comm_self))
521 [ - + ]: 28 : NCF_CHECK(cryst%ncwrite(ncid))
522 [ - + ]: 28 : NCF_CHECK(QP_bst%ncwrite(ncid))
523 : 28 : call mdfs_ncwrite(ncid, Bsp, green, eps_rpanlf, eps_gwnlf)
524 [ - + ]: 56 : NCF_CHECK(nf90_close(ncid))
525 : : end if
526 : :
527 [ - + ]: 32 : ABI_FREE(green)
528 [ - + ]: 32 : ABI_FREE(eps_rpanlf)
529 [ - + ]: 32 : ABI_FREE(eps_gwnlf)
530 [ - + ]: 32 : ABI_FREE(dos_ks)
531 [ - + ]: 32 : ABI_FREE(dos_gw)
532 : :
533 : 32 : call EPBSt%free()
534 : 55 : call EP_QPBst%free()
535 : : end do ! itemp loop
536 : :
537 [ - + ]: 23 : ABI_FREE(opt_cvk)
538 [ - + ]: 23 : ABI_FREE(kets)
539 : :
540 : 23 : call timab(694,2,tsec) ! exc_haydock_driver(apply
541 : 23 : call timab(695,1,tsec) ! exc_haydock_driver(end)
542 : :
543 : : !YG2014
544 : 23 : call hexc%free()
545 : 23 : call hexc_i%free()
546 : :
547 [ + + ]: 23 : ABI_SFREE(ep_renorms)
548 [ + + ]: 23 : ABI_SFREE(bs2eph)
549 : :
550 : 23 : call timab(695,2,tsec) ! exc_haydock_driver(end)
551 : 23 : call timab(690,2,tsec) ! exc_haydock_driver
552 : :
553 [ - + - + : 23 : end subroutine exc_haydock_driver
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + ]
554 : : !!***
555 : :
556 : : !----------------------------------------------------------------------
557 : :
558 : : !!****f* m_haydock/haydock_herm
559 : : !! NAME
560 : : !! haydock_herm
561 : : !!
562 : : !! FUNCTION
563 : : !! Reads the excitonic Hamiltonian from file and construct the Lanczos set of vectors
564 : : !! by iterative matrix-vector multiplications.
565 : : !!
566 : : !! INPUTS
567 : : !! BSp<excparam>=Parameters for the Bethe-Salpeter calculation.
568 : : !! BS_files<excparam>=Files associated to the bethe_salpeter code.
569 : : !! Cryst<crystal_t>=Info on the crystalline structure.
570 : : !! Pawtab(Cryst%ntypat*usepaw)<pawtab_type>=PAW tabulated starting data.
571 : : !! hize=Size of the excitonic matrix.
572 : : !! my_t1,my_t2=First and last columns treated by this node.
573 : : !! nkets=Number of starting vectors for Haydock method.
574 : : !! kets(hsize,nkets)=The kets in the eh representation.
575 : : !! comm=MPI communicator.
576 : : !!
577 : : !! OUTPUT
578 : : !! green(BSp%nomega,nkets)=
579 : : !!
580 : : !! SOURCE
581 : :
582 : 22 : subroutine haydock_herm(BSp,BS_files,Cryst,Hdr_bse,my_t1,my_t2,nkets,kets,green,hexc,hexc_i,comm)
583 : :
584 : : !Arguments ------------------------------------
585 : : !scalars
586 : : integer,intent(in) :: my_t1,my_t2,nkets,comm
587 : : type(crystal_t),intent(in) :: Cryst
588 : : type(excparam),intent(in) :: BSp
589 : : type(excfiles),intent(in) :: BS_files
590 : : type(Hdr_type),intent(in) :: Hdr_bse
591 : : type(hexc_t),intent(inout) :: hexc
592 : : type(hexc_interp_t),intent(inout) :: hexc_i
593 : : !arrays
594 : : complex(dp),intent(out) :: green(BSp%nomega,nkets)
595 : : complex(dp),intent(in) :: kets(hexc%hsize,nkets)
596 : :
597 : : !Local variables ------------------------------
598 : : !scalars
599 : : integer,parameter :: master=0
600 : : integer :: inn,nproc,my_rank,ierr
601 : : integer :: niter_file,niter_max,niter_done,nsppol,iq,my_nt,term_type,n_all_omegas
602 : : real(dp) :: norm,nfact
603 : : logical :: can_restart,is_converged
604 : : complex(dp) :: factor
605 : : character(len=500) :: msg
606 : : character(len=fnlen), parameter :: tag_file="_HAYDR_SAVE"
607 : : character(len=fnlen) :: restart_file,out_file
608 : 22 : type(haydock_type) :: haydock_file
609 : : !arrays
610 : 22 : real(dp),allocatable :: bb_file(:)
611 : 22 : real(dp),allocatable :: bb(:)
612 : 22 : complex(dp),allocatable :: aa(:),phi_nm1(:),phi_n(:),hphi_n(:),hphi_nm1(:)
613 : 22 : complex(dp),allocatable :: aa_file(:),phi_n_file(:),phi_nm1_file(:)
614 : 22 : complex(dp),allocatable :: ket0(:),all_omegas(:),green_temp(:,:)
615 : : ! complex(dp),allocatable :: diag_dense(:)
616 : : logical :: check(2)
617 : : !************************************************************************
618 : :
619 [ - + ]: 22 : ABI_CHECK(Bsp%nsppol==1,"nsppol > 1 not implemented yet")
620 : :
621 : 22 : nproc = xmpi_comm_size(comm); my_rank= xmpi_comm_rank(comm)
622 : 22 : nsppol = Hdr_bse%nsppol
623 : :
624 [ + + ]: 22 : if (BSp%use_interp) then
625 : 4 : ABI_COMMENT("No parallelization in Interpolation")
626 [ + + ]: 8 : my_nt = SUM(Bsp%nreh_interp)
627 : : else
628 : 18 : my_nt = my_t2-my_t1+1
629 : : end if
630 : :
631 [ - + ]: 22 : ABI_CHECK(my_nt>0,"One of the processors has zero columns")
632 : :
633 : 22 : write(msg,'(a,i0)')' Haydock algorithm with MAX number of iterations: ',BSp%niter
634 : 22 : call wrtout(std_out,msg)
635 : : !
636 : : ! Select the terminator for the continued fraction.
637 [ + + ]: 22 : term_type=0; if (Bsp%hayd_term>0) term_type=1
638 : 22 : call wrtout(std_out,sjoin("Using terminator type: ",itoa(term_type)))
639 : : !
640 : : ! Check for presence of the restart file.
641 : 22 : can_restart=.FALSE.
642 [ - + ]: 22 : if ( BS_files%in_haydock_basename /= BSE_NOFILE) then
643 [ # # ]: 0 : restart_file = TRIM(BS_files%in_haydock_basename)//TRIM(tag_file)
644 [ # # ]: 0 : if (file_exists(restart_file) ) then
645 : 0 : can_restart=.TRUE.
646 : 0 : msg = " Restarting Haydock calculation from file: "//TRIM(restart_file)
647 : 0 : call wrtout(std_out,msg)
648 : 0 : call wrtout(ab_out,msg)
649 : : else
650 : 0 : can_restart=.FALSE.
651 [ # # ]: 0 : call wrtout(ab_out," WARNING: cannot find restart file: "//TRIM(restart_file))
652 : : end if
653 : : end if
654 : 0 : ABI_CHECK(.not.can_restart,"restart not yet implemented")
655 : :
656 : : ! Open the file and write basic dimensions and info.
657 [ + + ]: 22 : if (my_rank==master) then
658 [ + - ]: 18 : out_file = TRIM(BS_files%out_basename)//TRIM(tag_file)
659 : 18 : call haydock_file%open(out_file)
660 : 18 : haydock_file%hsize = hexc%hsize
661 : 18 : haydock_file%use_coupling = Bsp%use_coupling
662 : 18 : haydock_file%op = BSE_HAYD_IMEPS
663 : 18 : haydock_file%nq = nkets
664 : 18 : haydock_file%broad = Bsp%broad
665 : 18 : call haydock_file%write_dim()
666 : : end if
667 : :
668 : : !
669 : : ! Calculate green(w) for the different starting points.
670 [ + + + + ]: 27244 : green=czero
671 [ + + ]: 154 : do iq=1,nkets
672 [ - + + - : 396 : ABI_MALLOC(ket0,(hexc%hsize))
+ - + - -
+ - + ]
673 [ + - - + : 63816 : ket0=kets(:,iq)
- - + + ]
674 : :
675 : 132 : niter_file=0
676 [ - + ]: 132 : if (can_restart) then
677 : : call haydock_restart(BSp,restart_file,BSE_HAYD_IMEPS,iq,hexc%hsize,&
678 [ # # # # : 0 : niter_file,aa_file,bb_file,phi_nm1_file,phi_n_file,comm)
# # # # ]
679 : : end if
680 : : !
681 : : ! For n>1, we have:
682 : : ! 1) a_n = <n|H|n>
683 : : ! 2) b_n = || H|n> - a_n|n> -b_{n-1}|n-1> ||
684 : : ! 3) |n+1> = [H|n> -a_n|n> -b_{n-1}|n-1>]/b_n
685 : : !
686 : : ! The sequences starts with |1> normalized to 1 and b_0 =0, therefore:
687 : : ! a_1 = <1|H|1>
688 : : ! b_1 = || H|1> - a_1|1> ||
689 : : ! |2> = [H|1> - a_1|1>]/b_1
690 : : !
691 [ - + + - : 396 : ABI_MALLOC(hphi_n,(hexc%hsize))
+ - + - -
+ - + -
+ ]
692 [ - + + - : 264 : ABI_MALLOC(hphi_nm1,(hexc%hsize))
- + - + -
+ ]
693 [ - + + - : 396 : ABI_MALLOC(phi_nm1,(my_nt))
+ - + - -
+ - + ]
694 [ - + + - : 264 : ABI_MALLOC(phi_n,(my_nt))
- + - + ]
695 : :
696 : 132 : niter_max = niter_file + Bsp%niter
697 [ - + + - : 396 : ABI_MALLOC(aa,(niter_max))
+ - + - -
+ - + ]
698 [ - + + - : 396 : ABI_MALLOC(bb,(niter_max))
+ - + - -
+ - + ]
699 [ + + + + ]: 34212 : aa=czero; bb=zero
700 : :
701 [ + - ]: 132 : if (niter_file==0) then ! Calculation from scratch.
702 [ + - - + : 54600 : phi_nm1=ket0(my_t1:my_t2) ! Select the slice treated by this node.
- - + + ]
703 : 132 : norm = DZNRM2(hexc%hsize,ket0,1) ! Normalization
704 [ + + ]: 54468 : phi_nm1=phi_nm1/norm
705 : :
706 : 132 : call hexc%matmul_tda(hexc_i,phi_nm1,hphi_n)
707 : :
708 : 132 : aa(1)=xdotc(my_nt,phi_nm1,1,hphi_n(my_t1:),1)
709 : 132 : call xmpi_sum(aa(1:1),comm,ierr)
710 : :
711 [ + - - + : 54600 : phi_n = hphi_n(my_t1:my_t2) - aa(1)*phi_nm1
- - + + ]
712 : :
713 : 132 : bb(1) = xdotc(my_nt,phi_n,1,phi_n,1)
714 : 132 : call xmpi_sum(bb(1:1),comm,ierr)
715 : 132 : bb(1) = SQRT(bb(1))
716 : :
717 [ + + ]: 54468 : phi_n = phi_n/bb(1)
718 : 132 : niter_done=1
719 : :
720 : : else ! Use the previous a and b.
721 : 0 : niter_done=niter_file
722 [ # # ]: 0 : aa(1:niter_done) = aa_file
723 [ # # ]: 0 : bb(1:niter_done) = bb_file
724 [ # # # # : 0 : phi_nm1=phi_nm1_file(my_t1:my_t2) ! Select the slice treated by this node.
# # # # ]
725 [ # # # # : 0 : phi_n =phi_n_file (my_t1:my_t2)
# # # # ]
726 : : end if
727 : :
728 [ - + ]: 132 : if (can_restart) then
729 [ # # ]: 0 : ABI_FREE(aa_file)
730 [ # # ]: 0 : ABI_FREE(bb_file)
731 [ # # ]: 0 : ABI_FREE(phi_nm1_file)
732 [ # # ]: 0 : ABI_FREE(phi_n_file)
733 : : end if
734 : :
735 : : ! Multiplicative factor (k-point sampling and unit cell volume)
736 : : ! TODO be careful with the spin here
737 : : ! TODO four_pi comes from the coulomb term 1/|q| is already included in the
738 : : ! oscillators hence the present approach wont work if a cutoff interaction is used.
739 : 132 : nfact = -four_pi/(Cryst%ucvol*hexc%nbz)
740 [ + - ]: 132 : if (nsppol==1) nfact=two*nfact
741 : :
742 : 132 : factor = nfact*(DZNRM2(hexc%hsize,ket0,1)**2)
743 : :
744 : : ! Which quantity should be checked for convergence?
745 : 132 : check = (/.TRUE.,.TRUE./)
746 [ - + ]: 132 : if (ABS(Bsp%haydock_tol(2)-one)<tol6) check = (/.TRUE. ,.FALSE./)
747 [ - + ]: 132 : if (ABS(Bsp%haydock_tol(2)-two)<tol6) check = (/.FALSE.,.TRUE./)
748 : :
749 : : ! Create new frequencies "mirror" in negative range to add
750 : : ! their contributions. Can be improved by computing only once
751 : : ! zero frequency, but loosing clearness
752 : 132 : n_all_omegas = 2*BSp%nomega
753 : :
754 [ - + + - : 396 : ABI_MALLOC(all_omegas,(n_all_omegas))
+ - + - -
+ - + ]
755 : : ! Put all omegas with frequency > 0 in table
756 [ + + ]: 27222 : all_omegas(BSp%nomega+1:n_all_omegas) = BSp%omega
757 : : ! Put all omegas with frequency < 0
758 : : ! Warning, the broadening must be kept positive
759 [ + + ]: 27222 : all_omegas(1:BSp%nomega) = -DBLE(BSp%omega(BSp%nomega:1:-1)) + j_dpc*AIMAG(BSp%omega(BSp%nomega:1:-1))
760 : :
761 [ - + + - : 528 : ABI_MALLOC(green_temp,(n_all_omegas,nkets))
- + + - +
- + - - +
- + ]
762 : :
763 : : call haydock_herm_algo(niter_done,niter_max,n_all_omegas,all_omegas,BSp%haydock_tol(1),check,&
764 : : my_t1,my_t2,factor,term_type,aa,bb,phi_nm1,phi_n,&
765 : : green_temp(:,iq),inn,is_converged,&
766 : 132 : hexc, hexc_i, comm)
767 : :
768 : : ! Computing result from two ranges of frequencies
769 : : ! The real part is added, the imaginary part is substracted
770 [ + + ]: 27222 : green(:,iq) = green_temp(BSp%nomega+1:n_all_omegas,iq)+CONJG(green_temp(BSp%nomega:1:-1,iq))
771 : :
772 [ - + ]: 132 : ABI_FREE(all_omegas)
773 [ - + ]: 132 : ABI_FREE(green_temp)
774 : : !
775 : : ! Save the a"s and the b"s for possible restarting.
776 : : ! 1) Info on the Q.
777 : : ! 2) Number of iterations performed.
778 : : ! 3) do iter=1,niter_performed
779 : : ! aa(iter),bb(iter)
780 : : ! end do
781 : : ! 4) |n-1>
782 : : ! |n>
783 : : !
784 [ + + ]: 63684 : hphi_nm1 = czero
785 [ + + ]: 54468 : hphi_nm1(my_t1:my_t2) = phi_nm1
786 : 132 : call xmpi_sum_master(hphi_nm1,master,comm,ierr)
787 : :
788 [ + + ]: 63684 : hphi_n = czero
789 [ + + ]: 54468 : hphi_n(my_t1:my_t2) = phi_n
790 : 132 : call xmpi_sum_master(hphi_n,master,comm,ierr)
791 : :
792 [ + + ]: 132 : if (my_rank==master) then
793 : : ! Write data for restarting
794 : 108 : call haydock_file%write(hexc%hsize, Bsp%q(:,iq), aa, bb, hphi_n, hphi_nm1, MIN(inn,niter_max), factor)
795 : : end if
796 : :
797 [ - + ]: 132 : ABI_FREE(hphi_n)
798 [ - + ]: 132 : ABI_FREE(hphi_nm1)
799 [ - + ]: 132 : ABI_FREE(phi_nm1)
800 [ - + ]: 132 : ABI_FREE(phi_n)
801 [ - + ]: 132 : ABI_FREE(aa)
802 [ - + ]: 132 : ABI_FREE(bb)
803 [ - + ]: 154 : ABI_FREE(ket0)
804 : : end do ! iq
805 : :
806 [ + + ]: 22 : if (my_rank==master) call haydock_file%close()
807 : :
808 : 22 : call xmpi_barrier(comm)
809 : :
810 [ - + - + : 22 : end subroutine haydock_herm
- + - + -
+ - + - +
- + ]
811 : : !!***
812 : :
813 : : !----------------------------------------------------------------------
814 : :
815 : : !!****f* m_haydock/haydock_herm_algo
816 : : !! NAME
817 : : !! haydock_herm_algo
818 : : !!
819 : : !! FUNCTION
820 : : !!
821 : : !! INPUTS
822 : : !! niter_done=Number of iterations already performed (0 if the run starts from scratch).
823 : : !! niter_max=Max number of iterations. Always > niter_done
824 : : !! nomega=Number of Frequency points for the evaluation of the matrix element.
825 : : !! omega(nomega)=Frequency set (imaginary part is already included).
826 : : !! tol_iter=Tolerance used to stop the algorithm.
827 : : !! check(2)=Logical flags to specify where both the real and the imaginary part of the
828 : : !! matrix elements of the Green functions have to be checked for convergence.
829 : : !! hsize=Size of the blocks.
830 : : !! my_t1,my_t2=Indices of the first and last column stored treated by this done.
831 : : !! term_type=0 if no terminator is used, 1 otherwise.
832 : : !! hmat(hsize,my_t1:my_t2)=The columns of the block.
833 : : !! factor
834 : : !! ntrans = Number of transitions
835 : : !! corresp = mapping between coarse points and neighbours
836 : : !! overlaps = overlaps of wavefunctions between dense k-point coarse neighbours and bands
837 : : !! comm=MPI communicator.
838 : : !!
839 : : !! OUTPUT
840 : : !! green(nomega)=Output matrix elements.
841 : : !! inn=Last iteration performed.
842 : : !! is_converged=.TRUE. of the algorithm converged.
843 : : !!
844 : : !! SIDE EFFECTS
845 : : !! phi_nm1(my_t2-my_t1+1), phi_n(my_t2-my_t1+1)
846 : : !! input: vectors used to initialize the iteration
847 : : !! output: the vectors obtained in the last iteration
848 : : !! aa(niter_max) and bb(niter_max)
849 : : !! if niter_done>0: aa(1:niter_done), bb(1:niter_done) store the coefficients of the previous run.
850 : : !! when the routine returns aa(1:inn) and bb(1:inn) contain the matrix elements of the tridiagonal form.
851 : : !!
852 : : !! SOURCE
853 : :
854 : 132 : subroutine haydock_herm_algo(niter_done,niter_max,nomega,omega,tol_iter,check,&
855 : 132 : my_t1,my_t2,factor,term_type,aa,bb,phi_nm1,phi_n,&
856 : 132 : green,inn,is_converged, hexc, hexc_i, comm)
857 : :
858 : : !Arguments ------------------------------------
859 : : !scalars
860 : : integer,intent(in) :: niter_max,niter_done,nomega
861 : : integer,intent(in) :: my_t1,my_t2,term_type
862 : : integer,intent(in) :: comm
863 : : integer,intent(out) :: inn
864 : : logical,intent(out) :: is_converged
865 : : real(dp),intent(in) :: tol_iter
866 : : complex(dp),intent(in) :: factor
867 : : type(hexc_t),intent(in) :: hexc
868 : : type(hexc_interp_t),intent(in) :: hexc_i
869 : : !arrays
870 : : real(dp),intent(inout) :: bb(niter_max)
871 : : complex(dp),intent(out) :: green(nomega)
872 : : complex(dp),intent(in) :: omega(nomega)
873 : : complex(dp),intent(inout) :: aa(niter_max)
874 : : complex(dp),intent(inout) :: phi_nm1(my_t2-my_t1+1)
875 : : complex(dp),intent(inout) :: phi_n (my_t2-my_t1+1)
876 : : logical,intent(in) :: check(2)
877 : :
878 : : !Local variables ------------------------------
879 : : !scalars
880 : : integer :: ierr,my_nt,niter_min,nconv
881 : : character(len=500) :: msg
882 : : logical,parameter :: force_real=.TRUE.
883 : : !arrays
884 : : integer :: units(2)
885 : 264 : real(dp) :: abs_err(nomega,2) !,rel_err(nomega,2)
886 : 132 : complex(dp),allocatable :: oldg(:),newg(:), phi_np1(:),hphi_n(:),cfact(:)
887 : : logical :: test(2)
888 : : !************************************************************************
889 : :
890 : : ! The sequences starts with |1> normalized to 1 and b_0 =0, therefore:
891 : : ! a_1 = <1|H|1>
892 : : ! b_1 = || H|1> - a_1|1> ||
893 : : ! |2> = [H|1> - a_1|1>]/b_1
894 : : !
895 : : ! For n>1 we have
896 : : ! 1) a_n = <n|H|n>
897 : : ! 2) b_n = || H|n> - a_n|n> -b_{n-1}|n-1> ||
898 : : ! 3) |n+1> = [H|n> -a_n|n> -b_{n-1}|n-1>]/b_n
899 : : !
900 : 132 : my_nt = my_t2-my_t1+1
901 [ + + ]: 396 : units = [std_out, ab_out]
902 : :
903 [ - + + - : 396 : ABI_MALLOC_OR_DIE(hphi_n,(hexc%hsize), ierr)
+ - + - +
- - + -
+ ]
904 [ - + + - : 396 : ABI_MALLOC(phi_np1,(my_nt))
+ - + - -
+ - + ]
905 [ - + + - : 54576 : ABI_CALLOC(oldg,(nomega))
+ - + - -
+ - + +
+ ]
906 [ - + + - : 54444 : ABI_CALLOC(newg,(nomega))
- + - + +
+ ]
907 [ - + + - : 54444 : ABI_CALLOC(cfact,(nomega))
- + - + +
+ ]
908 : :
909 : 132 : nconv=0
910 [ + - ]: 9661 : do inn=niter_done+1,niter_max
911 : :
912 : : !YG2014
913 : 9661 : call hexc%matmul_tda(hexc_i,phi_n,hphi_n)
914 : :
915 : 9661 : aa(inn) = xdotc(my_nt,phi_n,1,hphi_n(my_t1:),1)
916 : 9661 : call xmpi_sum(aa(inn:inn),comm,ierr)
917 : 9661 : if (force_real) aa(inn) = DBLE(aa(inn)) ! Matrix is Hermitian.
918 : :
919 : : ! |n+1> = H|n> - A(n)|n> - B(n-1)|n-1>
920 [ + - - + : 5417978 : phi_np1 = hphi_n(my_t1:my_t2) - aa(inn)*phi_n - bb(inn-1)*phi_nm1
- - + + ]
921 : :
922 : 9661 : bb(inn) = xdotc(my_nt,phi_np1,1,phi_np1,1)
923 : 9661 : call xmpi_sum(bb(inn),comm,ierr)
924 : 9661 : bb(inn) = SQRT(bb(inn))
925 : :
926 [ + + ]: 5408317 : phi_np1 = phi_np1/bb(inn)
927 [ + + ]: 5408317 : phi_nm1 = phi_n
928 [ + + ]: 5408317 : phi_n = phi_np1
929 : :
930 : 9661 : write(msg,'(a,i0,a,3es12.4)')' Iteration number ',inn,', b_i RE(a_i) IM(a_i) ',bb(inn),REAL(aa(inn)),AIMAG(aa(inn))
931 : 9661 : call wrtout(std_out, msg)
932 : :
933 : 9661 : call continued_fract(inn,term_type,aa,bb,nomega,omega,cfact)
934 [ - + + + ]: 2846830 : newg= factor*cfact
935 : :
936 : : ! Avoid spurious convergence.
937 [ - + ]: 9661 : niter_min=4; if (niter_done>1) niter_min=niter_done+1
938 [ + + ]: 9661 : if (inn>niter_min) then
939 [ + + ]: 27795 : test=.TRUE.
940 [ + + ]: 2674233 : abs_err(:,1) = ABS(DBLE (newg-oldg))
941 [ + + ]: 2674233 : abs_err(:,2) = ABS(AIMAG(newg-oldg))
942 : : !
943 [ + + ]: 9265 : if (tol_iter>zero) then
944 : : ! Test on the L1 norm.
945 [ + - + + : 3938750 : if (check(1)) test(1) = SUM(abs_err(:,1)) < tol_iter*SUM(ABS(DBLE (newg)))
+ + ]
946 [ + - + + : 3938750 : if (check(2)) test(2) = SUM(abs_err(:,2)) < tol_iter*SUM(ABS(AIMAG(newg)))
+ + ]
947 : : else
948 : : ! Stringent test for each point.
949 [ + - + + : 590552 : if (check(1)) test(1) = ALL( abs_err(:,1) < -tol_iter*ABS(DBLE (newg)))
+ + ]
950 [ + - + + : 593729 : if (check(2)) test(2) = ALL( abs_err(:,2) < -tol_iter*ABS(AIMAG(newg)))
+ + ]
951 : : end if
952 : : !
953 [ + + + + : 11028 : if (ALL(test)) then
+ + ]
954 : 371 : nconv = nconv+1
955 : : else
956 : : nconv = 0
957 : : end if
958 [ + + ]: 371 : if (nconv==2) then
959 [ + + ]: 132 : if(inn<100)then
960 : : write(msg,'(a,es10.2,a)')&
961 : 102 : " >>> Haydock algorithm converged twice within haydock_tol= ",tol_iter," after less than 100 iterations."
962 : : else
963 : : write(msg,'(a,es10.2,a,i0,a)')&
964 : 30 : " >>> Haydock algorithm converged twice within haydock_tol= ",tol_iter," after ",inn," iterations."
965 : : endif
966 : 132 : call wrtout(units, msg)
967 : 132 : EXIT
968 : : end if
969 : : end if
970 : :
971 [ - + + + ]: 2802047 : oldg = newg
972 : : end do ! inn
973 : :
974 [ + + ]: 54312 : green = newg
975 [ - + ]: 132 : if (nconv/=2) then
976 : : write(msg,'(a,es10.2,a,i0,a)')&
977 : 0 : " WARNING: Haydock algorithm did not converge within ",tol_iter," after ",niter_max," iterations."
978 : 0 : call wrtout(units, msg)
979 : : end if
980 : :
981 : 132 : is_converged = (nconv==2)
982 : :
983 [ - + ]: 132 : ABI_FREE(oldg)
984 [ - + ]: 132 : ABI_FREE(newg)
985 [ - + ]: 132 : ABI_FREE(cfact)
986 [ - + ]: 132 : ABI_FREE(hphi_n)
987 [ - + ]: 132 : ABI_FREE(phi_np1)
988 : :
989 : 132 : end subroutine haydock_herm_algo
990 : : !!***
991 : :
992 : : !----------------------------------------------------------------------
993 : :
994 : : !!****f* m_haydock/haydock_restart
995 : : !! NAME
996 : : !! haydock_restart
997 : : !!
998 : : !! FUNCTION
999 : : !! Restart the Haydock method from file reading the data produced in a previous run.
1000 : : !!
1001 : : !! INPUTS
1002 : : !! BSp<type(excparam)>=Parameters defining the Bethe-Salpeter calculation.
1003 : : !! omega(BSp%nomega)=Frequency mesh for the macroscopic dielectric function (broadening is already included).
1004 : : !! iq_search=The index of the q-point to be searched.
1005 : : !! hsize
1006 : : !! comm=MPI communicator.
1007 : : !! nsppol
1008 : : !! restart_file
1009 : : !!
1010 : : !! OUTPUT
1011 : : !! niter_file=Number of iterations already performed. 0 to signal that an error occurred during the reading
1012 : : !! bb_file(:)
1013 : : !! aa_file(:)
1014 : : !! phi_n_file(:)
1015 : : !! phi_nm1_file(:)
1016 : : !!
1017 : : !! SOURCE
1018 : :
1019 : 0 : subroutine haydock_restart(BSp,restart_file,ftype,iq_search,hsize,niter_file,aa_file,bb_file,phi_nm1_file,phi_n_file,comm)
1020 : :
1021 : : !Arguments ------------------------------------
1022 : : !scalars
1023 : : integer,intent(in) :: comm,hsize,iq_search,ftype
1024 : : integer,intent(out) :: niter_file
1025 : : character(len=*),intent(in) :: restart_file
1026 : : type(excparam),intent(in) :: BSp
1027 : : !arrays
1028 : : real(dp),allocatable,intent(out) :: bb_file(:)
1029 : : complex(dp),allocatable,intent(out) :: aa_file(:),phi_n_file(:),phi_nm1_file(:)
1030 : :
1031 : : !Local variables ------------------------------
1032 : : !scalars
1033 : : integer,parameter :: master=0
1034 : : integer :: nproc,my_rank,ierr,op_file, units(2)
1035 : : integer :: hsize_file,use_coupling_file
1036 : : complex(dp) :: factor_file
1037 : : character(len=500) :: msg
1038 : 0 : type(haydock_type) :: haydock_file
1039 : : !************************************************************************
1040 : :
1041 : 0 : nproc = xmpi_comm_size(comm); my_rank= xmpi_comm_rank(comm)
1042 [ # # ]: 0 : units = [std_out, ab_out]
1043 : :
1044 [ # # ]: 0 : if (my_rank==master) then
1045 : 0 : call haydock_file%open(restart_file)
1046 : 0 : call haydock_file%read_dim()
1047 : :
1048 [ # # ]: 0 : if (haydock_file%op/=ftype) then
1049 : 0 : write(msg,"(2(a,i0))")" Expecting restart file with filetype: ",ftype," but found ",op_file
1050 : 0 : ABI_ERROR(msg)
1051 : : end if
1052 : :
1053 [ # # ]: 0 : if (haydock_file%hsize/=hsize) then
1054 : : write(msg,"(2(a,i0))")&
1055 : 0 : " Rank of H_exc read from file: ",hsize_file," differs from the one used in this run: ",hsize
1056 : 0 : ABI_ERROR(msg)
1057 : : end if
1058 : :
1059 [ # # ]: 0 : if (haydock_file%use_coupling /= BSp%use_coupling) then
1060 : : write(msg,'(2(a,i0))')&
1061 : 0 : " use_coupling_file: ",use_coupling_file," differs from input file value: ",BSp%use_coupling
1062 : 0 : ABI_ERROR(msg)
1063 : : end if
1064 : :
1065 : : call haydock_file%read(Bsp%q(:,iq_search), aa_file, bb_file, &
1066 [ # # # # : 0 : phi_n_file, phi_nm1_file, niter_file, factor_file)
# # # # ]
1067 : :
1068 [ # # ]: 0 : if (niter_file == 0) then
1069 : : write(msg,"(a,3f8.4,3a)")&
1070 : 0 : " Could not find q-point: ",BSp%q(:,iq_search)," in file ",TRIM(restart_file),&
1071 : 0 : " Cannot restart Haydock iterations for this q-point"
1072 : 0 : ABI_COMMENT(msg)
1073 : : else
1074 : 0 : write(msg,'(a,i0)')" Number of iterations already performed: ",niter_file
1075 : 0 : call wrtout(units, msg)
1076 : :
1077 [ # # ]: 0 : if ( ABS(haydock_file%broad - BSp%broad) > tol6) then
1078 : : write(msg,'(2a,2(a,f8.4),a)')&
1079 : 0 : " Restart file has been produced with a different Lorentzian broadening: ",ch10,&
1080 : 0 : " broad_file: ",haydock_file%broad," input broadening: ",BSp%broad," Continuing anyway. "
1081 : 0 : ABI_WARNING(msg)
1082 : : end if
1083 : :
1084 : 0 : call haydock_file%close()
1085 : : end if
1086 : : end if
1087 : : !
1088 : : ! Master broadcasts the data.
1089 : 0 : call xmpi_bcast(niter_file,master,comm,ierr)
1090 : :
1091 [ # # ]: 0 : if (my_rank/=master) then
1092 [ # # # # : 0 : ABI_MALLOC(aa_file,(niter_file))
# # # # #
# # # #
# ]
1093 [ # # # # : 0 : ABI_MALLOC(bb_file,(niter_file))
# # # # #
# # # #
# ]
1094 [ # # # # : 0 : ABI_MALLOC(phi_nm1_file,(hsize))
# # # # #
# # # #
# ]
1095 [ # # # # : 0 : ABI_MALLOC(phi_n_file,(hsize))
# # # # #
# ]
1096 : : end if
1097 : :
1098 : 0 : call xmpi_bcast(aa_file,master,comm,ierr)
1099 : 0 : call xmpi_bcast(bb_file,master,comm,ierr)
1100 : 0 : call xmpi_bcast(phi_nm1_file,master,comm,ierr)
1101 : 0 : call xmpi_bcast(phi_n_file,master,comm,ierr)
1102 : :
1103 [ # # # # ]: 0 : end subroutine haydock_restart
1104 : : !!***
1105 : :
1106 : : !----------------------------------------------------------------------
1107 : :
1108 : : !!****f* m_haydock/haydock_mdf_to_tensor
1109 : : !! NAME
1110 : : !! haydock_mdf_to_tensor
1111 : : !!
1112 : : !! FUNCTION
1113 : : !! Transform macroscopic dielectric function from green function to each components of the tensor in red and cart coord.
1114 : : !!
1115 : : !! INPUTS
1116 : : !! BSp<type(excparam)>=Parameters defining the Bethe-Salpeter calculation.
1117 : : !! omega(BSp%nomega)=Frequency mesh for the macroscopic dielectric function (broadening is already included).
1118 : : !! Cryst=Parameters of the crystal
1119 : : !! eps(BSp%nomega,BSp%nq) = Macroscopic dielectric function to be written.
1120 : : !!
1121 : : !! OUTPUT
1122 : : !! tensor_cart(BSp%nomega,6) = dielectric tensor for each frequency, order (11,22,33,12,13,23) in cart. coord.
1123 : : !! tensor_red(BSp%nomega, 6) = idem in reduced coordinated
1124 : : !! ierr = 0 if the tensors have been successfully computed
1125 : : !! \= 0 if the system is ill-posed in terms of q-points (not enough or not independent q-points)
1126 : : !!
1127 : : !! SOURCE
1128 : :
1129 : 84 : subroutine haydock_mdf_to_tensor(BSp,Cryst,eps,tensor_cart,tensor_red,ierr)
1130 : :
1131 : : !Arguments ------------------------------------
1132 : : !scalars
1133 : : integer,intent(out) :: ierr
1134 : : type(excparam),intent(in) :: BSp
1135 : : type(crystal_t),intent(in) :: Cryst
1136 : : !arrays
1137 : : complex(dp),intent(in) :: eps(BSp%nomega,BSp%nq)
1138 : : complex(dp),intent(out) :: tensor_cart(BSp%nomega,6), tensor_red(BSp%nomega,6)
1139 : :
1140 : : !Local variables ------------------------------
1141 : : !scalars
1142 : : integer :: iq,info
1143 : : real(dp) :: normqcart, normqred
1144 : : !arrays
1145 : 84 : integer,allocatable :: ipiv(:)
1146 : : real(dp) :: qcart(3), qtmet(3)
1147 : : real(dp) :: qred2cart(3,3),qcart2red(3,3)
1148 : 168 : complex(dp) :: qqcart(BSp%nq,6), qqred(BSp%nq,6)
1149 : 0 : complex(dp) :: b(6,BSP%nomega)
1150 : : !************************************************************************
1151 : :
1152 : : ! Error flag
1153 : 84 : ierr = 0
1154 : :
1155 [ - + ]: 84 : if(BSp%nq /= 6) then
1156 : 0 : ierr = -1
1157 : 0 : return
1158 : : end if
1159 : :
1160 : : ! Transformation matrices from reduced coordinates to cartesian coordinates
1161 [ + + + + ]: 1092 : qred2cart = two_pi*Cryst%gprimd
1162 : 84 : qcart2red = qred2cart
1163 : 84 : call matrginv(qcart2red,3,3)
1164 [ + + ]: 588 : do iq = 1, 6
1165 : :
1166 : : ! Computing cartesian q-vector
1167 [ + + + + : 6552 : qcart = MATMUL(qred2cart, BSp%q(:,iq))
+ - ]
1168 : :
1169 : : ! Computing product 'metric - qred' to form quadratic form
1170 [ + + + + : 8568 : qtmet = (two_pi**2)*MATMUL(Cryst%gmet, BSp%q(:,iq))
+ - + + ]
1171 : :
1172 : : ! squared norms
1173 : 504 : normqcart = qcart(1)**2+qcart(2)**2+qcart(3)**2
1174 : 504 : normqred = (normv(BSp%q(:,iq),Cryst%gmet,"G"))**2
1175 : :
1176 : : ! Compute line 'iq' for matrix in cartesian coord
1177 : 504 : qqcart(iq,1) = (qcart(1))**2
1178 : 504 : qqcart(iq,2) = (qcart(2))**2
1179 : 504 : qqcart(iq,3) = (qcart(3))**2
1180 : 504 : qqcart(iq,4) = 2*(qcart(1)*qcart(2))
1181 : 504 : qqcart(iq,5) = 2*(qcart(1)*qcart(3))
1182 : 504 : qqcart(iq,6) = 2*(qcart(2)*qcart(3))
1183 : :
1184 : : ! Compute line 'iq' for matrix in reduced coord
1185 : 504 : qqred(iq,1) = (qtmet(1))**2
1186 : 504 : qqred(iq,2) = (qtmet(2))**2
1187 : 504 : qqred(iq,3) = (qtmet(3))**2
1188 : 504 : qqred(iq,4) = 2*(qtmet(1)*qtmet(2))
1189 : 504 : qqred(iq,5) = 2*(qtmet(1)*qtmet(3))
1190 : 504 : qqred(iq,6) = 2*(qtmet(2)*qtmet(3))
1191 : :
1192 : : ! Renormalize line
1193 [ + + ]: 3528 : qqcart(iq,:) = qqcart(iq,:)/normqcart
1194 [ + + ]: 3612 : qqred(iq,:) = qqred(iq,:)/normqred
1195 : : end do
1196 : :
1197 [ - + ]: 84 : ABI_MALLOC(ipiv,(6))
1198 : :
1199 : : ! Solving linear system
1200 [ + + + + ]: 279825 : b = TRANSPOSE(eps)
1201 : 84 : call ZGESV(6,BSp%nomega,qqcart,6,ipiv,b,6,info)
1202 [ + + + + ]: 240366 : tensor_cart = TRANSPOSE(b)
1203 : :
1204 [ - + ]: 84 : if(info /= 0) then
1205 : : ! Skipping the rest of the routine
1206 : 0 : ierr = info
1207 [ # # ]: 0 : ABI_FREE(ipiv)
1208 : 0 : return
1209 : : end if
1210 : :
1211 [ + + + + ]: 279825 : b = TRANSPOSE(eps)
1212 : 84 : call ZGESV(6,BSp%nomega,qqred,6,ipiv,b,6,info)
1213 [ + + + + ]: 240366 : tensor_red = TRANSPOSE(b)
1214 : :
1215 [ - + ]: 84 : if(info /= 0) ierr = info
1216 : :
1217 [ - + ]: 84 : ABI_FREE(ipiv)
1218 : :
1219 : : end subroutine haydock_mdf_to_tensor
1220 : : !!***
1221 : :
1222 : : !----------------------------------------------------------------------
1223 : :
1224 : : !!****f* m_haydock/haydock_psherm
1225 : : !! NAME
1226 : : !! haydock_psherm
1227 : : !!
1228 : : !! FUNCTION
1229 : : !! Reads the excitonic Hamiltonian from file and construct the Lanczos set of vectors
1230 : : !! by iterative matrix-vector multiplications.
1231 : : !!
1232 : : !! INPUTS
1233 : : !! BSp<type(excparam)>=Parameters defining the Bethe-Salpeter calculation.
1234 : : !! omega(BSp%nomega)=Frequency mesh for the macroscopic dielectric function (broadening is already included).
1235 : : !! hize
1236 : : !! my_t1,my_t2
1237 : : !! hreso(hsize,my_t1:my_t2)
1238 : : !! hcoup(hsize,my_t1:my_t2)
1239 : : !! nkets
1240 : : !! kets(hsize,nkets)
1241 : : !! comm=MPI communicator.
1242 : : !!
1243 : : !! OUTPUT
1244 : : !! green(BSp%nomega)=The imaginary part of the macroscopic dielectric function.
1245 : : !!
1246 : : !! SOURCE
1247 : :
1248 : 0 : subroutine haydock_psherm(BSp,BS_files,Cryst,Hdr_bse,hexc,hexc_i,hsize,my_t1,my_t2,nkets,kets,green,comm)
1249 : :
1250 : : !Arguments ------------------------------------
1251 : : !scalars
1252 : : integer,intent(in) :: hsize,my_t1,my_t2,nkets,comm
1253 : : type(crystal_t),intent(in) :: Cryst
1254 : : type(excparam),intent(in) :: BSp
1255 : : type(excfiles),intent(in) :: BS_files
1256 : : type(Hdr_type),intent(in) :: Hdr_bse
1257 : : type(hexc_t),intent(in) :: hexc
1258 : : type(hexc_interp_t),intent(in) :: hexc_i
1259 : : !arrays
1260 : : complex(dp),intent(out) :: green(BSp%nomega,BSp%nq)
1261 : : complex(dp),intent(in) :: kets(hsize,nkets)
1262 : :
1263 : : !Local variables ------------------------------
1264 : : !scalars
1265 : : integer,parameter :: master=0
1266 : : integer :: inn,itt,out_unt,nproc,my_rank,ierr
1267 : : integer :: niter_file,niter_max,niter_done,nsppol,iq,my_nt,term_type
1268 : : real(dp) :: ket0_hbar_norm,nfact
1269 : : logical :: can_restart,is_converged
1270 : : complex(dp) :: factor
1271 : : character(len=fnlen),parameter :: tag_file="_HAYDC_SAVE"
1272 : : character(len=500) :: msg
1273 : : character(len=fnlen) :: restart_file,out_file
1274 : : !arrays
1275 : : integer :: units(2)
1276 : 0 : real(dp),allocatable :: bb_file(:), bb(:)
1277 : 0 : complex(dp),allocatable :: aa(:),cc(:),phi_np1(:),phi_n(:),phi_nm1(:),cbuff(:)
1278 : 0 : complex(dp),allocatable :: aa_file(:),phi_n_file(:),phi_np1_file(:),cc_file(:), ket0(:)
1279 : : logical :: check(2)
1280 : : !************************************************************************
1281 : :
1282 : 0 : ABI_WARNING("Haydock + coupling is still under development")
1283 [ # # ]: 0 : units = [std_out, ab_out]
1284 : :
1285 [ # # ]: 0 : if(BSp%use_interp) then
1286 : 0 : ABI_ERROR("Coupling is not yet implemented with interpolation")
1287 : : end if
1288 : :
1289 : 0 : nproc = xmpi_comm_size(comm)
1290 : 0 : my_rank= xmpi_comm_rank(comm)
1291 : 0 : nsppol = Hdr_bse%nsppol
1292 : :
1293 : 0 : my_nt = my_t2-my_t1+1
1294 [ # # ]: 0 : ABI_CHECK(my_nt>0,"One of the processors has zero columns")
1295 : :
1296 : : ! Multiplicative factor (k-point sampling and unit cell volume)
1297 : : ! TODO be careful with the spin here
1298 : : ! TODO four_pi comes from the coulomb term 1/|q| is already included in the
1299 : : ! oscillators hence the present approach wont work if a cutoff interaction is used.
1300 : 0 : nfact = four_pi/(Cryst%ucvol*BSp%nkbz)
1301 [ # # ]: 0 : if (nsppol==1) nfact=two*nfact
1302 : :
1303 : 0 : write(msg,'(a,i0)')' Haydock algorithm with MAX number of iterations: ',BSp%niter
1304 : 0 : call wrtout(std_out,msg)
1305 : : !
1306 : : ! Check for presence of the restart file.
1307 : 0 : can_restart=.FALSE.
1308 : :
1309 [ # # ]: 0 : if ( BS_files%in_haydock_basename /= BSE_NOFILE) then
1310 [ # # ]: 0 : restart_file = strcat(BS_files%in_haydock_basename,tag_file)
1311 [ # # ]: 0 : if (file_exists(restart_file) ) then
1312 : 0 : can_restart=.TRUE.
1313 [ # # ]: 0 : msg = strcat(" Restarting Haydock calculation from file: ",restart_file)
1314 : 0 : call wrtout(units, msg)
1315 : 0 : ABI_ERROR("Restart is not tested")
1316 : : else
1317 : 0 : can_restart=.FALSE.
1318 : 0 : ABI_WARNING(strcat("Cannot find restart file: ",restart_file))
1319 : : end if
1320 : : end if
1321 : : !
1322 : : ! Open the file and writes basic dimensions and info.
1323 [ # # ]: 0 : if (my_rank==master) then
1324 [ # # ]: 0 : out_file = TRIM(BS_files%out_basename)//TRIM(tag_file)
1325 [ # # ]: 0 : if (open_file(out_file,msg,newunit=out_unt,form="unformatted") /= 0) then
1326 : 0 : ABI_ERROR(msg)
1327 : : end if
1328 : : ! write header TODO: standardize this part.
1329 : 0 : write(out_unt)hsize,Bsp%use_coupling,BSE_HAYD_IMEPS,nkets,Bsp%broad
1330 : : end if
1331 : : !
1332 : : ! Select the terminator for the continued fraction.
1333 : 0 : term_type=0 !; if (Bsp%hayd_term>0) term_type=2
1334 : 0 : call wrtout(std_out,sjoin("Using terminator type: ",itoa(term_type)))
1335 : : !
1336 : : ! Calculate green(w) for the different starting kets.
1337 [ # # # # ]: 0 : green=czero
1338 [ # # ]: 0 : do iq=1,nkets
1339 [ # # # # : 0 : ABI_MALLOC(ket0,(my_nt))
# # # # #
# # # ]
1340 [ # # # # : 0 : ket0 = kets(my_t1:my_t2,iq)
# # # # ]
1341 : : !
1342 : 0 : niter_file=0
1343 : :
1344 [ # # ]: 0 : if (can_restart) then
1345 : : call haydock_restart(BSp,restart_file,BSE_HAYD_IMEPS,iq,hsize,&
1346 [ # # # # : 0 : niter_file,aa_file,bb_file,phi_np1_file,phi_n_file,comm)
# # # # ]
1347 : : end if
1348 : :
1349 [ # # # # : 0 : ABI_MALLOC(phi_nm1,(my_nt))
# # # # ]
1350 [ # # # # : 0 : ABI_MALLOC(phi_n,(my_nt))
# # # # ]
1351 [ # # # # : 0 : ABI_MALLOC(phi_np1,(my_nt))
# # # # ]
1352 : : !
1353 : : ! TODO: Note the different convention used for the coefficients
1354 : : ! Should use the same convention in the Hermitian case.
1355 : 0 : niter_max = niter_file + Bsp%niter
1356 [ # # # # : 0 : ABI_MALLOC(aa,(niter_max))
# # # # #
# # # ]
1357 [ # # # # : 0 : ABI_MALLOC(bb,(niter_max+1))
# # # # #
# # # #
# ]
1358 [ # # # # : 0 : ABI_MALLOC(cc,(niter_max+1))
# # # # #
# # # ]
1359 [ # # # # : 0 : aa=czero; bb=czero; cc=czero
# # ]
1360 : :
1361 [ # # ]: 0 : if (niter_file==0) then ! Calculation from scratch.
1362 [ # # # # : 0 : phi_n = ket0
# # # # ]
1363 : 0 : call hexc%matmul_full(hexc_i, phi_n, phi_np1, -1)
1364 : : !phi_np1 = MATMUL(hreso,ket0) - MATMUL(hcoup,CONJG(ket0))
1365 [ # # ]: 0 : ket0_hbar_norm = SQRT(two*DBLE(DOT_PRODUCT(phi_n,phi_np1)))
1366 [ # # ]: 0 : phi_n = phi_n /ket0_hbar_norm
1367 [ # # ]: 0 : phi_np1 = phi_np1/ket0_hbar_norm
1368 : : !ket0 = ket0/ket0_hbar_norm
1369 : 0 : cc(1)=zero ! <P|F|P>
1370 : : !cc(1) = DOT_PRODUCT(ket0,phi_np1)
1371 : : !write(std_out,*)" cc(1), ket0_hbar_norm =",cc(1),ket0_hbar_norm
1372 : :
1373 [ # # ]: 0 : phi_nm1 = czero
1374 : 0 : niter_done=0 ! TODO Be careful here
1375 : :
1376 : : else ! Use the previously calculates a and b.
1377 : 0 : niter_done=niter_file
1378 : 0 : ABI_ERROR("Restart not coded")
1379 : : !aa(1:niter_done) = aa_file
1380 : : !bb(1:niter_done) = bb_file
1381 : : !phi_np1=phi_np1_file(my_t1:my_t2) ! Select the slice treated by this node.
1382 : : !phi_n =phi_n_file (my_t1:my_t2)
1383 : : end if
1384 : :
1385 [ # # ]: 0 : if (can_restart) then
1386 [ # # ]: 0 : ABI_FREE(aa_file)
1387 [ # # ]: 0 : ABI_FREE(bb_file)
1388 : 0 : ABI_FREE(cc_file)
1389 : : ABI_FREE(phi_np1_file)
1390 : : ABI_FREE(phi_n_file)
1391 : : end if
1392 : :
1393 : : ! This factor gives the correct results
1394 : 0 : factor = -nfact*ket0_hbar_norm / SQRT(two)
1395 : :
1396 : : ! Which quantity should be checked for convergence?
1397 : 0 : check = (/.TRUE.,.TRUE./)
1398 [ # # ]: 0 : if (ABS(Bsp%haydock_tol(2)-one)<tol6) check = (/.TRUE. ,.FALSE./)
1399 [ # # ]: 0 : if (ABS(Bsp%haydock_tol(2)-two)<tol6) check = (/.FALSE.,.TRUE./)
1400 : :
1401 : : call haydock_psherm_optalgo(niter_done,niter_max,BSp%nomega,BSp%omega,BSp%haydock_tol(1),check,hexc,hexc_i,&
1402 : : hsize,my_t1,my_t2,factor,term_type,aa,bb,cc,ket0,ket0_hbar_norm,phi_nm1,phi_n,phi_np1,&
1403 : 0 : green(:,iq),inn,is_converged,comm)
1404 : :
1405 : : ! Save the a"s and the b"s for possible restarting.
1406 : : ! 1) Info on the Q.
1407 : : ! 2) Number of iterations performed.
1408 : : ! 3) do iter=1,niter_performed
1409 : : ! aa(iter),bb(iter)
1410 : : ! end do
1411 : : ! 4) |n-1>
1412 : : ! |n>
1413 : : ! |n+1>
1414 : : !
1415 [ # # ]: 0 : if (my_rank==master) then ! Open the file and writes basic dimensions and info.
1416 : 0 : write(out_unt)Bsp%q(:,iq)
1417 : 0 : write(out_unt)MIN(inn,niter_max) ! NB: if the previous loop completed inn=niter_max+1
1418 [ # # ]: 0 : do itt=1,MIN(inn,niter_max) ! if we exited then inn is not incremented by one.
1419 : 0 : write(out_unt)itt,aa(itt),bb(itt)
1420 : : end do
1421 : : end if
1422 : : !
1423 : : ! cbuff is used as workspace to gather |n-1>, |n> and |n+1>.
1424 [ # # # # : 0 : ABI_MALLOC(cbuff,(hsize))
# # # # #
# # # #
# ]
1425 [ # # # # ]: 0 : cbuff=czero; cbuff(my_t1:my_t2) = phi_nm1
1426 : 0 : call xmpi_sum_master(cbuff,master,comm,ierr)
1427 [ # # ]: 0 : if (my_rank==master) write(out_unt) cbuff ! |n-1>
1428 : :
1429 [ # # # # ]: 0 : cbuff=czero; cbuff(my_t1:my_t2) = phi_n
1430 : 0 : call xmpi_sum_master(cbuff,master,comm,ierr)
1431 [ # # ]: 0 : if (my_rank==master) write(out_unt) cbuff ! |n>
1432 : :
1433 [ # # # # ]: 0 : cbuff=czero; cbuff(my_t1:my_t2) = phi_np1
1434 : 0 : call xmpi_sum_master(cbuff,master,comm,ierr)
1435 [ # # ]: 0 : if (my_rank==master) write(out_unt) cbuff ! |n+1>
1436 : :
1437 [ # # ]: 0 : ABI_FREE(phi_nm1)
1438 [ # # ]: 0 : ABI_FREE(phi_n)
1439 [ # # ]: 0 : ABI_FREE(phi_np1)
1440 [ # # ]: 0 : ABI_FREE(cbuff)
1441 [ # # ]: 0 : ABI_FREE(aa)
1442 [ # # ]: 0 : ABI_FREE(bb)
1443 [ # # ]: 0 : ABI_FREE(cc)
1444 [ # # ]: 0 : ABI_FREE(ket0)
1445 : : end do ! iq
1446 : :
1447 [ # # ]: 0 : if (my_rank==master) close(out_unt)
1448 : :
1449 : 0 : call xmpi_barrier(comm)
1450 : :
1451 [ # # # # : 0 : end subroutine haydock_psherm
# # # # #
# ]
1452 : : !!***
1453 : :
1454 : : !----------------------------------------------------------------------
1455 : :
1456 : : !!****f* m_haydock/haydock_psherm_optalgo
1457 : : !! NAME
1458 : : !! haydock_psherm_optalgo
1459 : : !!
1460 : : !! FUNCTION
1461 : : !! Haydock algorithm for pseudo-hermitian matrix
1462 : : !!
1463 : : !! INPUTS
1464 : : !! niter_done=Number of iterations already performed (0 if the run starts from scratch).
1465 : : !! niter_tot=Max number of iterations. Always > niter_done
1466 : : !! nomega=Number of Frequency points for the evaluation of the matrix element.
1467 : : !! omega(nomega)=Frequency set (imaginary part is already included).
1468 : : !! tol_iter=Tollerance used to stop the the algorithm.
1469 : : !! check(2)=Logical flags to specify where both the real and the imaginary part of the
1470 : : !! matrix elements of the Green functions have to be checked for convergence.
1471 : : !! hsize=Size of the blocks.
1472 : : !! my_t1,my_t2=Indeces of the first and last column stored treated by this done.
1473 : : !! term_type=0 if no terminator is used, 1 otherwise.
1474 : : !! hreso(hsize,my_t1:my_t2)=The columns of the resonant block.
1475 : : !! hcoup(hsize,my_t1:my_t2)=The columns of the coupling block.
1476 : : !! factor
1477 : : !! comm=MPI communicator.
1478 : : !!
1479 : : !! OUTPUT
1480 : : !! green(nomega)=Output matrix elements.
1481 : : !! inn=Last iteration performed.
1482 : : !! is_converged=.TRUE. of the algorithm converged.
1483 : : !!
1484 : : !! SIDE EFFECTS
1485 : : !! phi_nm1(my_t2-my_t1+1), phi_n(my_t2-my_t1+1)
1486 : : !! input: vectors used to initialize the iteration
1487 : : !! output: the vectors obtained in the last iteration
1488 : : !! aa(niter_tot) and bb(niter_tot+1)
1489 : : !! if niter_done>0: aa(1:niter_done), bb(1:niter_done) store the coefficients of the previous run.
1490 : : !! when the routine returns aa(1:inn) and bb(1:inn) contain the matrix elements of the tridiagonal form.
1491 : : !! cc(niter_tot+1)
1492 : : !!
1493 : : !! SOURCE
1494 : :
1495 : 0 : subroutine haydock_psherm_optalgo(niter_done,niter_tot,nomega,omega,tol_iter,check,hexc,hexc_i,hsize,my_t1,my_t2,&
1496 : 0 : & factor,term_type,aa,bb,cc,ket0,ket0_hbar_norm,phi_nm1,phi_n,phi_np1,green,inn,is_converged,comm)
1497 : :
1498 : : !Arguments ------------------------------------
1499 : : !scalars
1500 : : integer,intent(in) :: niter_tot,niter_done,nomega,comm,hsize,my_t1,my_t2,term_type
1501 : : integer,intent(out) :: inn
1502 : : logical,intent(out) :: is_converged
1503 : : real(dp),intent(in) :: tol_iter,ket0_hbar_norm
1504 : : complex(dp),intent(in) :: factor
1505 : : type(hexc_t),intent(in) :: hexc
1506 : : type(hexc_interp_t),intent(in) :: hexc_i
1507 : : !arrays
1508 : : real(dp),intent(inout) :: bb(niter_tot+1)
1509 : : complex(dp),intent(out) :: green(nomega)
1510 : : complex(dp),intent(in) :: omega(nomega)
1511 : : complex(dp),intent(inout) :: aa(niter_tot),cc(niter_tot+1)
1512 : : complex(dp),intent(in) :: ket0(my_t2-my_t1+1)
1513 : : complex(dp),intent(inout) :: phi_nm1(my_t2-my_t1+1)
1514 : : complex(dp),intent(inout) :: phi_n (my_t2-my_t1+1)
1515 : : complex(dp),intent(inout) :: phi_np1(my_t2-my_t1+1)
1516 : : logical,intent(in) :: check(2)
1517 : :
1518 : : !Local variables ------------------------------
1519 : : !scalars
1520 : : integer :: my_nt,niter_min,nconv,parity,ii,jj,tdim,ierr
1521 : : integer :: row_max,col_max,nlev
1522 : : character(len=500) :: msg
1523 : : real(dp) :: max_err,mean_err,mean_err2,std_dev,err
1524 : : logical :: keep_vectors=.TRUE.
1525 : : !arrays
1526 : : integer :: units(2)
1527 : 0 : real(dp) :: abs_err(nomega,2) !,ww_err(nomega,2)
1528 : 0 : complex(dp) :: gn0(nomega,niter_tot)
1529 : 0 : complex(dp),allocatable :: oldg(:),newg(:)
1530 : 0 : complex(dp),allocatable :: hphi_n(:),save_phi(:,:)
1531 : 0 : complex(dp),allocatable :: alpha(:,:),beta(:,:),ovlp(:,:)
1532 : 0 : complex(dp),allocatable :: phi_test(:),phi_test2(:),g00(:)
1533 : : logical :: test(2)
1534 : : !************************************************************************
1535 : :
1536 : : ABI_UNUSED(ket0_hbar_norm)
1537 [ # # ]: 0 : units = [std_out, ab_out]
1538 : :
1539 : 0 : my_nt = my_t2-my_t1+1
1540 : :
1541 [ # # # # : 0 : ABI_MALLOC(oldg,(nomega))
# # # # #
# # # ]
1542 [ # # # # : 0 : ABI_MALLOC(newg,(nomega))
# # # # ]
1543 [ # # # # : 0 : ABI_MALLOC(g00,(nomega))
# # # # ]
1544 [ # # # # : 0 : oldg=czero; newg=czero; g00=czero
# # ]
1545 : 0 : nconv=0
1546 : :
1547 [ # # # # ]: 0 : keep_vectors = (keep_vectors.and.xmpi_comm_size(comm)==1)
1548 [ # # ]: 0 : if (keep_vectors) then
1549 [ # # # # : 0 : ABI_MALLOC_OR_DIE(save_phi,(my_t2-my_t1+1,niter_tot), ierr)
# # # # #
# # # # #
# # # # #
# ]
1550 [ # # # # ]: 0 : save_phi=czero
1551 : : end if
1552 : :
1553 [ # # # # : 0 : ABI_MALLOC(hphi_n,(hsize))
# # # # #
# # # ]
1554 : :
1555 [ # # ]: 0 : do inn=niter_done+1,niter_tot
1556 : : !
1557 : : ! a(n) = <Vn+1|F|Vn+1> = <Vn|HFH|Vn>) = 0 by symmetry.
1558 : 0 : aa(inn)=zero
1559 : :
1560 : : ! |n+1> = |n+1> - a(n)|Vn> - a(n)|n-1>
1561 [ # # ]: 0 : phi_np1 = phi_np1 - bb(inn)*phi_nm1
1562 : : !
1563 : : ! |n-1> = |n>
1564 : : ! |n> = |n+1>
1565 [ # # ]: 0 : phi_nm1 = phi_n
1566 [ # # ]: 0 : phi_n = phi_np1
1567 : :
1568 : : !|n+1> = H |n> using all eh components.
1569 : 0 : parity = (-1)**(inn+1)
1570 : 0 : call hexc%matmul_full(hexc_i, phi_n, phi_np1, parity)
1571 : :
1572 : : !phi_np1 = MATMUL(hreso,phi_n) + parity * MATMUL(hcoup,CONJG(phi_n))
1573 : : !call xmpi_sum(hphi_np1,comm,ierr)
1574 : : !
1575 : : ! B(n+1)= <n|F|n+1>^(1/2) = <n|FH|n>^(1/2))= (2*Re(<n|V+1>))^(1/2)
1576 : : ! by symmetry, where the dot_product is done in the resonant eh sub-space.
1577 : : !
1578 [ # # ]: 0 : bb(inn+1)=SQRT(two*DBLE(DOT_PRODUCT(phi_n,phi_np1)))
1579 : : !bb(inn+1)=two*DBLE(DOT_PRODUCT(phi_n,phi_np1))
1580 : : !call xmpi_sum(bb(inn+1),comm,ierr)
1581 : : !bb(inn+1)=SQRT(bb(inn+1)
1582 : : !
1583 : : !|n+1> =|n+1>/B(n+1)
1584 [ # # ]: 0 : phi_n = phi_n /bb(inn+1)
1585 [ # # ]: 0 : phi_np1 = phi_np1/bb(inn+1)
1586 : :
1587 [ # # # # ]: 0 : if (keep_vectors) save_phi(:,inn) = phi_n
1588 : :
1589 : 0 : parity = (-1)**(inn+1)
1590 : : !if (parity==-1) then
1591 : : ! cc(inn+1)=czero
1592 : : !else
1593 [ # # # # ]: 0 : cc(inn+1)=DOT_PRODUCT(ket0,phi_n) + parity * DOT_PRODUCT(phi_n,ket0)
1594 : : !end if
1595 : : !call xmpi_sum(cc(inn+1),comm,ierr)
1596 : :
1597 : 0 : write(msg,'(a,i0,a,3es12.4)')' Iteration number ',inn,', b_i RE(c_i+1) IM(c_i+1) ',bb(inn),REAL(cc(inn+1)),AIMAG(cc(inn+1))
1598 : 0 : call wrtout(std_out,msg)
1599 : :
1600 : 0 : call continued_fract(inn,term_type,aa,bb(2:),nomega,omega,g00)
1601 [ # # ]: 0 : gn0(:,1) = g00
1602 : :
1603 : : if (.FALSE.) then
1604 : : gn0(:,2) = (one - omega(:)*g00(:))/bb(2)
1605 : : do ii=3,inn
1606 : : gn0(:,ii) = -(-bb(ii)*gn0(:,ii-2) -omega(:)*gn0(:,ii-1))/bb(ii+1)
1607 : : end do
1608 : : else
1609 [ # # ]: 0 : do ii=2,inn
1610 : 0 : nlev = inn-ii
1611 : 0 : call continued_fract(nlev,term_type,aa,bb(ii+1:),nomega,omega,g00)
1612 [ # # ]: 0 : gn0(:,ii) = +bb(ii+1) * g00 * gn0(:,ii-1)
1613 : : end do
1614 : : end if
1615 : :
1616 [ # # ]: 0 : newg=czero
1617 [ # # ]: 0 : do ii=1,inn
1618 [ # # ]: 0 : newg(:) = newg + cc(ii)* gn0(:,ii)
1619 : : end do
1620 [ # # ]: 0 : newg = factor*newg
1621 : : !
1622 : : ! Avoid spurious convergence.
1623 [ # # ]: 0 : niter_min=4; if (niter_done>1) niter_min=niter_done+1
1624 [ # # ]: 0 : if (inn>niter_min) then
1625 [ # # ]: 0 : test=.TRUE.
1626 [ # # ]: 0 : abs_err(:,1) = ABS(DBLE (newg-oldg))
1627 [ # # ]: 0 : abs_err(:,2) = ABS(AIMAG(newg-oldg))
1628 : : !
1629 [ # # ]: 0 : if (tol_iter>zero) then
1630 : : ! Test on the L1 norm.
1631 [ # # # # : 0 : if (check(1)) test(1) = SUM(abs_err(:,1)) < tol_iter*SUM(ABS(DBLE (newg)))
# # ]
1632 [ # # # # : 0 : if (check(2)) test(2) = SUM(abs_err(:,2)) < tol_iter*SUM(ABS(AIMAG(newg)))
# # ]
1633 : : else
1634 : : ! Stringent test for each point.
1635 [ # # # # : 0 : if (check(1)) test(1) = ALL( abs_err(:,1) < -tol_iter*ABS(DBLE (newg)))
# # ]
1636 [ # # # # : 0 : if (check(2)) test(2) = ALL( abs_err(:,2) < -tol_iter*ABS(AIMAG(newg)))
# # ]
1637 : : end if
1638 : : !
1639 [ # # # # : 0 : if (ALL(test)) then
# # ]
1640 : 0 : nconv = nconv+1
1641 : : else
1642 : : nconv = 0
1643 : : end if
1644 [ # # ]: 0 : if (nconv==2) then
1645 [ # # ]: 0 : if(inn<100)then
1646 : : write(msg,'(a,es10.2,a)')&
1647 : 0 : " >>> Haydock algorithm converged twice within haydock_tol= ",tol_iter," after less than 100 iterations."
1648 : : else
1649 : : write(msg,'(a,es10.2,a,i0,a)')&
1650 : 0 : " >>> Haydock algorithm converged twice within haydock_tol= ",tol_iter," after ",inn," iterations."
1651 : : endif
1652 : 0 : call wrtout(units, msg)
1653 : 0 : EXIT
1654 : : end if
1655 : : end if
1656 : : !
1657 [ # # # # ]: 0 : oldg = newg
1658 : : end do ! inn
1659 : :
1660 [ # # ]: 0 : green = newg
1661 [ # # ]: 0 : if (nconv/=2) then
1662 : : write(msg,'(a,es10.2,a,i0,a)')&
1663 : 0 : " WARNING: Haydock algorithm did not converge within ",tol_iter," after ",niter_tot," iterations."
1664 : 0 : call wrtout(units, msg)
1665 : : end if
1666 : :
1667 : 0 : is_converged = (nconv==2)
1668 : :
1669 [ # # ]: 0 : ABI_FREE(oldg)
1670 [ # # ]: 0 : ABI_FREE(newg)
1671 [ # # ]: 0 : ABI_FREE(g00)
1672 [ # # ]: 0 : ABI_FREE(hphi_n)
1673 : :
1674 [ # # ]: 0 : if (keep_vectors) then
1675 : 0 : tdim = MIN(inn,niter_tot)
1676 [ # # # # : 0 : ABI_MALLOC(ovlp,(tdim,tdim))
# # # # #
# # # # #
# # ]
1677 : :
1678 [ # # # # : 0 : ABI_MALLOC(phi_test,(hsize))
# # # # ]
1679 [ # # # # : 0 : ABI_MALLOC(phi_test2,(hsize))
# # # # ]
1680 : :
1681 : 0 : max_err=smallest_real; mean_err=zero; mean_err2=zero; row_max=-1
1682 [ # # ]: 0 : do ii=1,tdim
1683 : 0 : parity = (-1)**(ii+1)
1684 [ # # # # : 0 : phi_test = save_phi(:,ii)
# # # # ]
1685 : 0 : call hexc%matmul_full(hexc_i, phi_test, phi_test2, parity)
1686 : : !phi_test2 = MATMUL(hreso,phi_test) + parity * MATMUL(hcoup,CONJG(phi_test))
1687 [ # # # # ]: 0 : ovlp(ii,ii) = DOT_PRODUCT(phi_test,phi_test2) + DOT_PRODUCT(phi_test2,phi_test)
1688 : 0 : err = ABS(ovlp(ii,ii)-cone)
1689 : 0 : mean_err = mean_err + err
1690 : 0 : mean_err2 = mean_err2 + err**2
1691 [ # # ]: 0 : if (err > max_err) then
1692 : 0 : max_err = err
1693 : 0 : row_max = ii
1694 : : end if
1695 : : end do
1696 : 0 : mean_err = mean_err/tdim
1697 : 0 : std_dev = mean_err2/tdim -mean_err**2
1698 : : write(std_out,'(a,i0,1x,3es14.6)') &
1699 : 0 : " Error in normalization (ii, max_err,mean,std_dev): ",row_max,max_err,mean_err,std_dev
1700 : :
1701 [ # # ]: 0 : ABI_FREE(phi_test)
1702 [ # # ]: 0 : ABI_FREE(phi_test2)
1703 : :
1704 [ # # # # : 0 : ABI_MALLOC(alpha,(hsize,tdim))
# # # # #
# # # # #
# # ]
1705 : :
1706 : : ! Less efficient but for sake of simplicity with hexc_matmul
1707 : : ! TODO possibility to call hreso * phi, and hcoup * phi separately
1708 [ # # ]: 0 : do ii=1,tdim
1709 : 0 : parity = (-1)**(ii+1)
1710 : 0 : call hexc%matmul_full(hexc_i, save_phi(:,ii), alpha(:,ii), parity)
1711 : : end do
1712 : :
1713 : : !alpha = MATMUL(hreso,save_phi(:,1:tdim))
1714 : : !
1715 : : !do ii=1,tdim
1716 : : ! parity = (-1)**(ii+1)
1717 : : ! alpha(:,ii) = alpha(:,ii) + parity*MATMUL(hcoup,CONJG(save_phi(:,ii)))
1718 : : !end do
1719 : :
1720 [ # # # # : 0 : ovlp = MATMUL(TRANSPOSE(CONJG(save_phi(:,1:tdim))),alpha)
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # ]
1721 : :
1722 [ # # # # : 0 : ABI_MALLOC(beta,(hsize,tdim))
# # # # #
# # # ]
1723 [ # # ]: 0 : do ii=1,tdim
1724 : 0 : parity = (-1)**(ii+1)
1725 [ # # ]: 0 : beta(:,ii) = parity*save_phi(:,ii)
1726 [ # # ]: 0 : alpha(:,ii) = -parity*alpha(:,ii)
1727 : : end do
1728 : :
1729 [ # # # # : 0 : ovlp = ovlp - MATMUL(TRANSPOSE(CONJG(beta)),alpha)
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# ]
1730 : :
1731 : 0 : max_err=smallest_real; row_max=-1; col_max=-1
1732 : 0 : mean_err=zero; mean_err2=zero
1733 [ # # ]: 0 : do jj=1,tdim
1734 [ # # ]: 0 : do ii=1,jj
1735 : 0 : err = ABS(ovlp(ii,jj))
1736 [ # # ]: 0 : if (ii==jj) err = ABS(err - one)
1737 : 0 : mean_err = mean_err + err
1738 : 0 : mean_err2 = mean_err2 + err**2
1739 [ # # ]: 0 : if (err > max_err) then
1740 : 0 : max_err = err
1741 : 0 : row_max=ii
1742 : 0 : col_max=jj
1743 : : end if
1744 : : end do
1745 : : end do
1746 : :
1747 : 0 : mean_err = mean_err/(tdim*(tdim+1)/2)
1748 : 0 : std_dev = mean_err2/(tdim*(tdim+1)/2) - mean_err**2
1749 : : write(std_out,'(a,2(i0,1x),3es14.6)')&
1750 : 0 : " Error in Hbar-ortho (i,j), max_err, mean, std_dev ",row_max,col_max,max_err,mean_err,std_dev
1751 : : !call print_arr(ovlp,max_r=185,max_c=10,unit=std_out)
1752 : :
1753 [ # # ]: 0 : ABI_FREE(alpha)
1754 [ # # ]: 0 : ABI_FREE(beta)
1755 [ # # ]: 0 : ABI_FREE(ovlp)
1756 [ # # ]: 0 : ABI_FREE(save_phi)
1757 : : end if
1758 : :
1759 [ # # # # ]: 0 : end subroutine haydock_psherm_optalgo
1760 : : !!***
1761 : :
1762 : : !----------------------------------------------------------------------
1763 : :
1764 : : !!****f* m_haydock/haydock_bilanczos
1765 : : !! NAME
1766 : : !! haydock_bilanczos
1767 : : !!
1768 : : !! FUNCTION
1769 : : !! Reads the excitonic Hamiltonian from file and construct the Lanczos set of vectors
1770 : : !! by iterative matrix-vector multiplications for any general matrix.
1771 : : !!
1772 : : !! INPUTS
1773 : : !! BSp<type(excparam)>=Parameters defining the Bethe-Salpeter calculation.
1774 : : !! omega(BSp%nomega)=Frequency mesh for the macroscopic dielectric function (broadening is already included).
1775 : : !! hize
1776 : : !! my_t1,my_t2
1777 : : !! hreso(hsize,my_t1:my_t2)
1778 : : !! hcoup(hsize,my_t1:my_t2)
1779 : : !! nkets
1780 : : !! kets(hsize,nkets)
1781 : : !! comm=MPI communicator.
1782 : : !!
1783 : : !! OUTPUT
1784 : : !! green(BSp%nomega)=The imaginary part of the macroscopic dielectric function.
1785 : : !!
1786 : : !! SOURCE
1787 : :
1788 : 10 : subroutine haydock_bilanczos(BSp,BS_files,Cryst,Hdr_bse,hexc,hexc_i,hsize,my_t1,my_t2,nkets,kets,ep_renorms,green,comm)
1789 : :
1790 : : !Arguments ------------------------------------
1791 : : !scalars
1792 : : integer,intent(in) :: hsize,my_t1,my_t2,nkets,comm
1793 : : type(crystal_t),intent(in) :: Cryst
1794 : : type(excparam),intent(in) :: BSp
1795 : : type(excfiles),intent(in) :: BS_files
1796 : : type(Hdr_type),intent(in) :: Hdr_bse
1797 : : !arrays
1798 : : complex(dp),intent(out) :: green(BSp%nomega,BSp%nq)
1799 : : complex(dp),intent(in) :: kets(hsize,nkets)
1800 : : complex(dp),intent(in) :: ep_renorms(hsize)
1801 : :
1802 : : !Local variables ------------------------------
1803 : : !scalars
1804 : : integer,parameter :: master=0
1805 : : integer :: inn,itt,out_unt,nproc,my_rank,ierr
1806 : : integer :: niter_file,niter_max,niter_done,nsppol,iq,my_nt,term_type,n_all_omegas
1807 : : real(dp) :: ket0_hbar_norm,nfact,norm
1808 : : logical :: can_restart,is_converged
1809 : : complex(dp) :: factor
1810 : : character(len=fnlen),parameter :: tag_file="_HAYDC_SAVE"
1811 : : character(len=500) :: msg
1812 : : character(len=fnlen) :: restart_file,out_file
1813 : : type(hexc_t),intent(in) :: hexc
1814 : : type(hexc_interp_t),intent(in) :: hexc_i
1815 : : !arrays
1816 : 10 : complex(dp),allocatable :: aa_file(:),bb_file(:),cc_file(:)
1817 : 10 : complex(dp),allocatable :: aa(:),bb(:),cc(:)
1818 : 10 : complex(dp),allocatable :: phi_np1(:),phi_n(:),phi_nm1(:)
1819 : 10 : complex(dp),allocatable :: phit_np1(:),phit_n(:),phit_nm1(:)
1820 : 10 : complex(dp),allocatable :: cbuff(:), phi_n_file(:),phi_np1_file(:)
1821 : 10 : complex(dp),allocatable :: ket0(:), hphi_n(:), hphit_n(:)
1822 : 10 : complex(dp),allocatable :: all_omegas(:),green_temp(:,:)
1823 : : logical :: check(2)
1824 : : !************************************************************************
1825 : :
1826 : 10 : ABI_WARNING("Haydock with Bilanczos is still under development")
1827 : :
1828 [ - + ]: 10 : if(BSp%use_interp) then
1829 : 0 : ABI_ERROR("Bilanczos is not yet implemented with interpolation")
1830 : : end if
1831 : :
1832 : 10 : nproc = xmpi_comm_size(comm)
1833 : 10 : my_rank= xmpi_comm_rank(comm)
1834 : 10 : nsppol = Hdr_bse%nsppol
1835 : :
1836 : 10 : my_nt = my_t2-my_t1+1
1837 [ - + ]: 10 : ABI_CHECK(my_nt>0,"One of the processors has zero columns")
1838 : :
1839 : : ! Multiplicative factor (k-point sampling and unit cell volume)
1840 : : ! TODO be careful with the spin here
1841 : : ! TODO four_pi comes from the coulomb term 1/|q| is already included in the
1842 : : ! oscillators hence the present approach wont work if a cutoff interaction is used.
1843 : 10 : nfact = four_pi/(Cryst%ucvol*BSp%nkbz)
1844 [ + - ]: 10 : if (nsppol==1) nfact=two*nfact
1845 : :
1846 : 10 : write(msg,'(a,i0)')' Bi-Lanczos algorithm with MAX number of iterations: ',BSp%niter
1847 : 10 : call wrtout(std_out,msg)
1848 : : !
1849 : : ! Check for presence of the restart file.
1850 : 10 : can_restart=.FALSE.
1851 : :
1852 [ + - ]: 10 : if (BS_files%in_haydock_basename /= BSE_NOFILE) then
1853 [ # # ]: 0 : restart_file = strcat(BS_files%in_haydock_basename,tag_file)
1854 [ # # ]: 0 : if (file_exists(restart_file) ) then
1855 : 0 : can_restart=.TRUE.
1856 [ # # ]: 0 : msg = strcat(" Restarting Haydock calculation from file: ",restart_file)
1857 : 0 : call wrtout(std_out,msg)
1858 : 0 : call wrtout(ab_out,msg)
1859 : 0 : ABI_ERROR("Restart is not implemented")
1860 : : else
1861 : 0 : can_restart=.FALSE.
1862 : 0 : ABI_WARNING(strcat("Cannot find restart file: ",restart_file))
1863 : : end if
1864 : : end if
1865 : : !
1866 : : ! Open the file and writes basic dimensions and info.
1867 [ + - ]: 10 : if (my_rank==master) then
1868 [ + - ]: 10 : out_file = TRIM(BS_files%out_basename)//TRIM(tag_file)
1869 [ - + ]: 10 : if (open_file(out_file,msg,newunit=out_unt,form="unformatted") /= 0) then
1870 : 0 : ABI_ERROR(msg)
1871 : : end if
1872 : : ! write header TODO: standardize this part.
1873 : 10 : write(out_unt)hsize,Bsp%use_coupling,BSE_HAYD_IMEPS,nkets,Bsp%broad
1874 : : end if
1875 : : !
1876 : : ! Select the terminator for the continued fraction.
1877 : 10 : term_type=0 !; if (Bsp%hayd_term>0) term_type=2
1878 : 10 : call wrtout(std_out,sjoin("Using terminator type: ",itoa(term_type)))
1879 : : !
1880 : : ! Calculate green(w) for the different starting kets.
1881 [ + + + + ]: 60130 : green=czero
1882 [ + + ]: 70 : do iq=1,nkets
1883 [ - + + - : 180 : ABI_MALLOC(ket0,(hexc%hsize))
+ - + - -
+ - + ]
1884 [ + - - + : 5880 : ket0 = kets(:,iq)
- - + + ]
1885 : : !
1886 : 60 : niter_file=0
1887 : :
1888 : : if (can_restart) then
1889 : : ! call haydock_restart(BSp,restart_file,BSE_HAYD_IMEPS,iq,hsize,&
1890 : : !& niter_file,aa_file,bb_file,phi_np1_file,phi_n_file,comm)
1891 : : end if
1892 : : !
1893 [ - + + - : 180 : ABI_MALLOC(phi_nm1,(my_nt))
+ - + - -
+ - + ]
1894 [ - + + - : 120 : ABI_MALLOC(phi_n,(my_nt))
- + - + ]
1895 [ - + + - : 120 : ABI_MALLOC(phi_np1,(my_nt))
- + - + ]
1896 [ - + + - : 120 : ABI_MALLOC(phit_nm1,(my_nt))
- + - + ]
1897 [ - + + - : 120 : ABI_MALLOC(phit_n,(my_nt))
- + - + ]
1898 [ - + + - : 120 : ABI_MALLOC(phit_np1,(my_nt))
- + - + ]
1899 [ - + + - : 120 : ABI_MALLOC(hphi_n,(hexc%hsize))
- + - + ]
1900 [ - + + - : 120 : ABI_MALLOC(hphit_n,(hexc%hsize))
- + - + ]
1901 : : !
1902 : : ! TODO: Note the different convention used for the coefficients
1903 : : ! Should use the same convention in the Hermitian case.
1904 : 60 : niter_max = niter_file + Bsp%niter
1905 [ - + + - : 60180 : ABI_CALLOC(aa,(niter_max))
+ - + - -
+ - + +
+ ]
1906 [ - + + - : 60120 : ABI_CALLOC(bb,(niter_max))
- + - + +
+ ]
1907 [ - + + - : 60120 : ABI_CALLOC(cc,(niter_max))
- + - + +
+ ]
1908 : :
1909 : : if (niter_file==0) then ! Calculation from scratch.
1910 [ + - - + : 5880 : phi_nm1 = ket0(my_t1:my_t2)
- - + + ]
1911 [ + - - + : 5880 : phit_nm1 = ket0(my_t1:my_t2)
- - + + ]
1912 : 60 : norm = DZNRM2(hexc%hsize,ket0,1)
1913 [ + + ]: 5820 : phi_nm1=phi_nm1/norm
1914 [ + + ]: 5820 : phit_nm1=phit_nm1/norm
1915 : :
1916 : 60 : call hexc%matmul_elphon(phi_nm1,hphi_n,'N',ep_renorms)
1917 : 60 : call hexc%matmul_elphon(phit_nm1,hphit_n,'C',ep_renorms)
1918 : :
1919 : 60 : aa(1)=xdotc(my_nt,phit_nm1,1,hphi_n(my_t1:),1)
1920 : 60 : call xmpi_sum(aa(1:1),comm,ierr)
1921 : :
1922 [ + - - + : 5880 : phi_n = hphi_n(my_t1:my_t2) - aa(1)*phi_nm1
- - + + ]
1923 [ + - - + : 5880 : phit_n = hphit_n(my_t1:my_t2) - CONJG(aa(1))*phit_nm1
- - + + ]
1924 : :
1925 : 60 : bb(1)=xdotc(my_nt,phi_n,1,phi_n,1)
1926 : 60 : call xmpi_sum(bb(1:1),comm,ierr)
1927 : 60 : bb(1) = SQRT(bb(1))
1928 : :
1929 : 60 : cc(1)=xdotc(my_nt,phit_n,1,phi_n,1)
1930 : 60 : call xmpi_sum(cc(1:1),comm,ierr)
1931 : 60 : cc(1) = cc(1)/bb(1)
1932 : :
1933 [ + + ]: 5820 : phi_n = phi_n /bb(1)
1934 [ + + ]: 5820 : phit_n = phit_n /CONJG(cc(1))
1935 : 60 : niter_done=1 ! TODO Be careful here
1936 : :
1937 : : else ! Use the previously calculates a and b.
1938 : : niter_done=niter_file
1939 : : ABI_ERROR("Restart not coded")
1940 : : !aa(1:niter_done) = aa_file
1941 : : !bb(1:niter_done) = bb_file
1942 : : !phi_np1=phi_np1_file(my_t1:my_t2) ! Select the slice treated by this node.
1943 : : !phi_n =phi_n_file (my_t1:my_t2)
1944 : : end if
1945 : :
1946 [ - + ]: 60 : if (can_restart) then
1947 : 0 : ABI_FREE(aa_file)
1948 : : ABI_FREE(bb_file)
1949 : : ABI_FREE(cc_file)
1950 : : ABI_FREE(phi_np1_file)
1951 : : ABI_FREE(phi_n_file)
1952 : : end if
1953 : :
1954 : : ! This factor gives the correct results
1955 : 60 : factor = -nfact*(DZNRM2(hexc%hsize,ket0,1)**2)
1956 : :
1957 : : ! Which quantity should be checked for convergence?
1958 : 60 : check = (/.TRUE.,.TRUE./)
1959 [ - + ]: 60 : if (ABS(Bsp%haydock_tol(2)-one)<tol6) check = (/.TRUE. ,.FALSE./)
1960 [ - + ]: 60 : if (ABS(Bsp%haydock_tol(2)-two)<tol6) check = (/.FALSE.,.TRUE./)
1961 : : ! Create new frequencies "mirror" in negative range to add
1962 : : ! their contributions. Can be improved by computing only once
1963 : : ! zero frequency, but loosing clearness
1964 : 60 : n_all_omegas = 2*BSp%nomega
1965 : :
1966 [ - + + - : 180 : ABI_MALLOC(all_omegas,(n_all_omegas))
+ - + - -
+ - + ]
1967 : : ! Put all omegas with frequency > 0 in table
1968 [ + + ]: 60120 : all_omegas(BSp%nomega+1:n_all_omegas) = BSp%omega
1969 : : ! Put all omegas with frequency < 0
1970 : : ! Warning, the broadening must be kept positive
1971 [ + + ]: 60120 : all_omegas(1:BSp%nomega) = -DBLE(BSp%omega(BSp%nomega:1:-1)) + j_dpc*AIMAG(BSp%omega(BSp%nomega:1:-1))
1972 : :
1973 [ - + + - : 240 : ABI_MALLOC(green_temp,(n_all_omegas,nkets))
- + + - +
- + - - +
- + ]
1974 : :
1975 : : call haydock_bilanczos_optalgo(niter_done,niter_max,n_all_omegas,all_omegas,BSp%haydock_tol(1),check,hexc,hexc_i,&
1976 : : hsize,my_t1,my_t2,factor,term_type,ep_renorms,aa,bb,cc,ket0,ket0_hbar_norm,phi_nm1,phi_n,phi_np1,&
1977 : 60 : phit_nm1,phit_n,phit_np1,green_temp(:,iq),inn,is_converged,comm)
1978 : :
1979 : : ! Computing result from two ranges of frequencies
1980 : : ! The real part is added, the imaginary part is substracted
1981 [ + + ]: 60120 : green(:,iq) = green_temp(BSp%nomega+1:n_all_omegas,iq)+CONJG(green_temp(BSp%nomega:1:-1,iq))
1982 : :
1983 [ - + ]: 60 : ABI_FREE(all_omegas)
1984 [ - + ]: 60 : ABI_FREE(green_temp)
1985 : :
1986 : : ! Save the a"s and the b"s for possible restarting.
1987 : : ! 1) Info on the Q.
1988 : : ! 2) Number of iterations performed.
1989 : : ! 3) do iter=1,niter_performed
1990 : : ! aa(iter),bb(iter)
1991 : : ! end do
1992 : : ! 4) |n-1>
1993 : : ! |n>
1994 : : ! |n+1>
1995 : : !
1996 [ + - ]: 60 : if (my_rank==master) then ! Open the file and writes basic dimensions and info.
1997 : 60 : write(out_unt)Bsp%q(:,iq)
1998 : 60 : write(out_unt)MIN(inn,niter_max) ! NB: if the previous loop completed inn=niter_max+1
1999 [ + + ]: 3276 : do itt=1,MIN(inn,niter_max) ! if we exited then inn is not incremented by one.
2000 : 3276 : write(out_unt)itt,aa(itt),bb(itt)
2001 : : end do
2002 : : end if
2003 : : !
2004 : : ! cbuff is used as workspace to gather |n-1>, |n> and |n+1>.
2005 [ - + + - : 180 : ABI_MALLOC(cbuff,(hsize))
+ - + - -
+ - + -
+ ]
2006 [ + + + + ]: 11580 : cbuff=czero; cbuff(my_t1:my_t2) = phi_nm1
2007 : 60 : call xmpi_sum_master(cbuff,master,comm,ierr)
2008 [ + - ]: 60 : if (my_rank==master) write(out_unt) cbuff ! |n-1>
2009 : :
2010 [ + + + + ]: 11580 : cbuff=czero; cbuff(my_t1:my_t2) = phi_n
2011 : 60 : call xmpi_sum_master(cbuff,master,comm,ierr)
2012 [ + - ]: 60 : if (my_rank==master) write(out_unt) cbuff ! |n>
2013 : :
2014 [ + + + + ]: 11580 : cbuff=czero; cbuff(my_t1:my_t2) = phi_np1
2015 : 60 : call xmpi_sum_master(cbuff,master,comm,ierr)
2016 [ + - ]: 60 : if (my_rank==master) write(out_unt) cbuff ! |n+1>
2017 : :
2018 [ - + ]: 60 : ABI_FREE(phi_nm1)
2019 [ - + ]: 60 : ABI_FREE(phi_n)
2020 [ - + ]: 60 : ABI_FREE(phi_np1)
2021 [ - + ]: 60 : ABI_FREE(phit_nm1)
2022 [ - + ]: 60 : ABI_FREE(phit_n)
2023 [ - + ]: 60 : ABI_FREE(phit_np1)
2024 [ - + ]: 60 : ABI_FREE(hphi_n)
2025 [ - + ]: 60 : ABI_FREE(hphit_n)
2026 [ - + ]: 60 : ABI_FREE(cbuff)
2027 : 60 : ABI_FREE(aa)
2028 : 60 : ABI_FREE(bb)
2029 : 60 : ABI_FREE(cc)
2030 [ - + ]: 70 : ABI_FREE(ket0)
2031 : : end do ! iq
2032 : :
2033 [ + - ]: 10 : if (my_rank==master) close(out_unt)
2034 : 10 : call xmpi_barrier(comm)
2035 : :
2036 [ - + ]: 10 : end subroutine haydock_bilanczos
2037 : : !!***
2038 : :
2039 : : !----------------------------------------------------------------------
2040 : :
2041 : : !!****f* m_haydock/haydock_bilanczos_optalgo
2042 : : !! NAME
2043 : : !! haydock_bilanczos_optalgo
2044 : : !!
2045 : : !! FUNCTION
2046 : : !! Haydock algorithm for general matrix
2047 : : !!
2048 : : !! INPUTS
2049 : : !! niter_done=Number of iterations already performed (0 if the run starts from scratch).
2050 : : !! niter_tot=Max number of iterations. Always > niter_done
2051 : : !! nomega=Number of Frequency points for the evaluation of the matrix element.
2052 : : !! omega(nomega)=Frequency set (imaginary part is already included).
2053 : : !! tol_iter=Tollerance used to stop the the algorithm.
2054 : : !! check(2)=Logical flags to specify where both the real and the imaginary part of the
2055 : : !! matrix elements of the Green functions have to be checked for convergence.
2056 : : !! hsize=Size of the blocks.
2057 : : !! my_t1,my_t2=Indeces of the first and last column stored treated by this done.
2058 : : !! term_type=0 if no terminator is used, 1 otherwise.
2059 : : !! hreso(hsize,my_t1:my_t2)=The columns of the resonant block.
2060 : : !! hcoup(hsize,my_t1:my_t2)=The columns of the coupling block.
2061 : : !! factor
2062 : : !! comm=MPI communicator.
2063 : : !!
2064 : : !! OUTPUT
2065 : : !! green(nomega)=Output matrix elements.
2066 : : !! inn=Last iteration performed.
2067 : : !! is_converged=.TRUE. of the algorithm converged.
2068 : : !!
2069 : : !! SIDE EFFECTS
2070 : : !! phi_nm1(my_t2-my_t1+1), phi_n(my_t2-my_t1+1)
2071 : : !! input: vectors used to initialize the iteration
2072 : : !! output: the vectors obtained in the last iteration
2073 : : !! aa(niter_tot) and bb(niter_tot+1)
2074 : : !! if niter_done>0: aa(1:niter_done), bb(1:niter_done) store the coefficients of the previous run.
2075 : : !! when the routine returns aa(1:inn) and bb(1:inn) contain the matrix elements of the tridiagonal form.
2076 : : !! cc(niter_tot+1)
2077 : : !!
2078 : : !! SOURCE
2079 : :
2080 : 60 : subroutine haydock_bilanczos_optalgo(niter_done,niter_tot,nomega,omega,tol_iter,check,hexc,hexc_i,hsize,my_t1,my_t2,&
2081 : 60 : & factor,term_type,ep_renorms,aa,bb,cc,ket0,ket0_hbar_norm,phi_nm1,phi_n,phi_np1,phit_nm1,phit_n,phit_np1,&
2082 : 60 : & green,inn,is_converged,comm)
2083 : :
2084 : : !Arguments ------------------------------------
2085 : : !scalars
2086 : : integer,intent(in) :: niter_tot,niter_done,nomega,comm,hsize,my_t1,my_t2,term_type
2087 : : integer,intent(out) :: inn
2088 : : logical,intent(out) :: is_converged
2089 : : real(dp),intent(in) :: tol_iter,ket0_hbar_norm
2090 : : complex(dp),intent(in) :: factor
2091 : : type(hexc_t),intent(in) :: hexc
2092 : : type(hexc_interp_t),intent(in) :: hexc_i
2093 : : !arrays
2094 : : complex(dp),intent(inout) :: bb(niter_tot+1)
2095 : : complex(dp),intent(out) :: green(nomega)
2096 : : complex(dp),intent(in) :: omega(nomega)
2097 : : complex(dp),intent(inout) :: aa(niter_tot),cc(niter_tot+1)
2098 : : complex(dp),intent(in) :: ket0(my_t2-my_t1+1)
2099 : : complex(dp),intent(in) :: ep_renorms(hsize)
2100 : : complex(dp),intent(inout) :: phi_nm1(my_t2-my_t1+1)
2101 : : complex(dp),intent(inout) :: phi_n (my_t2-my_t1+1)
2102 : : complex(dp),intent(inout) :: phi_np1(my_t2-my_t1+1)
2103 : : complex(dp),intent(inout) :: phit_nm1(my_t2-my_t1+1)
2104 : : complex(dp),intent(inout) :: phit_n (my_t2-my_t1+1)
2105 : : complex(dp),intent(inout) :: phit_np1(my_t2-my_t1+1)
2106 : : logical,intent(in) :: check(2)
2107 : :
2108 : : !Local variables ------------------------------
2109 : : !scalars
2110 : : integer :: my_nt,niter_min,nconv !,ierr
2111 : : character(len=500) :: msg
2112 : : logical :: keep_vectors=.TRUE.
2113 : : !arrays
2114 : : integer :: units(2)
2115 : 120 : real(dp) :: abs_err(nomega,2) !,ww_err(nomega,2)
2116 : 60 : complex(dp),allocatable :: oldg(:),newg(:)
2117 : 60 : complex(dp),allocatable :: hphi_np1(:),hphit_np1(:),save_phi(:,:),save_phit(:,:)
2118 : 60 : complex(dp),allocatable :: g00(:)
2119 : : logical :: test(2)
2120 : : integer :: ierr
2121 : : !************************************************************************
2122 : :
2123 : : ABI_UNUSED(ket0_hbar_norm)
2124 : : ABI_UNUSED(ket0(1))
2125 : : ABI_UNUSED(hexc_i%hsize_dense)
2126 [ + + ]: 180 : units = [std_out, ab_out]
2127 : :
2128 : 60 : my_nt = my_t2-my_t1+1
2129 : :
2130 [ - + + - : 180 : ABI_MALLOC(oldg,(nomega))
+ - + - -
+ - + ]
2131 [ - + + - : 120 : ABI_MALLOC(newg,(nomega))
- + - + ]
2132 [ - + + - : 120 : ABI_MALLOC(g00,(nomega))
- + - + ]
2133 [ + + + + : 360420 : oldg=czero; newg=czero; g00=czero
+ + ]
2134 : 60 : nconv=0
2135 : :
2136 [ + - - + ]: 60 : keep_vectors = (keep_vectors.and.xmpi_comm_size(comm)==1)
2137 [ + - ]: 60 : if (keep_vectors) then
2138 [ + - - + : 300 : ABI_MALLOC(save_phi,(my_t2-my_t1+1,niter_tot))
+ - - + +
- + - + -
- + - + ]
2139 [ - + + - : 180 : ABI_MALLOC_OR_DIE(save_phit,(my_t2-my_t1+1,niter_tot),ierr)
- + + - +
- - + -
+ ]
2140 [ + + + + ]: 5820060 : save_phi=czero
2141 [ + + + + ]: 5820060 : save_phit=czero
2142 : : end if
2143 : :
2144 [ - + + - : 180 : ABI_MALLOC_OR_DIE(hphi_np1,(hexc%hsize),ierr)
+ - + - +
- - + -
+ ]
2145 [ - + + - : 180 : ABI_MALLOC_OR_DIE(hphit_np1,(hexc%hsize),ierr)
+ - + - +
- - + -
+ ]
2146 : :
2147 [ + - ]: 3156 : do inn=niter_done+1,niter_tot
2148 : :
2149 : : !|n+1> = H |n> using all eh components.
2150 : 3156 : call hexc%matmul_elphon(phi_n, hphi_np1, 'N', ep_renorms)
2151 : 3156 : call hexc%matmul_elphon(phit_n, hphit_np1, 'C', ep_renorms)
2152 : :
2153 : : ! a(n) = < phit_n | H | phi_n >
2154 : 3156 : aa(inn)=xdotc(my_nt,phit_n,1,hphi_np1(my_t1:),1)
2155 : 3156 : call xmpi_sum(aa(inn),comm,ierr)
2156 : :
2157 : : ! |n+1> = |n+1> - a(n)|Vn> - c(n)|n-1>
2158 [ + + ]: 306132 : phi_np1 = hphi_np1(my_t1:my_t2) - aa(inn)*phi_n - cc(inn-1)*phi_nm1
2159 [ + + ]: 306132 : phit_np1 = hphit_np1(my_t1:my_t2) - CONJG(aa(inn))*phit_n - CONJG(bb(inn-1))*phit_nm1
2160 : :
2161 : 3156 : bb(inn) = xdotc(my_nt,phi_np1,1,phi_np1,1)
2162 : 3156 : call xmpi_sum(bb(inn),comm,ierr)
2163 : 3156 : bb(inn) = SQRT(bb(inn))
2164 : :
2165 : 3156 : cc(inn) = xdotc(my_nt,phit_np1,1,phi_np1,1)
2166 : 3156 : call xmpi_sum(cc(inn),comm,ierr)
2167 : 3156 : cc(inn) = cc(inn)/bb(inn)
2168 : :
2169 [ + + ]: 306132 : phi_np1 = phi_np1 / bb(inn)
2170 [ + + ]: 306132 : phit_np1 = phit_np1 / CONJG(cc(inn))
2171 : :
2172 : : ! |n-1> = |n>
2173 : : ! |n> = |n+1>
2174 [ + + ]: 306132 : phi_nm1 = phi_n
2175 [ + + ]: 306132 : phi_n = phi_np1
2176 [ + + ]: 306132 : phit_nm1 = phit_n
2177 [ + + ]: 306132 : phit_n = phit_np1
2178 : :
2179 [ + - ]: 3156 : if (keep_vectors) then
2180 [ + + ]: 306132 : save_phi(:,inn) = phi_n
2181 [ + + ]: 306132 : save_phit(:,inn) = phit_n
2182 : : end if
2183 : 3156 : write(msg,'(a,i0,a,3es12.4)')' Iteration number ',inn,', b_i RE(c_i) IM(c_i) ',REAL(bb(inn)),REAL(cc(inn)),AIMAG(cc(inn))
2184 : 3156 : call wrtout(std_out,msg)
2185 : :
2186 : 3156 : call continued_fract_general(inn,term_type,aa,bb,cc,nomega,omega,g00)
2187 [ - + + + ]: 6324624 : newg = factor*g00
2188 : : !
2189 : : ! Avoid spurious convergence.
2190 [ - + ]: 3156 : niter_min=4; if (niter_done>1) niter_min=niter_done+1
2191 [ + + ]: 3156 : if (inn>niter_min) then
2192 [ + + ]: 8928 : test=.TRUE.
2193 [ + + ]: 5960928 : abs_err(:,1) = ABS(DBLE (newg-oldg))
2194 [ + + ]: 5960928 : abs_err(:,2) = ABS(AIMAG(newg-oldg))
2195 : :
2196 [ - + ]: 2976 : if (tol_iter>zero) then
2197 : : ! Test on the L1 norm.
2198 [ # # # # : 0 : if (check(1)) test(1) = SUM(abs_err(:,1)) < tol_iter*SUM(ABS(DBLE (newg)))
# # ]
2199 [ # # # # : 0 : if (check(2)) test(2) = SUM(abs_err(:,2)) < tol_iter*SUM(ABS(AIMAG(newg)))
# # ]
2200 : : else
2201 : : ! Stringent test for each point.
2202 [ + - + + : 4133190 : if (check(1)) test(1) = ALL( abs_err(:,1) < -tol_iter*ABS(DBLE (newg)))
+ + ]
2203 [ + - + + : 4102776 : if (check(2)) test(2) = ALL( abs_err(:,2) < -tol_iter*ABS(AIMAG(newg)))
+ + ]
2204 : : end if
2205 : : !
2206 [ + + + + : 3469 : if (ALL(test)) then
+ + ]
2207 : 186 : nconv = nconv+1
2208 : : else
2209 : : nconv = 0
2210 : : end if
2211 [ + + ]: 186 : if (nconv==2) then
2212 [ + - ]: 60 : if(inn<100)then
2213 : : write(msg,'(a,es10.2,a)')&
2214 : 60 : " >>> Haydock algorithm converged twice within haydock_tol= ",tol_iter," after less than 100 iterations."
2215 : : else
2216 : : write(msg,'(a,es10.2,a,i0,a)')&
2217 : 0 : " >>> Haydock algorithm converged twice within haydock_tol= ",tol_iter," after ",inn," iterations."
2218 : : endif
2219 : 60 : call wrtout(units, msg)
2220 : 60 : EXIT
2221 : : end if
2222 : : end if
2223 : : !
2224 [ - + + + ]: 6213852 : oldg = newg
2225 : : end do ! inn
2226 : :
2227 [ + + ]: 120180 : green = newg
2228 [ - + ]: 60 : if (nconv/=2) then
2229 : : write(msg,'(a,es10.2,a,i0,a)')&
2230 : 0 : " WARNING: Haydock algorithm did not converge within ",tol_iter," after ",niter_tot," iterations."
2231 : 0 : call wrtout(units, msg)
2232 : : end if
2233 : :
2234 : 60 : is_converged = (nconv==2)
2235 : :
2236 [ - + ]: 60 : ABI_FREE(oldg)
2237 [ - + ]: 60 : ABI_FREE(newg)
2238 [ - + ]: 60 : ABI_FREE(g00)
2239 [ - + ]: 60 : ABI_FREE(hphi_np1)
2240 [ - + ]: 60 : ABI_FREE(hphit_np1)
2241 : :
2242 [ + - ]: 60 : ABI_SFREE(save_phi)
2243 [ + - ]: 60 : ABI_SFREE(save_phit)
2244 : :
2245 : : !! if (keep_vectors) then
2246 : : !! tdim = MIN(inn,niter_tot)
2247 : : !! ABI_MALLOC(ovlp,(tdim,tdim))
2248 : :
2249 : : !! ABI_MALLOC(phi_test,(hsize))
2250 : : !! ABI_MALLOC(phi_test2,(hsize))
2251 : :
2252 : : !! max_err=smallest_real; mean_err=zero; mean_err2=zero; row_max=-1
2253 : : !! do ii=1,tdim
2254 : : !! parity = (-1)**(ii+1)
2255 : : !! phi_test = save_phi(:,ii)
2256 : : !! call hexc%matmul_full(hexc_i, phi_test, phi_test2, parity)
2257 : : !! !phi_test2 = MATMUL(hreso,phi_test) + parity * MATMUL(hcoup,CONJG(phi_test))
2258 : : !! ovlp(ii,ii) = DOT_PRODUCT(phi_test,phi_test2) + DOT_PRODUCT(phi_test2,phi_test)
2259 : : !! err = ABS(ovlp(ii,ii)-cone)
2260 : : !! mean_err = mean_err + err
2261 : : !! mean_err2 = mean_err2 + err**2
2262 : : !! if (err > max_err) then
2263 : : !! max_err = err
2264 : : !! row_max = ii
2265 : : !! end if
2266 : : !! end do
2267 : : !! mean_err = mean_err/tdim
2268 : : !! std_dev = mean_err2/tdim -mean_err**2
2269 : : !! write(std_out,'(a,i0,1x,3es14.6)')&
2270 : : !!& " Error in normalization (ii, max_err,mean,std_dev): ",row_max,max_err,mean_err,std_dev
2271 : :
2272 : : !! ABI_FREE(phi_test)
2273 : : !! ABI_FREE(phi_test2)
2274 : : !!
2275 : : !! ABI_MALLOC(alpha,(hsize,tdim))
2276 : :
2277 : : !! ! Less efficient but for sake of simplicity with hexc_matmul
2278 : : !! ! TODO possibility to call hreso * phi, and hcoup * phi separately
2279 : : !! do ii=1,tdim
2280 : : !! parity = (-1)**(ii+1)
2281 : : !! call hexc%matmul_full(hexc_i, save_phi(:,ii), alpha(:,ii), parity)
2282 : : !! end do
2283 : :
2284 : : !! !alpha = MATMUL(hreso,save_phi(:,1:tdim))
2285 : : !! !
2286 : : !! !do ii=1,tdim
2287 : : !! ! parity = (-1)**(ii+1)
2288 : : !! ! alpha(:,ii) = alpha(:,ii) + parity*MATMUL(hcoup,CONJG(save_phi(:,ii)))
2289 : : !! !end do
2290 : :
2291 : : !! ovlp = MATMUL(TRANSPOSE(CONJG(save_phi(:,1:tdim))),alpha)
2292 : :
2293 : : !! ABI_MALLOC(beta,(hsize,tdim))
2294 : : !! do ii=1,tdim
2295 : : !! parity = (-1)**(ii+1)
2296 : : !! beta(:,ii) = parity*save_phi(:,ii)
2297 : : !! alpha(:,ii) = -parity*alpha(:,ii)
2298 : : !! end do
2299 : :
2300 : : !! ovlp = ovlp - MATMUL(TRANSPOSE(CONJG(beta)),alpha)
2301 : :
2302 : : !! max_err=smallest_real; row_max=-1; col_max=-1
2303 : : !! mean_err=zero; mean_err2=zero
2304 : : !! do jj=1,tdim
2305 : : !! do ii=1,jj
2306 : : !! err = ABS(ovlp(ii,jj))
2307 : : !! if (ii==jj) err = ABS(err - one)
2308 : : !! mean_err = mean_err + err
2309 : : !! mean_err2 = mean_err2 + err**2
2310 : : !! if (err > max_err) then
2311 : : !! max_err = err
2312 : : !! row_max=ii
2313 : : !! col_max=jj
2314 : : !! end if
2315 : : !! end do
2316 : : !! end do
2317 : :
2318 : : !! mean_err = mean_err/(tdim*(tdim+1)/2)
2319 : : !! std_dev = mean_err2/(tdim*(tdim+1)/2) - mean_err**2
2320 : : !! write(std_out,'(a,2(i0,1x),3es14.6)')&
2321 : : !! " Error in Hbar-ortho (i,j), max_err, mean, std_dev ",row_max,col_max,max_err,mean_err,std_dev
2322 : : !! !call print_arr(ovlp,max_r=185,max_c=10,unit=std_out)
2323 : :
2324 : : !! ABI_FREE(alpha)
2325 : : !! ABI_FREE(beta)
2326 : : !! ABI_FREE(ovlp)
2327 : : !! ABI_FREE(save_phi)
2328 : : !! end if
2329 : :
2330 : 60 : end subroutine haydock_bilanczos_optalgo
2331 : : !!***
2332 : :
2333 : :
2334 : : !----------------------------------------------------------------------
2335 : :
2336 : : !!****f* m_numeric_tools/continued_fract_general
2337 : : !! NAME
2338 : : !! continued_fract
2339 : : !!
2340 : : !! FUNCTION
2341 : : !! This routine calculates the continued fraction:
2342 : : !!
2343 : : !! 1
2344 : : !! f(z) = _______________________________
2345 : : !! z - a1 - b1^2
2346 : : !! _____________________
2347 : : !! z - a2 - b2^2
2348 : : !! ___________
2349 : : !! z -a3 - ........
2350 : : !!
2351 : : !! INPUTS
2352 : : !! nlev=Number of "levels" in the continued fraction.
2353 : : !! term_type=Type of the terminator.
2354 : : !! 0 --> No terminator.
2355 : : !! -1 --> Assume constant coefficients for a_i and b_i for i>nlev with a_inf = a(nlev) and b_inf = b(nleb)
2356 : : !! 1 --> Same as above but a_inf and b_inf are obtained by averaging over the nlev values.
2357 : : !! aa(nlev)=Set of a_i coefficients.
2358 : : !! bb(nlev)=Set of b_i coefficients.
2359 : : !! nz=Number of points on the z-mesh.
2360 : : !! zpts(nz)=z-mesh.
2361 : : !!
2362 : : !! OUTPUT
2363 : : !! spectrum(nz)=Contains f(z) on the input mesh.
2364 : : !!
2365 : : !! SOURCE
2366 : :
2367 : 3156 : subroutine continued_fract_general(nlev,term_type,aa,bb,cc,nz,zpts,spectrum)
2368 : :
2369 : : !Arguments ------------------------------------
2370 : : !scalars
2371 : : integer,intent(in) :: nlev,term_type,nz
2372 : : !arrays
2373 : : complex(dp),intent(in) :: bb(nlev)
2374 : : complex(dp),intent(in) :: cc(nlev)
2375 : : complex(dp),intent(in) :: aa(nlev)
2376 : : complex(dp),intent(in) :: zpts(nz)
2377 : : complex(dp),intent(out) :: spectrum(nz)
2378 : :
2379 : : !Local variables ------------------------------
2380 : : !scalars
2381 : : integer :: it
2382 : : complex(dp) :: bb_inf,bg,bu,swap, aa_inf
2383 : : character(len=500) :: msg
2384 : : !arrays
2385 : 3156 : complex(dp),allocatable :: div(:),den(:)
2386 : : !************************************************************************
2387 : :
2388 [ - + + - : 9468 : ABI_MALLOC(div,(nz))
+ - + - -
+ - + ]
2389 [ - + + - : 6312 : ABI_MALLOC(den,(nz))
- + - + ]
2390 : :
2391 : 3156 : select case (term_type)
2392 : : case (0) ! No terminator.
2393 [ + + ]: 6321468 : div=czero
2394 : : case (-1,1)
2395 : 0 : ABI_ERROR("Not yet implemented")
2396 [ # # ]: 0 : if (term_type==-1) then
2397 : 0 : bb_inf=bb(nlev)
2398 : 0 : aa_inf=aa(nlev)
2399 : : else
2400 [ # # ]: 0 : bb_inf=SUM(bb)/nlev
2401 [ # # ]: 0 : aa_inf=SUM(aa)/nlev
2402 : : end if
2403 : : ! Be careful with the sign of the SQRT.
2404 [ # # ]: 0 : div(:) = half*(bb(nlev)/(bb_inf))**2 * ( zpts-aa_inf - SQRT((zpts-aa_inf)**2 - four*bb_inf**2) )
2405 : : case (2)
2406 : 0 : ABI_ERROR("Not yet implemented")
2407 [ # # ]: 0 : div = zero
2408 [ # # ]: 0 : if (nlev>4) then
2409 : 0 : bg=zero; bu=zero
2410 [ # # ]: 0 : do it=1,nlev,2
2411 [ # # ]: 0 : if (it+2<nlev) bg = bg + bb(it+2)
2412 [ # # ]: 0 : bu = bu + bb(it)
2413 : : end do
2414 : 0 : bg = bg/(nlev/2+MOD(nlev,2))
2415 : 0 : bu = bg/((nlev+1)/2)
2416 : : !if (iseven(nlev)) then
2417 [ # # ]: 0 : if (.not.iseven(nlev)) then
2418 : 0 : swap = bg
2419 : 0 : bg = bu
2420 : 0 : bu = bg
2421 : : end if
2422 : : !write(std_out,*)nlev,bg,bu
2423 : : !Here be careful with the sign of SQRT
2424 [ # # ]: 0 : do it=1,nz
2425 : : div(it) = half/zpts(it) * (bb(nlev)/bu)**2 * &
2426 : 0 : ( (zpts(it)**2 +bu**2 -bg**2) - SQRT( (zpts(it)**2+bu**2-bg**2)**2 -four*(zpts(it)*bu)**2) )
2427 : : end do
2428 : : end if
2429 : :
2430 : : case default
2431 : 0 : write(msg,'(a,i0)')" Wrong value for term_type: ",term_type
2432 [ + - - - ]: 3156 : ABI_ERROR(msg)
2433 : : end select
2434 : :
2435 [ + + ]: 89458 : do it=nlev,2,-1
2436 [ + + ]: 172862906 : den(:) = zpts(:) - aa(it) - div(:)
2437 [ + + ]: 172866062 : div(:) = (bb(it-1)*cc(it-1) )/ den(:)
2438 : : end do
2439 : :
2440 [ - + + + ]: 6324624 : den = zpts(:) - aa(1) - div(:)
2441 [ - + + + ]: 6324624 : div = one/den(:)
2442 : :
2443 [ + + ]: 6321468 : spectrum = div
2444 [ - + ]: 3156 : ABI_FREE(div)
2445 [ - + ]: 3156 : ABI_FREE(den)
2446 : :
2447 : 3156 : end subroutine continued_fract_general
2448 : : !!***
2449 : :
2450 : : !----------------------------------------------------------------------
2451 : :
2452 [ # # ]: 0 : end module m_haydock
2453 : : !!***
|
