TABLE OF CONTENTS


ABINIT/m_haydock [ Modules ]

[ Top ] [ Modules ]

NAME

 m_haydock

FUNCTION

COPYRIGHT

  Copyright (C) 2008-2018 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 .

PARENTS

CHILDREN

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 MODULE m_haydock
27 
28  use defs_basis
29  use m_abicore
30  use m_bs_defs
31  use m_xmpi
32  use m_errors
33  use m_nctk
34  use m_haydock_io
35  use m_linalg_interfaces
36  use m_ebands
37 #ifdef HAVE_NETCDF
38  use netcdf
39 #endif
40 
41  use m_time,              only : timab
42  use m_fstrings,          only : indent, strcat, sjoin, itoa, int2char4
43  use m_io_tools,          only : file_exists, open_file
44  use defs_abitypes,       only : Hdr_type
45  use defs_datatypes,      only : ebands_t, pseudopotential_type
46  use m_geometry,          only : normv
47  use m_hide_blas,         only : xdotc, xgemv
48  use m_hide_lapack,       only : matrginv
49  use m_numeric_tools,     only : print_arr, symmetrize, hermitianize, continued_fract, wrap2_pmhalf, iseven
50  use m_fft_mesh,          only : calc_ceigr
51  use m_kpts,              only : listkk
52  use m_crystal,           only : crystal_t
53  use m_crystal_io,        only : crystal_ncwrite
54  use m_bz_mesh,           only : kmesh_t, findqg0, get_bz_item
55  use m_double_grid,       only : double_grid_t, get_kpt_from_indices_coarse, compute_corresp
56  use m_paw_hr,            only : pawhur_t
57  use m_wfd,               only : wfd_t, wfd_sym_ur, wfd_get_ur, wfd_change_ngfft, wfd_wave_free
58  use m_bse_io,            only : exc_read_rcblock, exc_write_optme
59  use m_pawtab,            only : pawtab_type
60  use m_vcoul,             only : vcoul_t
61  use m_hexc,              only : hexc_init, hexc_interp_init, hexc_free, hexc_interp_free, &
62 &                                hexc_build_hinterp, hexc_matmul_tda, hexc_matmul_full, hexc_t, hexc_matmul_elphon, hexc_interp_t
63  use m_exc_spectra,       only : exc_write_data, exc_eps_rpa, exc_write_tensor, mdfs_ncwrite
64  use m_eprenorms,         only : eprenorms_t, renorm_bst
65  use m_wfd_optic,         only : calc_optical_mels
66 
67  implicit none
68 
69  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<wfd_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 LDA+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

PARENTS

      bethe_salpeter

CHILDREN

SOURCE

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

PARENTS

      m_haydock

CHILDREN

SOURCE

1941 subroutine haydock_bilanczos(BSp,BS_files,Cryst,Hdr_bse,hexc,hexc_i,hsize,my_t1,my_t2,nkets,kets,ep_renorms,green,comm)
1942 
1943 
1944 !This section has been created automatically by the script Abilint (TD).
1945 !Do not modify the following lines by hand.
1946 #undef ABI_FUNC
1947 #define ABI_FUNC 'haydock_bilanczos'
1948 !End of the abilint section
1949 
1950  implicit none
1951 
1952 !Arguments ------------------------------------
1953 !scalars
1954  integer,intent(in) :: hsize,my_t1,my_t2,nkets,comm
1955  type(crystal_t),intent(in) :: Cryst
1956  type(excparam),intent(in) :: BSp
1957  type(excfiles),intent(in) :: BS_files
1958  type(Hdr_type),intent(in) :: Hdr_bse
1959 !arrays
1960  complex(dp),intent(out) :: green(BSp%nomega,BSp%nq)
1961  complex(dpc),intent(in) :: kets(hsize,nkets)
1962  complex(dpc),intent(in) :: ep_renorms(hsize)
1963 
1964 !Local variables ------------------------------
1965 !scalars
1966  integer,parameter :: master=0
1967  integer :: inn,itt,out_unt,nproc,my_rank,ierr
1968  integer :: niter_file,niter_max,niter_done,nsppol,iq,my_nt,term_type,n_all_omegas
1969  real(dp) :: ket0_hbar_norm,nfact,norm
1970  logical :: can_restart,is_converged
1971  complex(dpc) :: factor
1972  character(len=fnlen),parameter :: tag_file="_HAYDC_SAVE"
1973  character(len=500) :: msg
1974  character(len=fnlen) :: restart_file,out_file
1975  type(hexc_t),intent(in) :: hexc
1976  type(hexc_interp_t),intent(in) :: hexc_i
1977 !arrays
1978  complex(dpc),allocatable :: aa_file(:),bb_file(:),cc_file(:)
1979  complex(dpc),allocatable :: aa(:),bb(:),cc(:)
1980  complex(dpc),allocatable :: phi_np1(:),phi_n(:),phi_nm1(:)
1981  complex(dpc),allocatable :: phit_np1(:),phit_n(:),phit_nm1(:)
1982  complex(dpc),allocatable :: cbuff(:), phi_n_file(:),phi_np1_file(:)
1983  complex(dpc),allocatable :: ket0(:)
1984  complex(dpc),allocatable :: hphi_n(:), hphit_n(:)
1985  complex(dpc),allocatable :: all_omegas(:),green_temp(:,:)
1986  logical :: check(2)
1987 
1988 !************************************************************************
1989 
1990  MSG_WARNING("Haydock with Bilanczos is still under development")
1991 
1992  if(BSp%use_interp) then
1993    MSG_ERROR("Bilanczos is not yet implemented with interpolation")
1994  end if
1995 
1996  nproc  = xmpi_comm_size(comm)
1997  my_rank= xmpi_comm_rank(comm)
1998  nsppol = Hdr_bse%nsppol
1999 
2000  my_nt = my_t2-my_t1+1
2001  ABI_CHECK(my_nt>0,"One of the processors has zero columns")
2002 
2003  ! Multiplicative factor (k-point sampling and unit cell volume)
2004  ! TODO be careful with the spin here
2005  ! TODO four_pi comes from the coulomb term 1/|q| is already included in the
2006  ! oscillators hence the present approach wont work if a cutoff interaction is used.
2007  nfact = four_pi/(Cryst%ucvol*BSp%nkbz)
2008  if (nsppol==1) nfact=two*nfact
2009 
2010  write(msg,'(a,i0)')' Bi-Lanczos algorithm with MAX number of iterations: ',BSp%niter
2011  call wrtout(std_out,msg,"COLL")
2012  !
2013  ! Check for presence of the restart file.
2014  can_restart=.FALSE.
2015 
2016  if ( BS_files%in_haydock_basename /= BSE_NOFILE) then
2017    restart_file = strcat(BS_files%in_haydock_basename,tag_file)
2018    if (file_exists(restart_file) ) then
2019      can_restart=.TRUE.
2020      msg = strcat(" Restarting Haydock calculation from file: ",restart_file)
2021      call wrtout(std_out,msg,"COLL")
2022      call wrtout(ab_out,msg,"COLL")
2023      MSG_ERROR("Restart is not implemented")
2024    else
2025      can_restart=.FALSE.
2026      MSG_WARNING(strcat("Cannot find restart file: ",restart_file))
2027    end if
2028  end if
2029  !
2030  ! Open the file and writes basic dimensions and info.
2031  if (my_rank==master) then
2032    out_file = TRIM(BS_files%out_basename)//TRIM(tag_file)
2033    if (open_file(out_file,msg,newunit=out_unt,form="unformatted") /= 0) then
2034      MSG_ERROR(msg)
2035    end if
2036    ! write header TODO: standardize this part.
2037    write(out_unt)hsize,Bsp%use_coupling,BSE_HAYD_IMEPS,nkets,Bsp%broad
2038  end if
2039  !
2040  ! Select the terminator for the continued fraction.
2041  term_type=0 !; if (Bsp%hayd_term>0) term_type=2
2042  call wrtout(std_out,sjoin("Using terminator type: ",itoa(term_type)),"COLL")
2043  !
2044  ! Calculate green(w) for the different starting kets.
2045  green=czero
2046  do iq=1,nkets
2047    ABI_MALLOC(ket0,(hexc%hsize))
2048    ket0 = kets(:,iq)
2049    !
2050    niter_file=0
2051 
2052    if (can_restart) then
2053 !     call haydock_restart(BSp,restart_file,BSE_HAYD_IMEPS,iq,hsize,&
2054 !&      niter_file,aa_file,bb_file,phi_np1_file,phi_n_file,comm)
2055    end if
2056    !
2057    ABI_MALLOC(phi_nm1,(my_nt))
2058    ABI_MALLOC(phi_n,(my_nt))
2059    ABI_MALLOC(phi_np1,(my_nt))
2060    ABI_MALLOC(phit_nm1,(my_nt))
2061    ABI_MALLOC(phit_n,(my_nt))
2062    ABI_MALLOC(phit_np1,(my_nt))
2063    ABI_MALLOC(hphi_n,(hexc%hsize))
2064    ABI_MALLOC(hphit_n,(hexc%hsize))
2065    !
2066    ! TODO: Note the different convention used for the coefficients
2067    ! Should use the same convention in the Hermitian case.
2068    niter_max = niter_file + Bsp%niter
2069    ABI_MALLOC(aa,(niter_max))
2070    ABI_MALLOC(bb,(niter_max))
2071    ABI_MALLOC(cc,(niter_max))
2072    aa=czero; bb=czero; cc=czero
2073 
2074    if (niter_file==0) then ! Calculation from scratch.
2075      phi_nm1 = ket0(my_t1:my_t2)
2076      phit_nm1 = ket0(my_t1:my_t2)
2077      norm = DZNRM2(hexc%hsize,ket0,1)
2078      phi_nm1=phi_nm1/norm
2079      phit_nm1=phit_nm1/norm
2080 
2081      call hexc_matmul_elphon(hexc,phi_nm1,hphi_n,'N',ep_renorms)
2082      call hexc_matmul_elphon(hexc,phit_nm1,hphit_n,'C',ep_renorms)
2083 
2084      aa(1)=xdotc(my_nt,phit_nm1,1,hphi_n(my_t1:),1)
2085      call xmpi_sum(aa(1:1),comm,ierr)
2086 
2087      phi_n = hphi_n(my_t1:my_t2) - aa(1)*phi_nm1
2088      phit_n = hphit_n(my_t1:my_t2) - CONJG(aa(1))*phit_nm1
2089 
2090      bb(1)=xdotc(my_nt,phi_n,1,phi_n,1)
2091      call xmpi_sum(bb(1:1),comm,ierr)
2092      bb(1) = SQRT(bb(1))
2093 
2094      cc(1)=xdotc(my_nt,phit_n,1,phi_n,1)
2095      call xmpi_sum(cc(1:1),comm,ierr)
2096      cc(1) = cc(1)/bb(1)
2097 
2098      phi_n   = phi_n  /bb(1)
2099      phit_n  = phit_n /CONJG(cc(1))
2100      niter_done=1  ! TODO Be careful here
2101 
2102    else ! Use the previously calculates a and b.
2103      niter_done=niter_file
2104      MSG_ERROR("Restart not coded")
2105      !aa(1:niter_done) = aa_file
2106      !bb(1:niter_done) = bb_file
2107      !phi_np1=phi_np1_file(my_t1:my_t2)   ! Select the slice treated by this node.
2108      !phi_n  =phi_n_file  (my_t1:my_t2)
2109    end if
2110 
2111    if (can_restart) then
2112      ABI_FREE(aa_file)
2113      ABI_FREE(bb_file)
2114      ABI_FREE(cc_file)
2115      ABI_FREE(phi_np1_file)
2116      ABI_FREE(phi_n_file)
2117    end if
2118 
2119    ! This factor gives the correct results
2120    factor = -nfact*(DZNRM2(hexc%hsize,ket0,1)**2)
2121 
2122    ! Which quantity should be checked for convergence?
2123    check = (/.TRUE.,.TRUE./)
2124    if (ABS(Bsp%haydock_tol(2)-one)<tol6) check = (/.TRUE. ,.FALSE./)
2125    if (ABS(Bsp%haydock_tol(2)-two)<tol6) check = (/.FALSE.,.TRUE./)
2126    ! Create new frequencies "mirror" in negative range to add
2127    ! their contributions. Can be improved by computing only once
2128    ! zero frequency, but loosing clearness
2129    n_all_omegas = 2*BSp%nomega
2130 
2131    ABI_MALLOC(all_omegas,(n_all_omegas))
2132    ! Put all omegas with frequency > 0 in table
2133    all_omegas(BSp%nomega+1:n_all_omegas) = BSp%omega
2134    ! Put all omegas with frequency < 0
2135    ! Warning, the broadening must be kept positive
2136    all_omegas(1:BSp%nomega) = -DBLE(BSp%omega(BSp%nomega:1:-1)) + j_dpc*AIMAG(BSp%omega(BSp%nomega:1:-1))
2137 
2138    ABI_MALLOC(green_temp,(n_all_omegas,nkets))
2139 
2140 
2141 
2142    call haydock_bilanczos_optalgo(niter_done,niter_max,n_all_omegas,all_omegas,BSp%haydock_tol(1),check,hexc,hexc_i,&
2143 &    hsize,my_t1,my_t2,factor,term_type,ep_renorms,aa,bb,cc,ket0,ket0_hbar_norm,phi_nm1,phi_n,phi_np1,&
2144 &    phit_nm1,phit_n,phit_np1,green_temp(:,iq),inn,is_converged,comm)
2145 
2146 
2147    ! Computing result from two ranges of frequencies
2148    ! The real part is added, the imaginary part is substracted
2149    green(:,iq) = green_temp(BSp%nomega+1:n_all_omegas,iq)+CONJG(green_temp(BSp%nomega:1:-1,iq))
2150 
2151    ABI_FREE(all_omegas)
2152    ABI_FREE(green_temp)
2153 
2154 
2155    !
2156    ! Save the a"s and the b"s for possible restarting.
2157    ! 1) Info on the Q.
2158    ! 2) Number of iterations performed.
2159    ! 3) do iter=1,niter_performed
2160    !      aa(iter),bb(iter)
2161    !    end do
2162    ! 4) |n-1>
2163    !    |n>
2164    !    |n+1>
2165    !
2166    if (my_rank==master) then ! Open the file and writes basic dimensions and info.
2167      write(out_unt)Bsp%q(:,iq)
2168      write(out_unt)MIN(inn,niter_max)  ! NB: if the previous loop completed inn=niter_max+1
2169      do itt=1,MIN(inn,niter_max)        !     if we exited then inn is not incremented by one.
2170        write(out_unt)itt,aa(itt),bb(itt)
2171      end do
2172    end if
2173    !
2174    ! cbuff is used as workspace to gather |n-1>, |n> and |n+1>.
2175    ABI_MALLOC(cbuff,(hsize))
2176    cbuff=czero; cbuff(my_t1:my_t2) = phi_nm1
2177    call xmpi_sum_master(cbuff,master,comm,ierr)
2178    if (my_rank==master) write(out_unt) cbuff ! |n-1>
2179 
2180    cbuff=czero; cbuff(my_t1:my_t2) = phi_n
2181    call xmpi_sum_master(cbuff,master,comm,ierr)
2182    if (my_rank==master) write(out_unt) cbuff ! |n>
2183 
2184    cbuff=czero; cbuff(my_t1:my_t2) = phi_np1
2185    call xmpi_sum_master(cbuff,master,comm,ierr)
2186    if (my_rank==master) write(out_unt) cbuff ! |n+1>
2187 
2188    ABI_FREE(phi_nm1)
2189    ABI_FREE(phi_n)
2190    ABI_FREE(phi_np1)
2191    ABI_FREE(phit_nm1)
2192    ABI_FREE(phit_n)
2193    ABI_FREE(phit_np1)
2194    ABI_FREE(hphi_n)
2195    ABI_FREE(hphit_n)
2196    ABI_FREE(cbuff)
2197    ABI_FREE(aa)
2198    ABI_FREE(bb)
2199    ABI_FREE(cc)
2200    ABI_FREE(ket0)
2201  end do ! iq
2202 
2203  if (my_rank==master) close(out_unt)
2204 
2205  call xmpi_barrier(comm)
2206 
2207 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)

PARENTS

      m_haydock

CHILDREN

SOURCE

2256 subroutine haydock_bilanczos_optalgo(niter_done,niter_tot,nomega,omega,tol_iter,check,hexc,hexc_i,hsize,my_t1,my_t2,&
2257 &  factor,term_type,ep_renorms,aa,bb,cc,ket0,ket0_hbar_norm,phi_nm1,phi_n,phi_np1,phit_nm1,phit_n,phit_np1,&
2258 &  green,inn,is_converged,comm)
2259 
2260 
2261 !This section has been created automatically by the script Abilint (TD).
2262 !Do not modify the following lines by hand.
2263 #undef ABI_FUNC
2264 #define ABI_FUNC 'haydock_bilanczos_optalgo'
2265 !End of the abilint section
2266 
2267  implicit none
2268 
2269 !Arguments ------------------------------------
2270 !scalars
2271  integer,intent(in) :: niter_tot,niter_done,nomega,comm,hsize,my_t1,my_t2,term_type
2272  integer,intent(out) :: inn
2273  logical,intent(out) :: is_converged
2274  real(dp),intent(in) :: tol_iter,ket0_hbar_norm
2275  complex(dpc),intent(in) :: factor
2276  type(hexc_t),intent(in) :: hexc
2277  type(hexc_interp_t),intent(in) :: hexc_i
2278 !arrays
2279  complex(dpc),intent(inout) :: bb(niter_tot+1)
2280  complex(dpc),intent(out) :: green(nomega)
2281  complex(dpc),intent(in) :: omega(nomega)
2282  complex(dpc),intent(inout) :: aa(niter_tot),cc(niter_tot+1)
2283  complex(dpc),intent(in) :: ket0(my_t2-my_t1+1)
2284  complex(dpc),intent(in) :: ep_renorms(hsize)
2285  complex(dpc),intent(inout) :: phi_nm1(my_t2-my_t1+1)
2286  complex(dpc),intent(inout) :: phi_n  (my_t2-my_t1+1)
2287  complex(dpc),intent(inout) :: phi_np1(my_t2-my_t1+1)
2288  complex(dpc),intent(inout) :: phit_nm1(my_t2-my_t1+1)
2289  complex(dpc),intent(inout) :: phit_n  (my_t2-my_t1+1)
2290  complex(dpc),intent(inout) :: phit_np1(my_t2-my_t1+1)
2291  logical,intent(in) :: check(2)
2292 
2293 !Local variables ------------------------------
2294 !scalars
2295  integer :: my_nt,niter_min,nconv !,ierr
2296  character(len=500) :: msg
2297  logical :: keep_vectors=.TRUE.
2298 !arrays
2299  real(dp) :: abs_err(nomega,2) !,ww_err(nomega,2)
2300  complex(dpc),allocatable :: oldg(:),newg(:)
2301  complex(dpc),allocatable :: hphi_np1(:),hphit_np1(:),save_phi(:,:),save_phit(:,:)
2302  complex(dpc),allocatable :: g00(:)
2303  logical :: test(2)
2304  integer :: ierr
2305 
2306 !************************************************************************
2307 
2308  ABI_UNUSED(ket0_hbar_norm)
2309  ABI_UNUSED(hexc_i%hsize_dense)
2310 
2311  my_nt = my_t2-my_t1+1
2312 
2313  ABI_MALLOC(oldg,(nomega))
2314  ABI_MALLOC(newg,(nomega))
2315  ABI_MALLOC(g00,(nomega))
2316  oldg=czero; newg=czero; g00=czero
2317  nconv=0
2318 
2319  keep_vectors = (keep_vectors.and.xmpi_comm_size(comm)==1)
2320  if (keep_vectors) then
2321    ABI_MALLOC(save_phi,(my_t2-my_t1+1,niter_tot))
2322    ABI_STAT_MALLOC(save_phit,(my_t2-my_t1+1,niter_tot),ierr)
2323    ABI_CHECK(ierr==0,"out of memory in save_phi")
2324    save_phi=czero
2325    save_phit=czero
2326  end if
2327 
2328  ABI_STAT_MALLOC(hphi_np1,(hexc%hsize),ierr)
2329  ABI_CHECK(ierr==0,"out-of-memory hphi_np1")
2330  ABI_STAT_MALLOC(hphit_np1,(hexc%hsize),ierr)
2331  ABI_CHECK(ierr==0,"out-of-memory hphit_np1")
2332 
2333  do inn=niter_done+1,niter_tot
2334 
2335    !|n+1> = H |n> using all eh components.
2336    call hexc_matmul_elphon(hexc, phi_n, hphi_np1, 'N', ep_renorms)
2337    call hexc_matmul_elphon(hexc, phit_n, hphit_np1, 'C', ep_renorms)
2338 
2339    ! a(n) = < phit_n | H  | phi_n >
2340    aa(inn)=xdotc(my_nt,phit_n,1,hphi_np1(my_t1:),1)
2341    call xmpi_sum(aa(inn),comm,ierr)
2342 
2343    ! |n+1> = |n+1> - a(n)|Vn> - c(n)|n-1>
2344    phi_np1 = hphi_np1(my_t1:my_t2) - aa(inn)*phi_n - cc(inn-1)*phi_nm1
2345    phit_np1 = hphit_np1(my_t1:my_t2) - CONJG(aa(inn))*phit_n - CONJG(bb(inn-1))*phit_nm1
2346 
2347    bb(inn) = xdotc(my_nt,phi_np1,1,phi_np1,1)
2348    call xmpi_sum(bb(inn),comm,ierr)
2349    bb(inn) = SQRT(bb(inn))
2350 
2351    cc(inn) = xdotc(my_nt,phit_np1,1,phi_np1,1)
2352    call xmpi_sum(cc(inn),comm,ierr)
2353    cc(inn) = cc(inn)/bb(inn)
2354 
2355    phi_np1 = phi_np1 / bb(inn)
2356    phit_np1 = phit_np1 / CONJG(cc(inn))
2357 
2358    !
2359    ! |n-1> = |n>
2360    ! |n>   = |n+1>
2361    phi_nm1 = phi_n
2362    phi_n   = phi_np1
2363    phit_nm1 = phit_n
2364    phit_n = phit_np1
2365 
2366    if (keep_vectors) then
2367      save_phi(:,inn) = phi_n
2368      save_phit(:,inn) = phit_n
2369    end if
2370    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))
2371    call wrtout(std_out,msg,"COLL")
2372 
2373    call continued_fract_general(inn,term_type,aa,bb,cc,nomega,omega,g00)
2374    newg = factor*g00
2375    !
2376    ! Avoid spurious convergence.
2377    niter_min=4; if (niter_done>1) niter_min=niter_done+1
2378    if (inn>niter_min) then
2379      test=.TRUE.
2380      abs_err(:,1) = ABS(DBLE (newg-oldg))
2381      abs_err(:,2) = ABS(AIMAG(newg-oldg))
2382      !
2383      if (tol_iter>zero) then
2384        ! Test on the L1 norm.
2385        if (check(1)) test(1) = SUM(abs_err(:,1)) < tol_iter*SUM(ABS(DBLE (newg)))
2386        if (check(2)) test(2) = SUM(abs_err(:,2)) < tol_iter*SUM(ABS(AIMAG(newg)))
2387      else
2388        ! Stringent test for each point.
2389        if (check(1)) test(1) = ALL( abs_err(:,1) < -tol_iter*ABS(DBLE (newg)))
2390        if (check(2)) test(2) = ALL( abs_err(:,2) < -tol_iter*ABS(AIMAG(newg)))
2391      end if
2392      !
2393      if (ALL(test)) then
2394        nconv = nconv+1
2395      else
2396        nconv = 0
2397      end if
2398      if (nconv==2) then
2399        write(msg,'(a,es10.2,a,i0,a)')&
2400 &        " >>> Haydock algorithm converged twice within haydock_tol= ",tol_iter," after ",inn," iterations."
2401        call wrtout(std_out,msg,'COLL')
2402        call wrtout(ab_out,msg,'COLL')
2403        EXIT
2404      end if
2405    end if
2406    !
2407    oldg = newg
2408  end do ! inn
2409 
2410  green = newg
2411  if (nconv/=2) then
2412    write(msg,'(a,es10.2,a,i0,a)')&
2413 &    " WARNING: Haydock algorithm did not converge within ",tol_iter," after ",niter_tot," iterations."
2414    call wrtout(std_out,msg,'COLL')
2415    call wrtout(ab_out,msg,'COLL')
2416  end if
2417 
2418  is_converged = (nconv==2)
2419 
2420  ABI_FREE(oldg)
2421  ABI_FREE(newg)
2422  ABI_FREE(g00)
2423  ABI_FREE(hphi_np1)
2424  ABI_FREE(hphit_np1)
2425 
2426  if (keep_vectors) then
2427    ABI_FREE(save_phi)
2428    ABI_FREE(save_phit)
2429  end if
2430 
2431  !! if (keep_vectors) then
2432  !!   tdim = MIN(inn,niter_tot)
2433  !!   ABI_MALLOC(ovlp,(tdim,tdim))
2434 
2435  !!   ABI_MALLOC(phi_test,(hsize))
2436  !!   ABI_MALLOC(phi_test2,(hsize))
2437 
2438  !!   max_err=smallest_real; mean_err=zero; mean_err2=zero; row_max=-1
2439  !!   do ii=1,tdim
2440  !!     parity = (-1)**(ii+1)
2441  !!     phi_test  = save_phi(:,ii)
2442  !!     call hexc_matmul_full(hexc, hexc_i, phi_test, phi_test2, parity)
2443  !!     !phi_test2 = MATMUL(hreso,phi_test) + parity * MATMUL(hcoup,CONJG(phi_test))
2444  !!     ovlp(ii,ii) = DOT_PRODUCT(phi_test,phi_test2) + DOT_PRODUCT(phi_test2,phi_test)
2445  !!     err = ABS(ovlp(ii,ii)-cone)
2446  !!     mean_err  = mean_err + err
2447  !!     mean_err2 = mean_err2 + err**2
2448  !!     if (err > max_err) then
2449  !!       max_err = err
2450  !!       row_max = ii
2451  !!     end if
2452  !!   end do
2453  !!   mean_err = mean_err/tdim
2454  !!   std_dev = mean_err2/tdim -mean_err**2
2455  !!   write(std_out,'(a,i0,1x,3es14.6)')&
2456  !!&    " Error in normalization (ii, max_err,mean,std_dev): ",row_max,max_err,mean_err,std_dev
2457 
2458  !!   ABI_FREE(phi_test)
2459  !!   ABI_FREE(phi_test2)
2460  !!
2461  !!   ABI_MALLOC(alpha,(hsize,tdim))
2462 
2463  !!   ! Less efficient but for sake of simplicity with hexc_matmul
2464  !!   ! TODO possibility to call hreso * phi, and hcoup * phi separately
2465  !!   do ii=1,tdim
2466  !!     parity = (-1)**(ii+1)
2467  !!     call hexc_matmul_full(hexc, hexc_i, save_phi(:,ii), alpha(:,ii), parity)
2468  !!   end do
2469 
2470  !!   !alpha = MATMUL(hreso,save_phi(:,1:tdim))
2471  !!   !
2472  !!   !do ii=1,tdim
2473  !!   !  parity = (-1)**(ii+1)
2474  !!   !  alpha(:,ii) =  alpha(:,ii) + parity*MATMUL(hcoup,CONJG(save_phi(:,ii)))
2475  !!   !end do
2476 
2477  !!   ovlp = MATMUL(TRANSPOSE(CONJG(save_phi(:,1:tdim))),alpha)
2478 
2479  !!   ABI_MALLOC(beta,(hsize,tdim))
2480  !!   do ii=1,tdim
2481  !!     parity = (-1)**(ii+1)
2482  !!     beta(:,ii)  =  parity*save_phi(:,ii)
2483  !!     alpha(:,ii) = -parity*alpha(:,ii)
2484  !!   end do
2485 
2486  !!   ovlp = ovlp - MATMUL(TRANSPOSE(CONJG(beta)),alpha)
2487 
2488  !!   max_err=smallest_real; row_max=-1; col_max=-1
2489  !!   mean_err=zero; mean_err2=zero
2490  !!   do jj=1,tdim
2491  !!     do ii=1,jj
2492  !!       err = ABS(ovlp(ii,jj))
2493  !!       if (ii==jj) err = ABS(err - one)
2494  !!       mean_err  = mean_err + err
2495  !!       mean_err2 = mean_err2 + err**2
2496  !!       if (err > max_err) then
2497  !!         max_err = err
2498  !!         row_max=ii
2499  !!         col_max=jj
2500  !!       end if
2501  !!     end do
2502  !!   end do
2503 
2504  !!   mean_err = mean_err/(tdim*(tdim+1)/2)
2505  !!   std_dev = mean_err2/(tdim*(tdim+1)/2) - mean_err**2
2506  !!   write(std_out,'(a,2(i0,1x),3es14.6)')&
2507  !!      " Error in Hbar-ortho (i,j), max_err, mean, std_dev ",row_max,col_max,max_err,mean_err,std_dev
2508  !!   !call print_arr(ovlp,max_r=185,max_c=10,unit=std_out)
2509 
2510  !!   ABI_FREE(alpha)
2511  !!   ABI_FREE(beta)
2512  !!   ABI_FREE(ovlp)
2513  !!   ABI_FREE(save_phi)
2514  !! end if
2515 
2516 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)=

PARENTS

      m_haydock

CHILDREN

SOURCE

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

PARENTS

      m_haydock

CHILDREN

SOURCE

 928 subroutine haydock_herm_algo(niter_done,niter_max,nomega,omega,tol_iter,check,&
 929 & my_t1,my_t2,factor,term_type,aa,bb,phi_nm1,phi_n,&
 930 & green,inn,is_converged,&
 931 & hexc, hexc_i, comm)
 932 
 933 
 934 !This section has been created automatically by the script Abilint (TD).
 935 !Do not modify the following lines by hand.
 936 #undef ABI_FUNC
 937 #define ABI_FUNC 'haydock_herm_algo'
 938 !End of the abilint section
 939 
 940  implicit none
 941 
 942 !Arguments ------------------------------------
 943 !scalars
 944  integer,intent(in) :: niter_max,niter_done,nomega
 945  integer,intent(in) :: my_t1,my_t2,term_type
 946  integer,intent(in) :: comm
 947  integer,intent(out) :: inn
 948  logical,intent(out) :: is_converged
 949  real(dp),intent(in) :: tol_iter
 950  complex(dpc),intent(in) :: factor
 951  type(hexc_t),intent(in) :: hexc
 952  type(hexc_interp_t),intent(in) :: hexc_i
 953 !arrays
 954  real(dp),intent(inout) :: bb(niter_max)
 955  complex(dpc),intent(out) :: green(nomega)
 956  complex(dpc),intent(in) :: omega(nomega)
 957  complex(dpc),intent(inout) :: aa(niter_max)
 958  complex(dpc),intent(inout) :: phi_nm1(my_t2-my_t1+1)
 959  complex(dpc),intent(inout) :: phi_n  (my_t2-my_t1+1)
 960  logical,intent(in) :: check(2)
 961 
 962 !Local variables ------------------------------
 963 !scalars
 964  integer :: ierr,my_nt,niter_min,nconv
 965  character(len=500) :: msg
 966  logical,parameter :: force_real=.TRUE.
 967 !arrays
 968  real(dp) :: abs_err(nomega,2) !,rel_err(nomega,2)
 969  complex(dpc),allocatable :: oldg(:),newg(:)
 970  complex(dpc),allocatable :: phi_np1(:),hphi_n(:),cfact(:)
 971  logical :: test(2)
 972 
 973 !************************************************************************
 974 
 975  ! The sequences starts with |1> normalized to 1 and b_0 =0, therefore:
 976  !  a_1 = <1|H|1>
 977  !  b_1 = || H|1> - a_1|1> ||
 978  !  |2> = [H|1> - a_1|1>]/b_1
 979  !
 980  ! For n>1 we have
 981  !  1) a_n = <n|H|n>
 982  !  2) b_n = || H|n> - a_n|n> -b_{n-1}|n-1> ||
 983  !  3) |n+1> = [H|n> -a_n|n> -b_{n-1}|n-1>]/b_n
 984  !
 985  my_nt = my_t2-my_t1+1
 986 
 987  ABI_STAT_MALLOC(hphi_n,(hexc%hsize), ierr)
 988  ABI_CHECK(ierr==0, "out-of-memory hphi_n")
 989 
 990  ABI_MALLOC(phi_np1,(my_nt))
 991 
 992  ABI_MALLOC(oldg,(nomega))
 993  oldg=czero
 994  ABI_MALLOC(newg,(nomega))
 995  newg=czero
 996  ABI_MALLOC(cfact,(nomega))
 997  cfact=czero
 998 
 999  nconv=0
1000  do inn=niter_done+1,niter_max
1001 
1002    !YG2014
1003    call hexc_matmul_tda(hexc,hexc_i,phi_n,hphi_n)
1004 
1005    aa(inn) = xdotc(my_nt,phi_n,1,hphi_n(my_t1:),1)
1006    call xmpi_sum(aa(inn:inn),comm,ierr)
1007    if (force_real) aa(inn) = DBLE(aa(inn)) ! Matrix is Hermitian.
1008 
1009    ! |n+1> = H|n> - A(n)|n> - B(n-1)|n-1>
1010    phi_np1 = hphi_n(my_t1:my_t2) - aa(inn)*phi_n - bb(inn-1)*phi_nm1
1011 
1012    bb(inn) = xdotc(my_nt,phi_np1,1,phi_np1,1)
1013    call xmpi_sum(bb(inn),comm,ierr)
1014    bb(inn) = SQRT(bb(inn))
1015 
1016    phi_np1 = phi_np1/bb(inn)
1017 
1018    phi_nm1 = phi_n
1019    phi_n   = phi_np1
1020 
1021    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))
1022    call wrtout(std_out,msg,"COLL")
1023 
1024    call continued_fract(inn,term_type,aa,bb,nomega,omega,cfact)
1025 
1026    newg= factor*cfact
1027    !
1028    ! Avoid spurious convergence.
1029    niter_min=4; if (niter_done>1) niter_min=niter_done+1
1030    if (inn>niter_min) then
1031      test=.TRUE.
1032      abs_err(:,1) = ABS(DBLE (newg-oldg))
1033      abs_err(:,2) = ABS(AIMAG(newg-oldg))
1034      !
1035      if (tol_iter>zero) then
1036        ! Test on the L1 norm.
1037        if (check(1)) test(1) = SUM(abs_err(:,1)) < tol_iter*SUM(ABS(DBLE (newg)))
1038        if (check(2)) test(2) = SUM(abs_err(:,2)) < tol_iter*SUM(ABS(AIMAG(newg)))
1039 
1040      else
1041        ! Stringent test for each point.
1042        if (check(1)) test(1) = ALL( abs_err(:,1) < -tol_iter*ABS(DBLE (newg)))
1043        if (check(2)) test(2) = ALL( abs_err(:,2) < -tol_iter*ABS(AIMAG(newg)))
1044      end if
1045      !
1046      if (ALL(test)) then
1047        nconv = nconv+1
1048      else
1049        nconv = 0
1050      end if
1051      if (nconv==2) then
1052        write(msg,'(a,es10.2,a,i0,a)')&
1053 &        " >>> Haydock algorithm converged twice within haydock_tol= ",tol_iter," after ",inn," iterations."
1054        call wrtout(std_out,msg,'COLL')
1055        call wrtout(ab_out,msg,'COLL')
1056        EXIT
1057      end if
1058    end if
1059 
1060    oldg = newg
1061  end do ! inn
1062 
1063  green = newg
1064  if (nconv/=2) then
1065    write(msg,'(a,es10.2,a,i0,a)')&
1066 &    " WARNING: Haydock algorithm did not converge within ",tol_iter," after ",niter_max," iterations."
1067    call wrtout(std_out,msg,'COLL')
1068    call wrtout(ab_out,msg,'COLL')
1069  end if
1070 
1071  is_converged = (nconv==2)
1072 
1073  ABI_FREE(oldg)
1074  ABI_FREE(newg)
1075  ABI_FREE(cfact)
1076  ABI_FREE(hphi_n)
1077  ABI_FREE(phi_np1)
1078 
1079 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)

PARENTS

      m_haydock

CHILDREN

SOURCE

1240 subroutine haydock_mdf_to_tensor(BSp,Cryst,eps,tensor_cart,tensor_red,ierr)
1241 
1242 
1243 !This section has been created automatically by the script Abilint (TD).
1244 !Do not modify the following lines by hand.
1245 #undef ABI_FUNC
1246 #define ABI_FUNC 'haydock_mdf_to_tensor'
1247 !End of the abilint section
1248 
1249  implicit none
1250 
1251 !Arguments ------------------------------------
1252 !scalars
1253  integer,intent(out) :: ierr
1254  type(excparam),intent(in) :: BSp
1255  type(crystal_t),intent(in) :: Cryst
1256 !arrays
1257  complex(dpc),intent(in) :: eps(BSp%nomega,BSp%nq)
1258  complex(dpc),intent(out) :: tensor_cart(BSp%nomega,6), tensor_red(BSp%nomega,6)
1259 
1260 !Local variables ------------------------------
1261 !scalars
1262  integer :: iq,info
1263  real(dp) :: normqcart, normqred
1264 !arrays
1265  integer,allocatable :: ipiv(:)
1266  real(dp) :: qcart(3), qtmet(3)
1267  real(dp) :: qred2cart(3,3),qcart2red(3,3)
1268  complex(dpc) :: qqcart(BSp%nq,6), qqred(BSp%nq,6)
1269  complex(dpc) :: b(6,BSP%nomega)
1270 
1271 !************************************************************************
1272 
1273  ! Error flag
1274  ierr = 0
1275 
1276  if(BSp%nq /= 6) then
1277     ierr = -1
1278     return
1279  end if
1280 
1281  ! Transformation matrices from reduced coordinates to cartesian coordinates
1282  qred2cart = two_pi*Cryst%gprimd
1283  qcart2red = qred2cart
1284  call matrginv(qcart2red,3,3)
1285  do iq = 1, 6
1286 
1287    ! Computing cartesian q-vector
1288    qcart = MATMUL(qred2cart, BSp%q(:,iq))
1289 
1290    ! Computing product 'metric - qred' to form quadratic form
1291    qtmet = (two_pi**2)*MATMUL(Cryst%gmet, BSp%q(:,iq))
1292 
1293    ! squared norms
1294    normqcart = qcart(1)**2+qcart(2)**2+qcart(3)**2
1295    normqred = (normv(BSp%q(:,iq),Cryst%gmet,"G"))**2
1296 
1297    ! Compute line 'iq' for matrix in cartesian coord
1298    qqcart(iq,1) = (qcart(1))**2
1299    qqcart(iq,2) = (qcart(2))**2
1300    qqcart(iq,3) = (qcart(3))**2
1301    qqcart(iq,4) = 2*(qcart(1)*qcart(2))
1302    qqcart(iq,5) = 2*(qcart(1)*qcart(3))
1303    qqcart(iq,6) = 2*(qcart(2)*qcart(3))
1304 
1305    ! Compute line 'iq' for matrix in reduced coord
1306    qqred(iq,1) = (qtmet(1))**2
1307    qqred(iq,2) = (qtmet(2))**2
1308    qqred(iq,3) = (qtmet(3))**2
1309    qqred(iq,4) = 2*(qtmet(1)*qtmet(2))
1310    qqred(iq,5) = 2*(qtmet(1)*qtmet(3))
1311    qqred(iq,6) = 2*(qtmet(2)*qtmet(3))
1312 
1313    ! Renormalize line
1314    qqcart(iq,:) = qqcart(iq,:)/normqcart
1315    qqred(iq,:) = qqred(iq,:)/normqred
1316  end do
1317 
1318  ABI_MALLOC(ipiv,(6))
1319 
1320  ! Solving linear system
1321  b = TRANSPOSE(eps)
1322  call ZGESV(6,BSp%nomega,qqcart,6,ipiv,b,6,info)
1323  tensor_cart = TRANSPOSE(b)
1324 
1325  if(info /= 0) then
1326    ! Skipping the rest of the routine
1327    ierr = info
1328    ABI_FREE(ipiv)
1329    return
1330  end if
1331 
1332  b = TRANSPOSE(eps)
1333  call ZGESV(6,BSp%nomega,qqred,6,ipiv,b,6,info)
1334  tensor_red = TRANSPOSE(b)
1335 
1336  if(info /= 0) ierr = info
1337 
1338  ABI_FREE(ipiv)
1339 
1340 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.

PARENTS

      m_haydock

CHILDREN

SOURCE

1374 subroutine haydock_psherm(BSp,BS_files,Cryst,Hdr_bse,hexc,hexc_i,hsize,my_t1,my_t2,nkets,kets,green,comm)
1375 
1376 
1377 !This section has been created automatically by the script Abilint (TD).
1378 !Do not modify the following lines by hand.
1379 #undef ABI_FUNC
1380 #define ABI_FUNC 'haydock_psherm'
1381 !End of the abilint section
1382 
1383  implicit none
1384 
1385 !Arguments ------------------------------------
1386 !scalars
1387  integer,intent(in) :: hsize,my_t1,my_t2,nkets,comm
1388  type(crystal_t),intent(in) :: Cryst
1389  type(excparam),intent(in) :: BSp
1390  type(excfiles),intent(in) :: BS_files
1391  type(Hdr_type),intent(in) :: Hdr_bse
1392 !arrays
1393  complex(dp),intent(out) :: green(BSp%nomega,BSp%nq)
1394  complex(dpc),intent(in) :: kets(hsize,nkets)
1395 
1396 !Local variables ------------------------------
1397 !scalars
1398  integer,parameter :: master=0
1399  integer :: inn,itt,out_unt,nproc,my_rank,ierr
1400  integer :: niter_file,niter_max,niter_done,nsppol,iq,my_nt,term_type
1401  real(dp) :: ket0_hbar_norm,nfact
1402  logical :: can_restart,is_converged
1403  complex(dpc) :: factor
1404  character(len=fnlen),parameter :: tag_file="_HAYDC_SAVE"
1405  character(len=500) :: msg
1406  character(len=fnlen) :: restart_file,out_file
1407  type(hexc_t),intent(in) :: hexc
1408  type(hexc_interp_t),intent(in) :: hexc_i
1409 !arrays
1410  real(dp),allocatable :: bb_file(:)
1411  real(dp),allocatable :: bb(:)
1412  complex(dpc),allocatable :: aa(:),cc(:),phi_np1(:),phi_n(:),phi_nm1(:),cbuff(:)
1413  complex(dpc),allocatable :: aa_file(:),phi_n_file(:),phi_np1_file(:),cc_file(:)
1414  complex(dpc),allocatable :: ket0(:)
1415  logical :: check(2)
1416 
1417 !************************************************************************
1418 
1419  MSG_WARNING("Haydock + coupling is still under development")
1420 
1421  if(BSp%use_interp) then
1422    MSG_ERROR("Coupling is not yet implemented with interpolation")
1423  end if
1424 
1425  nproc  = xmpi_comm_size(comm)
1426  my_rank= xmpi_comm_rank(comm)
1427  nsppol = Hdr_bse%nsppol
1428 
1429  my_nt = my_t2-my_t1+1
1430  ABI_CHECK(my_nt>0,"One of the processors has zero columns")
1431 
1432  ! Multiplicative factor (k-point sampling and unit cell volume)
1433  ! TODO be careful with the spin here
1434  ! TODO four_pi comes from the coulomb term 1/|q| is already included in the
1435  ! oscillators hence the present approach wont work if a cutoff interaction is used.
1436  nfact = four_pi/(Cryst%ucvol*BSp%nkbz)
1437  if (nsppol==1) nfact=two*nfact
1438 
1439  write(msg,'(a,i0)')' Haydock algorithm with MAX number of iterations: ',BSp%niter
1440  call wrtout(std_out,msg,"COLL")
1441  !
1442  ! Check for presence of the restart file.
1443  can_restart=.FALSE.
1444 
1445  if ( BS_files%in_haydock_basename /= BSE_NOFILE) then
1446    restart_file = strcat(BS_files%in_haydock_basename,tag_file)
1447    if (file_exists(restart_file) ) then
1448      can_restart=.TRUE.
1449      msg = strcat(" Restarting Haydock calculation from file: ",restart_file)
1450      call wrtout(std_out,msg,"COLL")
1451      call wrtout(ab_out,msg,"COLL")
1452      MSG_ERROR("Restart is not tested")
1453    else
1454      can_restart=.FALSE.
1455      MSG_WARNING(strcat("Cannot find restart file: ",restart_file))
1456    end if
1457  end if
1458  !
1459  ! Open the file and writes basic dimensions and info.
1460  if (my_rank==master) then
1461    out_file = TRIM(BS_files%out_basename)//TRIM(tag_file)
1462    if (open_file(out_file,msg,newunit=out_unt,form="unformatted") /= 0) then
1463      MSG_ERROR(msg)
1464    end if
1465    ! write header TODO: standardize this part.
1466    write(out_unt)hsize,Bsp%use_coupling,BSE_HAYD_IMEPS,nkets,Bsp%broad
1467  end if
1468  !
1469  ! Select the terminator for the continued fraction.
1470  term_type=0 !; if (Bsp%hayd_term>0) term_type=2
1471  call wrtout(std_out,sjoin("Using terminator type: ",itoa(term_type)),"COLL")
1472  !
1473  ! Calculate green(w) for the different starting kets.
1474  green=czero
1475  do iq=1,nkets
1476    ABI_MALLOC(ket0,(my_nt))
1477    ket0 = kets(my_t1:my_t2,iq)
1478    !
1479    niter_file=0
1480 
1481    if (can_restart) then
1482      call haydock_restart(BSp,restart_file,BSE_HAYD_IMEPS,iq,hsize,&
1483 &      niter_file,aa_file,bb_file,phi_np1_file,phi_n_file,comm)
1484    end if
1485    !
1486    ABI_MALLOC(phi_nm1,(my_nt))
1487    ABI_MALLOC(phi_n,(my_nt))
1488    ABI_MALLOC(phi_np1,(my_nt))
1489    !
1490    ! TODO: Note the different convention used for the coefficients
1491    ! Should use the same convention in the Hermitian case.
1492    niter_max = niter_file + Bsp%niter
1493    ABI_MALLOC(aa,(niter_max))
1494    ABI_MALLOC(bb,(niter_max+1))
1495    ABI_MALLOC(cc,(niter_max+1))
1496    aa=czero; bb=czero; cc=czero
1497 
1498    if (niter_file==0) then ! Calculation from scratch.
1499      phi_n   = ket0
1500      call hexc_matmul_full(hexc, hexc_i, phi_n, phi_np1, -1)
1501      !phi_np1 = MATMUL(hreso,ket0) - MATMUL(hcoup,CONJG(ket0))
1502      ket0_hbar_norm = SQRT(two*DBLE(DOT_PRODUCT(phi_n,phi_np1)))
1503      phi_n   = phi_n  /ket0_hbar_norm
1504      phi_np1 = phi_np1/ket0_hbar_norm
1505      !ket0    = ket0/ket0_hbar_norm
1506      cc(1)=zero ! <P|F|P>
1507      !cc(1) =  DOT_PRODUCT(ket0,phi_np1)
1508      write(std_out,*)" cc(1), ket0_hbar_norm =",cc(1),ket0_hbar_norm
1509 
1510      phi_nm1 = czero
1511      niter_done=0  ! TODO Be careful here
1512 
1513    else ! Use the previously calculates a and b.
1514      niter_done=niter_file
1515      MSG_ERROR("Restart not coded")
1516      !aa(1:niter_done) = aa_file
1517      !bb(1:niter_done) = bb_file
1518      !phi_np1=phi_np1_file(my_t1:my_t2)   ! Select the slice treated by this node.
1519      !phi_n  =phi_n_file  (my_t1:my_t2)
1520    end if
1521 
1522    if (can_restart) then
1523      ABI_FREE(aa_file)
1524      ABI_FREE(bb_file)
1525      ABI_FREE(cc_file)
1526      ABI_FREE(phi_np1_file)
1527      ABI_FREE(phi_n_file)
1528    end if
1529 
1530    ! This factor gives the correct results
1531    factor = -nfact*ket0_hbar_norm / SQRT(two)
1532 
1533    ! Which quantity should be checked for convergence?
1534    check = (/.TRUE.,.TRUE./)
1535    if (ABS(Bsp%haydock_tol(2)-one)<tol6) check = (/.TRUE. ,.FALSE./)
1536    if (ABS(Bsp%haydock_tol(2)-two)<tol6) check = (/.FALSE.,.TRUE./)
1537 
1538    call haydock_psherm_optalgo(niter_done,niter_max,BSp%nomega,BSp%omega,BSp%haydock_tol(1),check,hexc,hexc_i,&
1539 &    hsize,my_t1,my_t2,factor,term_type,aa,bb,cc,ket0,ket0_hbar_norm,phi_nm1,phi_n,phi_np1,&
1540 &    green(:,iq),inn,is_converged,comm)
1541    !
1542    ! Save the a"s and the b"s for possible restarting.
1543    ! 1) Info on the Q.
1544    ! 2) Number of iterations performed.
1545    ! 3) do iter=1,niter_performed
1546    !      aa(iter),bb(iter)
1547    !    end do
1548    ! 4) |n-1>
1549    !    |n>
1550    !    |n+1>
1551    !
1552    if (my_rank==master) then ! Open the file and writes basic dimensions and info.
1553      write(out_unt)Bsp%q(:,iq)
1554      write(out_unt)MIN(inn,niter_max)  ! NB: if the previous loop completed inn=niter_max+1
1555      do itt=1,MIN(inn,niter_max)        !     if we exited then inn is not incremented by one.
1556        write(out_unt)itt,aa(itt),bb(itt)
1557      end do
1558    end if
1559    !
1560    ! cbuff is used as workspace to gather |n-1>, |n> and |n+1>.
1561    ABI_MALLOC(cbuff,(hsize))
1562    cbuff=czero; cbuff(my_t1:my_t2) = phi_nm1
1563    call xmpi_sum_master(cbuff,master,comm,ierr)
1564    if (my_rank==master) write(out_unt) cbuff ! |n-1>
1565 
1566    cbuff=czero; cbuff(my_t1:my_t2) = phi_n
1567    call xmpi_sum_master(cbuff,master,comm,ierr)
1568    if (my_rank==master) write(out_unt) cbuff ! |n>
1569 
1570    cbuff=czero; cbuff(my_t1:my_t2) = phi_np1
1571    call xmpi_sum_master(cbuff,master,comm,ierr)
1572    if (my_rank==master) write(out_unt) cbuff ! |n+1>
1573 
1574    ABI_FREE(phi_nm1)
1575    ABI_FREE(phi_n)
1576    ABI_FREE(phi_np1)
1577    ABI_FREE(cbuff)
1578    ABI_FREE(aa)
1579    ABI_FREE(bb)
1580    ABI_FREE(cc)
1581    ABI_FREE(ket0)
1582  end do ! iq
1583 
1584  if (my_rank==master) close(out_unt)
1585 
1586  call xmpi_barrier(comm)
1587 
1588 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)

PARENTS

      m_haydock

CHILDREN

SOURCE

1637 subroutine haydock_psherm_optalgo(niter_done,niter_tot,nomega,omega,tol_iter,check,hexc,hexc_i,hsize,my_t1,my_t2,&
1638 &  factor,term_type,aa,bb,cc,ket0,ket0_hbar_norm,phi_nm1,phi_n,phi_np1,green,inn,is_converged,comm)
1639 
1640 
1641 !This section has been created automatically by the script Abilint (TD).
1642 !Do not modify the following lines by hand.
1643 #undef ABI_FUNC
1644 #define ABI_FUNC 'haydock_psherm_optalgo'
1645 !End of the abilint section
1646 
1647  implicit none
1648 
1649 !Arguments ------------------------------------
1650 !scalars
1651  integer,intent(in) :: niter_tot,niter_done,nomega,comm,hsize,my_t1,my_t2,term_type
1652  integer,intent(out) :: inn
1653  logical,intent(out) :: is_converged
1654  real(dp),intent(in) :: tol_iter,ket0_hbar_norm
1655  complex(dpc),intent(in) :: factor
1656  type(hexc_t),intent(in) :: hexc
1657  type(hexc_interp_t),intent(in) :: hexc_i
1658 !arrays
1659  real(dp),intent(inout) :: bb(niter_tot+1)
1660  complex(dpc),intent(out) :: green(nomega)
1661  complex(dpc),intent(in) :: omega(nomega)
1662  complex(dpc),intent(inout) :: aa(niter_tot),cc(niter_tot+1)
1663  complex(dpc),intent(in) :: ket0(my_t2-my_t1+1)
1664  complex(dpc),intent(inout) :: phi_nm1(my_t2-my_t1+1)
1665  complex(dpc),intent(inout) :: phi_n  (my_t2-my_t1+1)
1666  complex(dpc),intent(inout) :: phi_np1(my_t2-my_t1+1)
1667  logical,intent(in) :: check(2)
1668 
1669 !Local variables ------------------------------
1670 !scalars
1671  integer :: my_nt,niter_min,nconv,parity,ii,jj,tdim,ierr
1672  integer :: row_max,col_max,nlev
1673  character(len=500) :: msg
1674  real(dp) :: max_err,mean_err,mean_err2,std_dev,err
1675  logical :: keep_vectors=.TRUE.
1676 !arrays
1677  real(dp) :: abs_err(nomega,2) !,ww_err(nomega,2)
1678  complex(dpc) :: gn0(nomega,niter_tot)
1679  complex(dpc),allocatable :: oldg(:),newg(:)
1680  complex(dpc),allocatable :: hphi_n(:),save_phi(:,:)
1681  complex(dpc),allocatable ::  alpha(:,:),beta(:,:),ovlp(:,:)
1682  complex(dpc),allocatable :: phi_test(:),phi_test2(:),g00(:)
1683  logical :: test(2)
1684 
1685 !************************************************************************
1686 
1687  ABI_UNUSED(ket0_hbar_norm)
1688 
1689  my_nt = my_t2-my_t1+1
1690 
1691  ABI_MALLOC(oldg,(nomega))
1692  ABI_MALLOC(newg,(nomega))
1693  ABI_MALLOC(g00,(nomega))
1694  oldg=czero; newg=czero; g00=czero
1695  nconv=0
1696 
1697  keep_vectors = (keep_vectors.and.xmpi_comm_size(comm)==1)
1698  if (keep_vectors) then
1699    ABI_STAT_MALLOC(save_phi,(my_t2-my_t1+1,niter_tot), ierr)
1700    ABI_CHECK(ierr==0, "out of memory in save_phi")
1701    save_phi=czero
1702  end if
1703 
1704  ABI_MALLOC(hphi_n,(hsize))
1705 
1706  do inn=niter_done+1,niter_tot
1707    !
1708    ! a(n) = <Vn+1|F|Vn+1> = <Vn|HFH|Vn>) = 0 by symmetry.
1709    aa(inn)=zero
1710 
1711    ! |n+1> = |n+1> - a(n)|Vn> - a(n)|n-1>
1712    phi_np1 = phi_np1 - bb(inn)*phi_nm1
1713    !
1714    ! |n-1> = |n>
1715    ! |n>   = |n+1>
1716    phi_nm1 = phi_n
1717    phi_n   = phi_np1
1718    !
1719    !|n+1> = H |n> using all eh components.
1720    parity = (-1)**(inn+1)
1721    call hexc_matmul_full(hexc, hexc_i, phi_n, phi_np1, parity)
1722 
1723    !phi_np1 = MATMUL(hreso,phi_n) + parity * MATMUL(hcoup,CONJG(phi_n))
1724    !call xmpi_sum(hphi_np1,comm,ierr)
1725    !
1726    ! B(n+1)= <n|F|n+1>^(1/2) = <n|FH|n>^(1/2))= (2*Re(<n|V+1>))^(1/2)
1727    ! by symmetry, where the dot_product is done in the resonant eh sub-space.
1728    !
1729    bb(inn+1)=SQRT(two*DBLE(DOT_PRODUCT(phi_n,phi_np1)))
1730    !bb(inn+1)=two*DBLE(DOT_PRODUCT(phi_n,phi_np1))
1731    !call xmpi_sum(bb(inn+1),comm,ierr)
1732    !bb(inn+1)=SQRT(bb(inn+1)
1733    !
1734    !|n+1> =|n+1>/B(n+1)
1735    phi_n   = phi_n  /bb(inn+1)
1736    phi_np1 = phi_np1/bb(inn+1)
1737 
1738    if (keep_vectors) save_phi(:,inn) = phi_n
1739 
1740    parity = (-1)**(inn+1)
1741    !if (parity==-1) then
1742    !  cc(inn+1)=czero
1743    !else
1744      cc(inn+1)=DOT_PRODUCT(ket0,phi_n) + parity * DOT_PRODUCT(phi_n,ket0)
1745    !end if
1746    !call xmpi_sum(cc(inn+1),comm,ierr)
1747 
1748    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))
1749    call wrtout(std_out,msg,"COLL")
1750 
1751    call continued_fract(inn,term_type,aa,bb(2:),nomega,omega,g00)
1752    gn0(:,1) = g00
1753 
1754    if (.FALSE.) then
1755      gn0(:,2) = (one - omega(:)*g00(:))/bb(2)
1756      do ii=3,inn
1757        gn0(:,ii) = -(-bb(ii)*gn0(:,ii-2) -omega(:)*gn0(:,ii-1))/bb(ii+1)
1758      end do
1759    else
1760      do ii=2,inn
1761        nlev = inn-ii
1762        call continued_fract(nlev,term_type,aa,bb(ii+1:),nomega,omega,g00)
1763        gn0(:,ii) = +bb(ii+1) * g00 * gn0(:,ii-1)
1764      end do
1765    end if
1766 
1767    newg=czero
1768    do ii=1,inn
1769      newg(:) = newg + cc(ii)* gn0(:,ii)
1770    end do
1771    newg = factor*newg
1772    !
1773    ! Avoid spurious convergence.
1774    niter_min=4; if (niter_done>1) niter_min=niter_done+1
1775    if (inn>niter_min) then
1776      test=.TRUE.
1777      abs_err(:,1) = ABS(DBLE (newg-oldg))
1778      abs_err(:,2) = ABS(AIMAG(newg-oldg))
1779      !
1780      if (tol_iter>zero) then
1781        ! Test on the L1 norm.
1782        if (check(1)) test(1) = SUM(abs_err(:,1)) < tol_iter*SUM(ABS(DBLE (newg)))
1783        if (check(2)) test(2) = SUM(abs_err(:,2)) < tol_iter*SUM(ABS(AIMAG(newg)))
1784      else
1785        ! Stringent test for each point.
1786        if (check(1)) test(1) = ALL( abs_err(:,1) < -tol_iter*ABS(DBLE (newg)))
1787        if (check(2)) test(2) = ALL( abs_err(:,2) < -tol_iter*ABS(AIMAG(newg)))
1788      end if
1789      !
1790      if (ALL(test)) then
1791        nconv = nconv+1
1792      else
1793        nconv = 0
1794      end if
1795      if (nconv==2) then
1796        write(msg,'(a,es10.2,a,i0,a)')&
1797 &        " >>> Haydock algorithm converged twice within haydock_tol= ",tol_iter," after ",inn," iterations."
1798        call wrtout(std_out,msg,'COLL')
1799        call wrtout(ab_out,msg,'COLL')
1800        EXIT
1801      end if
1802    end if
1803    !
1804    oldg = newg
1805  end do ! inn
1806 
1807  green = newg
1808  if (nconv/=2) then
1809    write(msg,'(a,es10.2,a,i0,a)')&
1810 &    " WARNING: Haydock algorithm did not converge within ",tol_iter," after ",niter_tot," iterations."
1811    call wrtout(std_out,msg,'COLL')
1812    call wrtout(ab_out,msg,'COLL')
1813  end if
1814 
1815  is_converged = (nconv==2)
1816 
1817  ABI_FREE(oldg)
1818  ABI_FREE(newg)
1819  ABI_FREE(g00)
1820  ABI_FREE(hphi_n)
1821 
1822  if (keep_vectors) then
1823    tdim = MIN(inn,niter_tot)
1824    ABI_MALLOC(ovlp,(tdim,tdim))
1825 
1826    ABI_MALLOC(phi_test,(hsize))
1827    ABI_MALLOC(phi_test2,(hsize))
1828 
1829    max_err=smallest_real; mean_err=zero; mean_err2=zero; row_max=-1
1830    do ii=1,tdim
1831      parity = (-1)**(ii+1)
1832      phi_test  = save_phi(:,ii)
1833      call hexc_matmul_full(hexc, hexc_i, phi_test, phi_test2, parity)
1834      !phi_test2 = MATMUL(hreso,phi_test) + parity * MATMUL(hcoup,CONJG(phi_test))
1835      ovlp(ii,ii) = DOT_PRODUCT(phi_test,phi_test2) + DOT_PRODUCT(phi_test2,phi_test)
1836      err = ABS(ovlp(ii,ii)-cone)
1837      mean_err  = mean_err + err
1838      mean_err2 = mean_err2 + err**2
1839      if (err > max_err) then
1840        max_err = err
1841        row_max = ii
1842      end if
1843    end do
1844    mean_err = mean_err/tdim
1845    std_dev = mean_err2/tdim -mean_err**2
1846    write(std_out,'(a,i0,1x,3es14.6)')&
1847 &   " Error in normalization (ii, max_err,mean,std_dev): ",row_max,max_err,mean_err,std_dev
1848 
1849    ABI_FREE(phi_test)
1850    ABI_FREE(phi_test2)
1851 
1852    ABI_MALLOC(alpha,(hsize,tdim))
1853 
1854    ! Less efficient but for sake of simplicity with hexc_matmul
1855    ! TODO possibility to call hreso * phi, and hcoup * phi separately
1856    do ii=1,tdim
1857      parity = (-1)**(ii+1)
1858      call hexc_matmul_full(hexc, hexc_i, save_phi(:,ii), alpha(:,ii), parity)
1859    end do
1860 
1861    !alpha = MATMUL(hreso,save_phi(:,1:tdim))
1862    !
1863    !do ii=1,tdim
1864    !  parity = (-1)**(ii+1)
1865    !  alpha(:,ii) =  alpha(:,ii) + parity*MATMUL(hcoup,CONJG(save_phi(:,ii)))
1866    !end do
1867 
1868    ovlp = MATMUL(TRANSPOSE(CONJG(save_phi(:,1:tdim))),alpha)
1869 
1870    ABI_MALLOC(beta,(hsize,tdim))
1871    do ii=1,tdim
1872      parity = (-1)**(ii+1)
1873      beta(:,ii)  =  parity*save_phi(:,ii)
1874      alpha(:,ii) = -parity*alpha(:,ii)
1875    end do
1876 
1877    ovlp = ovlp - MATMUL(TRANSPOSE(CONJG(beta)),alpha)
1878 
1879    max_err=smallest_real; row_max=-1; col_max=-1
1880    mean_err=zero; mean_err2=zero
1881    do jj=1,tdim
1882      do ii=1,jj
1883        err = ABS(ovlp(ii,jj))
1884        if (ii==jj) err = ABS(err - one)
1885        mean_err  = mean_err + err
1886        mean_err2 = mean_err2 + err**2
1887        if (err > max_err) then
1888          max_err = err
1889          row_max=ii
1890          col_max=jj
1891        end if
1892      end do
1893    end do
1894 
1895    mean_err = mean_err/(tdim*(tdim+1)/2)
1896    std_dev = mean_err2/(tdim*(tdim+1)/2) - mean_err**2
1897    write(std_out,'(a,2(i0,1x),3es14.6)')&
1898 &     " Error in Hbar-ortho (i,j), max_err, mean, std_dev ",row_max,col_max,max_err,mean_err,std_dev
1899    !call print_arr(ovlp,max_r=185,max_c=10,unit=std_out)
1900 
1901    ABI_FREE(alpha)
1902    ABI_FREE(beta)
1903    ABI_FREE(ovlp)
1904    ABI_FREE(save_phi)
1905  end if
1906 
1907 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(:)

PARENTS

      m_haydock

CHILDREN

SOURCE

1114 subroutine haydock_restart(BSp,restart_file,ftype,iq_search,hsize,niter_file,aa_file,bb_file,phi_nm1_file,phi_n_file,comm)
1115 
1116 
1117 !This section has been created automatically by the script Abilint (TD).
1118 !Do not modify the following lines by hand.
1119 #undef ABI_FUNC
1120 #define ABI_FUNC 'haydock_restart'
1121 !End of the abilint section
1122 
1123  implicit none
1124 
1125 !Arguments ------------------------------------
1126 !scalars
1127  integer,intent(in) :: comm,hsize,iq_search,ftype
1128  integer,intent(out) :: niter_file
1129  character(len=*),intent(in) :: restart_file
1130  type(excparam),intent(in) :: BSp
1131 !arrays
1132  real(dp),allocatable,intent(out) :: bb_file(:)
1133  complex(dpc),allocatable,intent(out) :: aa_file(:),phi_n_file(:),phi_nm1_file(:)
1134 
1135 !Local variables ------------------------------
1136 !scalars
1137  integer,parameter :: master=0
1138  integer :: nproc,my_rank,ierr,op_file
1139  integer :: hsize_file,use_coupling_file
1140  complex(dpc) :: factor_file
1141  character(len=500) :: msg
1142  type(haydock_type) :: haydock_file
1143 
1144 !************************************************************************
1145 
1146  nproc = xmpi_comm_size(comm); my_rank= xmpi_comm_rank(comm)
1147 
1148  if (my_rank==master) then
1149    call open_haydock(restart_file, haydock_file)
1150 
1151    call read_dim_haydock(haydock_file)
1152 
1153    if (haydock_file%op/=ftype) then
1154      write(msg,"(2(a,i0))")" Expecting restart file with filetype: ",ftype," but found ",op_file
1155      MSG_ERROR(msg)
1156    end if
1157 
1158    if (haydock_file%hsize/=hsize) then
1159      write(msg,"(2(a,i0))")&
1160 &      " Rank of H_exc read from file: ",hsize_file," differs from the one used in this run: ",hsize
1161      MSG_ERROR(msg)
1162    end if
1163 
1164    if (haydock_file%use_coupling /= BSp%use_coupling) then
1165      write(msg,'(2(a,i0))')&
1166 &      " use_coupling_file: ",use_coupling_file," differs from input file value: ",BSp%use_coupling
1167      MSG_ERROR(msg)
1168    end if
1169 
1170    call read_haydock(haydock_file, Bsp%q(:,iq_search), aa_file, bb_file, &
1171 &                   phi_n_file, phi_nm1_file, niter_file, factor_file)
1172 
1173    if (niter_file == 0) then
1174      write(msg,"(a,3f8.4,3a)")&
1175 &      " Could not find q-point: ",BSp%q(:,iq_search)," in file ",TRIM(restart_file),&
1176 &      " Cannot restart Haydock iterations for this q-point"
1177      MSG_COMMENT(msg)
1178    else
1179      write(msg,'(a,i0)')" Number of iterations already performed: ",niter_file
1180      call wrtout(std_out,msg,"COLL")
1181      call wrtout(ab_out,msg,"COLL")
1182 
1183      if ( ABS(haydock_file%broad - BSp%broad) > tol6) then
1184        write(msg,'(2a,2(a,f8.4),a)')&
1185 &        " Restart file has been produced with a different Lorentzian broadening: ",ch10,&
1186 &        " broad_file: ",haydock_file%broad," input broadening: ",BSp%broad," Continuing anyway. "
1187        MSG_WARNING(msg)
1188      end if
1189 
1190      call close_haydock(haydock_file)
1191    end if
1192  end if
1193  !
1194  ! Master broadcasts the data.
1195  call xmpi_bcast(niter_file,master,comm,ierr)
1196 
1197  if (my_rank/=master) then
1198    ABI_MALLOC(aa_file,(niter_file))
1199    ABI_MALLOC(bb_file,(niter_file))
1200    ABI_MALLOC(phi_nm1_file,(hsize))
1201    ABI_MALLOC(phi_n_file,(hsize))
1202  end if
1203 
1204  call xmpi_bcast(aa_file,master,comm,ierr)
1205  call xmpi_bcast(bb_file,master,comm,ierr)
1206  call xmpi_bcast(phi_nm1_file,master,comm,ierr)
1207  call xmpi_bcast(phi_n_file,master,comm,ierr)
1208 
1209 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.

PARENTS

      bsepostproc,m_haydock

CHILDREN

SOURCE

2558 subroutine continued_fract_general(nlev,term_type,aa,bb,cc,nz,zpts,spectrum)
2559 
2560 
2561 !This section has been created automatically by the script Abilint (TD).
2562 !Do not modify the following lines by hand.
2563 #undef ABI_FUNC
2564 #define ABI_FUNC 'continued_fract_general'
2565 !End of the abilint section
2566 
2567  implicit none
2568 
2569 !Arguments ------------------------------------
2570 !scalars
2571  integer,intent(in) :: nlev,term_type,nz
2572 !arrays
2573  complex(dpc),intent(in) :: bb(nlev)
2574  complex(dpc),intent(in) :: cc(nlev)
2575  complex(dpc),intent(in) :: aa(nlev)
2576  complex(dpc),intent(in) :: zpts(nz)
2577  complex(dpc),intent(out) :: spectrum(nz)
2578 
2579 !Local variables ------------------------------
2580 !scalars
2581  integer :: it
2582  complex(dpc) ::  bb_inf,bg,bu,swap
2583  complex(dpc) :: aa_inf
2584  character(len=500) :: msg
2585 !arrays
2586  complex(dpc),allocatable :: div(:),den(:)
2587 
2588 !************************************************************************
2589 
2590  ABI_MALLOC(div,(nz))
2591  ABI_MALLOC(den,(nz))
2592 
2593  select case (term_type)
2594  case (0) ! No terminator.
2595    div=czero
2596  case (-1,1)
2597    MSG_ERROR("Not yet implemented")
2598    if (term_type==-1) then
2599      bb_inf=bb(nlev)
2600      aa_inf=aa(nlev)
2601    else
2602      bb_inf=SUM(bb)/nlev
2603      aa_inf=SUM(aa)/nlev
2604    end if
2605    ! Be careful with the sign of the SQRT.
2606    div(:) = half*(bb(nlev)/(bb_inf))**2 * ( zpts-aa_inf - SQRT((zpts-aa_inf)**2 - four*bb_inf**2) )
2607  case (2)
2608    MSG_ERROR("Not yet implemented")
2609    div = zero
2610    if (nlev>4) then
2611      bg=zero; bu=zero
2612      do it=1,nlev,2
2613        if (it+2<nlev) bg = bg + bb(it+2)
2614        bu = bu + bb(it)
2615      end do
2616      bg = bg/(nlev/2+MOD(nlev,2))
2617      bu = bg/((nlev+1)/2)
2618      !if (iseven(nlev)) then
2619      if (.not.iseven(nlev)) then
2620        swap = bg
2621        bg = bu
2622        bu = bg
2623      end if
2624      !write(std_out,*)nlev,bg,bu
2625      !Here be careful with the sign of SQRT
2626      do it=1,nz
2627        div(it) = half/zpts(it) * (bb(nlev)/bu)**2 * &
2628 &        ( (zpts(it)**2 +bu**2 -bg**2) - SQRT( (zpts(it)**2+bu**2-bg**2)**2 -four*(zpts(it)*bu)**2) )
2629      end do
2630    end if
2631 
2632  case default
2633    write(msg,'(a,i0)')" Wrong value for term_type : ",term_type
2634    MSG_ERROR(msg)
2635  end select
2636 
2637  do it=nlev,2,-1
2638    den(:) = zpts(:) - aa(it) - div(:)
2639    div(:) = (bb(it-1)*cc(it-1) )/ den(:)
2640  end do
2641 
2642  den = zpts(:) - aa(1) - div(:)
2643  div = one/den(:)
2644 
2645  spectrum = div
2646  ABI_FREE(div)
2647  ABI_FREE(den)
2648 
2649 end subroutine continued_fract_general