1! 2! Copyright (C) 2001-2015 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!-------------------------------------------------------------------- 9 10MODULE lr_exx_kernel 11 !------------------------------------------------------------------ 12 ! 13 ! This module handles the EXX part of the Liouvillian. It requires 14 ! all the correct initialization to have been done for the routines 15 ! in PW/src/exx.f90 (if you can't call vexx() then this module 16 ! won't work). It respects the divergence treatments etc set in 17 ! exx.f90 and also uses the custom grid defined there for gamma_only 18 ! calculations. 19 ! 20 ! lr_exx_alloc must be called first to allocate the appropriate 21 ! workspaces, lr_exx_dealloc can be used to dealloc them. Also 22 ! lr_exx_revc0_init must be called at the beginning to set up the 23 ! EXX operator. 24 ! 25 !------------------------------------------------------------------ 26 ! 27 28 USE kinds, ONLY : DP 29 USE lr_variables, ONLY : evc0, revc0 30 USE fft_base, ONLY : dffts 31 USE fft_interfaces, ONLY : invfft, fwfft 32 USE lsda_mod, ONLY : nspin 33 USE wvfct, ONLY : nbnd, npwx, wg 34 USE gvect, ONLY : g, ngm 35 USE klist, ONLY : xk, wk, nks 36 USE lr_variables, ONLY : gamma_only, lr_verbosity 37 USE realus, ONLY : invfft_orbital_gamma, fwfft_orbital_gamma,& 38 & invfft_orbital_k, fwfft_orbital_k 39 USE wavefunctions, ONLY : psic 40 USE cell_base, ONLY : omega 41 USE exx_base, ONLY : g2_convolution 42 USE exx, ONLY : exxalfa, npwt, gt, dfftt 43 44 45 REAL(kind=dp), PUBLIC, ALLOCATABLE :: revc_int(:,:) 46 COMPLEX(kind=dp), PUBLIC, ALLOCATABLE :: revc_int_c(:,:,:) 47 48 ! Workspaces used in the k*d_term functions. 49 50 ! Used for the density like \phi(r') \psi(r') products 51 COMPLEX(DP), PRIVATE, ALLOCATABLE :: pseudo_dens_c(:) 52 ! Used to hold the result of \int \phi(r') \psi(r') 1/|r_12| dr' 53 COMPLEX(DP), PRIVATE, ALLOCATABLE :: vhart(:,:) 54 55 ! This holds the groundstate orbitals, in real space on the 56 ! custom EXX grid. 57 COMPLEX(DP), PRIVATE, ALLOCATABLE :: red_revc0(:,:,:) 58 59 ! This provides the k -> q point correspondence for the 60 ! EXX operator. 61 INTEGER, PRIVATE, ALLOCATABLE :: k2q(:) 62 63CONTAINS 64 !------------------------------------------------------------------------ 65 SUBROUTINE lr_exx_restart( set_ace ) 66 !------------------------------------------------------------------------ 67 !This SUBROUTINE is called when restarting an exx calculation 68 USE funct, ONLY : get_exx_fraction, start_exx, & 69 exx_is_active, get_screening_parameter 70 USE cell_base, ONLY : at 71 USE exx_base, ONLY : exxdiv, erfc_scrlen, exx_divergence, exx_grid_init,& 72 exx_div_check 73 ! FIXME: are these variable useful? 74 USE exx, ONLY : fock0, exxenergy2, local_thr, use_ace 75 USE exx, ONLY : exxinit, aceinit, exx_gvec_reinit 76 77 IMPLICIT NONE 78 LOGICAL, INTENT(in) :: set_ace 79 ! 80 CALL exx_grid_init( reinit=.true. ) 81 CALL exx_gvec_reinit( at ) 82 CALL exx_div_check() 83 ! 84 use_ace = set_ace 85 erfc_scrlen = get_screening_parameter() 86 87 exxdiv = exx_divergence() 88 exxalfa = get_exx_fraction() 89 CALL start_exx() 90 CALL weights() 91 ! FIXME: is this useful ? 92 IF(local_thr.gt.0.0d0) CALL errore('exx_restart','SCDM with restart NYI',1) 93 CALL exxinit(.false.) 94 IF (use_ace) CALL aceinit ( DOLOC = .FALSE. ) 95 ! FIXME: are these variable useful? 96 fock0 = exxenergy2() 97 ! 98 RETURN 99 !------------------------------------------------------------------------ 100 END SUBROUTINE lr_exx_restart 101 !------------------------------------------------------------------------ 102 103 SUBROUTINE lr_exx_alloc() 104 105 USE exx_base, ONLY : nkqs 106 USE klist, ONLY : nks 107 108 IMPLICIT NONE 109 INTEGER :: nrxxs 110 111 IF (gamma_only) THEN 112 nrxxs= dfftt%nnr 113 ELSE 114 nrxxs= dffts%nnr 115 ENDIF 116 ! 117 ALLOCATE (vhart(nrxxs,nspin)) 118 ALLOCATE (pseudo_dens_c(nrxxs)) 119 ALLOCATE (red_revc0(nrxxs,nbnd,nkqs)) 120 red_revc0 = (0.0_dp, 0.0_dp) 121 ! 122 IF (gamma_only) THEN 123 ALLOCATE (revc_int(nrxxs,nbnd)) 124 ELSE 125 ALLOCATE(revc_int_c(nrxxs,nbnd,nks)) 126 ALLOCATE(k2q(nks)) 127 k2q=0 128 ENDIF 129 ! 130END SUBROUTINE lr_exx_alloc 131! 132SUBROUTINE lr_exx_dealloc() 133 134 IMPLICIT NONE 135 136 DEALLOCATE(pseudo_dens_c, vhart, red_revc0) 137 ! 138 IF (gamma_only) THEN 139 DEALLOCATE(revc_int) 140 ELSE 141 DEALLOCATE(revc_int_c, k2q) 142 ENDIF 143 ! 144END SUBROUTINE lr_exx_dealloc 145! 146SUBROUTINE lr_exx_sum_int() 147 USE mp_global, ONLY : inter_bgrp_comm 148 USE mp, ONLY : mp_sum 149 150 CALL start_clock('lr_exx_sum') 151 152 CALL mp_sum(revc_int, inter_bgrp_comm) 153 154 CALL stop_clock('lr_exx_sum') 155 156 RETURN 157 158END SUBROUTINE lr_exx_sum_int 159! 160SUBROUTINE lr_exx_apply_revc_int(psi, ibnd, nbnd, ik) 161 !------------------------------------------------------------------ 162 ! 163 ! This routine takes the interaction generated by lr_exx_kernel_int 164 ! and applies it to psi for a specified band, ibnd. 165 ! It handles the conversion to the smooth grid but assumes that 166 ! lr_exx_kernel_int has been called for every band and 167 ! lr_exx_sum_int has been called. 168 !------------------------------------------------------------------ 169 ! 170 171 USE klist, ONLY : ngk 172 173 IMPLICIT NONE 174 175 COMPLEX(DP), INTENT(INOUT) :: psi(:) 176 INTEGER, INTENT(IN) :: ibnd, nbnd, ik 177 178 COMPLEX(DP), ALLOCATABLE :: tempphic(:,:),temppsic(:,:), psi_t(:) 179 INTEGER :: nrxxs, j 180 INTEGER :: npw ! number of plane waves at point k 181 182 COMPLEX(DP) :: fp, fm 183 184 CALL start_clock('lr_exx_apply') 185 186 nrxxs= dfftt%nnr 187 188 IF (gamma_only) THEN 189 ! 190 npw = ngk(1) 191 ! 192 ALLOCATE ( tempphic(dffts%nnr,2),temppsic(dffts%nnr,2),& 193 & psi_t(dffts%nnr) ) 194 tempphic = (0.0_dp, 0.0_dp) 195 temppsic = (0.0_dp, 0.0_dp) 196 psi_t = (0.0_dp, 0.0_dp) 197 198 IF (ibnd < nbnd) THEN 199 ! 200 ! We need to do two bands at a time as per gamma tricks. 201 ! 202 tempphic(1:nrxxs,1)= 0.5d0 * CMPLX( revc_int(1:nrxxs,ibnd),& 203 & revc_int(1:nrxxs,ibnd+1), kind=dp ) 204 ! 205 ! To g-space 206 ! 207 CALL fwfft ('Wave', tempphic(:,1), dfftt) 208 ! 209 ! Now separate the two bands and apply the correct nl mapping 210 ! 211 DO j = 1, npw 212 fp = (tempphic(dfftt%nl(j),1) + & 213 tempphic(dfftt%nlm(j),1))*0.5d0 214 fm = (tempphic(dfftt%nl(j),1) - & 215 tempphic(dfftt%nlm(j),1))*0.5d0 216 temppsic( j, 1) = CMPLX( DBLE(fp), AIMAG(fm),kind=DP) 217 temppsic( j, 2) = CMPLX(AIMAG(fp),- DBLE(fm),kind=DP) 218 ENDDO 219 220 psi_t(dffts%nl(1:npw))= temppsic(1:npw,1)+ (0.0_dp,1.0_dp)& 221 &*temppsic(1:npw,2) 222 psi_t(dffts%nlm(1:npw))= CONJG(temppsic(1:npw,1)- (0.0_dp,1.0_dp)& 223 &*temppsic(1:npw,2)) 224 ELSE 225 tempphic(1:nrxxs,1) = 0.5d0 * CMPLX( revc_int(1:nrxxs,ibnd),& 226 & 0.0_dp, kind=dp ) 227 ! 228 ! To g-space 229 ! 230 CALL fwfft ('Wave', tempphic(:,1), dfftt) 231 ! 232 ! Correct the nl mapping for the two grids. 233 ! 234 temppsic(1:npw,1)=tempphic(dfftt%nl(1:npw),1) 235 psi_t(dffts%nl(1:npw))=temppsic(1:npw,1) 236 psi_t(dffts%nlm(1:npw))=CONJG(temppsic(1:npw,1)) 237 ! 238 ENDIF 239 ! 240 DEALLOCATE ( tempphic, temppsic ) 241 ! 242 CALL invfft ('Wave', psi_t, dffts) 243 ! 244 psi(:) = psi(:) + psi_t(:) 245 ! 246 DEALLOCATE ( psi_t ) 247 ! 248 ELSE 249 ! 250 psi(:) = psi(:) + revc_int_c(:,ibnd,ik) 251 ! 252 ENDIF 253 CALL stop_clock('lr_exx_apply') 254 RETURN 255 256END SUBROUTINE lr_exx_apply_revc_int 257 258 259SUBROUTINE lr_exx_revc0_init(orbital, ik) 260 !------------------------------------------------------------------ 261 ! 262 ! This routine should be called at the start of a TDDFT EXX calculation 263 ! and it take sthe ground state orbitals (passed into 'orbital'), 264 ! interpolates them onto the custom FFT grid used by exx.f90, and stores 265 ! the real-space result in red_revc0(). 266 ! 267 ! For the k-points version also it performs the necessary rotations to 268 ! obtain the q points from their respective k points. 269 !------------------------------------------------------------------ 270 ! 271 272 USE mp_global, ONLY : me_bgrp 273 USE exx_base, ONLY : rir, nkqs, index_sym, index_xk 274 USE scatter_mod, ONLY : gather_grid, scatter_grid 275 USE symm_base, ONLY : sname 276 277 IMPLICIT NONE 278 279 COMPLEX(DP), INTENT(IN) :: orbital(:,:,:) 280 INTEGER, INTENT(IN) :: ik 281 282 INTEGER :: ibnd, nnr_, ikq, isym, nxxs, nrxxs 283 COMPLEX(DP), ALLOCATABLE :: psic_all(:), temppsic_all(:), temppsic(:) 284 285 IF (gamma_only) THEN 286 ! 287 nnr_= dfftt%nnr 288 ! 289 DO ibnd=1,nbnd,2 290 ! 291 CALL invfft_orbital_custom_gamma(orbital(:,:,1), ibnd, nbnd, & 292 npwt, dfftt) 293 red_revc0(1:nnr_,ibnd,1)=psic(1:nnr_) 294 ! 295 ENDDO 296 ELSE 297 nxxs=dffts%nr1x * dffts%nr2x * dffts%nr3x 298 nrxxs=dffts%nnr 299 ! 300 ALLOCATE(temppsic_all(1:nxxs), psic_all(1:nxxs), temppsic(1:nrxxs)) 301 ! 302 DO ibnd=1,nbnd 303 ! 304 CALL invfft_orbital_k ( orbital(:,:,ik), ibnd, nbnd, ik) 305 ! 306 DO ikq=1,nkqs 307 ! 308 IF (index_xk(ikq) .NE. ik) CYCLE 309 ! 310 ! Now rotate k to q 311 isym = ABS(index_sym(ikq) ) 312 IF (index_sym(ikq) > 0) THEN 313 IF (TRIM(sname(index_sym(ikq))) .EQ. "identity" ) & 314 & k2q(ik) = ikq 315 ENDIF 316 ! 317#if defined(__MPI) 318 CALL gather_grid (dffts, psic, psic_all) 319 IF ( me_bgrp == 0 ) & 320 temppsic_all(1:nxxs) = psic_all(rir(1:nxxs, isym)) 321 CALL scatter_grid(dffts, temppsic_all, temppsic) 322#else 323 temppsic(1:nrxxs) = psic(rir(1:nrxxs, isym)) 324#endif 325 IF (index_sym(ikq) < 0 ) & 326 &temppsic(1:nrxxs) = CONJG(temppsic(1:nrxxs)) 327 ! 328 red_revc0(1:nrxxs, ibnd, ikq) = temppsic(:) 329 ! 330 ENDDO 331 ! 332 ENDDO 333 ! 334 DEALLOCATE(temppsic_all, psic_all) 335 ! 336 ENDIF 337 ! 338 RETURN 339 ! 340END SUBROUTINE lr_exx_revc0_init 341! 342SUBROUTINE lr_exx_kernel_noint ( evc, int_vect ) 343 !------------------------------------------------------------------ 344 ! 345 ! This routine computes the change of the self consistent potential 346 ! due to the perturbation. More specifically it combines the K^1d 347 ! and K^2d terms from Ref (1) in the 'non-interacting' (lower SBR 348 ! vectors) case. 349 ! (In the hybrid/exx CASE non-interacting is a misnomer as Hartree 350 ! like terms are applied even to the lower batch). 351 ! 352 ! (1) Rocca, Lu and Galli, J. Chem. Phys., 133, 164109 (2010) 353 !------------------------------------------------------------------ 354 ! 355 356 USE kinds, ONLY : DP 357 USE lr_variables, ONLY : gamma_only, lr_verbosity 358 USE exx_base, ONLY : g2_convolution, nqs, index_sym, & 359 index_xkq, index_xk, rir, nkqs 360 USE exx, ONLY : exxalfa 361 USE symm_base, ONLY : s 362 USE cell_base, ONLY : bg, at 363 USE funct, ONLY : exx_is_active 364 USE io_global, ONLY : stdout 365 USE mp_global, ONLY : inter_bgrp_comm, ibnd_start, ibnd_end,& 366 & me_bgrp 367 USE mp, ONLY : mp_sum 368 USE scatter_mod, ONLY : gather_grid, scatter_grid 369 370 IMPLICIT NONE 371 ! 372 COMPLEX(DP), INTENT(in) :: evc(npwx,nbnd,nks) 373 COMPLEX(DP), INTENT(out) :: int_vect(npwx,nbnd,nks) 374 ! 375 ! 376 REAL(kind=dp), ALLOCATABLE :: revc_int(:,:) 377 COMPLEX(kind=dp), ALLOCATABLE :: revc_int_c(:,:,:) 378 ! Container for the interaction 379 ! 380 INTEGER :: ibnd, ik 381 ! Counters 382 REAL(kind=dp) :: w1, w2, scale 383 ! 384 REAL(kind=DP), ALLOCATABLE :: fac(:) 385 ! Varible to hold to actual interaction (1/|r-r'| or similar in g-space). 386 INTEGER :: nxxs, nrxxs 387 ! 388 ! q-points related stuff 389 INTEGER :: iq, ikq, isym, ikk 390 REAL(DP) :: xk_cryst(3), sxk(3), xkq(3) 391 COMPLEX(DP), ALLOCATABLE :: psic_all(:), temppsic_all(:), temppsic(:) 392 ! 393 INTEGER :: ibnd_start_gamma, ibnd_end_gamma 394 395 LOGICAL :: exst 396 397 398 ! Offset ibnd_start and ibnd_ed to avoid conflicts with gamma_tricks for 399 ! even values 400 ! 401 CALL start_clock('lr_exx_noint') 402 ! 403 IF (lr_verbosity > 5 ) WRITE(stdout,'("<lr_exx_kernel_noint>")') 404 ! 405 ! Setup the variables that describe the FFT grid in use. 406 IF(gamma_only) THEN 407 ALLOCATE( fac(dfftt%ngm) ) 408 nrxxs= dfftt%nnr 409 ELSE 410 ALLOCATE( fac(ngm) ) 411 nxxs=dffts%nr1x * dffts%nr2x * dffts%nr3x 412 nrxxs = dffts%nnr 413 ALLOCATE(temppsic_all(1:nxxs), psic_all(1:nxxs), temppsic(1:nrxxs)) 414 ENDIF 415 ! 416 fac=0.d0 417 ! 418 IF (exx_is_active()) THEN 419 scale = exxalfa 420 ELSE 421 scale=1.d0 422 ENDIF 423 ! 424 ! 425 IF( gamma_only ) THEN 426 ! 427 ! Put the appropriate interaction in fac(). Note g2_convolution respects 428 ! the choice of divergence treatment etc set in the initial PWscf run. 429 ! 430 CALL g2_convolution(dfftt%ngm, gt, xk(:,1), xk(:,1), fac) 431 ! 432 ALLOCATE(revc_int(nrxxs,nbnd)) 433 ! 434 revc_int=0.d0 435 int_vect=(0.d0,0.d0) 436 ! 437 ! 438 ! If ibnd_start is even offset it so bands aren't skipped/double counted 439 ! when we use gamma_tricks. 440 ! 441 ibnd_start_gamma=ibnd_start 442 IF (MOD(ibnd_start, 2)==0) ibnd_start_gamma = ibnd_start+1 443 ibnd_end_gamma = MAX(ibnd_end, ibnd_start_gamma) 444 ! 445 DO ibnd=ibnd_start_gamma,ibnd_end_gamma,2 446 ! 447 CALL invfft_orbital_custom_gamma(evc(:,:,1), ibnd, nbnd, & 448 npwt, dfftt) 449 ! 450 w1=wg(ibnd,1)/omega 451 ! 452 IF (ibnd<nbnd) THEN 453 w2=wg(ibnd+1,1)/omega 454 ELSE 455 w2=0.0d0 456 ENDIF 457 ! Update the container with the actual interaction for this band(s). 458 revc_int(1:dfftt%nnr,:)= revc_int(1:dfftt%nnr,:) & 459 & -1.d0 * scale * k1d_term_gamma(w1,w2,psic,fac,ibnd,evc(:,:,1)) & 460 & +1.d0 * scale * k2d_vexx_term_gamma(w1,w2,psic,fac,ibnd,.false.) 461 ! 462 ENDDO 463 ! 464 CALL mp_sum( revc_int(:,:), inter_bgrp_comm ) 465 ! 466 ! Put the constructed interaction into the psic workspace and 467 ! then on to int_vec. 468 ! 469 DO ibnd=ibnd_start_gamma,ibnd_end_gamma,2!1,nbnd,2 470 ! 471 psic=(0.0_dp, 0.0_dp) 472 ! 473 IF (ibnd<nbnd) psic(1:nrxxs)=CMPLX(revc_int(1:nrxxs,ibnd),& 474 & revc_int(1:nrxxs,ibnd+1),dp) 475 IF (ibnd==nbnd) psic(1:nrxxs)=CMPLX(revc_int(1:nrxxs,ibnd)& 476 &,0.d0,dp) 477 ! 478 CALL fwfft_orbital_custom_gamma (int_vect(:,:,1), ibnd, nbnd, & 479 npwt, dfftt) 480 ! 481 ENDDO 482 ! 483 CALL mp_sum( int_vect(:,:,1), inter_bgrp_comm ) 484 ! 485 DEALLOCATE(revc_int, fac) 486 ! 487 int_vect=int_vect*0.5d0 488 ! 489 ELSE 490 ! 491 ! In the case of k-points... 492 ! 493 ALLOCATE(revc_int_c(dffts%nnr,nbnd,nks)) 494 ! 495 int_vect=(0.d0,0.d0) 496 revc_int_c=(0.d0,0.d0) 497 ! 498 ! Loop over all k-points 499 ! 500 DO ikk=1,nks 501 ! 502 DO ibnd=1,nbnd,1 503 ! 504 ! 505 CALL invfft_orbital_k (evc(:,:,ikk), ibnd, nbnd, ikk) 506 ! 507#if defined(__MPI) 508 psic_all=(0.d0,0.d0) 509 CALL gather_grid(dffts, psic, psic_all) 510#endif 511 ! 512 ! Now psic_all has the collected band ibnd at kpoint ikk. 513 ! We now search to see which q points are identical via symmetry 514 ! 515 DO ikq=1,nkqs 516 ! 517 IF (ikk /= index_xk(ikq)) CYCLE 518 ! 519 ! Now rotate k to q, scatter and put in temppsic 520 ! 521 isym = ABS(index_sym(ikq)) 522#if defined(__MPI) 523 IF ( me_bgrp == 0 ) & 524 temppsic_all(1:nxxs) = psic_all(rir(1:nxxs, isym)) 525 CALL scatter_grid(dffts, temppsic_all, temppsic) 526#else 527 528 temppsic(1:nrxxs) = psic(rir(1:nrxxs, isym)) 529#endif 530 IF (index_sym(ikq) < 0 ) & 531 &temppsic(1:nrxxs) = CONJG(temppsic(1:nrxxs)) 532 ! 533 ! We now search over all k+q to find which correspond 534 ! to the k+q point in temppsic 535 ! 536 DO ik=1,nks 537 DO iq=1,nqs 538 ! 539 ! 540 IF (ikq /= index_xkq(ik,iq)) CYCLE 541 ! 542 xk_cryst(:) = at(1,:)*xk(1,ikk) + at(2,:)*xk(2,ikk) & 543 + at(3,:)*xk(3,ikk) 544 ! 545 IF (index_sym(ikq) < 0 ) xk_cryst = - xk_cryst 546 ! 547 sxk(:) = s(:,1,isym)*xk_cryst(1) + & 548 s(:,2,isym)*xk_cryst(2) + & 549 s(:,3,isym)*xk_cryst(3) 550 xkq(:) = bg(:,1)*sxk(1) + bg(:,2)*sxk(2) + bg(:,3)*sxk(3) 551 ! 552 CALL g2_convolution(ngm, g, xk(:,ik), xkq, fac) 553 ! 554 ! 555 w1 = wg(ibnd,ik)/(wk(ik) * nqs) 556 ! 557 ! Update the container with the actual interaction for this band(s). 558 revc_int_c(:,:, ik)= revc_int_c(:,:, ik) & 559 &-1.d0 * scale * & 560 &k1d_term_k(w1,temppsic,fac,ibnd,ik,ikq) & 561 &+1.d0 * scale * & 562 &k2d_term_k(w1,temppsic,fac,ibnd,ik,ikq) 563 END DO 564 END DO 565 ENDDO 566 ENDDO 567 ENDDO 568 ! 569 ! 570 CALL mp_sum( revc_int_c, inter_bgrp_comm ) 571 ! 572 DO ik=1,nks 573 ! 574 DO ibnd=1,nbnd,1 575 ! 576 psic(:)=(0.0_dp,0.0_dp) 577 ! 578 ! Put the constructed interaction into the psic workspace 579 ! and then on to int_vec. 580 ! 581 psic(:)=revc_int_c(:,ibnd,ik) 582 ! 583 CALL fwfft_orbital_k (int_vect(:,:,ik), ibnd, nbnd, ik) 584 ! 585 ENDDO 586 ! 587 ENDDO 588 ! 589 DEALLOCATE(revc_int_c, fac) 590 DEALLOCATE(temppsic_all, psic_all, temppsic) 591 ! 592 ENDIF 593 ! 594 ! 595 CALL stop_clock('lr_exx_noint') 596 RETURN 597 ! 598END SUBROUTINE lr_exx_kernel_noint 599! 600SUBROUTINE lr_exx_kernel_int ( orbital, ibnd, nbnd, ikk ) 601 !------------------------------------------------------------------ 602 ! 603 ! This routine computes the change of the self consistent potential 604 ! due to the perturbation. More specifically it combines the K^1d 605 ! and K^2d terms from Ref (1) in the 'interacting' (upper SBR 606 ! vectors) case. 607 ! 608 ! (1) Rocca, Lu and Galli, J. Chem. Phys., 133, 164109 (2010) 609 ! 610 ! Unlike lr_exx_kernel_noint above this routine should be called 611 ! a band at a time. Also it does sum the interaction over the band 612 ! group (obviously) and therefore does not place the full interaction 613 ! in the revc_int container. See lr_bse_sum and lr_bse_apply_revc_int. 614 !------------------------------------------------------------------ 615 ! 616 617 USE kinds, ONLY : DP 618 USE klist, ONLY : xk 619 USE io_global, ONLY : stdout 620 USE exx, ONLY : exxalfa 621 USE exx_base, ONLY : g2_convolution, nqs, index_sym, & 622 index_xkq, index_xk, rir, nkqs 623 USE symm_base, ONLY : s 624 USE cell_base, ONLY : bg, at 625 USE funct, ONLY : exx_is_active 626 USE mp_global, ONLY : me_bgrp 627 USE scatter_mod, ONLY : gather_grid, scatter_grid 628 USE lr_variables, ONLY : ltammd 629 630 IMPLICIT NONE 631 632 COMPLEX(DP), INTENT(in) :: orbital(:,:) 633 INTEGER, INTENT(in) :: ibnd ! Current band under consideration 634 INTEGER, INTENT(in) :: ikk ! K-point of the current band 635 INTEGER, INTENT(in) :: nbnd ! Total number of bands 636 637 638 REAL(kind=DP), ALLOCATABLE :: fac(:) 639 ! Varible to hold to actual interaction (1/|r-r'| or similar in g-space). 640 641 REAL(kind=DP) :: scale 642 ! Variable to hold exxalfa 643 ! (or otherwise if interaction is used whilst exx_is_active == .false.) 644 645 INTEGER :: nxxs, nrxxs 646 REAL(kind=dp) :: w1, w2 647 ! Weights for the band(s) 648 649 ! q-points related stuff 650 INTEGER :: iq, ikq, isym, ik 651 REAL(DP) :: xk_cryst(3), sxk(3), xkq(3) 652 COMPLEX(DP), ALLOCATABLE :: psic_all(:), temppsic_all(:), temppsic(:) 653 654 CALL start_clock('lr_exx_int') 655 IF (lr_verbosity > 5 ) WRITE(stdout,'("<lr_exx_kernel_int>")') 656 657 IF(gamma_only) THEN 658 ALLOCATE( fac(dfftt%ngm) ) 659 nrxxs= dfftt%nnr 660 ELSE 661 ALLOCATE( fac(ngm) ) 662 nxxs=dffts%nr1x * dffts%nr2x * dffts%nr3x 663 nrxxs = dffts%nnr 664 ALLOCATE(temppsic_all(1:nxxs), psic_all(1:nxxs), temppsic(1:nrxxs)) 665 ENDIF 666 ! 667 fac=0.d0 668 ! 669 IF (exx_is_active()) THEN 670 scale = exxalfa 671 ELSE 672 scale=1.d0 673 ENDIF 674 ! 675 IF( gamma_only ) THEN 676 ! 677 CALL invfft_orbital_custom_gamma( orbital, ibnd, nbnd, npwt, dfftt ) 678 ! 679 w1=wg(ibnd,1)/omega 680 ! 681 IF (ibnd<nbnd) THEN 682 w2=wg(ibnd+1,1)/omega 683 ELSE 684 w2=0.0d0 685 ENDIF 686 ! 687 CALL g2_convolution(dfftt%ngm, gt,xk(:,1), xk(:,1), fac) 688 ! 689 IF (.NOT.ltammd) THEN 690 revc_int(1:dfftt%nnr,:)= revc_int(1:dfftt%nnr,:)& 691 & -1.d0 * scale * k1d_term_gamma(w1,w2,psic,fac,ibnd,orbital) & 692 & -1.d0 * scale * k2d_vexx_term_gamma(w1,w2,psic,fac,ibnd,.true.) 693 ELSE 694 ! 695 ! Slightly different interaction in the Tamm--Dancoff case. 696 ! Note the factor of 2 to account for the fact lr_apply_liovillian 697 ! scales the *whole interaction term* by a factor of 0.5 in the TD 698 ! case. 699 ! 700 revc_int(1:dfftt%nnr,:)= revc_int(1:dfftt%nnr,:)& 701 & -2.d0 * scale * k1d_term_gamma(w1,w2,psic,fac,ibnd,orbital) 702 ENDIF 703 ! 704 ELSE 705 ! 706 CALL invfft_orbital_k (orbital(:,:), ibnd, nbnd, ikk) 707 ! 708#if defined(__MPI) 709 CALL gather_grid(dffts, psic, psic_all) 710#endif 711 ! 712 ! Now psic_all has the collected band ibnd at kpoint ikk. 713 ! We now search to see which q points are identical via symmetry 714 ! 715 DO ikq=1,nkqs 716 ! 717 IF (ikk /= index_xk(ikq)) CYCLE 718 ! 719 ! Now rotate k to q, scatter and put in temppsic 720 ! 721 isym = ABS(index_sym(ikq)) 722#if defined(__MPI) 723 IF ( me_bgrp == 0 ) & 724 temppsic_all(1:nxxs) = psic_all(rir(1:nxxs, isym)) 725 CALL scatter_grid(dffts, temppsic_all, temppsic) 726#else 727 728 temppsic(1:nrxxs) = psic(rir(1:nrxxs, isym)) 729#endif 730 IF (index_sym(ikq) < 0 ) & 731 &temppsic(1:nrxxs) = CONJG(temppsic(1:nrxxs)) 732 ! 733 ! We now search over all k+q to find which correspond 734 ! to the k+q point in temppsic 735 ! 736 DO ik=1,nks 737 DO iq=1,nqs 738 ! 739 ! 740 IF (ikq /= index_xkq(ik,iq)) CYCLE 741 ! 742 xk_cryst(:) = at(1,:)*xk(1,ikk) + at(2,:)*xk(2,ikk) & 743 + at(3,:)*xk(3,ikk) 744 ! 745 IF (index_sym(ikq) < 0 ) xk_cryst = - xk_cryst 746 ! 747 sxk(:) = s(:,1,isym)*xk_cryst(1) + & 748 s(:,2,isym)*xk_cryst(2) + & 749 s(:,3,isym)*xk_cryst(3) 750 xkq(:) = bg(:,1)*sxk(1) + bg(:,2)*sxk(2) + bg(:,3)*sxk(3) 751 ! 752 CALL g2_convolution(ngm, g, xk(:,ik), xkq, fac) 753 ! 754 ! 755 w1=wg(ibnd,ik)/(wk(ik) * nqs) 756 ! 757 IF(.NOT.ltammd) THEN 758 revc_int_c(:,:,ik)= revc_int_c(:,:,ik) & 759 & -1.d0 * scale * & 760 &k1d_term_k(w1,temppsic,fac,ibnd,ik,ikq) & 761 & -1.d0 * scale * & 762 &k2d_term_k(w1,temppsic,fac,ibnd,ik,ikq) 763 ELSE 764 ! 765 ! Slightly different interaction in the Tamm--Dancoff case. 766 ! Note the factor of 2 to account for the fact 767 ! lr_apply_liovillian scales the *whole interaction term* 768 ! by a factor of 0.5 in the TD case. 769 ! 770 revc_int_c(:,:,ik)= revc_int_c(:,:,ik) & 771 & -2.d0 * scale * & 772 &k1d_term_k(w1,temppsic,fac,ibnd,ik,ikq) 773 ENDIF 774 ! 775 ENDDO 776 ENDDO 777 ENDDO 778 ! 779 DEALLOCATE(temppsic_all, psic_all, temppsic) 780 ! 781 ENDIF 782 IF (lr_verbosity > 5 ) WRITE(stdout,'("<lr_exx_kernel_int> exit")') 783 CALL stop_clock('lr_exx_int') 784 785 ! Note at this point revc_int(_c) may be incomplete 786 ! (depending if all the bands have been processed) and still needs to be summed 787 ! and converted to the smooth FFT grid. 788 789 RETURN 790 791END SUBROUTINE lr_exx_kernel_int 792! 793FUNCTION k1d_term_gamma(w1, w2, psi, fac_in, ibnd, orbital) RESULT (psi_int) 794 !------------------------------------------------------------------ 795 ! 796 ! This routine computes the K^1d term, Eq (21) from Eq Ref (1). 797 ! As the integral in this term remains the same throughout the 798 ! chain it can also be calculated once for each combination of 799 ! bands and stored in k1d_pot(). 800 ! 801 ! psi contains two bands of |a> with w1, w2 the associated weights 802 ! fac_in contains the interaction W(G) and ibnd the band index n'. 803 ! 804 ! in gamma, bands are real, so is not necessary computes vhart 805 ! integrals for each couple of bands (integral in eq(21), 806 ! because v, v' is the same of v', v, so only nbnd(nbnd+1)/2 807 ! couple are computed 808 ! 809 ! (1) Rocca, Lu and Galli, J. Chem. Phys., 133, 164109 (2010) 810 ! 811 !------------------------------------------------------------------ 812 ! 813 814 IMPLICIT NONE 815 816 COMPLEX(KIND=DP),DIMENSION(:), INTENT(IN) :: psi 817 REAL(KIND=DP),DIMENSION(:), INTENT(IN) :: fac_in 818 REAL(kind=DP), INTENT(IN) :: w1, w2 819 REAL(KIND=DP), ALLOCATABLE :: psi_int(:,:) 820 INTEGER, INTENT(IN) :: ibnd 821 COMPLEX(DP), INTENT(IN) :: orbital(:,:) 822 COMPLEX(DP) :: psitemp(dfftt%nnr) 823 ! 824 ! Workspaces 825 ! 826 INTEGER :: ibnd2, is, npw_, ngm_, nnr_ 827 INTEGER :: nrec 828 ! 829 npw_=npwt 830 ngm_=dfftt%ngm 831 nnr_=dfftt%nnr 832 ! 833 ALLOCATE(psi_int(nnr_, nbnd)) 834 psi_int = 0.d0 835 ! 836 ! 837 IF (ibnd < nbnd) then 838 ! 839 vhart(:,:) = (0.d0,0.d0) 840 pseudo_dens_c(:) = (0.d0,0.d0) 841 ! 842 ! calculate vhart for couples ibnd,ibnd and ibnd+1,ibnd 843 ! 844 pseudo_dens_c(1:nnr_) = CMPLX( w1*DBLE(red_revc0(1:nnr_,ibnd,1)) *& 845 & DBLE(red_revc0(1:nnr_,ibnd,1)), & 846 & w2*AIMAG(red_revc0(1:nnr_,ibnd, 1)) *& 847 & DBLE(red_revc0(1:nnr_,ibnd,1)), kind=DP ) 848 ! 849 CALL fwfft ('Rho', pseudo_dens_c, dfftt) 850 ! 851 DO is = 1, nspin 852 ! 853 vhart(dfftt%nl(1:ngm_),is) =& 854 & pseudo_dens_c(dfftt%nl(1:ngm_)) *& 855 & fac_in(1:ngm_) 856 IF (gamma_only) vhart(dfftt%nlm(1:ngm_),is) = & 857 & pseudo_dens_c(dfftt%nlm(1:ngm_)) *& 858 & fac_in(1:ngm_) 859 ! 860 ! and transformed back to real space 861 ! 862 CALL invfft ('Rho', vhart (:, is), dfftt) 863 ENDDO 864 ! 865 ! Finally return the interaction psi_int for terms: 866 ! ibnd, ibnd, ibnd 867 ! ibnd+1, ibnd+1, ibnd 868 ! ibnd, ibnd, ibnd+1 869 ! 870 psi_int(1:nnr_,ibnd) = psi_int(1:nnr_,ibnd) & 871 & + DBLE(vhart(1:nnr_, 1)) * DBLE(psi(1:nnr_)) & 872 & + AIMAG(vhart(1:nnr_,1)) * AIMAG(psi(1:nnr_)) 873 ! 874 psi_int(1:nnr_,ibnd+1) = psi_int(1:nnr_,ibnd+1) & 875 & + AIMAG(vhart(1:nnr_,1)) * DBLE(psi(1:nnr_)) 876 ! 877 ! calculate vhart for couple ibnd+1,ibnd+1 878 ! 879 vhart(:,:) = (0.d0,0.d0) 880 pseudo_dens_c(:) = (0.d0,0.d0) 881 ! 882 pseudo_dens_c(1:nnr_) = CMPLX( w2*AIMAG(red_revc0(1:nnr_,ibnd,1)) *& 883 & AIMAG(red_revc0(1:nnr_,ibnd,1)), 0.0d0, kind=DP ) 884 ! 885 CALL fwfft ('Rho', pseudo_dens_c, dfftt) 886 ! 887 DO is = 1, nspin 888 ! 889 vhart(dfftt%nl(1:ngm_),is) =& 890 & pseudo_dens_c(dfftt%nl(1:ngm_)) *& 891 & fac_in(1:ngm_) 892 IF (gamma_only) vhart(dfftt%nlm(1:ngm_),is) = & 893 & pseudo_dens_c(dfftt%nlm(1:ngm_)) *& 894 & fac_in(1:ngm_) 895 ! 896 ! and transformed back to real space 897 ! 898 CALL invfft ('Rho', vhart (:, is), dfftt) 899 ENDDO 900 ! 901 ! Finally return the interaction psi_int for term: 902 ! ibnd+1, ibnd+1, ibnd+1 903 ! 904 psi_int(1:nnr_,ibnd+1) = psi_int(1:nnr_,ibnd+1) & 905 & + DBLE(vhart(1:nnr_,1)) * AIMAG(psi(1:nnr_)) 906 ! 907 ! 908 ! start second loop over bands 909 ! 910 DO ibnd2=1,ibnd-1,1 911 ! 912 ! calculate vhart for couples ibnd,ibnd2 and ibnd+1,ibnd2 913 ! 914 vhart(:,:) = (0.d0,0.d0) 915 pseudo_dens_c(:) = (0.d0,0.d0) 916 ! 917 ! The following code is to check if the value of ibnd2 is even or odd 918 ! and therefore whether the REAL or IMAGINARY part of red_revc0 is to be 919 ! used. This is because red_revc0 is stored using gamma_tricks. 920 ! 921 IF (MOD(ibnd2,2)==1) THEN 922 pseudo_dens_c(1:nnr_) = CMPLX( w1*DBLE(red_revc0(1:nnr_,ibnd,1)) *& 923 & DBLE(red_revc0(1:nnr_,ibnd2,1)), & 924 & w2*AIMAG(red_revc0(1:nnr_,ibnd, 1)) *& 925 & DBLE(red_revc0(1:nnr_,ibnd2,1)), kind=DP ) 926 ELSE 927 pseudo_dens_c(1:nnr_) = CMPLX( w1*DBLE(red_revc0(1:nnr_,ibnd,1)) *& 928 &AIMAG(red_revc0(1:nnr_,ibnd2-1,1)),& 929 &w2*AIMAG(red_revc0(1:nnr_,ibnd,1)) *& 930 &AIMAG(red_revc0(1:nnr_,ibnd2-1,1)), kind=DP ) 931 ENDIF 932 ! 933 CALL fwfft ('Rho', pseudo_dens_c, dfftt) 934 ! 935 ! hartree contribution is computed in reciprocal space 936 ! 937 DO is = 1, nspin 938 ! 939 vhart(dfftt%nl(1:ngm_),is) =& 940 & pseudo_dens_c(dfftt%nl(1:ngm_)) *& 941 & fac_in(1:ngm_) 942 IF (gamma_only) vhart(dfftt%nlm(1:ngm_),is) = & 943 & pseudo_dens_c(dfftt%nlm(1:ngm_)) *& 944 & fac_in(1:ngm_) 945 ! 946 ! and transformed back to real space 947 ! 948 CALL invfft ('Rho', vhart (:, is), dfftt) 949 ! 950 ENDDO 951 ! 952 ! Finally return the interaction psi_int for terms: 953 ! ibnd, ibnd, ibnd2 954 ! ibnd+1, ibnd+1, ibnd2 955 ! ibnd2, ibnd2, ibnd 956 ! ibnd2, ibnd2, ibnd+1 957 ! 958 psi_int(1:nnr_,ibnd2) = psi_int(1:nnr_,ibnd2) & 959 & + DBLE(vhart(1:nnr_, 1)) * DBLE(psi(1:nnr_)) & 960 & + AIMAG(vhart(1:nnr_,1)) * AIMAG(psi(1:nnr_)) 961 ! 962 CALL invfft_orbital_ibnd2_gamma(orbital(:,:), psitemp, ibnd2, npw_, dfftt) 963 ! 964 ! 965 psi_int(1:nnr_,ibnd) = psi_int(1:nnr_,ibnd) & 966 & + DBLE(vhart(1:nnr_, 1)) * DBLE(psitemp(1:nnr_)) 967 ! 968 psi_int(1:nnr_,ibnd+1) = psi_int(1:nnr_,ibnd+1) & 969 & + AIMAG(vhart(1:nnr_, 1)) * DBLE(psitemp(1:nnr_)) 970 ! 971 ENDDO 972 ! 973 ELSE 974 ! 975 ! calculate vhart for couple ibnd,ibnd 976 ! 977 vhart(:,:) = (0.d0,0.d0) 978 pseudo_dens_c(:) = (0.d0,0.d0) 979 ! 980 pseudo_dens_c(1:nnr_) = CMPLX( w1*DBLE(red_revc0(1:nnr_,ibnd,1)) *& 981 & DBLE(red_revc0(1:nnr_,ibnd,1)), 0.0d0, kind=DP ) 982 ! 983 CALL fwfft ('Rho', pseudo_dens_c, dfftt) 984 ! 985 DO is = 1, nspin 986 ! 987 vhart(dfftt%nl(1:ngm_),is) =& 988 & pseudo_dens_c(dfftt%nl(1:ngm_)) *& 989 & fac_in(1:ngm_) 990 IF (gamma_only) vhart(dfftt%nlm(1:ngm_),is) = & 991 & pseudo_dens_c(dfftt%nlm(1:ngm_)) *& 992 & fac_in(1:ngm_) 993 ! 994 ! and transformed back to real space 995 ! 996 CALL invfft ('Rho', vhart (:, is), dfftt) 997 ENDDO 998 ! 999 ! Finally return the interaction psi_int for term: 1000 ! ibnd, ibnd, ibnd 1001 ! 1002 psi_int(1:nnr_,ibnd) = psi_int(1:nnr_,ibnd) & 1003 & + DBLE(vhart(1:nnr_,1)) * DBLE(psi(1:nnr_)) 1004 ! 1005 ! start second loop over bands 1006 ! 1007 DO ibnd2=1,ibnd-1,1 1008 ! 1009 ! 1010 ! Set up the pseudo density for this Hartree like interaction. 1011 ! calculate vhart for couple ibnd,ibnd2 1012 ! 1013 vhart(:,:) = (0.d0,0.d0) 1014 pseudo_dens_c(:) = (0.d0,0.d0) 1015 ! 1016 ! The following code is to check if the value of ibnd2 is even or odd 1017 ! and therefore whether the REAL or IMAGINARY part of red_revc0 is to be 1018 ! used. This is because red_revc0 is stored using gamma_tricks. 1019 ! 1020 IF (MOD(ibnd2,2)==1) THEN 1021 pseudo_dens_c(1:nnr_) = CMPLX( w1*DBLE(red_revc0(1:nnr_,ibnd,1)) *& 1022 & DBLE(red_revc0(1:nnr_,ibnd2,1)), 0.0d0, kind=DP ) 1023 ELSE 1024 pseudo_dens_c(1:nnr_) = CMPLX( w1*DBLE(red_revc0(1:nnr_,ibnd,1)) *& 1025 & AIMAG(red_revc0(1:nnr_,ibnd2-1,1)), 0.0d0, kind=DP ) 1026 ENDIF 1027 ! 1028 CALL fwfft ('Rho', pseudo_dens_c, dfftt) 1029 ! 1030 ! hartree contribution is computed in reciprocal space 1031 ! 1032 DO is = 1, nspin 1033 ! 1034 vhart(dfftt%nl(1:ngm_),is) =& 1035 & pseudo_dens_c(dfftt%nl(1:ngm_)) *& 1036 & fac_in(1:ngm_) 1037 IF (gamma_only) vhart(dfftt%nlm(1:ngm_),is) = & 1038 & pseudo_dens_c(dfftt%nlm(1:ngm_)) *& 1039 & fac_in(1:ngm_) 1040 ! 1041 ! and transformed back to real space 1042 ! 1043 CALL invfft ('Rho', vhart (:, is), dfftt) 1044 ! 1045 ENDDO 1046 ! 1047 ! Finally return the interaction psi_int for terms: 1048 ! ibnd, ibnd, ibnd2 1049 ! ibn2, ibnd2, ibnd 1050 ! 1051 psi_int(1:nnr_,ibnd2) = psi_int(1:nnr_,ibnd2) & 1052 & + DBLE(vhart(1:nnr_, 1)) * DBLE(psi(1:nnr_)) 1053 ! 1054 CALL invfft_orbital_ibnd2_gamma(orbital(:,:), psitemp, ibnd2, npw_, dfftt) 1055 ! 1056 ! 1057 psi_int(1:nnr_,ibnd) = psi_int(1:nnr_,ibnd) & 1058 & + DBLE(vhart(1:nnr_, 1)) * DBLE(psitemp(1:nnr_)) 1059 ! 1060 ENDDO 1061 ! 1062 ENDIF 1063 ! 1064 RETURN 1065 ! 1066END FUNCTION k1d_term_gamma 1067! 1068FUNCTION k1d_term_k(w1, psi, fac_in, ibnd, ik,ikq) RESULT (psi_int) 1069 !------------------------------------------------------------------ 1070 ! 1071 ! This routine computes the K^1d term, Eq (21) from Eq Ref (1). 1072 ! As the integral in this term remains the same throughout the 1073 ! chain it can also be calculated once for each combination of 1074 ! bands and stored in k1d_pot(). 1075 ! 1076 ! psi contains two bands of |a> with w1, w2 the associated weights 1077 ! fac_in contains the interaction W(G) and ibnd the band index n'. 1078 ! 1079 ! (1) Rocca, Lu and Galli, J. Chem. Phys., 133, 164109 (2010) 1080 !------------------------------------------------------------------ 1081 ! 1082 1083 IMPLICIT NONE 1084 1085 COMPLEX(KIND=DP),DIMENSION(:), INTENT(IN) :: psi 1086 REAL(KIND=DP),DIMENSION(:), INTENT(IN) :: fac_in 1087 REAL(kind=DP), INTENT(IN) :: w1 1088 COMPLEX(KIND=DP), ALLOCATABLE :: psi_int(:,:) 1089 INTEGER, INTENT(IN) :: ibnd, ik,ikq 1090 ! 1091 ! Workspaces 1092 ! 1093 INTEGER :: ibnd2, is 1094 1095 ALLOCATE(psi_int(dffts%nnr, nbnd)) 1096 psi_int = (0.d0,0.d0) 1097 ! 1098 DO ibnd2=1,nbnd,1 1099 ! 1100 ! 1101 ! Set up the pseudo density for this Hartree like interaction. 1102 ! 1103 vhart(:,:) = (0.d0,0.d0) 1104 pseudo_dens_c(:) = (0.d0,0.d0) 1105 ! 1106 pseudo_dens_c(:) = CONJG(red_revc0(:,ibnd,ikq))*& 1107 &red_revc0(:,ibnd2,k2q(ik))/omega 1108 ! 1109 CALL fwfft ('Rho', pseudo_dens_c, dffts) 1110 ! 1111 ! hartree contribution is computed in reciprocal space 1112 ! 1113 DO is = 1, nspin 1114 ! 1115 vhart(dffts%nl(1:ngm),is) = w1*pseudo_dens_c(dffts%nl(1:ngm)) *& 1116 & fac_in(1:ngm) 1117 ! 1118 ! and transformed back to real space 1119 ! 1120 CALL invfft ('Rho', vhart (:, is), dffts) 1121 ! 1122 ENDDO 1123 ! 1124 ! Finally return the interaction 1125 ! 1126 psi_int(:,ibnd2) = psi_int(:,ibnd2) + vhart(:,1) * psi(:) 1127 ! 1128 ENDDO 1129 ! 1130END FUNCTION k1d_term_k 1131! 1132FUNCTION k2d_vexx_term_gamma(w1, w2, psi, fac_in, ibnd, interaction) RESULT (psi_int) 1133 !------------------------------------------------------------------ 1134 ! 1135 ! This routine computes the K^2d term, Eq (22) from Eq Ref (1). 1136 ! psi contains two bands of |b> with w1, w2 the associated weights 1137 ! fac_in contains the interaction W(G) and ibnd the band index n'. 1138 ! 1139 ! Also computes the Vexx term, Eq (15) from Eq Ref (2). 1140 ! So in gamma calculations, h_psi soubrite doesn't call anymore 1141 ! vexx routine 1142 ! 1143 ! (1) Rocca, Lu and Galli, J. Chem. Phys., 133, 164109 (2010) 1144 ! (2) Ge, Binnie, Rocca, Gebauer, Baroni, Comp. Phys. Comm., 1145 ! 185, 2080 (2014) 1146 !------------------------------------------------------------------ 1147 ! 1148 1149 IMPLICIT NONE 1150 1151 COMPLEX(KIND=DP),DIMENSION(:), INTENT(IN) :: psi 1152 REAL(KIND=DP),DIMENSION(:), INTENT(IN) :: fac_in 1153 REAL(kind=DP), INTENT(IN) :: w1, w2 1154 INTEGER, INTENT(IN) :: ibnd 1155 REAL(KIND=DP), ALLOCATABLE :: psi_int(:,:) 1156 LOGICAL, INTENT(IN) :: interaction 1157 ! 1158 ! Workspaces 1159 ! 1160 INTEGER :: ibnd2, is, npw_,ngm_, nnr_ 1161 ! 1162 nnr_ = dfftt%nnr 1163 ngm_ = dfftt%ngm 1164 npw_ = npwt 1165 ! 1166 ALLOCATE(psi_int(nnr_, nbnd)) 1167 psi_int = 0.d0 1168 ! 1169 ! 1170 DO ibnd2=1,nbnd,1 1171 ! 1172 ! Set up the pseudo density for this Hartree like interaction. 1173 ! 1174 vhart(:,:) = (0.d0,0.d0) 1175 pseudo_dens_c(:) = (0.d0,0.d0) 1176 ! 1177 ! The following code is to check if the value of ibnd2 is even or odd 1178 ! and therefore whether the REAL or IMAGINARY part of red_revc0 is to be 1179 ! used. This is because red_revc0 is stored using gamma_tricks. 1180 ! 1181 IF (MOD(ibnd2,2)==1) THEN 1182 pseudo_dens_c(1:nnr_) = CMPLX( & 1183 & w1*DBLE(psi(1:nnr_))*DBLE(red_revc0(1:nnr_,ibnd2,1)),& 1184 & w2*AIMAG(psi(1:nnr_))*DBLE(red_revc0(1:nnr_,ibnd2,1)), kind=DP ) 1185 ! 1186 CALL fwfft ('Rho', pseudo_dens_c, dfftt) 1187 ! 1188 ! hartree contribution is computed in reciprocal space 1189 ! 1190 DO is = 1, nspin 1191 ! 1192 vhart(dfftt%nl(1:ngm_),is) = pseudo_dens_c(dfftt%nl(1:ngm_)) * fac_in(1:ngm_) 1193 IF (gamma_only) vhart(dfftt%nlm(1:ngm_),is) = & 1194 & pseudo_dens_c(dfftt%nlm(1:ngm_)) *& 1195 & fac_in(1:ngm_) 1196 ! 1197 ! and transformed back to real space 1198 ! 1199 CALL invfft ('Rho', vhart (:, is), dfftt) 1200 ! 1201 ENDDO 1202 ! 1203 ! 1204 ! Finally return the interaction 1205 ! 1206 psi_int(1:nnr_,ibnd2) = psi_int(1:nnr_,ibnd2) & 1207 & +DBLE(vhart(1:nnr_,1)) * DBLE(red_revc0(1:nnr_,ibnd,1)) & 1208 & +AIMAG(vhart(1:nnr_,1)) * AIMAG(red_revc0(1:nnr_,ibnd,1)) 1209 ! 1210 IF (interaction) then 1211 psi_int(1:nnr_,ibnd) = psi_int(1:nnr_,ibnd) & 1212 & +DBLE(vhart(1:nnr_,1)) * DBLE(red_revc0(1:nnr_,ibnd2,1)) 1213 IF (ibnd < nbnd) then 1214 psi_int(1:nnr_,ibnd+1) = psi_int(1:nnr_,ibnd+1) & 1215 & +AIMAG(vhart(1:nnr_,1)) * DBLE(red_revc0(1:nnr_,ibnd2,1)) 1216 ENDIF 1217 ELSE 1218 psi_int(1:nnr_,ibnd) = psi_int(1:nnr_,ibnd) & 1219 & -DBLE(vhart(1:nnr_,1)) * DBLE(red_revc0(1:nnr_,ibnd2,1)) 1220 IF (ibnd < nbnd) then 1221 psi_int(1:nnr_,ibnd+1) = psi_int(1:nnr_,ibnd+1) & 1222 & -AIMAG(vhart(1:nnr_,1)) * DBLE(red_revc0(1:nnr_,ibnd2,1)) 1223 ENDIF 1224 ENDIF 1225 ! 1226 ELSE 1227 pseudo_dens_c(1:nnr_) = CMPLX( & 1228 & w1*DBLE(psi(1:nnr_))*AIMAG(red_revc0(1:nnr_,ibnd2-1,1)),& 1229 & w2*AIMAG(psi(1:nnr_))*AIMAG(red_revc0(1:nnr_,ibnd2-1,1)), kind=DP ) 1230 ! 1231 CALL fwfft ('Rho', pseudo_dens_c, dfftt) 1232 ! 1233 ! hartree contribution is computed in reciprocal space 1234 ! 1235 DO is = 1, nspin 1236 ! 1237 vhart(dfftt%nl(1:ngm_),is) = pseudo_dens_c(dfftt%nl(1:ngm_)) * fac_in(1:ngm_) 1238 IF (gamma_only) vhart(dfftt%nlm(1:ngm_),is) = & 1239 & pseudo_dens_c(dfftt%nlm(1:ngm_)) *& 1240 & fac_in(1:ngm_) 1241 ! 1242 ! and transformed back to real space 1243 ! 1244 CALL invfft ('Rho', vhart (:, is), dfftt) 1245 ! 1246 ENDDO 1247 ! 1248 ! 1249 ! Finally return the interaction 1250 ! 1251 psi_int(1:nnr_,ibnd2) = psi_int(1:nnr_,ibnd2) & 1252 & +DBLE(vhart(1:nnr_,1)) * DBLE(red_revc0(1:nnr_,ibnd,1)) & 1253 & +AIMAG(vhart(1:nnr_,1)) * AIMAG(red_revc0(1:nnr_,ibnd,1)) 1254 ! 1255 ! 1256 ! and interaction for Vexx term 1257 ! 1258 IF (interaction) then 1259 psi_int(1:nnr_,ibnd) = psi_int(1:nnr_,ibnd) & 1260 & +DBLE(vhart(1:nnr_,1)) * AIMAG(red_revc0(1:nnr_,ibnd2-1,1)) 1261 IF (ibnd < nbnd) then 1262 psi_int(1:nnr_,ibnd+1) = psi_int(1:nnr_,ibnd+1) & 1263 & +AIMAG(vhart(1:nnr_,1)) * AIMAG(red_revc0(1:nnr_,ibnd2-1,1)) 1264 ENDIF 1265 ELSE 1266 psi_int(1:nnr_,ibnd) = psi_int(1:nnr_,ibnd) & 1267 & -DBLE(vhart(1:nnr_,1)) * AIMAG(red_revc0(1:nnr_,ibnd2-1,1)) 1268 IF (ibnd < nbnd) then 1269 psi_int(1:nnr_,ibnd+1) = psi_int(1:nnr_,ibnd+1) & 1270 & -AIMAG(vhart(1:nnr_,1)) * AIMAG(red_revc0(1:nnr_,ibnd2-1,1)) 1271 ENDIF 1272 ENDIF 1273 ! 1274 ENDIF 1275 ! 1276 ! 1277 ENDDO 1278 ! 1279 RETURN 1280 ! 1281END FUNCTION k2d_vexx_term_gamma 1282! 1283FUNCTION k2d_term_k(w1, psi, fac_in, ibnd, ik, ikq) RESULT (psi_int) 1284 !------------------------------------------------------------------ 1285 ! 1286 ! This routine computes the K^2d term, Eq (22) from Eq Ref (1). 1287 ! psi contains two bands of |b> with w1, w2 the associated weights 1288 ! fac_in contains the interaction W(G) and ibnd the band index n'. 1289 ! 1290 ! (1) Rocca, Lu and Galli, J. Chem. Phys., 133, 164109 (2010) 1291 !------------------------------------------------------------------ 1292 ! 1293 1294 IMPLICIT NONE 1295 1296 COMPLEX(KIND=DP),DIMENSION(:), INTENT(IN) :: psi 1297 REAL(KIND=DP),DIMENSION(:), INTENT(IN) :: fac_in 1298 REAL(kind=DP), INTENT(IN) :: w1 1299 INTEGER, INTENT(IN) :: ibnd, ik, ikq 1300 COMPLEX(KIND=DP), ALLOCATABLE :: psi_int(:,:) 1301 ! 1302 ! Workspaces 1303 ! 1304 INTEGER :: ibnd2, is 1305 ! 1306 ALLOCATE(psi_int(dffts%nnr, nbnd)) 1307 psi_int = (0.d0,0.d0) 1308 ! 1309 DO ibnd2=1,nbnd,1 1310 ! 1311 ! Set up the pseudo density for this Hartree like interaction. 1312 ! 1313 vhart(:,:) = (0.d0,0.d0) 1314 pseudo_dens_c(:) = (0.d0,0.d0) 1315 ! 1316 pseudo_dens_c(:) = CONJG(psi(:))*red_revc0(:,ibnd2,k2q(ik))/omega 1317 ! 1318 CALL fwfft ('Rho', pseudo_dens_c, dffts) 1319 ! 1320 ! hartree contribution is computed in reciprocal space 1321 ! 1322 DO is = 1, nspin 1323 ! 1324 vhart(dffts%nl(1:ngm),is) = w1*pseudo_dens_c(dffts%nl(1:ngm)) *& 1325 & fac_in(1:ngm) 1326 ! 1327 ! and transformed back to real space 1328 ! 1329 CALL invfft ('Rho', vhart (:, is), dffts) 1330 ! 1331 ENDDO 1332 ! 1333 ! 1334 ! Finally return the interaction 1335 ! 1336 psi_int(:,ibnd2) = psi_int(:,ibnd2) + vhart(:,1) *& 1337 &red_revc0(:,ibnd,ikq) 1338 ! 1339 ENDDO 1340 ! 1341END FUNCTION k2d_term_k 1342 1343!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1344!! 1345!! The following routines are wrapper functions for inv/fw fft handling both 1346!! the custom FFT grids and gamma tricks. They are the analogues of those found 1347!! in PW/src/realus.f90. Ideally both these and their counterparts should be 1348!! moved somewhere else but for now they live here. 1349!! 1350!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1351SUBROUTINE invfft_orbital_custom_gamma(orbital, ibnd, nbnd, npwt, dfftt) 1352 1353 USE kinds, ONLY : DP 1354 USE fft_types, ONLY : fft_type_descriptor 1355 1356 IMPLICIT NONE 1357 1358 COMPLEX(DP), INTENT(IN) :: orbital(:,:) 1359 INTEGER, INTENT(IN) :: ibnd, nbnd, npwt 1360 TYPE(fft_type_descriptor), INTENT(IN) :: dfftt 1361 ! 1362 psic=(0.0_dp, 0.0_dp) 1363 ! 1364 IF (ibnd < nbnd) THEN 1365 ! 1366 psic(dfftt%nl(1:npwt)) = orbital(1:npwt,ibnd) + & 1367 &(0.0_dp, 1.0_dp) * orbital(1:npwt,ibnd+1) 1368 psic(dfftt%nlm(1:npwt)) = CONJG(orbital(1:npwt,ibnd) - & 1369 &(0.0_dp, 1.0_dp) * orbital(1:npwt,ibnd+1)) 1370 ! 1371 ELSE 1372 ! 1373 psic(dfftt%nl(1:npwt)) = orbital(1:npwt,ibnd) 1374 psic(dfftt%nlm(1:npwt)) =CONJG(orbital(1:npwt,ibnd)) 1375 ! 1376 ENDIF 1377 ! 1378 CALL invfft ('Wave', psic, dfftt) 1379 ! 1380 RETURN 1381 ! 1382END SUBROUTINE invfft_orbital_custom_gamma 1383 1384SUBROUTINE fwfft_orbital_custom_gamma(orbital, ibnd, nbnd, npwt, dfftt) 1385 1386 USE kinds, ONLY : DP 1387 USE fft_types, ONLY : fft_type_descriptor 1388 1389 IMPLICIT NONE 1390 1391 COMPLEX(DP), INTENT(INOUT) :: orbital(:,:) 1392 INTEGER, INTENT(IN) :: ibnd, nbnd, npwt 1393 TYPE(fft_type_descriptor), INTENT(IN) :: dfftt 1394 1395 ! Workspaces 1396 COMPLEX(DP) :: fp, fm 1397 ! Counters 1398 INTEGER :: j 1399 ! 1400 CALL fwfft ('Wave', psic(:), dfftt) 1401 ! 1402 IF (ibnd < nbnd) THEN 1403 ! 1404 ! two ffts at the same time 1405 DO j = 1, npwt 1406 fp = (psic(dfftt%nl(j)) + psic(dfftt%nlm(j)))*0.5d0 1407 fm = (psic(dfftt%nl(j)) - psic(dfftt%nlm(j)))*0.5d0 1408 orbital( j, ibnd) = CMPLX( DBLE(fp), AIMAG(fm),kind=DP) 1409 orbital( j, ibnd+1) = CMPLX(AIMAG(fp),- DBLE(fm),kind=DP) 1410 ENDDO 1411 ! 1412 ELSE 1413 ! 1414 orbital(1:npwt,ibnd)=psic(dfftt%nl(1:npwt)) 1415 ! 1416 ENDIF 1417 ! 1418 RETURN 1419 ! 1420END SUBROUTINE fwfft_orbital_custom_gamma 1421 1422SUBROUTINE invfft_orbital_ibnd2_gamma(orbital, psitemp, ibnd2, npwt, dfftt) 1423 1424 USE kinds, ONLY : DP 1425 USE fft_types, ONLY : fft_type_descriptor 1426 1427 IMPLICIT NONE 1428 1429 COMPLEX(DP), INTENT(IN) :: orbital(:,:) 1430 INTEGER, INTENT(IN) :: ibnd2, npwt 1431 TYPE(fft_type_descriptor), INTENT(IN) :: dfftt 1432 COMPLEX(DP), INTENT(OUT) :: psitemp(:) 1433 ! 1434 psitemp=(0.0_dp, 0.0_dp) 1435 ! 1436 psitemp(dfftt%nl(1:npwt)) = orbital(1:npwt,ibnd2) 1437 psitemp(dfftt%nlm(1:npwt)) = CONJG(orbital(1:npwt,ibnd2)) 1438 ! 1439 CALL invfft ('Wave', psitemp, dfftt) 1440 ! 1441 RETURN 1442 ! 1443END SUBROUTINE invfft_orbital_ibnd2_gamma 1444 1445END MODULE lr_exx_kernel 1446