1! 2! Copyright (C) 2001-2016 Quantum ESPRESSO group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8MODULE lr_lanczos 9 10CONTAINS 11 12!-------------------------------------------------------------------------------------------! 13! 14! There are 2 flavors of the Lanczos algorithm implemented in TDDFPT: 15! - non-Hermitian Lanczos biorthogonalization algorithm 16! - pseudo-Hermitian Lanczos algorithm (twice faster) 17! For details see: 18! - O. Malcioglu, R. Gebauer, D. Rocca, S. Baroni, 19! Comput. Phys. Commun. 182, 1744 (2011). 20! - X. Ge, S. Binnie, D. Rocca, R. Gebauer, S. Baroni, 21! Comput. Phys. Commun. 185, 2080 (2014). 22! 23! Non-Hermitian algorithm: 24! This subroutine handles two interleaved non-Hermitian chains for x and y at the same time, 25! in the beginning 26! for even steps evc1(:,:,:,1) corresponds to q of y and evc1(:,:,:,2) corresponds to p of x 27! for odd steps evc1(:,:,:,1) corresponds to q of x and evc1(:,:,:,2) corresponds to p of y 28! 29! using the terminology of Lanczos chains by Y.Saad, this is effectively equal to 30! altering the A matrix correspondingly for even and odd steps correspondingly. 31! This change is controlled by the interaction parameter in lr_apply_liouvillian 32! 33! For further reference please refer to eq. (32) and (33) in 34! Ralph Gebauer, Brent Walker J. Chem. Phys., 127, 164106 (2007) 35! 36! Modified by Osman Baris Malcioglu (2009) 37! Modified by Xiaochuan Ge to introduce pseudo-Hermitian algorithm (2013) 38! Modified by Iurii Timrov to introduce EELS (2013) 39! 40!--------------------------------------------------------------------------------------------! 41SUBROUTINE one_lanczos_step() 42 43 USE io_global, ONLY : ionode, stdout 44 USE kinds, ONLY : dp 45 USE klist, ONLY : nks,xk 46 USE lr_variables, ONLY : n_ipol, ltammd, pseudo_hermitian, itermax, & 47 evc1, evc1_new, evc1_old, sevc1_new, sevc1, & 48 evc0, sevc0, d0psi, d0psi2, eels, test_case_no, lr_verbosity, & 49 alpha_store, beta_store, gamma_store, zeta_store, & 50 charge_response, size_evc, LR_polarization, LR_iteration 51 USE uspp, ONLY : vkb, nkb, okvan 52 USE wvfct, ONLY : nbnd, npwx 53 USE control_flags, ONLY : gamma_only 54 USE becmod, ONLY : bec_type, becp, calbec 55 USE realus, ONLY : invfft_orbital_gamma, initialisation_level, & 56 fwfft_orbital_gamma, calbec_rs_gamma, add_vuspsir_gamma, & 57 v_loc_psir, s_psir_gamma 58 USE charg_resp, ONLY : w_T_beta_store, w_T, lr_calc_F 59 USE lr_us, ONLY : lr_apply_s 60 USE noncollin_module, ONLY : npol 61 ! Debugging 62 USE iso_c_binding, ONLY : c_int 63 64 USE qpoint, ONLY : nksq 65 ! 66 IMPLICIT NONE 67 ! 68 COMPLEX(kind=dp),EXTERNAL :: lr_dot 69 INTEGER :: ik, ip, ibnd,ig, pol_index 70 CHARACTER(len=6), EXTERNAL :: int_to_char 71 ! 72 ! Local variables 73 ! 74 REAL(kind=dp) :: alpha, beta, gamma, angle 75 COMPLEX(kind=dp) :: zeta(n_ipol) 76 INTEGER(kind=c_int) :: kilobytes 77 ! 78 IF (lr_verbosity > 5) THEN 79 WRITE(stdout,'("<lr_lanczos_one_step>")') 80 ENDIF 81 ! 82 ! Memory usage 83 ! 84 !IF (eels) THEN 85 ! CALL memstat( kilobytes ) 86 ! IF ( kilobytes > 0 ) WRITE(stdout,'(5X,"lr_lanczos, & 87 ! & per-process dynamical memory:",f7.1,"Mb")' ) kilobytes/1000.0 88 !ENDIF 89 ! 90 CALL start_clock('one_step') 91 ! 92 pol_index = 1 93 ! 94 IF ( n_ipol /= 1 ) pol_index = LR_polarization 95 ! 96 ! Application of the Liouvillian superoperator 97 ! 98 IF (eels) THEN 99 ! 100 ! EELS case 101 ! 102 IF (mod(LR_iteration,2)==0) THEN 103 ! 104 CALL lr_apply_liouvillian_eels(evc1(1,1,1,1),evc1_new(1,1,1,1),.true.) 105 IF (.NOT.pseudo_hermitian) & 106 CALL lr_apply_liouvillian_eels(evc1(1,1,1,2),evc1_new(1,1,1,2),.false.) 107 ! 108 ELSE 109 ! 110 CALL lr_apply_liouvillian_eels(evc1(1,1,1,1),evc1_new(1,1,1,1),.false.) 111 IF (.not. pseudo_hermitian) & 112 CALL lr_apply_liouvillian_eels(evc1(1,1,1,2),evc1_new(1,1,1,2),.true.) 113 ! 114 ENDIF 115 ! 116 ELSE 117 ! 118 ! Optical case 119 ! 120 IF (.NOT.ltammd) THEN 121 ! 122 ! Normal regime 123 ! 124 IF (mod(LR_iteration,2)==0) THEN 125 ! 126 CALL lr_apply_liouvillian(evc1(1,1,1,1),evc1_new(1,1,1,1),.true.) 127 IF (.NOT.pseudo_hermitian) & 128 CALL lr_apply_liouvillian(evc1(1,1,1,2),evc1_new(1,1,1,2),.false.) 129 ! 130 ELSE 131 ! 132 CALL lr_apply_liouvillian(evc1(1,1,1,1),evc1_new(1,1,1,1),.false.) 133 IF (.not. pseudo_hermitian) & 134 CALL lr_apply_liouvillian(evc1(1,1,1,2),evc1_new(1,1,1,2),.true.) 135 ! 136 ENDIF 137 ! 138 ELSE 139 ! 140 ! Tamm-Dancoff approximation 141 ! 142 CALL lr_apply_liouvillian(evc1(1,1,1,1),evc1_new(1,1,1,1),.true.) 143 ! 144 IF (.not. pseudo_hermitian) THEN 145 CALL zcopy(size_evc,evc1(1,1,1,1),1,evc1(1,1,1,2),1) 146 CALL zcopy(size_evc,evc1_new(1,1,1,1),1,evc1_new(1,1,1,2),1) 147 ENDIF 148 ! 149 ENDIF 150 ! 151 ENDIF 152 ! 153 ! Apply S operator if USPP, otherwise just copy 154 ! one array into another. 155 ! 156 IF (pseudo_hermitian) THEN 157 CALL lr_apply_s(evc1_new(:,:,:,1), sevc1_new(:,:,:)) 158 ELSE 159 CALL lr_apply_s(evc1(:,:,:,2), sevc1(:,:,:)) 160 ENDIF 161 ! 162 ! call general lanczos iteration routines 163 ! O. Baseggio (2019) 164 ! 165 IF (pseudo_hermitian) THEN 166 IF (eels) THEN 167 CALL lanczos_pseudohermitian(LR_iteration,size(evc1,1), size(evc1,2), size(evc1,3),& 168 &evc1(:,:,:,1), evc1_new(:,:,:,1), sevc1_new(:,:,:), & 169 &evc1_old(:,:,:,1), n_ipol, d0psi2, alpha, beta, & 170 &gamma, zeta) 171 ELSE 172 CALL lanczos_pseudohermitian(LR_iteration,size(evc1,1), size(evc1,2), size(evc1,3),& 173 &evc1(:,:,:,1), evc1_new(:,:,:,1), sevc1_new(:,:,:), & 174 &evc1_old(:,:,:,1), n_ipol, d0psi(:,:,:,:), alpha, beta, & 175 &gamma, zeta) 176 ENDIF 177 ELSE 178 IF (eels) THEN 179 CALL lanczos_nonhermitian(LR_iteration,size(evc1,1), size(evc1,2), size(evc1,3),& 180 &evc1(:,:,:,:), evc1_new(:,:,:,:), sevc1(:,:,:), & 181 &evc1_old(:,:,:,1), n_ipol, d0psi2, alpha, beta, & 182 &gamma, zeta) 183 ELSE 184 CALL lanczos_nonhermitian(LR_iteration,size(evc1,1), size(evc1,2), size(evc1,3),& 185 &evc1(:,:,:,1), evc1_new(:,:,:,1), sevc1(:,:,:), & 186 &evc1_old(:,:,:,1), n_ipol, d0psi, alpha, beta, & 187 &gamma, zeta) 188 ENDIF 189 ENDIF 190 ! 191 alpha_store(pol_index,LR_iteration) = alpha 192 WRITE(stdout,'(5X,"alpha(",i8.8,")=",f10.6)') LR_iteration, alpha 193 ! 194 ! beta<0 is a serious error for the pseudo-Hermitian algorithm 195 ! 196 IF ( beta < 0.0d0 .and. pseudo_hermitian) CALL error_beta() 197 ! 198 IF ( abs(beta)<1.0d-12 ) THEN 199 ! 200 WRITE(stdout,'(5x,"lr_lanczos: Left and right Lanczos vectors are orthogonal, & 201 & this is a violation of oblique projection")') 202 ! 203 ENDIF 204 ! 205 beta_store (pol_index,LR_iteration) = beta 206 gamma_store(pol_index,LR_iteration) = gamma 207 WRITE(stdout,'(5X,"beta (",i8.8,")=",f10.6)') LR_iteration, beta 208 WRITE(stdout,'(5X,"gamma(",i8.8,")=",f10.6)') LR_iteration, gamma 209 ! 210 !IF (LR_iteration==1) WRITE(stdout,'(5X,"Norm of initial Lanczos vectors=",1x,f21.15)') beta 211 ! 212 ! Analysis of the evolution of the beta coefficients. 213 ! 214 IF (lr_verbosity > 3) THEN 215 ! 216 IF ( LR_iteration > 6 ) THEN 217 ! 218 WRITE(stdout,'(5X,"(oscillatory variation) : ",f6.2,"%")') & 219 & abs(100*(beta_store(pol_index,LR_iteration)- & 220 & (beta_store(pol_index,LR_iteration-6)+beta_store(pol_index,LR_iteration-4)+& 221 & beta_store(pol_index,LR_iteration-2))/3.0)/beta_store(pol_index,LR_iteration)) 222 ! 223 WRITE(stdout,'(5X,"(linear variation) : ",f6.2,"%")') & 224 & abs(100*(beta_store(pol_index,LR_iteration)- & 225 & (beta_store(pol_index,LR_iteration-3)+beta_store(pol_index,LR_iteration-2)+& 226 & beta_store(pol_index,LR_iteration-1))/3.0)/beta_store(pol_index,LR_iteration)) 227 ! 228 ENDIF 229 ! 230 ENDIF 231 ! 232 IF (mod(LR_iteration,2)==0) THEN 233 ! 234 DO ip = 1, n_ipol 235 ! 236 zeta_store (pol_index,ip,LR_iteration) = zeta(ip) 237 WRITE(stdout,'(5x,"z1= ",1x,i6,2(1x,e22.15))') ip,real(zeta(ip)),aimag(zeta(ip)) 238 ! 239 ENDDO 240 ! 241 ! OBM: evc1(:,:,:,1) contains the q of x for even steps, 242 ! lets calculate the response related observables. 243 ! 244 IF (charge_response == 1 .and. .not.eels) THEN 245 CALL lr_calc_dens(evc1_old(:,:,:,1), .true.) 246 CALL lr_calc_F(evc1_old(:,:,:,1)) 247 ENDIF 248 ! 249 ELSE 250 ! 251 DO ip = 1, n_ipol 252 ! 253 zeta_store (pol_index,ip,LR_iteration) = zeta(ip) 254 WRITE(stdout,'(5x,"z1= ",1x,i6,2(1x,e22.15))') ip,real(zeta(ip)),aimag(zeta(ip)) 255 ! 256 ENDDO 257 ! 258 ENDIF 259 ! 260 ! X. Ge: To increase the stability, apply lr_ortho. 261 ! I.Timrov: Actually, without this trick, it turns out that 262 ! the Lanczos chain is not stable when pseudo_hermitian=.false., 263 ! because there is a warning from lr_calc_dens that the integral of 264 ! the response charge density does not some to zero, i.e. the non-zero 265 ! value of the integral increases during the Lanczos chain. 266 ! 267 IF (.not.eels) THEN 268 ! 269 DO ik=1, nks 270 CALL lr_ortho(evc1(:,:,ik,1), evc0(:,:,ik), ik, ik, sevc0(:,:,ik),.true.) 271 IF (.not. pseudo_hermitian) & 272 CALL lr_ortho(evc1(:,:,ik,2), evc0(:,:,ik), ik, ik, sevc0(:,:,ik),.true.) 273 ENDDO 274 ! 275 ENDIF 276 ! 277 IF (ionode) THEN 278 IF ( charge_response == 1 .and. lr_verbosity > 0) THEN 279 WRITE (stdout,'(5x,"(calc=",e22.15," read=",e22.15,")")') & 280 & beta, w_T_beta_store(LR_iteration) 281 WRITE (stdout,'(5x,"Weight for this step=",2(e12.5,1x))') w_T(LR_iteration) 282 ENDIF 283 ENDIF 284 ! 285 CALL stop_clock('one_step') 286 ! 287 RETURN 288 ! 289CONTAINS 290 ! 291 SUBROUTINE error_beta() 292 ! 293 ! IT: This subroutine is called when beta<0 in the pseudo-Hermitian 294 ! algorithm. This subroutine prints the response orbitals to 295 ! the file for the analysis of the error. Note, in the parallel 296 ! execution only a part of the information about the response 297 ! orbitals will be written to the file. 298 ! 299 ! Written by X. Ge (2013) 300 ! 301 IMPLICIT NONE 302 ! 303 INTEGER :: ibnd, ig 304 ! 305 WRITE(stdout,'(5x,"Error: Beta is negative:",5x,e22.15)') beta 306 ! 307 OPEN(301, file="evc1_1.dat") 308 DO ibnd = 1, nbnd 309 DO ig = 1, npwx 310 WRITE(301,*) ibnd, ig, dble(evc1(ig,ibnd,1,1)), aimag(evc1(ig,ibnd,1,1)) 311 ENDDO 312 ENDDO 313 CLOSE(301) 314 ! 315 OPEN(302, file="evc1_2.dat") 316 DO ibnd = 1, nbnd 317 DO ig = 1, npwx 318 WRITE(302,*) ibnd, ig, dble(evc1(ig,ibnd,1,2)), aimag(evc1(ig,ibnd,1,2)) 319 ENDDO 320 ENDDO 321 CLOSE(302) 322 ! 323 WRITE (stdout,'(/5x,"!!!WARNING!!! The pseudo-Hermitian Lanczos algorithm is stopping...")') 324 WRITE (stdout,'(/5x,"Try to use the non-Hermitian Lanczos algorithm.")') 325 ! 326 itermax = LR_iteration-1 327 CALL clean_pw( .FALSE. ) 328 CALL stop_lr( .TRUE. ) 329 ! 330 END SUBROUTINE error_beta 331 332END SUBROUTINE one_lanczos_step 333 334END MODULE lr_lanczos 335