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 : !!***
|
