TABLE OF CONTENTS


ABINIT/m_haydock [ Modules ]

[ Top ] [ Modules ]

NAME

 m_haydock

FUNCTION

COPYRIGHT

  Copyright (C) 2008-2022 ABINIT group (M.Giantomassi, Y. Gillet, L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

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_nctk
29  use m_haydock_io
30  use m_linalg_interfaces
31  use m_ebands
32  use m_hdr
33  use netcdf
34 
35  use m_time,              only : timab
36  use m_fstrings,          only : strcat, sjoin, itoa, int2char4
37  use m_io_tools,          only : file_exists, open_file
38  use defs_datatypes,      only : ebands_t, pseudopotential_type
39  use m_geometry,          only : normv
40  use m_hide_blas,         only : xdotc, xgemv
41  use m_hide_lapack,       only : matrginv
42  use m_numeric_tools,     only : print_arr, symmetrize, hermitianize, continued_fract, wrap2_pmhalf, iseven
43  use m_kpts,              only : listkk
44  use m_crystal,           only : crystal_t
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_init, hexc_interp_init, hexc_free, hexc_interp_free, &
53 &                                hexc_build_hinterp, hexc_matmul_tda, hexc_matmul_full, hexc_t, hexc_matmul_elphon, hexc_interp_t
54  use m_exc_spectra,       only : exc_write_data, exc_eps_rpa, exc_write_tensor, mdfs_ncwrite
55  use m_eprenorms,         only : eprenorms_t, renorm_bst
56  use m_wfd_optic,         only : calc_optical_mels
57 
58  implicit none
59 
60  private

m_haydock/exc_haydock_driver [ Functions ]

[ Top ] [ m_haydock ] [ Functions ]

NAME

 exc_haydock_driver

FUNCTION

  Calculate the imaginary part of the macroscopic dielectric function with the Haydock recursive method.

INPUTS

 BSp<type(excparam)=The parameter for the Bethe-Salpeter run.
 BS_files<excparam>=Files associated to the bethe_salpeter code.
 Cryst<crystal_t>=Info on the crystalline structure.
 Kmesh<type(kmesh_t)>=The list of k-points in the BZ, IBZ and symmetry tables.
 Cryst<type(crystal_t)>=Info on the crystalline structure.
 Hdr_bse
 KS_BSt=The KS energies.
 QP_BSt=The QP energies.
 Wfd<wfdgw_t>=Wavefunction descriptor (input k-mesh)
 Psps <type(pseudopotential_type)>=variables related to pseudopotentials.
 Pawtab(Cryst%ntypat*usepaw)<pawtab_type>=PAW tabulated starting data.
 Hur(Cryst%natom*usepaw)<type(pawhur_t)>=Only for PAW and DFT+U, quantities used to evaluate the commutator [H_u,r].

OUTPUT

  The imaginary part of the macroscopic dielectric function is written on the external file _EXC_MDF

SOURCE

 94 subroutine exc_haydock_driver(BSp,BS_files,Cryst,Kmesh,Hdr_bse,KS_BSt,QP_Bst,Wfd,Psps,Pawtab,Hur,Epren, &
 95  Kmesh_dense, KS_BSt_dense, QP_BSt_dense, Wfd_dense, Vcp_dense, grid) ! Optional args
 96 
 97 !Arguments ------------------------------------
 98 !scalars
 99  type(excparam),intent(in) :: BSp
100  type(excfiles),intent(in) :: BS_files
101  type(kmesh_t),intent(in) :: Kmesh
102  type(crystal_t),intent(in) :: Cryst
103  type(Hdr_type),intent(in) :: Hdr_bse
104  type(wfdgw_t),intent(inout) :: Wfd
105  type(pseudopotential_type),intent(in) :: Psps
106  type(ebands_t),intent(in) :: KS_BSt,QP_Bst
107  type(double_grid_t),intent(in),optional :: grid
108  type(kmesh_t),intent(in),optional :: Kmesh_dense
109  type(wfdgw_t),intent(inout),optional :: Wfd_dense
110  type(ebands_t),intent(in),optional :: KS_BSt_dense, QP_Bst_dense
111  type(vcoul_t),intent(in),optional :: Vcp_dense
112  type(eprenorms_t),intent(in) :: Epren
113 !arrays
114  type(pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd%usepaw)
115  type(pawhur_t),intent(in) :: Hur(Cryst%natom*Wfd%usepaw)
116 
117 !Local variables ------------------------------
118 !scalars
119  integer,parameter :: master=0
120  integer :: io,my_rank,iq,itt,ierr
121  integer :: hsize,comm,my_t1,my_t2,nsppol,nkets,nproc,ncid
122  integer :: spin,spad,ik_bz,iv,ic,trans_idx,lomo_min,max_band
123  real(dp) :: omegaev,rand_phi !,norm
124  complex(dpc) :: ks_avg,gw_avg,exc_avg
125  logical :: use_mpio,prtdos
126  character(len=500) :: msg
127  type(hexc_t) :: hexc
128  type(hexc_interp_t) :: hexc_i
129 !arrays
130  real(dp) :: tsec(2)
131  real(dp),allocatable :: dos(:),dos_gw(:),dos_ks(:)
132  complex(dpc),allocatable :: green(:,:)
133  complex(dpc),allocatable :: opt_cvk(:,:,:,:,:),kets(:,:)
134  complex(dpc),allocatable :: eps_rpanlf(:,:),eps_gwnlf(:,:)
135  complex(dpc),allocatable :: tensor_cart(:,:),tensor_cart_rpanlf(:,:),tensor_cart_gwnlf(:,:)
136  complex(dpc),allocatable :: tensor_red(:,:),tensor_red_rpanlf(:,:),tensor_red_gwnlf(:,:)
137 
138  !Temperature
139  real(dp) :: dksqmax, en
140  integer,allocatable :: bs2eph(:,:)
141  integer :: sppoldbl, timrev
142  logical :: do_ep_renorm, do_ep_lifetime
143  integer :: ntemp
144  character(len=4) :: ts
145  character(len=fnlen) :: prefix
146  character(len=fnlen) :: path
147  complex(dpc),allocatable :: ep_renorms(:)
148  integer :: ep_ik, ik, ireh, isppol
149  integer :: itemp
150  type(ebands_t) :: EPBSt, EP_QPBSt
151 
152 !************************************************************************
153 
154  call timab(690,1,tsec) ! exc_haydock_driver
155  call timab(691,1,tsec) ! exc_haydock_driver(read)
156 
157  if (BSp%have_complex_ene) then
158    ABI_ERROR("Complex energies are not supported yet")
159  end if
160 
161  my_rank = Wfd%my_rank
162  comm    = Wfd%comm
163  nsppol  = Wfd%nsppol
164  nproc   = Wfd%nproc
165 
166  use_mpio=.FALSE.
167 #ifdef HAVE_MPI_IO
168  use_mpio = (nproc > 1)
169  !use_mpio = .TRUE.
170 #endif
171  use_mpio=.FALSE.
172  !use_mpio = .TRUE.
173 
174  ! Hsize refers to the size of the individual blocks (resonant and coupling).
175  ! Thanks to the symmetry property of the starting vector, the Haydock method
176  ! can be reformulated in terms of matrix-vector multiplication involving the
177  ! blocks thus avoiding to allocation of the full matrix ( R   C )
178  !                                                        -C* -R*)
179  hsize=SUM(BSp%nreh)
180 
181  !YG2014
182  call hexc_init(hexc, BSp, BS_files, Cryst, Kmesh, Wfd, KS_BSt, QP_BSt, comm)
183 
184  !YG2014
185  if(BSp%use_interp) then
186    call hexc_interp_init(hexc_i, hexc, BSp%interp_m3_width, BSp%interp_method,&
187 &     Kmesh_dense, Vcp_dense, grid, Wfd_dense, &
188 &     KS_BSt_dense, QP_BSt_dense, Psps, Pawtab)
189 
190  end if
191 
192  call timab(691,2,tsec) ! exc_haydock_driver(read)
193  call timab(692,1,tsec) ! exc_haydock_driver(prep)
194  !
195  ! Prepare the starting vectors for the Lanczos chain.
196  nkets=Bsp%nq
197 
198  prtdos=.FALSE. !prtdos=.TRUE.
199  if (prtdos) then
200    nkets=nkets+1
201    if (Bsp%use_coupling>0) then
202      ABI_ERROR("DOS with coupling not coded")
203      nkets=nkets+1
204    end if
205  end if
206 
207  !YG2014
208  ABI_MALLOC_OR_DIE(kets,(hexc%hsize,nkets), ierr)
209  kets=czero
210  !
211  ! Prepare the kets for the macroscopic dielectric function.
212  lomo_min=Bsp%lomo_min; max_band=Bsp%nbnds
213 
214  !YG2014
215  ABI_MALLOC_OR_DIE(opt_cvk,(lomo_min:max_band,lomo_min:max_band,hexc%nbz,Wfd%nsppol,BSp%nq), ierr)
216 
217  do iq=1,Bsp%nq
218  ! Note KS_BSt is used here to calculate the commutator.
219    call calc_optical_mels(hexc%Wfd,hexc%Kmesh,hexc%KS_BSt,Cryst,Psps,Pawtab,Hur, &
220 &     BSp%inclvkb,BSp%lomo_spin,lomo_min,max_band,hexc%nbz,BSp%q(:,iq),opt_cvk(:,:,:,:,iq))
221 
222  ! Fill ket0 using the same ordering for the indeces as the one used for the excitonic Hamiltonian.
223  ! Note that only the resonant part is used here.
224    do spin=1,nsppol
225 
226      if(BSp%use_interp) then
227        spad=(spin-1)*BSp%nreh_interp(spin)
228      else
229        spad=(spin-1)*BSp%nreh(spin)
230      end if
231 
232      do ik_bz=1,hexc%nbz
233        do iv=BSp%lomo_spin(spin),BSp%homo_spin(spin)
234          do ic=BSp%lumo_spin(spin),BSp%nbnds
235 
236            if(BSp%use_interp) then
237              trans_idx = BSp%vcks2t_interp(iv,ic,ik_bz,spin)
238            else
239              trans_idx = BSp%vcks2t(iv,ic,ik_bz,spin)
240            end if
241 
242            if (trans_idx>0) kets(trans_idx+spad,iq)=opt_cvk(ic,iv,ik_bz,spin,iq)
243          end do
244        end do
245      end do
246    end do
247 
248  end do
249 
250  !
251  ! ========================================================
252  ! === Write the Optical Matrix Elements to NetCDF file ===
253  ! ========================================================
254 
255  !if (.false.) then
256  !  ome_fname='test_OME.nc'
257  !  call exc_write_optme(ome_fname,minb,maxb,BSp%nkbz,Wfd%nsppol,BSp%nq,opt_cvk,ierr)
258  !end if
259 
260  ! Free WFD descriptor, we don't need ur and ug anymore !
261  ! We make space for interpolated hamiltonian
262  call wfd%wave_free("All")
263  if(BSp%use_interp) call wfd_dense%wave_free("All")
264 
265  ! Build interpolated hamiltonian
266  if(BSp%use_interp) then
267    if (any(BSp%interp_mode == [2,3,4])) then
268      call hexc_build_hinterp(hexc, hexc_i)
269    end if
270  end if
271 
272  call timab(692,2,tsec) ! exc_haydock_driver(prep)
273  call timab(693,1,tsec) ! exc_haydock_driver(wo lf) - that is, without local field
274 
275  do_ep_renorm = .FALSE.
276  ntemp = 1
277  do_ep_lifetime = .FALSE.
278 
279  if(BSp%do_ep_renorm) then
280    if (BSp%nsppol == 2) then
281      ABI_ERROR('Elphon renorm with nsppol == 2 not yet coded !')
282    end if
283    do_ep_renorm = .TRUE.
284    ntemp = Epren%ntemp
285    if(BSp%do_lifetime) then
286      do_ep_lifetime = .TRUE.
287    end if
288 
289    ! Force elphon linewidth
290    do_ep_lifetime = .TRUE.
291 
292    ! Map points from BSE to elphon kpoints
293    sppoldbl = 1 !; if (any(Cryst%symafm == -1) .and. Epren%nsppol == 1) nsppoldbl=2
294    ABI_MALLOC(bs2eph, (Kmesh%nbz*sppoldbl, 6))
295    ABI_MALLOC(ep_renorms, (hsize))
296    timrev = 1
297    call listkk(dksqmax, Cryst%gmet, bs2eph, Epren%kpts, Kmesh%bz, Epren%nkpt, Kmesh%nbz, Cryst%nsym, &
298       sppoldbl, Cryst%symafm, Cryst%symrel, timrev, comm, use_symrec=.False.)
299  end if
300 
301  call timab(693,2,tsec) ! exc_haydock_driver(wo lf    - that is, without local field
302  call timab(694,1,tsec) ! exc_haydock_driver(apply
303 
304  prefix = ""
305  do itemp = 1, ntemp
306 
307    call ebands_copy(hexc%KS_BSt, EPBSt)
308    call ebands_copy(hexc%QP_BSt, EP_QPBSt)
309 
310    ! =================================================
311    ! == Calculate elphon vector in transition space ==
312    ! =================================================
313    if (do_ep_renorm) then
314 
315      ! Will perform elphon renormalization for itemp
316 
317      call int2char4(itemp,ts)
318      prefix = TRIM("_T") // ts
319 
320      ! No scissor with KSBST
321      call renorm_bst(Epren, EPBSt, Cryst, itemp, do_lifetime=.TRUE.,do_check=.TRUE.)
322 
323      call renorm_bst(Epren, EP_QPBSt, Cryst, itemp, do_lifetime=.TRUE.,do_check=.FALSE.)
324 
325      do isppol = 1, BSp%nsppol
326        do ireh = 1, BSp%nreh(isppol)
327          ic = BSp%Trans(ireh,isppol)%c
328          iv = BSp%Trans(ireh,isppol)%v
329          ik = BSp%Trans(ireh,isppol)%k ! In the full bz
330          en = BSp%Trans(ireh,isppol)%en
331 
332          ep_ik = bs2eph(ik,1)
333 
334          !TODO support multiple spins !
335          if(ABS(en - (Epren%eigens(ic,ep_ik,isppol)-Epren%eigens(iv,ep_ik,isppol)+BSp%mbpt_sciss)) > tol3) then
336            ABI_ERROR("Eigen from the transition does not correspond to the EP file !")
337          end if
338          ep_renorms(ireh) = (Epren%renorms(1,ic,ik,isppol,itemp) - Epren%renorms(1,iv,ik,isppol,itemp))
339 
340          ! Add linewith
341          if(do_ep_lifetime) then
342            ep_renorms(ireh) = ep_renorms(ireh) - j_dpc*(Epren%linewidth(1,ic,ik,isppol,itemp) +&
343 &                                                       Epren%linewidth(1,iv,ik,isppol,itemp))
344          end if
345 
346        end do
347      end do
348    end if
349 
350 
351    ! =======================================================
352    ! === Make EPS RPA and GW without local-field effects ===
353    ! =======================================================
354    ABI_MALLOC(eps_rpanlf,(BSp%nomega,BSp%nq))
355    ABI_MALLOC(dos_ks,(BSp%nomega))
356    ABI_MALLOC(eps_gwnlf ,(BSp%nomega,BSp%nq))
357    ABI_MALLOC(dos_gw,(BSp%nomega))
358 
359    call wrtout(std_out," Calculating RPA NLF and QP NLF epsilon","COLL")
360 
361    call exc_eps_rpa(BSp%nbnds,BSp%lomo_spin,BSp%lomo_min,BSp%homo_spin,hexc%Kmesh,EPBSt,BSp%nq,nsppol,&
362 &    opt_cvk,Cryst%ucvol,BSp%broad,BSp%nomega,BSp%omega,eps_rpanlf,dos_ks)
363 
364    call exc_eps_rpa(BSp%nbnds,BSp%lomo_spin,BSp%lomo_min,BSp%homo_spin,hexc%Kmesh,EP_QPBSt,BSp%nq,nsppol,&
365 &    opt_cvk,Cryst%ucvol,Bsp%broad,BSp%nomega,BSp%omega,eps_gwnlf,dos_gw)
366 
367    if (my_rank==master) then ! Only master works.
368      !
369      ! Master node writes final results on file.
370      call exc_write_data(BSp,BS_files,"RPA_NLF_MDF",eps_rpanlf,prefix=prefix,dos=dos_ks)
371 
372      call exc_write_data(BSp,BS_files,"GW_NLF_MDF",eps_gwnlf,prefix=prefix,dos=dos_gw)
373 
374      ! Computing and writing tensor in files
375 
376      ! RPA_NLF
377      ABI_MALLOC(tensor_cart_rpanlf,(BSp%nomega,6))
378      ABI_MALLOC(tensor_red_rpanlf,(BSp%nomega,6))
379 
380      call wrtout(std_out," Calculating RPA NLF dielectric tensor","COLL")
381      call haydock_mdf_to_tensor(BSp,Cryst,eps_rpanlf,tensor_cart_rpanlf, tensor_red_rpanlf, ierr)
382 
383      if(ierr == 0) then
384         ! Writing tensor
385         call exc_write_tensor(BSp,BS_files,"RPA_NLF_TSR_CART",tensor_cart_rpanlf)
386         call exc_write_tensor(BSp,BS_files,"RPA_NLF_TSR_RED",tensor_red_rpanlf)
387      else
388         write(msg,'(3a)')&
389 &         'The RPA_NLF dielectric complex tensor cannot be computed',ch10,&
390 &         'There must be 6 different q-points in long wavelength limit (see gw_nqlwl)'
391         ABI_COMMENT(msg)
392      end if
393 
394      ABI_FREE(tensor_cart_rpanlf)
395      ABI_FREE(tensor_red_rpanlf)
396 
397      ! GW_NLF
398      ABI_MALLOC(tensor_cart_gwnlf,(BSp%nomega,6))
399      ABI_MALLOC(tensor_red_gwnlf,(BSp%nomega,6))
400 
401      call wrtout(std_out," Calculating GW NLF dielectric tensor","COLL")
402 
403      call haydock_mdf_to_tensor(BSp,Cryst,eps_gwnlf,tensor_cart_gwnlf, tensor_red_gwnlf, ierr)
404 
405      if(ierr == 0) then
406         ! Writing tensor
407         call exc_write_tensor(BSp,BS_files,"GW_NLF_TSR_CART",tensor_cart_gwnlf)
408         call exc_write_tensor(BSp,BS_files,"GW_NLF_TSR_RED",tensor_red_gwnlf)
409      else
410         write(msg,'(3a)')&
411 &         'The GW_NLF dielectric complex tensor cannot be computed',ch10,&
412 &         'There must be 6 different q-points in long wavelength limit (see gw_nqlwl)'
413         ABI_COMMENT(msg)
414      end if
415 
416      ABI_FREE(tensor_cart_gwnlf)
417      ABI_FREE(tensor_red_gwnlf)
418 
419      !call wrtout(std_out," Checking Kramers Kronig on Excitonic Macroscopic Epsilon","COLL")
420      !call check_kramerskronig(BSp%nomega,REAL(BSp%omega),eps_exc(:,1))
421 
422      !call wrtout(std_out," Checking Kramers Kronig on RPA NLF Macroscopic Epsilon","COLL")
423      !call check_kramerskronig(BSp%nomega,REAL(BSp%omega),eps_rpanlf(:,1))
424 
425      !call wrtout(std_out," Checking Kramers Kronig on GW NLF Macroscopic Epsilon","COLL")
426      !call check_kramerskronig(BSp%nomega,REAL(BSp%omega),eps_gwnlf(:,1))
427 
428      !call wrtout(std_out," Checking f-sum rule on Excitonic Macroscopic Epsilon","COLL")
429 
430      !if (BSp%exchange_term>0) then
431      !  ABI_COMMENT(' f-sum rule should be checked without LF')
432      !end if
433      !call check_fsumrule(BSp%nomega,REAL(BSp%omega),AIMAG(eps_exc(:,1)),drude_plsmf)
434 
435      !call wrtout(std_out," Checking f-sum rule on RPA NLF Macroscopic Epsilon","COLL")
436      !call check_fsumrule(BSp%nomega,REAL(BSp%omega),AIMAG(eps_rpanlf(:,1)),drude_plsmf)
437 
438      !call wrtout(std_out," Checking f-sum rule on GW NLF Macroscopic Epsilon","COLL")
439      !call check_fsumrule(BSp%nomega,REAL(BSp%omega),AIMAG(eps_gwnlf(:,1)),drude_plsmf)
440    end if ! my_rank==master
441 
442    !call xmpi_barrier(comm)
443    !
444    ! The ket for the approximated DOS.
445    if (prtdos) then
446      ABI_WARNING("Calculating DOS with Haydock method")
447      ABI_CHECK(BSp%use_coupling==0,"DOS with coupling not coded")
448      iq = BSp%nq + 1
449      if (my_rank==master) then
450        !call random_seed()
451        do itt=1,SUM(Bsp%nreh)
452          call RANDOM_NUMBER(rand_phi)
453          rand_phi = two_pi*rand_phi
454          kets(itt,iq) = CMPLX( COS(rand_phi), SIN(rand_phi) )
455        end do
456        ! Normalize the vector.
457        !norm = SQRT( DOT_PRODUCT(kets(:,iq), kets(:,iq)) )
458        !kets(:,iq) = kets(:,iq)/norm
459      end if
460      call xmpi_bcast(kets(:,iq),master,comm,ierr)
461    end if
462 
463    ABI_MALLOC(green,(BSp%nomega,nkets))
464 
465    if (BSp%use_coupling==0) then
466      if(do_ep_renorm) then
467        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)
468      else
469        !YG2014
470        call haydock_herm(BSp,BS_files,Cryst,Hdr_bse,hexc%my_t1,hexc%my_t2,&
471 &         nkets,kets,green,hexc,hexc_i,comm)
472      end if
473    else
474      if (BSp%use_interp) then
475        ABI_ERROR("BSE Interpolation with coupling is not supported")
476      else
477        call haydock_psherm(BSp,BS_files,Cryst,Hdr_bse,hexc,hexc_i,hsize,my_t1,my_t2,nkets,kets,green,comm)
478      end if
479    end if
480 
481    !
482    ! Add 1 to have the real part right.
483    green = one + green
484 
485    if (my_rank==master) then ! Master writes the final results.
486      !
487      if (prtdos) then
488        ABI_MALLOC(dos,(BSp%nomega))
489        dos = -AIMAG(green(:,BSp%nq+1))
490        call exc_write_data(BSp,BS_files,"EXC_MDF",green,prefix=prefix,dos=dos)
491        ABI_FREE(dos)
492      else
493        call exc_write_data(BSp,BS_files,"EXC_MDF",green,prefix=prefix)
494      end if
495      !
496      ! =========================
497      ! === Write out Epsilon ===
498      ! =========================
499 
500      ABI_MALLOC(tensor_cart,(BSp%nomega,6))
501      ABI_MALLOC(tensor_red,(BSp%nomega,6))
502 
503      call wrtout(std_out," Calculating EXC dielectric tensor","COLL")
504      call haydock_mdf_to_tensor(BSp,Cryst,green,tensor_cart,tensor_red,ierr)
505 
506      if (ierr == 0) then
507          ! Writing tensor
508          call exc_write_tensor(BSp,BS_files,"EXC_TSR_CART",tensor_cart)
509          call exc_write_tensor(BSp,BS_files,"EXC_TSR_RED",tensor_red)
510      else
511          write(msg,'(3a)')&
512 &          'The EXC dielectric complex tensor cannot be computed',ch10,&
513 &          'There must be 6 different q-points in long wavelength limit (see gw_nqlwl)'
514          ABI_COMMENT(msg)
515      end if
516 
517      ABI_FREE(tensor_cart)
518      ABI_FREE(tensor_red)
519      !
520      ! This part will be removed when fldiff will be able to compare two mdf files.
521      write(ab_out,*)" "
522      write(ab_out,*)"Macroscopic dielectric function:"
523      write(ab_out,*)"omega [eV] <KS_RPA_nlf>  <GW_RPA_nlf>  <BSE> "
524      do io=1,MIN(BSp%nomega,10)
525        omegaev = REAL(BSp%omega(io))*Ha_eV
526        ks_avg  = SUM( eps_rpanlf(io,:)) / Bsp%nq
527        gw_avg  = SUM( eps_gwnlf (io,:)) / Bsp%nq
528        exc_avg = SUM( green     (io,:)) / BSp%nq
529        write(ab_out,'(7f9.4)')omegaev,ks_avg,gw_avg,exc_avg
530      end do
531      write(ab_out,*)" "
532 
533      ! Write MDF file with the final results.
534      ! FIXME: It won't work if prtdos == True
535      path = strcat(BS_files%out_basename,strcat(prefix,"_MDF.nc"))
536      NCF_CHECK(nctk_open_create(ncid, path, xmpi_comm_self))
537      NCF_CHECK(cryst%ncwrite(ncid))
538      NCF_CHECK(ebands_ncwrite(QP_bst, ncid))
539      call mdfs_ncwrite(ncid, Bsp, green, eps_rpanlf, eps_gwnlf)
540      NCF_CHECK(nf90_close(ncid))
541    end if
542 
543    ABI_FREE(green)
544    ABI_FREE(eps_rpanlf)
545    ABI_FREE(eps_gwnlf)
546    ABI_FREE(dos_ks)
547    ABI_FREE(dos_gw)
548 
549    call ebands_free(EPBSt)
550    call ebands_free(EP_QPBst)
551 
552  end do ! itemp loop
553 
554  ABI_FREE(opt_cvk)
555 
556  ABI_FREE(kets)
557 
558  call timab(694,2,tsec) ! exc_haydock_driver(apply
559  call timab(695,1,tsec) ! exc_haydock_driver(end)
560 
561  !YG2014
562  call hexc_free(hexc)
563  call hexc_interp_free(hexc_i)
564 
565  if (do_ep_renorm) then
566    ABI_FREE(ep_renorms)
567    ABI_FREE(bs2eph)
568  end if
569 
570  call timab(695,2,tsec) ! exc_haydock_driver(end)
571  call timab(690,2,tsec) ! exc_haydock_driver
572 
573 end subroutine exc_haydock_driver

m_haydock/haydock_bilanczos [ Functions ]

[ Top ] [ m_haydock ] [ Functions ]

NAME

 haydock_bilanczos

FUNCTION

  Reads the excitonic Hamiltonian from file and construct the Lanczos set of vectors
  by iterative matrix-vector multiplications for any general matrix.

INPUTS

  BSp<type(excparam)>=Parameters defining the Bethe-Salpeter calculation.
    omega(BSp%nomega)=Frequency mesh for the macroscopic dielectric function (broadening is already included).
 hize
 my_t1,my_t2
 hreso(hsize,my_t1:my_t2)
 hcoup(hsize,my_t1:my_t2)
 nkets
 kets(hsize,nkets)
 comm=MPI communicator.

OUTPUT

  green(BSp%nomega)=The imaginary part of the macroscopic dielectric function.

SOURCE

1827 subroutine haydock_bilanczos(BSp,BS_files,Cryst,Hdr_bse,hexc,hexc_i,hsize,my_t1,my_t2,nkets,kets,ep_renorms,green,comm)
1828 
1829 !Arguments ------------------------------------
1830 !scalars
1831  integer,intent(in) :: hsize,my_t1,my_t2,nkets,comm
1832  type(crystal_t),intent(in) :: Cryst
1833  type(excparam),intent(in) :: BSp
1834  type(excfiles),intent(in) :: BS_files
1835  type(Hdr_type),intent(in) :: Hdr_bse
1836 !arrays
1837  complex(dp),intent(out) :: green(BSp%nomega,BSp%nq)
1838  complex(dpc),intent(in) :: kets(hsize,nkets)
1839  complex(dpc),intent(in) :: ep_renorms(hsize)
1840 
1841 !Local variables ------------------------------
1842 !scalars
1843  integer,parameter :: master=0
1844  integer :: inn,itt,out_unt,nproc,my_rank,ierr
1845  integer :: niter_file,niter_max,niter_done,nsppol,iq,my_nt,term_type,n_all_omegas
1846  real(dp) :: ket0_hbar_norm,nfact,norm
1847  logical :: can_restart,is_converged
1848  complex(dpc) :: factor
1849  character(len=fnlen),parameter :: tag_file="_HAYDC_SAVE"
1850  character(len=500) :: msg
1851  character(len=fnlen) :: restart_file,out_file
1852  type(hexc_t),intent(in) :: hexc
1853  type(hexc_interp_t),intent(in) :: hexc_i
1854 !arrays
1855  complex(dpc),allocatable :: aa_file(:),bb_file(:),cc_file(:)
1856  complex(dpc),allocatable :: aa(:),bb(:),cc(:)
1857  complex(dpc),allocatable :: phi_np1(:),phi_n(:),phi_nm1(:)
1858  complex(dpc),allocatable :: phit_np1(:),phit_n(:),phit_nm1(:)
1859  complex(dpc),allocatable :: cbuff(:), phi_n_file(:),phi_np1_file(:)
1860  complex(dpc),allocatable :: ket0(:)
1861  complex(dpc),allocatable :: hphi_n(:), hphit_n(:)
1862  complex(dpc),allocatable :: all_omegas(:),green_temp(:,:)
1863  logical :: check(2)
1864 
1865 !************************************************************************
1866 
1867  ABI_WARNING("Haydock with Bilanczos is still under development")
1868 
1869  if(BSp%use_interp) then
1870    ABI_ERROR("Bilanczos is not yet implemented with interpolation")
1871  end if
1872 
1873  nproc  = xmpi_comm_size(comm)
1874  my_rank= xmpi_comm_rank(comm)
1875  nsppol = Hdr_bse%nsppol
1876 
1877  my_nt = my_t2-my_t1+1
1878  ABI_CHECK(my_nt>0,"One of the processors has zero columns")
1879 
1880  ! Multiplicative factor (k-point sampling and unit cell volume)
1881  ! TODO be careful with the spin here
1882  ! TODO four_pi comes from the coulomb term 1/|q| is already included in the
1883  ! oscillators hence the present approach wont work if a cutoff interaction is used.
1884  nfact = four_pi/(Cryst%ucvol*BSp%nkbz)
1885  if (nsppol==1) nfact=two*nfact
1886 
1887  write(msg,'(a,i0)')' Bi-Lanczos algorithm with MAX number of iterations: ',BSp%niter
1888  call wrtout(std_out,msg,"COLL")
1889  !
1890  ! Check for presence of the restart file.
1891  can_restart=.FALSE.
1892 
1893  if ( BS_files%in_haydock_basename /= BSE_NOFILE) then
1894    restart_file = strcat(BS_files%in_haydock_basename,tag_file)
1895    if (file_exists(restart_file) ) then
1896      can_restart=.TRUE.
1897      msg = strcat(" Restarting Haydock calculation from file: ",restart_file)
1898      call wrtout(std_out,msg,"COLL")
1899      call wrtout(ab_out,msg,"COLL")
1900      ABI_ERROR("Restart is not implemented")
1901    else
1902      can_restart=.FALSE.
1903      ABI_WARNING(strcat("Cannot find restart file: ",restart_file))
1904    end if
1905  end if
1906  !
1907  ! Open the file and writes basic dimensions and info.
1908  if (my_rank==master) then
1909    out_file = TRIM(BS_files%out_basename)//TRIM(tag_file)
1910    if (open_file(out_file,msg,newunit=out_unt,form="unformatted") /= 0) then
1911      ABI_ERROR(msg)
1912    end if
1913    ! write header TODO: standardize this part.
1914    write(out_unt)hsize,Bsp%use_coupling,BSE_HAYD_IMEPS,nkets,Bsp%broad
1915  end if
1916  !
1917  ! Select the terminator for the continued fraction.
1918  term_type=0 !; if (Bsp%hayd_term>0) term_type=2
1919  call wrtout(std_out,sjoin("Using terminator type: ",itoa(term_type)),"COLL")
1920  !
1921  ! Calculate green(w) for the different starting kets.
1922  green=czero
1923  do iq=1,nkets
1924    ABI_MALLOC(ket0,(hexc%hsize))
1925    ket0 = kets(:,iq)
1926    !
1927    niter_file=0
1928 
1929    if (can_restart) then
1930 !     call haydock_restart(BSp,restart_file,BSE_HAYD_IMEPS,iq,hsize,&
1931 !&      niter_file,aa_file,bb_file,phi_np1_file,phi_n_file,comm)
1932    end if
1933    !
1934    ABI_MALLOC(phi_nm1,(my_nt))
1935    ABI_MALLOC(phi_n,(my_nt))
1936    ABI_MALLOC(phi_np1,(my_nt))
1937    ABI_MALLOC(phit_nm1,(my_nt))
1938    ABI_MALLOC(phit_n,(my_nt))
1939    ABI_MALLOC(phit_np1,(my_nt))
1940    ABI_MALLOC(hphi_n,(hexc%hsize))
1941    ABI_MALLOC(hphit_n,(hexc%hsize))
1942    !
1943    ! TODO: Note the different convention used for the coefficients
1944    ! Should use the same convention in the Hermitian case.
1945    niter_max = niter_file + Bsp%niter
1946    ABI_MALLOC(aa,(niter_max))
1947    ABI_MALLOC(bb,(niter_max))
1948    ABI_MALLOC(cc,(niter_max))
1949    aa=czero; bb=czero; cc=czero
1950 
1951    if (niter_file==0) then ! Calculation from scratch.
1952      phi_nm1 = ket0(my_t1:my_t2)
1953      phit_nm1 = ket0(my_t1:my_t2)
1954      norm = DZNRM2(hexc%hsize,ket0,1)
1955      phi_nm1=phi_nm1/norm
1956      phit_nm1=phit_nm1/norm
1957 
1958      call hexc_matmul_elphon(hexc,phi_nm1,hphi_n,'N',ep_renorms)
1959      call hexc_matmul_elphon(hexc,phit_nm1,hphit_n,'C',ep_renorms)
1960 
1961      aa(1)=xdotc(my_nt,phit_nm1,1,hphi_n(my_t1:),1)
1962      call xmpi_sum(aa(1:1),comm,ierr)
1963 
1964      phi_n = hphi_n(my_t1:my_t2) - aa(1)*phi_nm1
1965      phit_n = hphit_n(my_t1:my_t2) - CONJG(aa(1))*phit_nm1
1966 
1967      bb(1)=xdotc(my_nt,phi_n,1,phi_n,1)
1968      call xmpi_sum(bb(1:1),comm,ierr)
1969      bb(1) = SQRT(bb(1))
1970 
1971      cc(1)=xdotc(my_nt,phit_n,1,phi_n,1)
1972      call xmpi_sum(cc(1:1),comm,ierr)
1973      cc(1) = cc(1)/bb(1)
1974 
1975      phi_n   = phi_n  /bb(1)
1976      phit_n  = phit_n /CONJG(cc(1))
1977      niter_done=1  ! TODO Be careful here
1978 
1979    else ! Use the previously calculates a and b.
1980      niter_done=niter_file
1981      ABI_ERROR("Restart not coded")
1982      !aa(1:niter_done) = aa_file
1983      !bb(1:niter_done) = bb_file
1984      !phi_np1=phi_np1_file(my_t1:my_t2)   ! Select the slice treated by this node.
1985      !phi_n  =phi_n_file  (my_t1:my_t2)
1986    end if
1987 
1988    if (can_restart) then
1989      ABI_FREE(aa_file)
1990      ABI_FREE(bb_file)
1991      ABI_FREE(cc_file)
1992      ABI_FREE(phi_np1_file)
1993      ABI_FREE(phi_n_file)
1994    end if
1995 
1996    ! This factor gives the correct results
1997    factor = -nfact*(DZNRM2(hexc%hsize,ket0,1)**2)
1998 
1999    ! Which quantity should be checked for convergence?
2000    check = (/.TRUE.,.TRUE./)
2001    if (ABS(Bsp%haydock_tol(2)-one)<tol6) check = (/.TRUE. ,.FALSE./)
2002    if (ABS(Bsp%haydock_tol(2)-two)<tol6) check = (/.FALSE.,.TRUE./)
2003    ! Create new frequencies "mirror" in negative range to add
2004    ! their contributions. Can be improved by computing only once
2005    ! zero frequency, but loosing clearness
2006    n_all_omegas = 2*BSp%nomega
2007 
2008    ABI_MALLOC(all_omegas,(n_all_omegas))
2009    ! Put all omegas with frequency > 0 in table
2010    all_omegas(BSp%nomega+1:n_all_omegas) = BSp%omega
2011    ! Put all omegas with frequency < 0
2012    ! Warning, the broadening must be kept positive
2013    all_omegas(1:BSp%nomega) = -DBLE(BSp%omega(BSp%nomega:1:-1)) + j_dpc*AIMAG(BSp%omega(BSp%nomega:1:-1))
2014 
2015    ABI_MALLOC(green_temp,(n_all_omegas,nkets))
2016 
2017 
2018 
2019    call haydock_bilanczos_optalgo(niter_done,niter_max,n_all_omegas,all_omegas,BSp%haydock_tol(1),check,hexc,hexc_i,&
2020 &    hsize,my_t1,my_t2,factor,term_type,ep_renorms,aa,bb,cc,ket0,ket0_hbar_norm,phi_nm1,phi_n,phi_np1,&
2021 &    phit_nm1,phit_n,phit_np1,green_temp(:,iq),inn,is_converged,comm)
2022 
2023 
2024    ! Computing result from two ranges of frequencies
2025    ! The real part is added, the imaginary part is substracted
2026    green(:,iq) = green_temp(BSp%nomega+1:n_all_omegas,iq)+CONJG(green_temp(BSp%nomega:1:-1,iq))
2027 
2028    ABI_FREE(all_omegas)
2029    ABI_FREE(green_temp)
2030 
2031 
2032    !
2033    ! Save the a"s and the b"s for possible restarting.
2034    ! 1) Info on the Q.
2035    ! 2) Number of iterations performed.
2036    ! 3) do iter=1,niter_performed
2037    !      aa(iter),bb(iter)
2038    !    end do
2039    ! 4) |n-1>
2040    !    |n>
2041    !    |n+1>
2042    !
2043    if (my_rank==master) then ! Open the file and writes basic dimensions and info.
2044      write(out_unt)Bsp%q(:,iq)
2045      write(out_unt)MIN(inn,niter_max)  ! NB: if the previous loop completed inn=niter_max+1
2046      do itt=1,MIN(inn,niter_max)        !     if we exited then inn is not incremented by one.
2047        write(out_unt)itt,aa(itt),bb(itt)
2048      end do
2049    end if
2050    !
2051    ! cbuff is used as workspace to gather |n-1>, |n> and |n+1>.
2052    ABI_MALLOC(cbuff,(hsize))
2053    cbuff=czero; cbuff(my_t1:my_t2) = phi_nm1
2054    call xmpi_sum_master(cbuff,master,comm,ierr)
2055    if (my_rank==master) write(out_unt) cbuff ! |n-1>
2056 
2057    cbuff=czero; cbuff(my_t1:my_t2) = phi_n
2058    call xmpi_sum_master(cbuff,master,comm,ierr)
2059    if (my_rank==master) write(out_unt) cbuff ! |n>
2060 
2061    cbuff=czero; cbuff(my_t1:my_t2) = phi_np1
2062    call xmpi_sum_master(cbuff,master,comm,ierr)
2063    if (my_rank==master) write(out_unt) cbuff ! |n+1>
2064 
2065    ABI_FREE(phi_nm1)
2066    ABI_FREE(phi_n)
2067    ABI_FREE(phi_np1)
2068    ABI_FREE(phit_nm1)
2069    ABI_FREE(phit_n)
2070    ABI_FREE(phit_np1)
2071    ABI_FREE(hphi_n)
2072    ABI_FREE(hphit_n)
2073    ABI_FREE(cbuff)
2074    ABI_FREE(aa)
2075    ABI_FREE(bb)
2076    ABI_FREE(cc)
2077    ABI_FREE(ket0)
2078  end do ! iq
2079 
2080  if (my_rank==master) close(out_unt)
2081 
2082  call xmpi_barrier(comm)
2083 
2084 end subroutine haydock_bilanczos

m_haydock/haydock_bilanczos_optalgo [ Functions ]

[ Top ] [ m_haydock ] [ Functions ]

NAME

 haydock_bilanczos_optalgo

FUNCTION

  Haydock algorithm for general matrix

INPUTS

  niter_done=Number of iterations already performed (0 if the run starts from scratch).
  niter_tot=Max number of iterations. Always > niter_done
  nomega=Number of Frequency points for the evaluation of the matrix element.
  omega(nomega)=Frequency set (imaginary part is already included).
  tol_iter=Tollerance used to stop the the algorithm.
  check(2)=Logical flags to specify where both the real and the imaginary part of the
    matrix elements of the Green functions have to be checked for convergence.
  hsize=Size of the blocks.
  my_t1,my_t2=Indeces of the first and last column stored treated by this done.
  term_type=0 if no terminator is used, 1 otherwise.
  hreso(hsize,my_t1:my_t2)=The columns of the resonant block.
  hcoup(hsize,my_t1:my_t2)=The columns of the coupling block.
  factor
  comm=MPI communicator.

OUTPUT

  green(nomega)=Output matrix elements.
  inn=Last iteration performed.
  is_converged=.TRUE. of the algorithm converged.

SIDE EFFECTS

  phi_nm1(my_t2-my_t1+1), phi_n(my_t2-my_t1+1)
    input: vectors used to initialize the iteration
    output: the vectors obtained in the last iteration
  aa(niter_tot) and bb(niter_tot+1)
    if niter_done>0: aa(1:niter_done), bb(1:niter_done) store the coefficients of the previous run.
    when the routine returns aa(1:inn) and bb(1:inn) contain the matrix elements of the tridiagonal form.
  cc(niter_tot+1)

SOURCE

2128 subroutine haydock_bilanczos_optalgo(niter_done,niter_tot,nomega,omega,tol_iter,check,hexc,hexc_i,hsize,my_t1,my_t2,&
2129 &  factor,term_type,ep_renorms,aa,bb,cc,ket0,ket0_hbar_norm,phi_nm1,phi_n,phi_np1,phit_nm1,phit_n,phit_np1,&
2130 &  green,inn,is_converged,comm)
2131 
2132 !Arguments ------------------------------------
2133 !scalars
2134  integer,intent(in) :: niter_tot,niter_done,nomega,comm,hsize,my_t1,my_t2,term_type
2135  integer,intent(out) :: inn
2136  logical,intent(out) :: is_converged
2137  real(dp),intent(in) :: tol_iter,ket0_hbar_norm
2138  complex(dpc),intent(in) :: factor
2139  type(hexc_t),intent(in) :: hexc
2140  type(hexc_interp_t),intent(in) :: hexc_i
2141 !arrays
2142  complex(dpc),intent(inout) :: bb(niter_tot+1)
2143  complex(dpc),intent(out) :: green(nomega)
2144  complex(dpc),intent(in) :: omega(nomega)
2145  complex(dpc),intent(inout) :: aa(niter_tot),cc(niter_tot+1)
2146  complex(dpc),intent(in) :: ket0(my_t2-my_t1+1)
2147  complex(dpc),intent(in) :: ep_renorms(hsize)
2148  complex(dpc),intent(inout) :: phi_nm1(my_t2-my_t1+1)
2149  complex(dpc),intent(inout) :: phi_n  (my_t2-my_t1+1)
2150  complex(dpc),intent(inout) :: phi_np1(my_t2-my_t1+1)
2151  complex(dpc),intent(inout) :: phit_nm1(my_t2-my_t1+1)
2152  complex(dpc),intent(inout) :: phit_n  (my_t2-my_t1+1)
2153  complex(dpc),intent(inout) :: phit_np1(my_t2-my_t1+1)
2154  logical,intent(in) :: check(2)
2155 
2156 !Local variables ------------------------------
2157 !scalars
2158  integer :: my_nt,niter_min,nconv !,ierr
2159  character(len=500) :: msg
2160  logical :: keep_vectors=.TRUE.
2161 !arrays
2162  real(dp) :: abs_err(nomega,2) !,ww_err(nomega,2)
2163  complex(dpc),allocatable :: oldg(:),newg(:)
2164  complex(dpc),allocatable :: hphi_np1(:),hphit_np1(:),save_phi(:,:),save_phit(:,:)
2165  complex(dpc),allocatable :: g00(:)
2166  logical :: test(2)
2167  integer :: ierr
2168 
2169 !************************************************************************
2170 
2171  ABI_UNUSED(ket0_hbar_norm)
2172  ABI_UNUSED(ket0(1))
2173  ABI_UNUSED(hexc_i%hsize_dense)
2174 
2175  my_nt = my_t2-my_t1+1
2176 
2177  ABI_MALLOC(oldg,(nomega))
2178  ABI_MALLOC(newg,(nomega))
2179  ABI_MALLOC(g00,(nomega))
2180  oldg=czero; newg=czero; g00=czero
2181  nconv=0
2182 
2183  keep_vectors = (keep_vectors.and.xmpi_comm_size(comm)==1)
2184  if (keep_vectors) then
2185    ABI_MALLOC(save_phi,(my_t2-my_t1+1,niter_tot))
2186    ABI_MALLOC_OR_DIE(save_phit,(my_t2-my_t1+1,niter_tot),ierr)
2187    save_phi=czero
2188    save_phit=czero
2189  end if
2190 
2191  ABI_MALLOC_OR_DIE(hphi_np1,(hexc%hsize),ierr)
2192  ABI_MALLOC_OR_DIE(hphit_np1,(hexc%hsize),ierr)
2193 
2194  do inn=niter_done+1,niter_tot
2195 
2196    !|n+1> = H |n> using all eh components.
2197    call hexc_matmul_elphon(hexc, phi_n, hphi_np1, 'N', ep_renorms)
2198    call hexc_matmul_elphon(hexc, phit_n, hphit_np1, 'C', ep_renorms)
2199 
2200    ! a(n) = < phit_n | H  | phi_n >
2201    aa(inn)=xdotc(my_nt,phit_n,1,hphi_np1(my_t1:),1)
2202    call xmpi_sum(aa(inn),comm,ierr)
2203 
2204    ! |n+1> = |n+1> - a(n)|Vn> - c(n)|n-1>
2205    phi_np1 = hphi_np1(my_t1:my_t2) - aa(inn)*phi_n - cc(inn-1)*phi_nm1
2206    phit_np1 = hphit_np1(my_t1:my_t2) - CONJG(aa(inn))*phit_n - CONJG(bb(inn-1))*phit_nm1
2207 
2208    bb(inn) = xdotc(my_nt,phi_np1,1,phi_np1,1)
2209    call xmpi_sum(bb(inn),comm,ierr)
2210    bb(inn) = SQRT(bb(inn))
2211 
2212    cc(inn) = xdotc(my_nt,phit_np1,1,phi_np1,1)
2213    call xmpi_sum(cc(inn),comm,ierr)
2214    cc(inn) = cc(inn)/bb(inn)
2215 
2216    phi_np1 = phi_np1 / bb(inn)
2217    phit_np1 = phit_np1 / CONJG(cc(inn))
2218 
2219    !
2220    ! |n-1> = |n>
2221    ! |n>   = |n+1>
2222    phi_nm1 = phi_n
2223    phi_n   = phi_np1
2224    phit_nm1 = phit_n
2225    phit_n = phit_np1
2226 
2227    if (keep_vectors) then
2228      save_phi(:,inn) = phi_n
2229      save_phit(:,inn) = phit_n
2230    end if
2231    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))
2232    call wrtout(std_out,msg,"COLL")
2233 
2234    call continued_fract_general(inn,term_type,aa,bb,cc,nomega,omega,g00)
2235    newg = factor*g00
2236    !
2237    ! Avoid spurious convergence.
2238    niter_min=4; if (niter_done>1) niter_min=niter_done+1
2239    if (inn>niter_min) then
2240      test=.TRUE.
2241      abs_err(:,1) = ABS(DBLE (newg-oldg))
2242      abs_err(:,2) = ABS(AIMAG(newg-oldg))
2243      !
2244      if (tol_iter>zero) then
2245        ! Test on the L1 norm.
2246        if (check(1)) test(1) = SUM(abs_err(:,1)) < tol_iter*SUM(ABS(DBLE (newg)))
2247        if (check(2)) test(2) = SUM(abs_err(:,2)) < tol_iter*SUM(ABS(AIMAG(newg)))
2248      else
2249        ! Stringent test for each point.
2250        if (check(1)) test(1) = ALL( abs_err(:,1) < -tol_iter*ABS(DBLE (newg)))
2251        if (check(2)) test(2) = ALL( abs_err(:,2) < -tol_iter*ABS(AIMAG(newg)))
2252      end if
2253      !
2254      if (ALL(test)) then
2255        nconv = nconv+1
2256      else
2257        nconv = 0
2258      end if
2259      if (nconv==2) then
2260        if(inn<100)then
2261          write(msg,'(a,es10.2,a)')&
2262 &          " >>> Haydock algorithm converged twice within haydock_tol= ",tol_iter," after less than 100 iterations."
2263        else
2264          write(msg,'(a,es10.2,a,i0,a)')&
2265 &          " >>> Haydock algorithm converged twice within haydock_tol= ",tol_iter," after ",inn," iterations."
2266        endif
2267        call wrtout(std_out,msg,'COLL')
2268        call wrtout(ab_out,msg,'COLL')
2269        EXIT
2270      end if
2271    end if
2272    !
2273    oldg = newg
2274  end do ! inn
2275 
2276  green = newg
2277  if (nconv/=2) then
2278    write(msg,'(a,es10.2,a,i0,a)')&
2279 &    " WARNING: Haydock algorithm did not converge within ",tol_iter," after ",niter_tot," iterations."
2280    call wrtout(std_out,msg,'COLL')
2281    call wrtout(ab_out,msg,'COLL')
2282  end if
2283 
2284  is_converged = (nconv==2)
2285 
2286  ABI_FREE(oldg)
2287  ABI_FREE(newg)
2288  ABI_FREE(g00)
2289  ABI_FREE(hphi_np1)
2290  ABI_FREE(hphit_np1)
2291 
2292  if (keep_vectors) then
2293    ABI_FREE(save_phi)
2294    ABI_FREE(save_phit)
2295  end if
2296 
2297  !! if (keep_vectors) then
2298  !!   tdim = MIN(inn,niter_tot)
2299  !!   ABI_MALLOC(ovlp,(tdim,tdim))
2300 
2301  !!   ABI_MALLOC(phi_test,(hsize))
2302  !!   ABI_MALLOC(phi_test2,(hsize))
2303 
2304  !!   max_err=smallest_real; mean_err=zero; mean_err2=zero; row_max=-1
2305  !!   do ii=1,tdim
2306  !!     parity = (-1)**(ii+1)
2307  !!     phi_test  = save_phi(:,ii)
2308  !!     call hexc_matmul_full(hexc, hexc_i, phi_test, phi_test2, parity)
2309  !!     !phi_test2 = MATMUL(hreso,phi_test) + parity * MATMUL(hcoup,CONJG(phi_test))
2310  !!     ovlp(ii,ii) = DOT_PRODUCT(phi_test,phi_test2) + DOT_PRODUCT(phi_test2,phi_test)
2311  !!     err = ABS(ovlp(ii,ii)-cone)
2312  !!     mean_err  = mean_err + err
2313  !!     mean_err2 = mean_err2 + err**2
2314  !!     if (err > max_err) then
2315  !!       max_err = err
2316  !!       row_max = ii
2317  !!     end if
2318  !!   end do
2319  !!   mean_err = mean_err/tdim
2320  !!   std_dev = mean_err2/tdim -mean_err**2
2321  !!   write(std_out,'(a,i0,1x,3es14.6)')&
2322  !!&    " Error in normalization (ii, max_err,mean,std_dev): ",row_max,max_err,mean_err,std_dev
2323 
2324  !!   ABI_FREE(phi_test)
2325  !!   ABI_FREE(phi_test2)
2326  !!
2327  !!   ABI_MALLOC(alpha,(hsize,tdim))
2328 
2329  !!   ! Less efficient but for sake of simplicity with hexc_matmul
2330  !!   ! TODO possibility to call hreso * phi, and hcoup * phi separately
2331  !!   do ii=1,tdim
2332  !!     parity = (-1)**(ii+1)
2333  !!     call hexc_matmul_full(hexc, hexc_i, save_phi(:,ii), alpha(:,ii), parity)
2334  !!   end do
2335 
2336  !!   !alpha = MATMUL(hreso,save_phi(:,1:tdim))
2337  !!   !
2338  !!   !do ii=1,tdim
2339  !!   !  parity = (-1)**(ii+1)
2340  !!   !  alpha(:,ii) =  alpha(:,ii) + parity*MATMUL(hcoup,CONJG(save_phi(:,ii)))
2341  !!   !end do
2342 
2343  !!   ovlp = MATMUL(TRANSPOSE(CONJG(save_phi(:,1:tdim))),alpha)
2344 
2345  !!   ABI_MALLOC(beta,(hsize,tdim))
2346  !!   do ii=1,tdim
2347  !!     parity = (-1)**(ii+1)
2348  !!     beta(:,ii)  =  parity*save_phi(:,ii)
2349  !!     alpha(:,ii) = -parity*alpha(:,ii)
2350  !!   end do
2351 
2352  !!   ovlp = ovlp - MATMUL(TRANSPOSE(CONJG(beta)),alpha)
2353 
2354  !!   max_err=smallest_real; row_max=-1; col_max=-1
2355  !!   mean_err=zero; mean_err2=zero
2356  !!   do jj=1,tdim
2357  !!     do ii=1,jj
2358  !!       err = ABS(ovlp(ii,jj))
2359  !!       if (ii==jj) err = ABS(err - one)
2360  !!       mean_err  = mean_err + err
2361  !!       mean_err2 = mean_err2 + err**2
2362  !!       if (err > max_err) then
2363  !!         max_err = err
2364  !!         row_max=ii
2365  !!         col_max=jj
2366  !!       end if
2367  !!     end do
2368  !!   end do
2369 
2370  !!   mean_err = mean_err/(tdim*(tdim+1)/2)
2371  !!   std_dev = mean_err2/(tdim*(tdim+1)/2) - mean_err**2
2372  !!   write(std_out,'(a,2(i0,1x),3es14.6)')&
2373  !!      " Error in Hbar-ortho (i,j), max_err, mean, std_dev ",row_max,col_max,max_err,mean_err,std_dev
2374  !!   !call print_arr(ovlp,max_r=185,max_c=10,unit=std_out)
2375 
2376  !!   ABI_FREE(alpha)
2377  !!   ABI_FREE(beta)
2378  !!   ABI_FREE(ovlp)
2379  !!   ABI_FREE(save_phi)
2380  !! end if
2381 
2382 end subroutine haydock_bilanczos_optalgo

m_haydock/haydock_herm [ Functions ]

[ Top ] [ m_haydock ] [ Functions ]

NAME

 haydock_herm

FUNCTION

  Reads the excitonic Hamiltonian from file and construct the Lanczos set of vectors
  by iterative matrix-vector multiplications.

INPUTS

 BSp<excparam>=Parameters for the Bethe-Salpeter calculation.
 BS_files<excparam>=Files associated to the bethe_salpeter code.
 Cryst<crystal_t>=Info on the crystalline structure.
 Pawtab(Cryst%ntypat*usepaw)<pawtab_type>=PAW tabulated starting data.
 hize=Size of the excitonic matrix.
 my_t1,my_t2=First and last columns treated by this node.
 nkets=Number of starting vectors for Haydock method.
 kets(hsize,nkets)=The kets in the eh representation.
 comm=MPI communicator.

OUTPUT

  green(BSp%nomega,nkets)=

SOURCE

602 subroutine haydock_herm(BSp,BS_files,Cryst,Hdr_bse,my_t1,my_t2,&
603 & nkets,kets,green,hexc,hexc_i,comm)
604 
605 !Arguments ------------------------------------
606 !scalars
607  integer,intent(in) :: my_t1,my_t2,nkets,comm
608  type(crystal_t),intent(in) :: Cryst
609  type(excparam),intent(in) :: BSp
610  type(excfiles),intent(in) :: BS_files
611  type(Hdr_type),intent(in) :: Hdr_bse
612  type(hexc_t),intent(inout) :: hexc
613  type(hexc_interp_t),intent(inout) :: hexc_i
614 !arrays
615  complex(dp),intent(out) :: green(BSp%nomega,nkets)
616  complex(dpc),intent(in) :: kets(hexc%hsize,nkets)
617 
618 !Local variables ------------------------------
619 !scalars
620  integer,parameter :: master=0
621  integer :: inn,nproc,my_rank,ierr
622  integer :: niter_file,niter_max,niter_done,nsppol,iq,my_nt,term_type,n_all_omegas
623  real(dp) :: norm,nfact
624  logical :: can_restart,is_converged
625  complex(dpc) :: factor
626  character(len=500) :: msg
627  character(len=fnlen),parameter :: tag_file="_HAYDR_SAVE"
628  character(len=fnlen) :: restart_file,out_file
629  type(haydock_type) :: haydock_file
630 !arrays
631  real(dp),allocatable :: bb_file(:)
632  real(dp),allocatable :: bb(:)
633  complex(dpc),allocatable :: aa(:),phi_nm1(:),phi_n(:),hphi_n(:),hphi_nm1(:)
634  complex(dpc),allocatable :: aa_file(:),phi_n_file(:),phi_nm1_file(:)
635  complex(dpc),allocatable :: ket0(:),all_omegas(:),green_temp(:,:)
636 ! complex(dpc),allocatable :: diag_dense(:)
637  logical :: check(2)
638 
639 !************************************************************************
640 
641  ABI_CHECK(Bsp%nsppol==1,"nsppol > 1 not implemented yet")
642 
643  nproc  = xmpi_comm_size(comm); my_rank= xmpi_comm_rank(comm)
644  nsppol = Hdr_bse%nsppol
645 
646  if (BSp%use_interp) then
647    ABI_COMMENT("No parallelization in Interpolation")
648    my_nt = SUM(Bsp%nreh_interp)
649  else
650    my_nt = my_t2-my_t1+1
651  end if
652 
653  ABI_CHECK(my_nt>0,"One of the processors has zero columns")
654 
655  write(msg,'(a,i0)')' Haydock algorithm with MAX number of iterations: ',BSp%niter
656  call wrtout(std_out,msg,"COLL")
657  !
658  ! Select the terminator for the continued fraction.
659  term_type=0; if (Bsp%hayd_term>0) term_type=1
660  call wrtout(std_out,sjoin("Using terminator type: ",itoa(term_type)),"COLL")
661  !
662  ! Check for presence of the restart file.
663  can_restart=.FALSE.
664  if ( BS_files%in_haydock_basename /= BSE_NOFILE) then
665    restart_file = TRIM(BS_files%in_haydock_basename)//TRIM(tag_file)
666    if (file_exists(restart_file) ) then
667      can_restart=.TRUE.
668      msg = " Restarting Haydock calculation from file: "//TRIM(restart_file)
669      call wrtout(std_out,msg,"COLL")
670      call wrtout(ab_out,msg,"COLL")
671    else
672      can_restart=.FALSE.
673      call wrtout(ab_out," WARNING: cannot find restart file: "//TRIM(restart_file),"COLL")
674    end if
675  end if
676  ABI_CHECK(.not.can_restart,"restart not yet implemented")
677 
678  ! Open the file and write basic dimensions and info.
679  if (my_rank==master) then
680    out_file = TRIM(BS_files%out_basename)//TRIM(tag_file)
681    call open_haydock(out_file,haydock_file)
682    haydock_file%hsize = hexc%hsize
683    haydock_file%use_coupling = Bsp%use_coupling
684    haydock_file%op = BSE_HAYD_IMEPS
685    haydock_file%nq = nkets
686    haydock_file%broad = Bsp%broad
687    call write_dim_haydock(haydock_file)
688  end if
689 
690  !
691  ! Calculate green(w) for the different starting points.
692  green=czero
693  do iq=1,nkets
694    ABI_MALLOC(ket0,(hexc%hsize))
695    ket0=kets(:,iq)
696 
697    !
698    niter_file=0
699    if (can_restart) then
700      call haydock_restart(BSp,restart_file,BSE_HAYD_IMEPS,iq,hexc%hsize,&
701 &      niter_file,aa_file,bb_file,phi_nm1_file,phi_n_file,comm)
702    end if
703    !
704    ! For n>1, we have:
705    !  1) a_n = <n|H|n>
706    !  2) b_n = || H|n> - a_n|n> -b_{n-1}|n-1> ||
707    !  3) |n+1> = [H|n> -a_n|n> -b_{n-1}|n-1>]/b_n
708    !
709    ! The sequences starts with |1> normalized to 1 and b_0 =0, therefore:
710    !  a_1 = <1|H|1>
711    !  b_1 = || H|1> - a_1|1> ||
712    !  |2> = [H|1> - a_1|1>]/b_1
713    !
714    ABI_MALLOC(hphi_n,(hexc%hsize))
715    ABI_MALLOC(hphi_nm1,(hexc%hsize))
716    ABI_MALLOC(phi_nm1,(my_nt))
717    ABI_MALLOC(phi_n,(my_nt))
718 
719    niter_max = niter_file + Bsp%niter
720    ABI_MALLOC(aa,(niter_max))
721    ABI_MALLOC(bb,(niter_max))
722    aa=czero; bb=zero
723 
724    if (niter_file==0) then       ! Calculation from scratch.
725      phi_nm1=ket0(my_t1:my_t2)   ! Select the slice treated by this node.
726      norm = DZNRM2(hexc%hsize,ket0,1) ! Normalization
727      phi_nm1=phi_nm1/norm
728 
729      call hexc_matmul_tda(hexc,hexc_i,phi_nm1,hphi_n)
730 
731      aa(1)=xdotc(my_nt,phi_nm1,1,hphi_n(my_t1:),1)
732      call xmpi_sum(aa(1:1),comm,ierr)
733 
734      phi_n = hphi_n(my_t1:my_t2) - aa(1)*phi_nm1
735 
736      bb(1) = xdotc(my_nt,phi_n,1,phi_n,1)
737      call xmpi_sum(bb(1:1),comm,ierr)
738      bb(1) = SQRT(bb(1))
739 
740      phi_n = phi_n/bb(1)
741      niter_done=1
742 
743    else ! Use the previous a and b.
744      niter_done=niter_file
745      aa(1:niter_done) = aa_file
746      bb(1:niter_done) = bb_file
747      phi_nm1=phi_nm1_file(my_t1:my_t2)   ! Select the slice treated by this node.
748      phi_n  =phi_n_file  (my_t1:my_t2)
749    end if
750 
751    if (can_restart) then
752      ABI_FREE(aa_file)
753      ABI_FREE(bb_file)
754      ABI_FREE(phi_nm1_file)
755      ABI_FREE(phi_n_file)
756    end if
757 
758    ! Multiplicative factor (k-point sampling and unit cell volume)
759    ! TODO be careful with the spin here
760    ! TODO four_pi comes from the coulomb term 1/|q| is already included in the
761    ! oscillators hence the present approach wont work if a cutoff interaction is used.
762    nfact = -four_pi/(Cryst%ucvol*hexc%nbz)
763    if (nsppol==1) nfact=two*nfact
764 
765    factor = nfact*(DZNRM2(hexc%hsize,ket0,1)**2)
766 
767    ! Which quantity should be checked for convergence?
768    check = (/.TRUE.,.TRUE./)
769    if (ABS(Bsp%haydock_tol(2)-one)<tol6) check = (/.TRUE. ,.FALSE./)
770    if (ABS(Bsp%haydock_tol(2)-two)<tol6) check = (/.FALSE.,.TRUE./)
771 
772    ! Create new frequencies "mirror" in negative range to add
773    ! their contributions. Can be improved by computing only once
774    ! zero frequency, but loosing clearness
775    n_all_omegas = 2*BSp%nomega
776 
777    ABI_MALLOC(all_omegas,(n_all_omegas))
778    ! Put all omegas with frequency > 0 in table
779    all_omegas(BSp%nomega+1:n_all_omegas) = BSp%omega
780    ! Put all omegas with frequency < 0
781    ! Warning, the broadening must be kept positive
782    all_omegas(1:BSp%nomega) = -DBLE(BSp%omega(BSp%nomega:1:-1)) + j_dpc*AIMAG(BSp%omega(BSp%nomega:1:-1))
783 
784    ABI_MALLOC(green_temp,(n_all_omegas,nkets))
785 
786    call haydock_herm_algo(niter_done,niter_max,n_all_omegas,all_omegas,BSp%haydock_tol(1),check,&
787 &    my_t1,my_t2,factor,term_type,aa,bb,phi_nm1,phi_n,&
788 &    green_temp(:,iq),inn,is_converged,&
789 &    hexc, hexc_i, comm)
790 
791    ! Computing result from two ranges of frequencies
792    ! The real part is added, the imaginary part is substracted
793    green(:,iq) = green_temp(BSp%nomega+1:n_all_omegas,iq)+CONJG(green_temp(BSp%nomega:1:-1,iq))
794 
795    ABI_FREE(all_omegas)
796    ABI_FREE(green_temp)
797    !
798    ! Save the a"s and the b"s for possible restarting.
799    ! 1) Info on the Q.
800    ! 2) Number of iterations performed.
801    ! 3) do iter=1,niter_performed
802    !      aa(iter),bb(iter)
803    !    end do
804    ! 4) |n-1>
805    !    |n>
806    !
807    hphi_nm1 = czero
808    hphi_nm1(my_t1:my_t2) = phi_nm1
809    call xmpi_sum_master(hphi_nm1,master,comm,ierr)
810 
811    hphi_n = czero
812    hphi_n(my_t1:my_t2) = phi_n
813    call xmpi_sum_master(hphi_n,master,comm,ierr)
814 
815    if (my_rank==master) then
816      ! Write data for restarting
817      call write_haydock(haydock_file, hexc%hsize, Bsp%q(:,iq), aa, bb, hphi_n, hphi_nm1, MIN(inn,niter_max), factor)
818    end if
819 
820    ABI_FREE(hphi_n)
821    ABI_FREE(hphi_nm1)
822    ABI_FREE(phi_nm1)
823    ABI_FREE(phi_n)
824    ABI_FREE(aa)
825    ABI_FREE(bb)
826    ABI_FREE(ket0)
827  end do ! iq
828 
829  if (my_rank==master) call close_haydock(haydock_file)
830 
831  call xmpi_barrier(comm)
832 
833 end subroutine haydock_herm

m_haydock/haydock_herm_algo [ Functions ]

[ Top ] [ m_haydock ] [ Functions ]

NAME

 haydock_herm_algo

FUNCTION

INPUTS

  niter_done=Number of iterations already performed (0 if the run starts from scratch).
  niter_max=Max number of iterations. Always > niter_done
  nomega=Number of Frequency points for the evaluation of the matrix element.
  omega(nomega)=Frequency set (imaginary part is already included).
  tol_iter=Tolerance used to stop the algorithm.
  check(2)=Logical flags to specify where both the real and the imaginary part of the
    matrix elements of the Green functions have to be checked for convergence.
  hsize=Size of the blocks.
  my_t1,my_t2=Indices of the first and last column stored treated by this done.
  term_type=0 if no terminator is used, 1 otherwise.
  hmat(hsize,my_t1:my_t2)=The columns of the block.
  factor
  ntrans = Number of transitions
  corresp = mapping between coarse points and neighbours
  overlaps = overlaps of wavefunctions between dense k-point coarse neighbours and bands
  comm=MPI communicator.

OUTPUT

  green(nomega)=Output matrix elements.
  inn=Last iteration performed.
  is_converged=.TRUE. of the algorithm converged.

SIDE EFFECTS

  phi_nm1(my_t2-my_t1+1), phi_n(my_t2-my_t1+1)
    input: vectors used to initialize the iteration
    output: the vectors obtained in the last iteration
  aa(niter_max) and bb(niter_max)
    if niter_done>0: aa(1:niter_done), bb(1:niter_done) store the coefficients of the previous run.
    when the routine returns aa(1:inn) and bb(1:inn) contain the matrix elements of the tridiagonal form.

SOURCE

 877 subroutine haydock_herm_algo(niter_done,niter_max,nomega,omega,tol_iter,check,&
 878 & my_t1,my_t2,factor,term_type,aa,bb,phi_nm1,phi_n,&
 879 & green,inn,is_converged,&
 880 & hexc, hexc_i, comm)
 881 
 882 !Arguments ------------------------------------
 883 !scalars
 884  integer,intent(in) :: niter_max,niter_done,nomega
 885  integer,intent(in) :: my_t1,my_t2,term_type
 886  integer,intent(in) :: comm
 887  integer,intent(out) :: inn
 888  logical,intent(out) :: is_converged
 889  real(dp),intent(in) :: tol_iter
 890  complex(dpc),intent(in) :: factor
 891  type(hexc_t),intent(in) :: hexc
 892  type(hexc_interp_t),intent(in) :: hexc_i
 893 !arrays
 894  real(dp),intent(inout) :: bb(niter_max)
 895  complex(dpc),intent(out) :: green(nomega)
 896  complex(dpc),intent(in) :: omega(nomega)
 897  complex(dpc),intent(inout) :: aa(niter_max)
 898  complex(dpc),intent(inout) :: phi_nm1(my_t2-my_t1+1)
 899  complex(dpc),intent(inout) :: phi_n  (my_t2-my_t1+1)
 900  logical,intent(in) :: check(2)
 901 
 902 !Local variables ------------------------------
 903 !scalars
 904  integer :: ierr,my_nt,niter_min,nconv
 905  character(len=500) :: msg
 906  logical,parameter :: force_real=.TRUE.
 907 !arrays
 908  real(dp) :: abs_err(nomega,2) !,rel_err(nomega,2)
 909  complex(dpc),allocatable :: oldg(:),newg(:)
 910  complex(dpc),allocatable :: phi_np1(:),hphi_n(:),cfact(:)
 911  logical :: test(2)
 912 
 913 !************************************************************************
 914 
 915  ! The sequences starts with |1> normalized to 1 and b_0 =0, therefore:
 916  !  a_1 = <1|H|1>
 917  !  b_1 = || H|1> - a_1|1> ||
 918  !  |2> = [H|1> - a_1|1>]/b_1
 919  !
 920  ! For n>1 we have
 921  !  1) a_n = <n|H|n>
 922  !  2) b_n = || H|n> - a_n|n> -b_{n-1}|n-1> ||
 923  !  3) |n+1> = [H|n> -a_n|n> -b_{n-1}|n-1>]/b_n
 924  !
 925  my_nt = my_t2-my_t1+1
 926 
 927  ABI_MALLOC_OR_DIE(hphi_n,(hexc%hsize), ierr)
 928 
 929  ABI_MALLOC(phi_np1,(my_nt))
 930 
 931  ABI_MALLOC(oldg,(nomega))
 932  oldg=czero
 933  ABI_MALLOC(newg,(nomega))
 934  newg=czero
 935  ABI_MALLOC(cfact,(nomega))
 936  cfact=czero
 937 
 938  nconv=0
 939  do inn=niter_done+1,niter_max
 940 
 941    !YG2014
 942    call hexc_matmul_tda(hexc,hexc_i,phi_n,hphi_n)
 943 
 944    aa(inn) = xdotc(my_nt,phi_n,1,hphi_n(my_t1:),1)
 945    call xmpi_sum(aa(inn:inn),comm,ierr)
 946    if (force_real) aa(inn) = DBLE(aa(inn)) ! Matrix is Hermitian.
 947 
 948    ! |n+1> = H|n> - A(n)|n> - B(n-1)|n-1>
 949    phi_np1 = hphi_n(my_t1:my_t2) - aa(inn)*phi_n - bb(inn-1)*phi_nm1
 950 
 951    bb(inn) = xdotc(my_nt,phi_np1,1,phi_np1,1)
 952    call xmpi_sum(bb(inn),comm,ierr)
 953    bb(inn) = SQRT(bb(inn))
 954 
 955    phi_np1 = phi_np1/bb(inn)
 956 
 957    phi_nm1 = phi_n
 958    phi_n   = phi_np1
 959 
 960    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))
 961    call wrtout(std_out,msg,"COLL")
 962 
 963    call continued_fract(inn,term_type,aa,bb,nomega,omega,cfact)
 964 
 965    newg= factor*cfact
 966    !
 967    ! Avoid spurious convergence.
 968    niter_min=4; if (niter_done>1) niter_min=niter_done+1
 969    if (inn>niter_min) then
 970      test=.TRUE.
 971      abs_err(:,1) = ABS(DBLE (newg-oldg))
 972      abs_err(:,2) = ABS(AIMAG(newg-oldg))
 973      !
 974      if (tol_iter>zero) then
 975        ! Test on the L1 norm.
 976        if (check(1)) test(1) = SUM(abs_err(:,1)) < tol_iter*SUM(ABS(DBLE (newg)))
 977        if (check(2)) test(2) = SUM(abs_err(:,2)) < tol_iter*SUM(ABS(AIMAG(newg)))
 978      else
 979        ! Stringent test for each point.
 980        if (check(1)) test(1) = ALL( abs_err(:,1) < -tol_iter*ABS(DBLE (newg)))
 981        if (check(2)) test(2) = ALL( abs_err(:,2) < -tol_iter*ABS(AIMAG(newg)))
 982      end if
 983      !
 984      if (ALL(test)) then
 985        nconv = nconv+1
 986      else
 987        nconv = 0
 988      end if
 989      if (nconv==2) then
 990        if(inn<100)then
 991          write(msg,'(a,es10.2,a)')&
 992 &          " >>> Haydock algorithm converged twice within haydock_tol= ",tol_iter," after less than 100 iterations."
 993        else
 994          write(msg,'(a,es10.2,a,i0,a)')&
 995 &          " >>> Haydock algorithm converged twice within haydock_tol= ",tol_iter," after ",inn," iterations."
 996        endif
 997        call wrtout(std_out,msg,'COLL')
 998        call wrtout(ab_out,msg,'COLL')
 999        EXIT
1000      end if
1001    end if
1002 
1003    oldg = newg
1004  end do ! inn
1005 
1006  green = newg
1007  if (nconv/=2) then
1008    write(msg,'(a,es10.2,a,i0,a)')&
1009 &    " WARNING: Haydock algorithm did not converge within ",tol_iter," after ",niter_max," iterations."
1010    call wrtout(std_out,msg,'COLL')
1011    call wrtout(ab_out,msg,'COLL')
1012  end if
1013 
1014  is_converged = (nconv==2)
1015 
1016  ABI_FREE(oldg)
1017  ABI_FREE(newg)
1018  ABI_FREE(cfact)
1019  ABI_FREE(hphi_n)
1020  ABI_FREE(phi_np1)
1021 
1022 end subroutine haydock_herm_algo

m_haydock/haydock_mdf_to_tensor [ Functions ]

[ Top ] [ m_haydock ] [ Functions ]

NAME

 haydock_mdf_to_tensor

FUNCTION

 Transform macroscopic dielectric function from green function to each components of the tensor in red and cart coord.

INPUTS

  BSp<type(excparam)>=Parameters defining the Bethe-Salpeter calculation.
    omega(BSp%nomega)=Frequency mesh for the macroscopic dielectric function (broadening is already included).
  Cryst=Parameters of the crystal
  eps(BSp%nomega,BSp%nq) = Macroscopic dielectric function to be written.

OUTPUT

  tensor_cart(BSp%nomega,6) = dielectric tensor for each frequency, order (11,22,33,12,13,23) in cart. coord.
  tensor_red(BSp%nomega, 6) = idem in reduced coordinated
  ierr = 0 if the tensors have been successfully computed
      \= 0 if the system is ill-posed in terms of q-points (not enough or not independent q-points)

SOURCE

1164 subroutine haydock_mdf_to_tensor(BSp,Cryst,eps,tensor_cart,tensor_red,ierr)
1165 
1166 !Arguments ------------------------------------
1167 !scalars
1168  integer,intent(out) :: ierr
1169  type(excparam),intent(in) :: BSp
1170  type(crystal_t),intent(in) :: Cryst
1171 !arrays
1172  complex(dpc),intent(in) :: eps(BSp%nomega,BSp%nq)
1173  complex(dpc),intent(out) :: tensor_cart(BSp%nomega,6), tensor_red(BSp%nomega,6)
1174 
1175 !Local variables ------------------------------
1176 !scalars
1177  integer :: iq,info
1178  real(dp) :: normqcart, normqred
1179 !arrays
1180  integer,allocatable :: ipiv(:)
1181  real(dp) :: qcart(3), qtmet(3)
1182  real(dp) :: qred2cart(3,3),qcart2red(3,3)
1183  complex(dpc) :: qqcart(BSp%nq,6), qqred(BSp%nq,6)
1184  complex(dpc) :: b(6,BSP%nomega)
1185 
1186 !************************************************************************
1187 
1188  ! Error flag
1189  ierr = 0
1190 
1191  if(BSp%nq /= 6) then
1192     ierr = -1
1193     return
1194  end if
1195 
1196  ! Transformation matrices from reduced coordinates to cartesian coordinates
1197  qred2cart = two_pi*Cryst%gprimd
1198  qcart2red = qred2cart
1199  call matrginv(qcart2red,3,3)
1200  do iq = 1, 6
1201 
1202    ! Computing cartesian q-vector
1203    qcart = MATMUL(qred2cart, BSp%q(:,iq))
1204 
1205    ! Computing product 'metric - qred' to form quadratic form
1206    qtmet = (two_pi**2)*MATMUL(Cryst%gmet, BSp%q(:,iq))
1207 
1208    ! squared norms
1209    normqcart = qcart(1)**2+qcart(2)**2+qcart(3)**2
1210    normqred = (normv(BSp%q(:,iq),Cryst%gmet,"G"))**2
1211 
1212    ! Compute line 'iq' for matrix in cartesian coord
1213    qqcart(iq,1) = (qcart(1))**2
1214    qqcart(iq,2) = (qcart(2))**2
1215    qqcart(iq,3) = (qcart(3))**2
1216    qqcart(iq,4) = 2*(qcart(1)*qcart(2))
1217    qqcart(iq,5) = 2*(qcart(1)*qcart(3))
1218    qqcart(iq,6) = 2*(qcart(2)*qcart(3))
1219 
1220    ! Compute line 'iq' for matrix in reduced coord
1221    qqred(iq,1) = (qtmet(1))**2
1222    qqred(iq,2) = (qtmet(2))**2
1223    qqred(iq,3) = (qtmet(3))**2
1224    qqred(iq,4) = 2*(qtmet(1)*qtmet(2))
1225    qqred(iq,5) = 2*(qtmet(1)*qtmet(3))
1226    qqred(iq,6) = 2*(qtmet(2)*qtmet(3))
1227 
1228    ! Renormalize line
1229    qqcart(iq,:) = qqcart(iq,:)/normqcart
1230    qqred(iq,:) = qqred(iq,:)/normqred
1231  end do
1232 
1233  ABI_MALLOC(ipiv,(6))
1234 
1235  ! Solving linear system
1236  b = TRANSPOSE(eps)
1237  call ZGESV(6,BSp%nomega,qqcart,6,ipiv,b,6,info)
1238  tensor_cart = TRANSPOSE(b)
1239 
1240  if(info /= 0) then
1241    ! Skipping the rest of the routine
1242    ierr = info
1243    ABI_FREE(ipiv)
1244    return
1245  end if
1246 
1247  b = TRANSPOSE(eps)
1248  call ZGESV(6,BSp%nomega,qqred,6,ipiv,b,6,info)
1249  tensor_red = TRANSPOSE(b)
1250 
1251  if(info /= 0) ierr = info
1252 
1253  ABI_FREE(ipiv)
1254 
1255 end subroutine haydock_mdf_to_tensor

m_haydock/haydock_psherm [ Functions ]

[ Top ] [ m_haydock ] [ Functions ]

NAME

 haydock_psherm

FUNCTION

  Reads the excitonic Hamiltonian from file and construct the Lanczos set of vectors
  by iterative matrix-vector multiplications.

INPUTS

  BSp<type(excparam)>=Parameters defining the Bethe-Salpeter calculation.
    omega(BSp%nomega)=Frequency mesh for the macroscopic dielectric function (broadening is already included).
 hize
 my_t1,my_t2
 hreso(hsize,my_t1:my_t2)
 hcoup(hsize,my_t1:my_t2)
 nkets
 kets(hsize,nkets)
 comm=MPI communicator.

OUTPUT

  green(BSp%nomega)=The imaginary part of the macroscopic dielectric function.

SOURCE

1284 subroutine haydock_psherm(BSp,BS_files,Cryst,Hdr_bse,hexc,hexc_i,hsize,my_t1,my_t2,nkets,kets,green,comm)
1285 
1286 !Arguments ------------------------------------
1287 !scalars
1288  integer,intent(in) :: hsize,my_t1,my_t2,nkets,comm
1289  type(crystal_t),intent(in) :: Cryst
1290  type(excparam),intent(in) :: BSp
1291  type(excfiles),intent(in) :: BS_files
1292  type(Hdr_type),intent(in) :: Hdr_bse
1293 !arrays
1294  complex(dp),intent(out) :: green(BSp%nomega,BSp%nq)
1295  complex(dpc),intent(in) :: kets(hsize,nkets)
1296 
1297 !Local variables ------------------------------
1298 !scalars
1299  integer,parameter :: master=0
1300  integer :: inn,itt,out_unt,nproc,my_rank,ierr
1301  integer :: niter_file,niter_max,niter_done,nsppol,iq,my_nt,term_type
1302  real(dp) :: ket0_hbar_norm,nfact
1303  logical :: can_restart,is_converged
1304  complex(dpc) :: factor
1305  character(len=fnlen),parameter :: tag_file="_HAYDC_SAVE"
1306  character(len=500) :: msg
1307  character(len=fnlen) :: restart_file,out_file
1308  type(hexc_t),intent(in) :: hexc
1309  type(hexc_interp_t),intent(in) :: hexc_i
1310 !arrays
1311  real(dp),allocatable :: bb_file(:)
1312  real(dp),allocatable :: bb(:)
1313  complex(dpc),allocatable :: aa(:),cc(:),phi_np1(:),phi_n(:),phi_nm1(:),cbuff(:)
1314  complex(dpc),allocatable :: aa_file(:),phi_n_file(:),phi_np1_file(:),cc_file(:)
1315  complex(dpc),allocatable :: ket0(:)
1316  logical :: check(2)
1317 
1318 !************************************************************************
1319 
1320  ABI_WARNING("Haydock + coupling is still under development")
1321 
1322  if(BSp%use_interp) then
1323    ABI_ERROR("Coupling is not yet implemented with interpolation")
1324  end if
1325 
1326  nproc  = xmpi_comm_size(comm)
1327  my_rank= xmpi_comm_rank(comm)
1328  nsppol = Hdr_bse%nsppol
1329 
1330  my_nt = my_t2-my_t1+1
1331  ABI_CHECK(my_nt>0,"One of the processors has zero columns")
1332 
1333  ! Multiplicative factor (k-point sampling and unit cell volume)
1334  ! TODO be careful with the spin here
1335  ! TODO four_pi comes from the coulomb term 1/|q| is already included in the
1336  ! oscillators hence the present approach wont work if a cutoff interaction is used.
1337  nfact = four_pi/(Cryst%ucvol*BSp%nkbz)
1338  if (nsppol==1) nfact=two*nfact
1339 
1340  write(msg,'(a,i0)')' Haydock algorithm with MAX number of iterations: ',BSp%niter
1341  call wrtout(std_out,msg,"COLL")
1342  !
1343  ! Check for presence of the restart file.
1344  can_restart=.FALSE.
1345 
1346  if ( BS_files%in_haydock_basename /= BSE_NOFILE) then
1347    restart_file = strcat(BS_files%in_haydock_basename,tag_file)
1348    if (file_exists(restart_file) ) then
1349      can_restart=.TRUE.
1350      msg = strcat(" Restarting Haydock calculation from file: ",restart_file)
1351      call wrtout(std_out,msg,"COLL")
1352      call wrtout(ab_out,msg,"COLL")
1353      ABI_ERROR("Restart is not tested")
1354    else
1355      can_restart=.FALSE.
1356      ABI_WARNING(strcat("Cannot find restart file: ",restart_file))
1357    end if
1358  end if
1359  !
1360  ! Open the file and writes basic dimensions and info.
1361  if (my_rank==master) then
1362    out_file = TRIM(BS_files%out_basename)//TRIM(tag_file)
1363    if (open_file(out_file,msg,newunit=out_unt,form="unformatted") /= 0) then
1364      ABI_ERROR(msg)
1365    end if
1366    ! write header TODO: standardize this part.
1367    write(out_unt)hsize,Bsp%use_coupling,BSE_HAYD_IMEPS,nkets,Bsp%broad
1368  end if
1369  !
1370  ! Select the terminator for the continued fraction.
1371  term_type=0 !; if (Bsp%hayd_term>0) term_type=2
1372  call wrtout(std_out,sjoin("Using terminator type: ",itoa(term_type)),"COLL")
1373  !
1374  ! Calculate green(w) for the different starting kets.
1375  green=czero
1376  do iq=1,nkets
1377    ABI_MALLOC(ket0,(my_nt))
1378    ket0 = kets(my_t1:my_t2,iq)
1379    !
1380    niter_file=0
1381 
1382    if (can_restart) then
1383      call haydock_restart(BSp,restart_file,BSE_HAYD_IMEPS,iq,hsize,&
1384 &      niter_file,aa_file,bb_file,phi_np1_file,phi_n_file,comm)
1385    end if
1386    !
1387    ABI_MALLOC(phi_nm1,(my_nt))
1388    ABI_MALLOC(phi_n,(my_nt))
1389    ABI_MALLOC(phi_np1,(my_nt))
1390    !
1391    ! TODO: Note the different convention used for the coefficients
1392    ! Should use the same convention in the Hermitian case.
1393    niter_max = niter_file + Bsp%niter
1394    ABI_MALLOC(aa,(niter_max))
1395    ABI_MALLOC(bb,(niter_max+1))
1396    ABI_MALLOC(cc,(niter_max+1))
1397    aa=czero; bb=czero; cc=czero
1398 
1399    if (niter_file==0) then ! Calculation from scratch.
1400      phi_n   = ket0
1401      call hexc_matmul_full(hexc, hexc_i, phi_n, phi_np1, -1)
1402      !phi_np1 = MATMUL(hreso,ket0) - MATMUL(hcoup,CONJG(ket0))
1403      ket0_hbar_norm = SQRT(two*DBLE(DOT_PRODUCT(phi_n,phi_np1)))
1404      phi_n   = phi_n  /ket0_hbar_norm
1405      phi_np1 = phi_np1/ket0_hbar_norm
1406      !ket0    = ket0/ket0_hbar_norm
1407      cc(1)=zero ! <P|F|P>
1408      !cc(1) =  DOT_PRODUCT(ket0,phi_np1)
1409      write(std_out,*)" cc(1), ket0_hbar_norm =",cc(1),ket0_hbar_norm
1410 
1411      phi_nm1 = czero
1412      niter_done=0  ! TODO Be careful here
1413 
1414    else ! Use the previously calculates a and b.
1415      niter_done=niter_file
1416      ABI_ERROR("Restart not coded")
1417      !aa(1:niter_done) = aa_file
1418      !bb(1:niter_done) = bb_file
1419      !phi_np1=phi_np1_file(my_t1:my_t2)   ! Select the slice treated by this node.
1420      !phi_n  =phi_n_file  (my_t1:my_t2)
1421    end if
1422 
1423    if (can_restart) then
1424      ABI_FREE(aa_file)
1425      ABI_FREE(bb_file)
1426      ABI_FREE(cc_file)
1427      ABI_FREE(phi_np1_file)
1428      ABI_FREE(phi_n_file)
1429    end if
1430 
1431    ! This factor gives the correct results
1432    factor = -nfact*ket0_hbar_norm / SQRT(two)
1433 
1434    ! Which quantity should be checked for convergence?
1435    check = (/.TRUE.,.TRUE./)
1436    if (ABS(Bsp%haydock_tol(2)-one)<tol6) check = (/.TRUE. ,.FALSE./)
1437    if (ABS(Bsp%haydock_tol(2)-two)<tol6) check = (/.FALSE.,.TRUE./)
1438 
1439    call haydock_psherm_optalgo(niter_done,niter_max,BSp%nomega,BSp%omega,BSp%haydock_tol(1),check,hexc,hexc_i,&
1440 &    hsize,my_t1,my_t2,factor,term_type,aa,bb,cc,ket0,ket0_hbar_norm,phi_nm1,phi_n,phi_np1,&
1441 &    green(:,iq),inn,is_converged,comm)
1442    !
1443    ! Save the a"s and the b"s for possible restarting.
1444    ! 1) Info on the Q.
1445    ! 2) Number of iterations performed.
1446    ! 3) do iter=1,niter_performed
1447    !      aa(iter),bb(iter)
1448    !    end do
1449    ! 4) |n-1>
1450    !    |n>
1451    !    |n+1>
1452    !
1453    if (my_rank==master) then ! Open the file and writes basic dimensions and info.
1454      write(out_unt)Bsp%q(:,iq)
1455      write(out_unt)MIN(inn,niter_max)  ! NB: if the previous loop completed inn=niter_max+1
1456      do itt=1,MIN(inn,niter_max)        !     if we exited then inn is not incremented by one.
1457        write(out_unt)itt,aa(itt),bb(itt)
1458      end do
1459    end if
1460    !
1461    ! cbuff is used as workspace to gather |n-1>, |n> and |n+1>.
1462    ABI_MALLOC(cbuff,(hsize))
1463    cbuff=czero; cbuff(my_t1:my_t2) = phi_nm1
1464    call xmpi_sum_master(cbuff,master,comm,ierr)
1465    if (my_rank==master) write(out_unt) cbuff ! |n-1>
1466 
1467    cbuff=czero; cbuff(my_t1:my_t2) = phi_n
1468    call xmpi_sum_master(cbuff,master,comm,ierr)
1469    if (my_rank==master) write(out_unt) cbuff ! |n>
1470 
1471    cbuff=czero; cbuff(my_t1:my_t2) = phi_np1
1472    call xmpi_sum_master(cbuff,master,comm,ierr)
1473    if (my_rank==master) write(out_unt) cbuff ! |n+1>
1474 
1475    ABI_FREE(phi_nm1)
1476    ABI_FREE(phi_n)
1477    ABI_FREE(phi_np1)
1478    ABI_FREE(cbuff)
1479    ABI_FREE(aa)
1480    ABI_FREE(bb)
1481    ABI_FREE(cc)
1482    ABI_FREE(ket0)
1483  end do ! iq
1484 
1485  if (my_rank==master) close(out_unt)
1486 
1487  call xmpi_barrier(comm)
1488 
1489 end subroutine haydock_psherm

m_haydock/haydock_psherm_optalgo [ Functions ]

[ Top ] [ m_haydock ] [ Functions ]

NAME

 haydock_psherm_optalgo

FUNCTION

  Haydock algorithm for pseudo-hermitian matrix

INPUTS

  niter_done=Number of iterations already performed (0 if the run starts from scratch).
  niter_tot=Max number of iterations. Always > niter_done
  nomega=Number of Frequency points for the evaluation of the matrix element.
  omega(nomega)=Frequency set (imaginary part is already included).
  tol_iter=Tollerance used to stop the the algorithm.
  check(2)=Logical flags to specify where both the real and the imaginary part of the
    matrix elements of the Green functions have to be checked for convergence.
  hsize=Size of the blocks.
  my_t1,my_t2=Indeces of the first and last column stored treated by this done.
  term_type=0 if no terminator is used, 1 otherwise.
  hreso(hsize,my_t1:my_t2)=The columns of the resonant block.
  hcoup(hsize,my_t1:my_t2)=The columns of the coupling block.
  factor
  comm=MPI communicator.

OUTPUT

  green(nomega)=Output matrix elements.
  inn=Last iteration performed.
  is_converged=.TRUE. of the algorithm converged.

SIDE EFFECTS

  phi_nm1(my_t2-my_t1+1), phi_n(my_t2-my_t1+1)
    input: vectors used to initialize the iteration
    output: the vectors obtained in the last iteration
  aa(niter_tot) and bb(niter_tot+1)
    if niter_done>0: aa(1:niter_done), bb(1:niter_done) store the coefficients of the previous run.
    when the routine returns aa(1:inn) and bb(1:inn) contain the matrix elements of the tridiagonal form.
  cc(niter_tot+1)

SOURCE

1533 subroutine haydock_psherm_optalgo(niter_done,niter_tot,nomega,omega,tol_iter,check,hexc,hexc_i,hsize,my_t1,my_t2,&
1534 &  factor,term_type,aa,bb,cc,ket0,ket0_hbar_norm,phi_nm1,phi_n,phi_np1,green,inn,is_converged,comm)
1535 
1536 !Arguments ------------------------------------
1537 !scalars
1538  integer,intent(in) :: niter_tot,niter_done,nomega,comm,hsize,my_t1,my_t2,term_type
1539  integer,intent(out) :: inn
1540  logical,intent(out) :: is_converged
1541  real(dp),intent(in) :: tol_iter,ket0_hbar_norm
1542  complex(dpc),intent(in) :: factor
1543  type(hexc_t),intent(in) :: hexc
1544  type(hexc_interp_t),intent(in) :: hexc_i
1545 !arrays
1546  real(dp),intent(inout) :: bb(niter_tot+1)
1547  complex(dpc),intent(out) :: green(nomega)
1548  complex(dpc),intent(in) :: omega(nomega)
1549  complex(dpc),intent(inout) :: aa(niter_tot),cc(niter_tot+1)
1550  complex(dpc),intent(in) :: ket0(my_t2-my_t1+1)
1551  complex(dpc),intent(inout) :: phi_nm1(my_t2-my_t1+1)
1552  complex(dpc),intent(inout) :: phi_n  (my_t2-my_t1+1)
1553  complex(dpc),intent(inout) :: phi_np1(my_t2-my_t1+1)
1554  logical,intent(in) :: check(2)
1555 
1556 !Local variables ------------------------------
1557 !scalars
1558  integer :: my_nt,niter_min,nconv,parity,ii,jj,tdim,ierr
1559  integer :: row_max,col_max,nlev
1560  character(len=500) :: msg
1561  real(dp) :: max_err,mean_err,mean_err2,std_dev,err
1562  logical :: keep_vectors=.TRUE.
1563 !arrays
1564  real(dp) :: abs_err(nomega,2) !,ww_err(nomega,2)
1565  complex(dpc) :: gn0(nomega,niter_tot)
1566  complex(dpc),allocatable :: oldg(:),newg(:)
1567  complex(dpc),allocatable :: hphi_n(:),save_phi(:,:)
1568  complex(dpc),allocatable ::  alpha(:,:),beta(:,:),ovlp(:,:)
1569  complex(dpc),allocatable :: phi_test(:),phi_test2(:),g00(:)
1570  logical :: test(2)
1571 
1572 !************************************************************************
1573 
1574  ABI_UNUSED(ket0_hbar_norm)
1575 
1576  my_nt = my_t2-my_t1+1
1577 
1578  ABI_MALLOC(oldg,(nomega))
1579  ABI_MALLOC(newg,(nomega))
1580  ABI_MALLOC(g00,(nomega))
1581  oldg=czero; newg=czero; g00=czero
1582  nconv=0
1583 
1584  keep_vectors = (keep_vectors.and.xmpi_comm_size(comm)==1)
1585  if (keep_vectors) then
1586    ABI_MALLOC_OR_DIE(save_phi,(my_t2-my_t1+1,niter_tot), ierr)
1587    save_phi=czero
1588  end if
1589 
1590  ABI_MALLOC(hphi_n,(hsize))
1591 
1592  do inn=niter_done+1,niter_tot
1593    !
1594    ! a(n) = <Vn+1|F|Vn+1> = <Vn|HFH|Vn>) = 0 by symmetry.
1595    aa(inn)=zero
1596 
1597    ! |n+1> = |n+1> - a(n)|Vn> - a(n)|n-1>
1598    phi_np1 = phi_np1 - bb(inn)*phi_nm1
1599    !
1600    ! |n-1> = |n>
1601    ! |n>   = |n+1>
1602    phi_nm1 = phi_n
1603    phi_n   = phi_np1
1604    !
1605    !|n+1> = H |n> using all eh components.
1606    parity = (-1)**(inn+1)
1607    call hexc_matmul_full(hexc, hexc_i, phi_n, phi_np1, parity)
1608 
1609    !phi_np1 = MATMUL(hreso,phi_n) + parity * MATMUL(hcoup,CONJG(phi_n))
1610    !call xmpi_sum(hphi_np1,comm,ierr)
1611    !
1612    ! B(n+1)= <n|F|n+1>^(1/2) = <n|FH|n>^(1/2))= (2*Re(<n|V+1>))^(1/2)
1613    ! by symmetry, where the dot_product is done in the resonant eh sub-space.
1614    !
1615    bb(inn+1)=SQRT(two*DBLE(DOT_PRODUCT(phi_n,phi_np1)))
1616    !bb(inn+1)=two*DBLE(DOT_PRODUCT(phi_n,phi_np1))
1617    !call xmpi_sum(bb(inn+1),comm,ierr)
1618    !bb(inn+1)=SQRT(bb(inn+1)
1619    !
1620    !|n+1> =|n+1>/B(n+1)
1621    phi_n   = phi_n  /bb(inn+1)
1622    phi_np1 = phi_np1/bb(inn+1)
1623 
1624    if (keep_vectors) save_phi(:,inn) = phi_n
1625 
1626    parity = (-1)**(inn+1)
1627    !if (parity==-1) then
1628    !  cc(inn+1)=czero
1629    !else
1630      cc(inn+1)=DOT_PRODUCT(ket0,phi_n) + parity * DOT_PRODUCT(phi_n,ket0)
1631    !end if
1632    !call xmpi_sum(cc(inn+1),comm,ierr)
1633 
1634    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))
1635    call wrtout(std_out,msg,"COLL")
1636 
1637    call continued_fract(inn,term_type,aa,bb(2:),nomega,omega,g00)
1638    gn0(:,1) = g00
1639 
1640    if (.FALSE.) then
1641      gn0(:,2) = (one - omega(:)*g00(:))/bb(2)
1642      do ii=3,inn
1643        gn0(:,ii) = -(-bb(ii)*gn0(:,ii-2) -omega(:)*gn0(:,ii-1))/bb(ii+1)
1644      end do
1645    else
1646      do ii=2,inn
1647        nlev = inn-ii
1648        call continued_fract(nlev,term_type,aa,bb(ii+1:),nomega,omega,g00)
1649        gn0(:,ii) = +bb(ii+1) * g00 * gn0(:,ii-1)
1650      end do
1651    end if
1652 
1653    newg=czero
1654    do ii=1,inn
1655      newg(:) = newg + cc(ii)* gn0(:,ii)
1656    end do
1657    newg = factor*newg
1658    !
1659    ! Avoid spurious convergence.
1660    niter_min=4; if (niter_done>1) niter_min=niter_done+1
1661    if (inn>niter_min) then
1662      test=.TRUE.
1663      abs_err(:,1) = ABS(DBLE (newg-oldg))
1664      abs_err(:,2) = ABS(AIMAG(newg-oldg))
1665      !
1666      if (tol_iter>zero) then
1667        ! Test on the L1 norm.
1668        if (check(1)) test(1) = SUM(abs_err(:,1)) < tol_iter*SUM(ABS(DBLE (newg)))
1669        if (check(2)) test(2) = SUM(abs_err(:,2)) < tol_iter*SUM(ABS(AIMAG(newg)))
1670      else
1671        ! Stringent test for each point.
1672        if (check(1)) test(1) = ALL( abs_err(:,1) < -tol_iter*ABS(DBLE (newg)))
1673        if (check(2)) test(2) = ALL( abs_err(:,2) < -tol_iter*ABS(AIMAG(newg)))
1674      end if
1675      !
1676      if (ALL(test)) then
1677        nconv = nconv+1
1678      else
1679        nconv = 0
1680      end if
1681      if (nconv==2) then
1682        if(inn<100)then
1683          write(msg,'(a,es10.2,a)')&
1684 &          " >>> Haydock algorithm converged twice within haydock_tol= ",tol_iter," after less than 100 iterations."
1685        else
1686          write(msg,'(a,es10.2,a,i0,a)')&
1687 &          " >>> Haydock algorithm converged twice within haydock_tol= ",tol_iter," after ",inn," iterations."
1688        endif
1689        call wrtout(std_out,msg,'COLL')
1690        call wrtout(ab_out,msg,'COLL')
1691        EXIT
1692      end if
1693    end if
1694    !
1695    oldg = newg
1696  end do ! inn
1697 
1698  green = newg
1699  if (nconv/=2) then
1700    write(msg,'(a,es10.2,a,i0,a)')&
1701 &    " WARNING: Haydock algorithm did not converge within ",tol_iter," after ",niter_tot," iterations."
1702    call wrtout(std_out,msg,'COLL')
1703    call wrtout(ab_out,msg,'COLL')
1704  end if
1705 
1706  is_converged = (nconv==2)
1707 
1708  ABI_FREE(oldg)
1709  ABI_FREE(newg)
1710  ABI_FREE(g00)
1711  ABI_FREE(hphi_n)
1712 
1713  if (keep_vectors) then
1714    tdim = MIN(inn,niter_tot)
1715    ABI_MALLOC(ovlp,(tdim,tdim))
1716 
1717    ABI_MALLOC(phi_test,(hsize))
1718    ABI_MALLOC(phi_test2,(hsize))
1719 
1720    max_err=smallest_real; mean_err=zero; mean_err2=zero; row_max=-1
1721    do ii=1,tdim
1722      parity = (-1)**(ii+1)
1723      phi_test  = save_phi(:,ii)
1724      call hexc_matmul_full(hexc, hexc_i, phi_test, phi_test2, parity)
1725      !phi_test2 = MATMUL(hreso,phi_test) + parity * MATMUL(hcoup,CONJG(phi_test))
1726      ovlp(ii,ii) = DOT_PRODUCT(phi_test,phi_test2) + DOT_PRODUCT(phi_test2,phi_test)
1727      err = ABS(ovlp(ii,ii)-cone)
1728      mean_err  = mean_err + err
1729      mean_err2 = mean_err2 + err**2
1730      if (err > max_err) then
1731        max_err = err
1732        row_max = ii
1733      end if
1734    end do
1735    mean_err = mean_err/tdim
1736    std_dev = mean_err2/tdim -mean_err**2
1737    write(std_out,'(a,i0,1x,3es14.6)')&
1738 &   " Error in normalization (ii, max_err,mean,std_dev): ",row_max,max_err,mean_err,std_dev
1739 
1740    ABI_FREE(phi_test)
1741    ABI_FREE(phi_test2)
1742 
1743    ABI_MALLOC(alpha,(hsize,tdim))
1744 
1745    ! Less efficient but for sake of simplicity with hexc_matmul
1746    ! TODO possibility to call hreso * phi, and hcoup * phi separately
1747    do ii=1,tdim
1748      parity = (-1)**(ii+1)
1749      call hexc_matmul_full(hexc, hexc_i, save_phi(:,ii), alpha(:,ii), parity)
1750    end do
1751 
1752    !alpha = MATMUL(hreso,save_phi(:,1:tdim))
1753    !
1754    !do ii=1,tdim
1755    !  parity = (-1)**(ii+1)
1756    !  alpha(:,ii) =  alpha(:,ii) + parity*MATMUL(hcoup,CONJG(save_phi(:,ii)))
1757    !end do
1758 
1759    ovlp = MATMUL(TRANSPOSE(CONJG(save_phi(:,1:tdim))),alpha)
1760 
1761    ABI_MALLOC(beta,(hsize,tdim))
1762    do ii=1,tdim
1763      parity = (-1)**(ii+1)
1764      beta(:,ii)  =  parity*save_phi(:,ii)
1765      alpha(:,ii) = -parity*alpha(:,ii)
1766    end do
1767 
1768    ovlp = ovlp - MATMUL(TRANSPOSE(CONJG(beta)),alpha)
1769 
1770    max_err=smallest_real; row_max=-1; col_max=-1
1771    mean_err=zero; mean_err2=zero
1772    do jj=1,tdim
1773      do ii=1,jj
1774        err = ABS(ovlp(ii,jj))
1775        if (ii==jj) err = ABS(err - one)
1776        mean_err  = mean_err + err
1777        mean_err2 = mean_err2 + err**2
1778        if (err > max_err) then
1779          max_err = err
1780          row_max=ii
1781          col_max=jj
1782        end if
1783      end do
1784    end do
1785 
1786    mean_err = mean_err/(tdim*(tdim+1)/2)
1787    std_dev = mean_err2/(tdim*(tdim+1)/2) - mean_err**2
1788    write(std_out,'(a,2(i0,1x),3es14.6)')&
1789 &     " Error in Hbar-ortho (i,j), max_err, mean, std_dev ",row_max,col_max,max_err,mean_err,std_dev
1790    !call print_arr(ovlp,max_r=185,max_c=10,unit=std_out)
1791 
1792    ABI_FREE(alpha)
1793    ABI_FREE(beta)
1794    ABI_FREE(ovlp)
1795    ABI_FREE(save_phi)
1796  end if
1797 
1798 end subroutine haydock_psherm_optalgo

m_haydock/haydock_restart [ Functions ]

[ Top ] [ m_haydock ] [ Functions ]

NAME

 haydock_restart

FUNCTION

 Restart the Haydock method from file reading the data produced in a previous run.

INPUTS

  BSp<type(excparam)>=Parameters defining the Bethe-Salpeter calculation.
    omega(BSp%nomega)=Frequency mesh for the macroscopic dielectric function (broadening is already included).
  iq_search=The index of the q-point to be searched.
  hsize
  comm=MPI communicator.
  nsppol
  restart_file

OUTPUT

  niter_file=Number of iterations already performed. 0 to signal that an error occurred during the reading
  bb_file(:)
  aa_file(:)
  phi_n_file(:)
  phi_nm1_file(:)

SOURCE

1052 subroutine haydock_restart(BSp,restart_file,ftype,iq_search,hsize,niter_file,aa_file,bb_file,phi_nm1_file,phi_n_file,comm)
1053 
1054 !Arguments ------------------------------------
1055 !scalars
1056  integer,intent(in) :: comm,hsize,iq_search,ftype
1057  integer,intent(out) :: niter_file
1058  character(len=*),intent(in) :: restart_file
1059  type(excparam),intent(in) :: BSp
1060 !arrays
1061  real(dp),allocatable,intent(out) :: bb_file(:)
1062  complex(dpc),allocatable,intent(out) :: aa_file(:),phi_n_file(:),phi_nm1_file(:)
1063 
1064 !Local variables ------------------------------
1065 !scalars
1066  integer,parameter :: master=0
1067  integer :: nproc,my_rank,ierr,op_file
1068  integer :: hsize_file,use_coupling_file
1069  complex(dpc) :: factor_file
1070  character(len=500) :: msg
1071  type(haydock_type) :: haydock_file
1072 
1073 !************************************************************************
1074 
1075  nproc = xmpi_comm_size(comm); my_rank= xmpi_comm_rank(comm)
1076 
1077  if (my_rank==master) then
1078    call open_haydock(restart_file, haydock_file)
1079 
1080    call read_dim_haydock(haydock_file)
1081 
1082    if (haydock_file%op/=ftype) then
1083      write(msg,"(2(a,i0))")" Expecting restart file with filetype: ",ftype," but found ",op_file
1084      ABI_ERROR(msg)
1085    end if
1086 
1087    if (haydock_file%hsize/=hsize) then
1088      write(msg,"(2(a,i0))")&
1089 &      " Rank of H_exc read from file: ",hsize_file," differs from the one used in this run: ",hsize
1090      ABI_ERROR(msg)
1091    end if
1092 
1093    if (haydock_file%use_coupling /= BSp%use_coupling) then
1094      write(msg,'(2(a,i0))')&
1095 &      " use_coupling_file: ",use_coupling_file," differs from input file value: ",BSp%use_coupling
1096      ABI_ERROR(msg)
1097    end if
1098 
1099    call read_haydock(haydock_file, Bsp%q(:,iq_search), aa_file, bb_file, &
1100 &                   phi_n_file, phi_nm1_file, niter_file, factor_file)
1101 
1102    if (niter_file == 0) then
1103      write(msg,"(a,3f8.4,3a)")&
1104 &      " Could not find q-point: ",BSp%q(:,iq_search)," in file ",TRIM(restart_file),&
1105 &      " Cannot restart Haydock iterations for this q-point"
1106      ABI_COMMENT(msg)
1107    else
1108      write(msg,'(a,i0)')" Number of iterations already performed: ",niter_file
1109      call wrtout(std_out,msg,"COLL")
1110      call wrtout(ab_out,msg,"COLL")
1111 
1112      if ( ABS(haydock_file%broad - BSp%broad) > tol6) then
1113        write(msg,'(2a,2(a,f8.4),a)')&
1114 &        " Restart file has been produced with a different Lorentzian broadening: ",ch10,&
1115 &        " broad_file: ",haydock_file%broad," input broadening: ",BSp%broad," Continuing anyway. "
1116        ABI_WARNING(msg)
1117      end if
1118 
1119      call close_haydock(haydock_file)
1120    end if
1121  end if
1122  !
1123  ! Master broadcasts the data.
1124  call xmpi_bcast(niter_file,master,comm,ierr)
1125 
1126  if (my_rank/=master) then
1127    ABI_MALLOC(aa_file,(niter_file))
1128    ABI_MALLOC(bb_file,(niter_file))
1129    ABI_MALLOC(phi_nm1_file,(hsize))
1130    ABI_MALLOC(phi_n_file,(hsize))
1131  end if
1132 
1133  call xmpi_bcast(aa_file,master,comm,ierr)
1134  call xmpi_bcast(bb_file,master,comm,ierr)
1135  call xmpi_bcast(phi_nm1_file,master,comm,ierr)
1136  call xmpi_bcast(phi_n_file,master,comm,ierr)
1137 
1138 end subroutine haydock_restart

m_numeric_tools/continued_fract_general [ Functions ]

[ Top ] [ m_numeric_tools ] [ Functions ]

NAME

  continued_fract

FUNCTION

  This routine calculates the continued fraction:

                        1
 f(z) =  _______________________________
           z - a1 -        b1^2
                   _____________________
                     z - a2 -    b2^2
                             ___________
                                z -a3 -    ........

INPUTS

  nlev=Number of "levels" in the continued fraction.
  term_type=Type of the terminator.
    0 --> No terminator.
   -1 --> Assume constant coefficients for a_i and b_i for i>nlev with a_inf = a(nlev) and b_inf = b(nleb)
    1 --> Same as above but a_inf and b_inf are obtained by averaging over the nlev values.
  aa(nlev)=Set of a_i coefficients.
  bb(nlev)=Set of b_i coefficients.
  nz=Number of points on the z-mesh.
  zpts(nz)=z-mesh.

OUTPUT

  spectrum(nz)=Contains f(z) on the input mesh.

SOURCE

2419 subroutine continued_fract_general(nlev,term_type,aa,bb,cc,nz,zpts,spectrum)
2420 
2421 !Arguments ------------------------------------
2422 !scalars
2423  integer,intent(in) :: nlev,term_type,nz
2424 !arrays
2425  complex(dpc),intent(in) :: bb(nlev)
2426  complex(dpc),intent(in) :: cc(nlev)
2427  complex(dpc),intent(in) :: aa(nlev)
2428  complex(dpc),intent(in) :: zpts(nz)
2429  complex(dpc),intent(out) :: spectrum(nz)
2430 
2431 !Local variables ------------------------------
2432 !scalars
2433  integer :: it
2434  complex(dpc) ::  bb_inf,bg,bu,swap
2435  complex(dpc) :: aa_inf
2436  character(len=500) :: msg
2437 !arrays
2438  complex(dpc),allocatable :: div(:),den(:)
2439 
2440 !************************************************************************
2441 
2442  ABI_MALLOC(div,(nz))
2443  ABI_MALLOC(den,(nz))
2444 
2445  select case (term_type)
2446  case (0) ! No terminator.
2447    div=czero
2448  case (-1,1)
2449    ABI_ERROR("Not yet implemented")
2450    if (term_type==-1) then
2451      bb_inf=bb(nlev)
2452      aa_inf=aa(nlev)
2453    else
2454      bb_inf=SUM(bb)/nlev
2455      aa_inf=SUM(aa)/nlev
2456    end if
2457    ! Be careful with the sign of the SQRT.
2458    div(:) = half*(bb(nlev)/(bb_inf))**2 * ( zpts-aa_inf - SQRT((zpts-aa_inf)**2 - four*bb_inf**2) )
2459  case (2)
2460    ABI_ERROR("Not yet implemented")
2461    div = zero
2462    if (nlev>4) then
2463      bg=zero; bu=zero
2464      do it=1,nlev,2
2465        if (it+2<nlev) bg = bg + bb(it+2)
2466        bu = bu + bb(it)
2467      end do
2468      bg = bg/(nlev/2+MOD(nlev,2))
2469      bu = bg/((nlev+1)/2)
2470      !if (iseven(nlev)) then
2471      if (.not.iseven(nlev)) then
2472        swap = bg
2473        bg = bu
2474        bu = bg
2475      end if
2476      !write(std_out,*)nlev,bg,bu
2477      !Here be careful with the sign of SQRT
2478      do it=1,nz
2479        div(it) = half/zpts(it) * (bb(nlev)/bu)**2 * &
2480 &        ( (zpts(it)**2 +bu**2 -bg**2) - SQRT( (zpts(it)**2+bu**2-bg**2)**2 -four*(zpts(it)*bu)**2) )
2481      end do
2482    end if
2483 
2484  case default
2485    write(msg,'(a,i0)')" Wrong value for term_type : ",term_type
2486    ABI_ERROR(msg)
2487  end select
2488 
2489  do it=nlev,2,-1
2490    den(:) = zpts(:) - aa(it) - div(:)
2491    div(:) = (bb(it-1)*cc(it-1) )/ den(:)
2492  end do
2493 
2494  den = zpts(:) - aa(1) - div(:)
2495  div = one/den(:)
2496 
2497  spectrum = div
2498  ABI_FREE(div)
2499  ABI_FREE(den)
2500 
2501 end subroutine continued_fract_general