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