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! 8!----------------------------------------------------------------------- 9SUBROUTINE lr_read_wf() 10 !--------------------------------------------------------------------- 11 ! 12 ! ... reads in and stores the ground state wavefunctions 13 ! ... for use in Lanczos linear response calculation 14 ! 15 ! Modified by Osman Baris Malcioglu (2009) 16 ! Modified by Xiaochuan Ge (2013) to fix some bugs of virt_read and include Davidson 17 ! 18 USE kinds, ONLY : dp 19 USE cell_base, ONLY : at 20 USE io_global, ONLY : stdout 21 USE klist, ONLY : nks, xk, ngk, igk_k 22 USE gvect, ONLY : ngm, g 23 USE io_files, ONLY : nwordwfc, iunwfc, prefix, diropn,& 24 & tmp_dir, wfc_dir 25 USE lr_variables, ONLY : evc0, sevc0 ,revc0, evc0_virt, & 26 & sevc0_virt, nbnd_total, becp1_virt, & 27 & becp1_c_virt, no_hxc, becp_1, becp1_c, & 28 & test_case_no, size_evc, project, & 29 & lr_verbosity, lr_exx, davidson, eels 30 USE wvfct, ONLY : nbnd, npwx 31 USE control_flags, ONLY : gamma_only,io_level 32 USE fft_base, ONLY : dffts 33 USE fft_interfaces, ONLY : invfft 34 USE uspp, ONLY : vkb, nkb, okvan 35 USE becmod, ONLY : bec_type, becp, calbec 36 USE realus, ONLY : real_space, invfft_orbital_gamma,& 37 initialisation_level,& 38 fwfft_orbital_gamma, calbec_rs_gamma,& 39 add_vuspsir_gamma, v_loc_psir,& 40 s_psir_gamma 41 USE funct, ONLY : dft_is_hybrid 42 USE lr_exx_kernel, ONLY : lr_exx_revc0_init, lr_exx_alloc, & 43 lr_exx_restart 44 USE wavefunctions, ONLY : evc 45 USE buffers, ONLY : open_buffer 46 USE qpoint, ONLY : nksq 47 USE noncollin_module, ONLY : npol 48 USE symm_base, ONLY : fft_fact 49 USE fft_helper_subroutines 50 ! 51 IMPLICIT NONE 52 ! 53 ! local variables 54 ! 55 INTEGER :: ik, ibnd, ig, itmp1,itmp2,itmp3 56 LOGICAL :: exst 57 CHARACTER(len=256) :: filename, tmp_dir_saved 58 ! 59 IF (lr_verbosity > 5) THEN 60 WRITE(stdout,'("<lr_read_wf>")') 61 ENDIF 62 ! 63 CALL start_clock("read_wf") 64 ! 65 IF ((nbnd_total>nbnd .and. davidson) .OR. project) THEN 66 CALL virt_read() 67 ELSE 68 CALL normal_read() 69 ENDIF 70 ! 71 IF (.NOT.eels) evc(:,:) = evc0(:,:,1) 72 ! 73 IF ( dft_is_hybrid() ) THEN 74 ! 75 ! Initialize fft_fact 76 ! Warning: If there are fractional translations and 77 ! they are not commensurate with the dfftt grid, then 78 ! fft_fact is different from (1,1,1) and it must be 79 ! properly initialized. This matters when a symmetrization 80 ! is real space is used. 81 ! 82 fft_fact(:) = 1 83 ! 84 CALL open_buffer ( iunwfc, 'wfc', nwordwfc, io_level, exst ) 85 ! 86 ! set_ace=.false. disables Lin Lin's ACE for TD-DFPT 87 ! 88 CALL lr_exx_restart( set_ace=.false.) 89 ! 90 IF (.NOT. no_hxc) THEN 91 ! 92 lr_exx = .TRUE. 93 ! 94 CALL lr_exx_alloc() 95 ! 96 DO ik=1,nks 97 CALL lr_exx_revc0_init(evc0,ik) 98 ENDDO 99 ! 100 ENDIF 101 ! 102 WRITE(stdout,'(5x,"Finished exx setting.")') 103 ! 104 ENDIF 105 ! 106 CALL stop_clock("read_wf") 107 ! 108 RETURN 109 ! 110 CONTAINS 111 112SUBROUTINE normal_read() 113 ! 114 ! The usual way of reading wavefunctions 115 ! 116 USE lr_variables, ONLY : tg_revc0 117 USE wavefunctions, ONLY : psic 118 USE realus, ONLY : tg_psic 119 USE mp_global, ONLY : me_bgrp 120 ! 121 IMPLICIT NONE 122 ! 123 INTEGER :: v_siz, incr, ioff, j 124 ! 125 WRITE( stdout, '(/5x,"Normal read")' ) 126 ! 127 incr = 2 128 ! 129 size_evc = nbnd * npwx * npol * nksq 130 nwordwfc = nbnd * npwx * npol 131 ! 132 ! Read in the ground state wavefunctions. 133 ! This is a parallel read, done in wfc_dir. 134 ! 135 tmp_dir_saved = tmp_dir 136 ! 137 !IF ( wfc_dir /= 'undefined' ) tmp_dir = wfc_dir 138 ! 139 wfc_dir = tmp_dir 140 ! 141 CALL diropn ( iunwfc, 'wfc', 2*nwordwfc, exst) 142 ! 143 IF (.NOT.exst .AND. wfc_dir == 'undefined') & 144 & CALL errore('lr_read_wfc', TRIM( prefix )//'.wfc'//' not found',1) 145 ! 146 IF (.NOT.exst .AND. wfc_dir /= 'undefined') THEN 147 ! 148 WRITE( stdout, '(/5x,"Attempting to read wfc from outdir instead of wfcdir")' ) 149 CLOSE( UNIT = iunwfc) 150 tmp_dir = tmp_dir_saved 151 CALL diropn ( iunwfc, 'wfc', 2*nwordwfc, exst) 152 IF (.NOT.exst) CALL errore('lr_read_wfc', TRIM( prefix )//'.wfc'//' not found',1) 153 ! 154 ENDIF 155 ! 156 IF (gamma_only) THEN 157 WRITE( stdout, '(/5x,"Gamma point algorithm")' ) 158 ELSE 159 WRITE( stdout, '(/5x,"WARNING: Generalised k-points algorithm")' ) 160 ENDIF 161 ! 162 DO ik = 1, nks 163 ! 164 CALL davcio(evc0(:,:,ik),2*nwordwfc,iunwfc,ik,-1) 165 ! 166 ENDDO 167 ! 168 CLOSE( UNIT = iunwfc) 169 ! 170 ! End of file reading 171 ! 172 tmp_dir = tmp_dir_saved 173 ! 174 ! vkb * evc0, and initialization of sevc0. 175 ! 176 IF ( okvan ) THEN 177 ! 178 IF ( gamma_only ) THEN 179 ! 180 ! Following line is to be removed when real space 181 ! implementation is complete. 182 ! 183 CALL init_us_2(ngk(1),igk_k(:,1),xk(1,1),vkb) 184 ! 185 IF (real_space) THEN 186 ! 187 DO ibnd = 1, nbnd, 2 188 ! 189 CALL invfft_orbital_gamma(evc0(:,:,1),ibnd,nbnd) 190 CALL calbec_rs_gamma(ibnd,nbnd,becp_1) 191 becp%r(:,ibnd) = becp_1(:,ibnd) 192 IF (ibnd + 1 <= nbnd) becp%r(:,ibnd+1) = becp_1(:,ibnd+1) 193 CALL s_psir_gamma(ibnd,nbnd) 194 CALL fwfft_orbital_gamma(sevc0(:,:,1),ibnd,nbnd) 195 ! 196 ENDDO 197 ! 198 ELSE 199 ! 200 CALL calbec(ngk(1),vkb,evc0(:,:,1),becp_1) 201 becp%r = becp_1 202 CALL s_psi(npwx, ngk(1), nbnd, evc0(:,:,1), sevc0(:,:,1)) 203 ! 204 ENDIF 205 ! 206 ELSE 207 ! 208 ! K point generalized stuff starts here 209 ! 210 DO ik = 1, nks 211 ! 212 CALL init_us_2(ngk(ik),igk_k(1,ik),xk(1,ik),vkb) 213 CALL calbec(ngk(ik),vkb,evc0(:,:,ik),becp1_c(:,:,ik)) 214 becp%k = becp1_c(:,:,ik) 215 CALL s_psi (npwx, ngk(ik), nbnd, evc0(:,:,ik), sevc0(:,:,ik)) 216 ! 217 ENDDO 218 ! 219 ENDIF 220 ! 221 ELSE 222 ! 223 sevc0 = evc0 224 ! 225 ENDIF 226 ! 227 ! Calculation of the unperturbed wavefunctions in R-space revc0. 228 ! Inverse Fourier transform of evc0. 229 ! 230 IF ( dffts%has_task_groups ) THEN 231 ! 232 v_siz = dffts%nnr_tg 233 incr = 2 * fftx_ntgrp(dffts) 234 tg_revc0 = (0.0d0,0.0d0) 235 ! 236 ELSE 237 ! 238 revc0 = (0.0d0,0.0d0) 239 ! 240 ENDIF 241 ! 242 IF ( gamma_only ) THEN 243 ! 244 DO ibnd = 1, nbnd, incr 245 ! 246 CALL invfft_orbital_gamma ( evc0(:,:,1), ibnd, nbnd) 247 ! 248 IF (dffts%has_task_groups) THEN 249 ! 250 DO j = 1, dffts%nr1x*dffts%nr2x*dffts%my_nr3p 251 ! 252 tg_revc0(j,ibnd,1) = tg_psic(j) 253 ! 254 ENDDO 255 ! 256 ELSE 257 ! 258 revc0(1:dffts%nnr,ibnd,1) = psic(1:dffts%nnr) 259 ! 260 ENDIF 261 ! 262 ENDDO 263 ! 264 ELSE 265 ! 266 ! The FFT is done in the same way as in invfft_orbital_k 267 ! (where also the task groups is implemented but must be checked). 268 ! 269 DO ik = 1, nks 270 DO ibnd = 1, nbnd 271 DO ig = 1, ngk(ik) 272 ! 273 revc0(dffts%nl(igk_k(ig,ik)),ibnd,ik) = evc0(ig,ibnd,ik) 274 ! 275 ENDDO 276 ! 277 CALL invfft ('Wave', revc0(:,ibnd,ik), dffts) 278 ! 279 ENDDO 280 ENDDO 281 ! 282 ENDIF 283 ! 284 IF ( real_space .AND. .NOT. gamma_only ) & 285 CALL errore( ' iosys ', ' Linear response calculation ' // & 286 & 'real space algorithms with k-points not implemented', 1 ) 287 ! 288 RETURN 289 ! 290END SUBROUTINE normal_read 291 292 293SUBROUTINE virt_read() 294 ! 295 ! The modifications to read also the virtual orbitals. 296 ! 297 USE control_lr, ONLY : nbnd_occ 298 USE becmod, ONLY : allocate_bec_type, deallocate_bec_type 299 ! 300 IMPLICIT NONE 301 ! 302 REAL(kind=dp), ALLOCATABLE :: becp1_all(:,:) 303 COMPLEX(kind=dp), ALLOCATABLE :: evc_all(:,:,:), sevc_all(:,:,:), & 304 becp1_c_all(:,:,:), revc_all(:,:,:) 305 ! 306 ! Check for task groups 307 ! 308 WRITE( stdout, '(/5x,"Virt read")' ) 309 ! 310 IF (dffts%has_task_groups) CALL errore ( 'virt_read', 'Task & 311 & groups not supported when there are virtual states in the & 312 & input.', 1 ) 313 ! 314 ! First pretend everything is normal 315 ! 316 nbnd = nbnd_total 317 ! 318 ALLOCATE(revc_all(dffts%nnr,nbnd,nks)) 319 ALLOCATE(evc_all(npwx,nbnd,nks)) 320 ALLOCATE(sevc_all(npwx,nbnd,nks)) 321 ALLOCATE(evc0_virt(npwx,(nbnd-nbnd_occ(1)),nks)) 322 ALLOCATE(sevc0_virt(npwx,(nbnd-nbnd_occ(1)),nks)) 323 ! 324 ! Xiaochun Ge: It's kind of messy to deallocate and alloacte again becp. This is mainly due to the 325 ! change of nbnd. Later this operation need to be done again when we change nbnd back to nbnd_occ. 326 ! Again, I suggest please, in the future try the best not to change the meaning of global variables. 327 ! 328 CALL deallocate_bec_type ( becp ) 329 CALL allocate_bec_type ( nkb, nbnd, becp ) 330 ! 331 IF (nkb > 0) THEN 332 ! 333 IF (gamma_only) THEN 334 ALLOCATE(becp1_all(nkb,nbnd)) 335 becp1_all(:,:)=0.0d0 336 ELSE 337 ALLOCATE(becp1_c_all(nkb,nbnd,nks)) 338 becp1_c_all(:,:,:)=(0.0d0,0.0d0) 339 ENDIF 340 ! 341 ENDIF 342 ! 343 size_evc = nbnd_occ(1) * npwx * npol * nksq 344 nwordwfc = nbnd * npwx * npol ! nbnd > nbnd_occ(1) 345 ! 346 ! Read in the ground state wavefunctions 347 ! This is a parallel read, done in wfc_dir 348 ! 349 tmp_dir_saved = tmp_dir 350 ! 351 !IF ( wfc_dir /= 'undefined' ) tmp_dir = wfc_dir 352 ! 353 wfc_dir = tmp_dir 354 ! 355 CALL diropn ( iunwfc, 'wfc', 2*nwordwfc, exst) 356 ! 357 IF (.NOT.exst .AND. wfc_dir == 'undefined') & 358 & CALL errore('lr_read_wfc', TRIM( prefix )//'.wfc'//' not found',1) 359 ! 360 IF (.NOT.exst .AND. wfc_dir /= 'undefined') THEN 361 WRITE( stdout, '(/5x,"Attempting to read from outdir instead of wfcdir")' ) 362 CLOSE( UNIT = iunwfc) 363 tmp_dir = tmp_dir_saved 364 CALL diropn ( iunwfc, 'wfc', 2*nwordwfc, exst) 365 IF (.NOT.exst) CALL errore('lr_read_wfc', TRIM( prefix )//'.wfc'//' not found',1) 366 ENDIF 367 ! 368 IF (gamma_only) THEN 369 WRITE( stdout, '(/5x,"Gamma point algorithm")' ) 370 ELSE 371 WRITE( stdout, '(/5x,"WARNING: Generalised k-points algorithm")' ) 372 ENDIF 373 ! 374 ! Read in the ground state wavefunctions. 375 ! This is a parallel read, done in wfc_dir. 376 ! 377 DO ik = 1, nks 378 CALL davcio(evc_all(:,:,ik),2*nwordwfc,iunwfc,ik,-1) 379 ENDDO 380 ! 381 CLOSE( UNIT = iunwfc) 382 ! 383 ! End of file reading 384 ! 385 tmp_dir = tmp_dir_saved 386 ! 387 ! vkb * evc_all and initialization of sevc_all 388 ! 389 IF ( okvan ) THEN 390 ! 391 IF ( gamma_only ) THEN 392 ! 393 ! Following line is to be removed when real space 394 ! implementation is complete. 395 ! 396 CALL init_us_2(ngk(1),igk_k(:,1),xk(1,1),vkb) 397 ! 398 IF (real_space) THEN 399 ! 400 DO ibnd=1,nbnd,2 401 CALL invfft_orbital_gamma(evc_all(:,:,1),ibnd,nbnd) 402 CALL calbec_rs_gamma(ibnd,nbnd,becp1_all) 403 becp%r(:,ibnd) = becp1_all(:,ibnd) 404 IF (ibnd + 1 <= nbnd) becp%r(:,ibnd+1) = becp1_all(:,ibnd+1) 405 CALL s_psir_gamma(ibnd,nbnd) 406 CALL fwfft_orbital_gamma(sevc_all(:,:,1),ibnd,nbnd) 407 ENDDO 408 ! 409 ELSE 410 CALL calbec(ngk(1),vkb,evc_all(:,:,1),becp1_all) 411 becp%r=becp1_all 412 CALL s_psi(npwx, ngk(1), nbnd, evc_all(:,:,1), sevc_all(:,:,1)) 413 ENDIF 414 ! 415 ELSE 416 ! 417 ! K point generalized case 418 ! 419 DO ik = 1, nks 420 ! 421 CALL init_us_2(ngk(ik),igk_k(1,ik),xk(1,ik),vkb) 422 CALL calbec(ngk(ik),vkb,evc_all(:,:,ik),becp1_c_all(:,:,ik),nbnd) 423 becp%k=becp1_c_all(:,:,ik) 424 CALL s_psi (npwx, ngk(ik), nbnd, evc_all(:,:,ik), sevc_all(:,:,ik)) 425 ! 426 ENDDO 427 ! 428 ENDIF 429 ! 430 ELSE 431 ! 432 sevc_all = evc_all 433 ! 434 ENDIF 435 ! 436 ! Calculation of the unperturbed wavefunctions in R-space revc0_all. 437 ! Inverse fourier transform of evc_all 438 ! 439 revc_all = (0.0d0,0.0d0) 440 ! 441 ! X. Ge: Very important, otherwise there will be bugs. 442 ! 443 nbnd = nbnd_occ(1) 444 ! 445 nwordwfc = nbnd * npwx * npol ! needed for EXX 446 ! 447 CALL deallocate_bec_type(becp) 448 CALL allocate_bec_type ( nkb, nbnd, becp ) 449 ! 450 IF ( gamma_only ) THEN 451 ! 452 ! The FFT is done in the same way as in invfft_orbital_gamma 453 ! 454 DO ibnd=1,nbnd,2 455 IF (ibnd<nbnd) THEN 456 DO ig=1,ngk(1) 457 ! 458 revc_all(dffts%nl(igk_k(ig,1)),ibnd,1) = evc_all(ig,ibnd,1)& 459 &+(0.0d0,1.0d0)*evc_all(ig,ibnd+1,1) 460 revc_all(dffts%nlm(igk_k(ig,1)),ibnd,1) = & 461 &CONJG(evc_all(ig,ibnd,1)& 462 &-(0.0d0,1.0d0)*evc_all(ig,ibnd+1,1)) 463 ! 464 ENDDO 465 ELSE 466 DO ig=1,ngk(1) 467 ! 468 revc_all(dffts%nl(igk_k(ig,1)),ibnd,1) = evc_all(ig,ibnd,1) 469 revc_all(dffts%nlm(igk_k(ig,1)),ibnd,1) = CONJG(evc_all(ig,ibnd,1)) 470 ! 471 ENDDO 472 ENDIF 473 ! 474 CALL invfft ('Wave', revc_all(:,ibnd,1), dffts) 475 ! 476 ENDDO 477 ! 478 ELSE 479 ! 480 ! The FFT is done in the same way as in invfft_orbital_k 481 ! 482 DO ik=1,nks 483 DO ibnd=1,nbnd 484 DO ig=1,ngk(ik) 485 ! 486 revc_all(dffts%nl(igk_k(ig,ik)),ibnd,ik) = evc_all(ig,ibnd,ik) 487 ! 488 ENDDO 489 ! 490 CALL invfft ('Wave', revc_all(:,ibnd,ik), dffts) 491 ! 492 ENDDO 493 ENDDO 494 ! 495 ENDIF 496 ! 497 ! Now everything goes into right place 498 ! 499 evc0 = (0.0d0,0.0d0) 500 evc0(:,:,:) = evc_all(:,1:nbnd,:) 501 ! 502 sevc0 = (0.0d0,0.0d0) 503 sevc0(:,:,:) = sevc_all(:,1:nbnd,:) 504 ! 505 revc0 = (0.0d0,0.0d0) 506 revc0(:,:,:) = revc_all(:,1:nbnd,:) 507 ! 508 IF (nkb>0) THEN 509 ! 510 IF (gamma_only) THEN 511 becp_1(:,:)=becp1_all(:,1:nbnd) 512 becp%r=0.0d0 513 becp%r=becp_1 514 ELSE 515 becp1_c(:,:,:)=becp1_c_all(:,1:nbnd,:) 516 becp%k=(0.0d0,0.0d0) 517 becp%k=becp1_c(:,:,1) 518 ENDIF 519 ! 520 ENDIF 521 ! 522 ! Finally retain the conduction bands if needed for projection 523 ! 524 evc0_virt(:,:,:) = evc_all(:,nbnd+1:nbnd_total,:) 525 sevc0_virt(:,:,:) = sevc_all(:,nbnd+1:nbnd_total,:) 526 ! 527 IF (nkb>0) THEN 528 ! 529 IF (gamma_only) THEN 530 becp1_virt(:,:) = becp1_all(:,nbnd+1:nbnd_total) 531 DEALLOCATE(becp1_all) 532 ELSE 533 becp1_c_virt(:,:,:) = becp1_c_all(:,nbnd+1:nbnd_total,:) 534 DEALLOCATE(becp1_c_all) 535 ENDIF 536 ! 537 ENDIF 538 ! 539 DEALLOCATE(evc_all) 540 DEALLOCATE(sevc_all) 541 DEALLOCATE(revc_all) 542 ! 543 IF ( real_space .and. .not. gamma_only ) & 544 & CALL errore( ' iosys ', ' Linear response calculation ' // & 545 & 'real space algorithms with k-points not implemented', 1 ) 546 ! 547END SUBROUTINE virt_read 548 549END SUBROUTINE lr_read_wf 550