1! 2! Copyright (C) 2001-2013 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 10! 11 12SUBROUTINE o_rcgdiagg( npwx, npw, nbnd, psi, e, precondition, & 13 ethr, maxter, reorder, notconv, avg_iter ,& 14 numv, v_states,hdiag,ptype,fcw_number,fcw_state,fcw_mat) 15 !---------------------------------------------------------------------------- 16 ! 17 ! ... "poor man" iterative diagonalization of a complex hermitian matrix 18 ! ... through preconditioned conjugate gradient algorithm 19 ! ... Band-by-band algorithm with minimal use of memory 20 21 ! 22 USE constants, ONLY : pi 23 USE kinds, ONLY : DP 24 USE gvect, ONLY : gstart 25 USE mp, ONLY : mp_sum 26 USE mp_world, ONLY : world_comm 27 USE fft_base, ONLY : dffts 28 USE io_global, ONLY :stdout 29 30 ! 31 IMPLICIT NONE 32 ! 33 ! ... I/O variables 34 ! 35 INTEGER, INTENT(IN) :: npwx, npw, nbnd, maxter 36 REAL (DP), INTENT(IN) :: precondition(npw), ethr 37 COMPLEX (DP), INTENT(INOUT) :: psi(npwx,nbnd) 38 REAL (DP), INTENT(INOUT) :: e(nbnd) 39 INTEGER, INTENT(OUT) :: notconv 40 REAL (DP), INTENT(OUT) :: avg_iter 41 42 INTEGER, INTENT(in) :: numv!number of valence states 43 REAL(kind=DP), INTENT(in) :: v_states(dffts%nnr,numv)!valence states in real space 44 REAL(kind=DP), INTENT(in) :: hdiag(npw)!inverse of estimation of diagonal part of hamiltonian 45 INTEGER, INTENT(in) :: ptype!type of approximation for O operator 46 INTEGER, INTENT(in) :: fcw_number!number of "fake conduction" states for O matrix method 47 COMPLEX(kind=DP) :: fcw_state(npw,fcw_number)! "fake conduction" states for O matrix method 48 REAL(kind=DP) :: fcw_mat(fcw_number,fcw_number)! "fake conduction" matrix 49 50 51 ! 52 ! ... local variables 53 ! 54 INTEGER :: i, j, m, iter, moved,ig 55 REAL (DP), ALLOCATABLE :: lagrange(:) 56 COMPLEX (DP), ALLOCATABLE :: hpsi(:), spsi(:), g(:), cg(:), & 57 scg(:), ppsi(:), g0(:) 58 REAL (DP) :: psi_norm, a0, b0, gg0, gamma, gg, gg1, & 59 cg0, e0, es(2) 60 REAL (DP) :: theta, cost, sint, cos2t, sin2t 61 LOGICAL :: reorder 62 INTEGER :: npw2, npwx2 63 REAL (DP) :: empty_ethr,sca 64 ! 65 ! ... external functions 66 ! 67 REAL (DP), EXTERNAL :: ddot 68 ! 69 ! 70 CALL start_clock( 'rcgdiagg' ) 71 ! 72 empty_ethr = MAX( ( ethr * 5.D0 ), 1.D-5 ) 73 ! 74 npw2 = 2 * npw 75 npwx2 = 2 * npwx 76 ! 77 ALLOCATE( spsi( npwx ) ) 78 ALLOCATE( scg( npwx ) ) 79 ALLOCATE( hpsi( npwx ) ) 80 ALLOCATE( g( npwx ) ) 81 ALLOCATE( cg( npwx ) ) 82 ALLOCATE( g0( npwx ) ) 83 ALLOCATE( ppsi( npwx ) ) 84 ! 85 ALLOCATE( lagrange( nbnd ) ) 86 ! 87 avg_iter = 0.D0 88 notconv = 0 89 moved = 0 90 ! 91 ! ... every eigenfunction is calculated separately 92 ! 93 DO m = 1, nbnd 94 ! 95 ! ... calculate S|psi> 96 ! 97 spsi(1:npw)=psi(1:npw,m) 98 99 ! 100 ! ... orthogonalize starting eigenfunction to those already calculated 101 ! 102 CALL DGEMV( 'T', npw2, m, 2.D0, psi, npwx2, spsi, 1, 0.D0, lagrange, 1 ) 103 ! 104 IF ( gstart == 2 ) lagrange(1:m) = lagrange(1:m) - psi(1,1:m) * spsi(1) 105 ! 106 CALL mp_sum( lagrange( 1:m ), world_comm) 107 ! 108 psi_norm = lagrange(m) 109 ! 110 DO j = 1, m - 1 111 ! 112 psi(:,m) = psi(:,m) - lagrange(j) * psi(:,j) 113 ! 114 psi_norm = psi_norm - lagrange(j)**2 115 ! 116 END DO 117 ! 118 psi_norm = SQRT( psi_norm ) 119 ! 120 psi(:,m) = psi(:,m) / psi_norm 121 ! ... set Im[ psi(G=0) ] - needed for numerical stability 122 IF ( gstart == 2 ) psi(1,m) = CMPLX( DBLE(psi(1,m)), 0.D0 ,kind=DP) 123 ! 124 ! ... calculate starting gradient (|hpsi> = H|psi>) ... 125 ! 126 !CALL hs_1psi( npwx, npw, psi(1,m), hpsi, spsi ) 127 call o_1psi_gamma( numv, v_states, psi(1,m), hpsi,.false.,hdiag, ptype,fcw_number,fcw_state,fcw_mat,ethr) 128 spsi(1:npw)=psi(1:npw,m) 129 130 131 ! 132 ! ... and starting eigenvalue (e = <y|PHP|y> = <psi|H|psi>) 133 ! 134 ! ... NB: ddot(2*npw,a,1,b,1) = DBLE( zdotc(npw,a,1,b,1) ) 135 ! 136 e(m) = 2.D0 * ddot( npw2, psi(1,m), 1, hpsi, 1 ) 137 ! 138 IF ( gstart == 2 ) e(m) = e(m) - psi(1,m) * hpsi(1) 139 ! 140 CALL mp_sum( e(m), world_comm) 141 ! 142 ! ... start iteration for this band 143 ! 144 iterate: DO iter = 1, maxter 145 ! 146 ! ... calculate P (PHP)|y> 147 ! ... ( P = preconditioning matrix, assumed diagonal ) 148 ! 149 g(1:npw) = hpsi(1:npw)! / precondition(:) 150 ppsi(1:npw) = spsi(1:npw)! / precondition(:) 151 ! 152 ! ... ppsi is now S P(P^2)|y> = S P^2|psi>) 153 ! 154 es(1) = 2.D0 * ddot( npw2, spsi(1), 1, g(1), 1 ) 155 es(2) = 2.D0 * ddot( npw2, spsi(1), 1, ppsi(1), 1 ) 156 ! 157 IF ( gstart == 2 ) THEN 158 ! 159 es(1) = es(1) - spsi(1) * g(1) 160 es(2) = es(2) - spsi(1) * ppsi(1) 161 ! 162 END IF 163 ! 164 CALL mp_sum( es, world_comm ) 165 ! 166 es(1) = es(1) / es(2) 167 ! 168 g(:) = g(:) - es(1) * ppsi(:) 169 ! 170 ! ... e1 = <y| S P^2 PHP|y> / <y| S S P^2|y> ensures that 171 ! ... <g| S P^2|y> = 0 172 ! 173 ! ... orthogonalize to lowest eigenfunctions (already calculated) 174 ! 175 ! ... scg is used as workspace 176 ! 177 !CALL s_1psi( npwx, npw, g(1), scg(1) ) 178 scg(1:npw)=g(1:npw) 179 ! 180 CALL DGEMV( 'T', npw2, ( m - 1 ), 2.D0, & 181 psi, npwx2, scg, 1, 0.D0, lagrange, 1 ) 182 ! 183 IF ( gstart == 2 ) & 184 lagrange(1:m-1) = lagrange(1:m-1) - psi(1,1:m-1) * scg(1) 185 ! 186 CALL mp_sum( lagrange( 1 : m-1 ), world_comm) 187 ! 188 DO j = 1, ( m - 1 ) 189 ! 190 g(:) = g(:) - lagrange(j) * psi(:,j) 191 scg(:) = scg(:) - lagrange(j) * psi(:,j) 192 ! 193 END DO 194 ! 195 IF ( iter /= 1 ) THEN 196 ! 197 ! ... gg1 is <g(n+1)|S|g(n)> (used in Polak-Ribiere formula) 198 ! 199 gg1 = 2.D0 * ddot( npw2, g(1), 1, g0(1), 1 ) 200 ! 201 IF ( gstart == 2 ) gg1 = gg1 - g(1) * g0(1) 202 ! 203 CALL mp_sum( gg1, world_comm ) 204 ! 205 END IF 206 ! 207 ! ... gg is <g(n+1)|S|g(n+1)> 208 ! 209 g0(:) = scg(:) 210 ! 211 g0(1:npw) = g0(1:npw)! * precondition(:) 212 ! 213 gg = 2.D0 * ddot( npw2, g(1), 1, g0(1), 1 ) 214 ! 215 IF ( gstart == 2 ) gg = gg - g(1) * g0(1) 216 ! 217 CALL mp_sum( gg, world_comm ) 218 ! 219 IF ( iter == 1 ) THEN 220 ! 221 ! ... starting iteration, the conjugate gradient |cg> = |g> 222 ! 223 gg0 = gg 224 ! 225 cg(:) = g(:) 226 ! 227 ELSE 228 ! 229 ! ... |cg(n+1)> = |g(n+1)> + gamma(n) * |cg(n)> 230 ! 231 ! ... Polak-Ribiere formula : 232 ! 233 gamma = ( gg - gg1 ) / gg0 234 gg0 = gg 235 ! 236 cg(:) = cg(:) * gamma 237 cg(:) = g + cg(:) 238 ! 239 ! ... The following is needed because <y(n+1)| S P^2 |cg(n+1)> 240 ! ... is not 0. In fact : 241 ! ... <y(n+1)| S P^2 |cg(n)> = sin(theta)*<cg(n)|S|cg(n)> 242 ! 243 psi_norm = gamma * cg0 * sint 244 ! 245 cg(:) = cg(:) - psi_norm * psi(:,m) 246 ! 247 END IF 248 ! 249 ! ... |cg> contains now the conjugate gradient 250 ! ... set Im[ cg(G=0) ] - needed for numerical stability 251 IF ( gstart == 2 ) cg(1) = CMPLX( DBLE(cg(1)), 0.D0 ,kind=DP) 252 ! 253 ! ... |scg> is S|cg> 254 ! 255 !CALL hs_1psi( npwx, npw, cg(1), ppsi(1), scg(1) ) 256 call o_1psi_gamma( numv, v_states, cg, ppsi,.false.,hdiag, ptype,fcw_number,fcw_state,fcw_mat,ethr) 257 sca=0.d0 258 do ig=1,npw 259 sca=sca+2.d0*dble(conjg(cg(ig))*ppsi(ig)) 260 enddo 261 if(gstart==2) sca=sca-dble(conjg(cg(1))*ppsi(1)) 262 call mp_sum(sca, world_comm) 263 264 265 scg(1:npw)=cg(1:npw) 266 ! 267 cg0 = 2.D0 * ddot( npw2, cg(1), 1, scg(1), 1 ) 268 ! 269 IF ( gstart == 2 ) cg0 = cg0 - cg(1) * scg(1) 270 ! 271 CALL mp_sum( cg0, world_comm ) 272 ! 273 cg0 = SQRT( cg0 ) 274 ! 275 ! ... |ppsi> contains now HP|cg> 276 ! ... minimize <y(t)|PHP|y(t)> , where : 277 ! ... |y(t)> = cos(t)|y> + sin(t)/cg0 |cg> 278 ! ... Note that <y|P^2S|y> = 1, <y|P^2S|cg> = 0 , 279 ! ... <cg|P^2S|cg> = cg0^2 280 ! ... so that the result is correctly normalized : 281 ! ... <y(t)|P^2S|y(t)> = 1 282 ! 283 a0 = 4.D0 * ddot( npw2, psi(1,m), 1, ppsi(1), 1 ) 284 ! 285 IF ( gstart == 2 ) a0 = a0 - 2.D0 * psi(1,m) * ppsi(1) 286 ! 287 a0 = a0 / cg0 288 ! 289 CALL mp_sum( a0, world_comm ) 290 ! 291 b0 = 2.D0 * ddot( npw2, cg(1), 1, ppsi(1), 1 ) 292 ! 293 IF ( gstart == 2 ) b0 = b0 - cg(1) * ppsi(1) 294 ! 295 b0 = b0 / cg0**2 296 ! 297 CALL mp_sum( b0, world_comm ) 298 ! 299 e0 = e(m) 300 ! 301 theta = 0.5D0 * ATAN( a0 / ( e0 - b0 ) ) 302 ! 303 cost = COS( theta ) 304 sint = SIN( theta ) 305 ! 306 cos2t = cost*cost - sint*sint 307 sin2t = 2.D0*cost*sint 308 ! 309 es(1) = 0.5D0 * ( ( e0 - b0 ) * cos2t + a0 * sin2t + e0 + b0 ) 310 es(2) = 0.5D0 * ( - ( e0 - b0 ) * cos2t - a0 * sin2t + e0 + b0 ) 311 ! 312 ! ... there are two possible solutions, choose the minimum 313 ! 314 IF ( es(2) < es(1) ) THEN 315 ! 316 theta = theta + 0.5D0 * pi 317 ! 318 cost = COS( theta ) 319 sint = SIN( theta ) 320 ! 321 END IF 322 ! 323 ! ... new estimate of the eigenvalue 324 ! 325 e(m) = MIN( es(1), es(2) ) 326 327 328 ! 329 ! ... upgrade |psi> 330 ! 331 psi(:,m) = cost * psi(:,m) + sint / cg0 * cg(:) 332 ! 333 ! ... here one could test convergence on the energy 334 ! 335 336 ! 337 IF ( ABS( e(m) - e0 ) < ethr ) EXIT iterate 338 ! 339 340 ! 341 ! ... upgrade H|psi> and S|psi> 342 ! 343 spsi(:) = cost * spsi(:) + sint / cg0 * scg(:) 344 ! 345 hpsi(:) = cost * hpsi(:) + sint / cg0 * ppsi(:) 346 ! 347 END DO iterate 348 ! 349 IF ( iter >= maxter ) notconv = notconv + 1 350 ! 351 avg_iter = avg_iter + iter + 1 352 ! 353 ! ... reorder eigenvalues if they are not in the right order 354 ! ... ( this CAN and WILL happen in not-so-special cases ) 355 ! 356 IF ( m > 1 .AND. reorder ) THEN 357 ! 358 IF ( e(m) - e(m-1) < - 2.D0 * ethr ) THEN 359 ! 360 ! ... if the last calculated eigenvalue is not the largest... 361 ! 362 DO i = m - 2, 1, - 1 363 ! 364 IF ( e(m) - e(i) > 2.D0 * ethr ) EXIT 365 ! 366 END DO 367 ! 368 i = i + 1 369 ! 370 moved = moved + 1 371 ! 372 ! ... last calculated eigenvalue should be in the 373 ! ... i-th position: reorder 374 ! 375 e0 = e(m) 376 ! 377 ppsi(:) = psi(:,m) 378 ! 379 DO j = m, i + 1, - 1 380 ! 381 e(j) = e(j-1) 382 ! 383 psi(:,j) = psi(:,j-1) 384 ! 385 END DO 386 ! 387 e(i) = e0 388 ! 389 psi(:,i) = ppsi(:) 390 ! 391 ! ... this procedure should be good if only a few inversions occur, 392 ! ... extremely inefficient if eigenvectors are often in bad order 393 ! ... ( but this should not happen ) 394 ! 395 END IF 396 ! 397 END IF 398 ! 399 END DO 400 ! 401 avg_iter = avg_iter / DBLE( nbnd ) 402 ! 403 DEALLOCATE( lagrange ) 404 DEALLOCATE( ppsi ) 405 DEALLOCATE( g0 ) 406 DEALLOCATE( cg ) 407 DEALLOCATE( g ) 408 DEALLOCATE( hpsi ) 409 DEALLOCATE( scg ) 410 DEALLOCATE( spsi ) 411 ! 412 CALL stop_clock( 'rcgdiagg' ) 413 ! 414 RETURN 415 ! 416END SUBROUTINE o_rcgdiagg 417 418 419 420! 421!---------------------------------------------------------------------------- 422SUBROUTINE o_1psi_gamma( numv, v_states, psi, opsi,l_freq,hdiag, ptype,fcw_number,fcw_state,fcw_mat,ethr) 423 !---------------------------------------------------------------------------- 424 ! 425 ! 426 !this subroutines applies the O oprator to a state psi 427 !IT WORKS ONLY FOR NORMCONSERVING PSEUDOPOTENTIALS 428 !the valence states in G space must be in evc 429 ! Gamma point version 430 431 USE io_global, ONLY : stdout, ionode, ionode_id 432 USE kinds, ONLY : DP 433 USE wannier_gw 434 USE gvect 435 USE constants, ONLY : e2, pi, tpi, fpi 436 USE cell_base, ONLY: at, alat, tpiba, omega, tpiba2 437 USE wvfct, ONLY : g2kin, npwx, npw, nbnd, et 438 USE wavefunctions, ONLY : evc, psic 439 USE mp, ONLY : mp_sum, mp_barrier, mp_bcast 440 USE mp_world, ONLY : world_comm, mpime, nproc 441 USE gvecs, ONLY : doublegrid 442 USE kinds, ONLY : DP 443 USE fft_base, ONLY : dfftp, dffts 444 USE fft_interfaces, ONLY : fwfft, invfft 445 USE becmod, ONLY : becp,allocate_bec_type,deallocate_bec_type 446 USE uspp, ONLY : vkb, nkb, okvan 447 USE g_psi_mod, ONLY : h_diag, s_diag 448 USE klist, ONLY : xk,igk_k 449 450 ! 451 IMPLICIT NONE 452 453 INTEGER, INTENT(in) :: numv!number of valence states 454 REAL(kind=DP), INTENT(in) :: v_states(dffts%nnr,numv)!valence states in real space 455 COMPLEX(kind=DP), INTENT(in) :: psi(npw)!input wavefunction 456 COMPLEX(kind=DP), INTENT(out) :: opsi(npw)!O|\psi> 457 LOGICAL, INTENT(in) :: l_freq!if true estimates the operator a 0 frequency 458 REAL(kind=DP), INTENT(in) :: hdiag(npw)!inverse of estimation of diagonal part of hamiltonian 459 INTEGER, INTENT(in) :: ptype!type of approximation for O operator 460 INTEGER, INTENT(in) :: fcw_number!number of "fake conduction" states for O matrix method 461 COMPLEX(kind=DP) :: fcw_state(npw,fcw_number)! "fake conduction" states for O matrix method 462 REAL(kind=DP) :: fcw_mat(fcw_number,fcw_number)! "fake conduction" matrix 463 REAL(kind=DP), INTENT(in) :: ethr!threshold on (H-e) inversion 464 465 466 REAL(kind=DP), ALLOCATABLE :: psi_r(:,:), psi_v(:) 467 COMPLEX(kind=DP), ALLOCATABLE :: psi_g(:,:), h_psi_g(:,:),s_psi_g(:,:) 468 REAL(kind=DP) :: ec 469 INTEGER :: iv,ii,jj,ig 470 471 REAL(kind=DP),ALLOCATABLE :: p_terms(:),s_terms(:) 472 INTEGER :: l_blk,nbegin,nend,nsize 473 !REAL(kind=DP), ALLOCATABLE :: h_diag (:,:) 474 COMPLEX(kind=DP), ALLOCATABLE :: psi_g2(:,:) 475 INTEGER :: kter 476 LOGICAL :: lconv_root,lfirst 477 REAL(kind=DP) :: anorm 478 EXTERNAL :: hpsi_pw4gww,cg_psi_pw4gww 479 COMPLEX(kind=DP), POINTER, SAVE :: tmp_psi(:,:) 480 481 482 allocate(psi_r(dffts%nnr,2),psi_v(dffts%nnr)) 483 if(pmat_type/=0) then 484 allocate(h_psi_g(npw,2),s_psi_g(npw,2),psi_g(npw,2)) 485 else 486 allocate(h_psi_g(npw,numv),s_psi_g(npw,numv),psi_g(npw,numv)) 487 endif 488 489!fourier transform psi to R space 490 491 opsi(1:npw)=(0.d0,0.d0) 492 493 494 if(pmat_type==0 .or. pmat_type==1 .or. pmat_type==2) then 495 psic(:)=(0.d0,0.d0) 496 psic(dffts%nl(1:npw)) = psi(1:npw) 497 psic(dffts%nlm(1:npw)) = CONJG( psi(1:npw) ) 498 CALL invfft ('Wave', psic, dffts) 499 psi_v(:)= DBLE(psic(:)) 500 endif 501 502 if(pmat_type==0) then 503 call start_clock('opsi_total') 504 call allocate_bec_type ( nkb, numv, becp) 505 IF ( nkb > 0 ) CALL init_us_2( npw, igk_k(1,1), xk(1,1), vkb ) 506 if(.not.associated(tmp_psi)) then 507 allocate( tmp_psi(npw,num_nbndv(1))) 508 lfirst=.true. 509 else 510 lfirst=.false. 511 endif 512 allocate (h_diag(npw, numv),s_diag(npw,numv)) 513 allocate(psi_g2(npw,numv)) 514 ! 515 ! compute the kinetic energy 516 ! 517 do ig = 1, npw 518 g2kin (ig) = ( g (1,ig)**2 + g (2,ig)**2 + g (3,ig)**2 ) * tpiba2 519 enddo 520 h_diag=0.d0 521 do iv = 1, numv 522 do ig = 1, npw 523 !h_diag(ig,ibnd)=1.d0/max(1.0d0,g2kin(ig)/eprec(ibnd,ik)) 524 !h_diag(ig,iv) = 1.D0 + g2kin(ig) + & 525 ! SQRT( 1.D0 + ( g2kin(ig) - 1.D0 )**2 ) 526 h_diag(ig,iv)=g2kin(ig) 527 !h_diag(ig,iv)=1.d0/max(1.0d0,g2kin(ig)/8.d0) 528 enddo 529 enddo 530 do iv=1,numv,2 531!!product with psi_v 532 if(iv/=numv) then 533 psi_r(1:dffts%nnr,1)=psi_v(1:dffts%nnr)*v_states(1:dffts%nnr, iv) 534 psi_r(1:dffts%nnr,2)=psi_v(1:dffts%nnr)*v_states(1:dffts%nnr, iv+1) 535 else 536 psi_r(1:dffts%nnr,1)=psi_v(1:dffts%nnr)*v_states(1:dffts%nnr, iv) 537 endif 538!!fourier transfrm to G 539 if(iv/=numv) then 540 psic(1:dffts%nnr)=cmplx(psi_r(1:dffts%nnr,1),psi_r(1:dffts%nnr,2)) 541 else 542 psic(1:dffts%nnr)=cmplx(psi_r(1:dffts%nnr,1),0.d0) 543 endif 544 call start_clock('opsi_fft') 545 CALL fwfft ('Wave', psic, dffts) 546 call stop_clock('opsi_fft') 547 if(iv/=numv) then 548 psi_g(1:npw,iv)=0.5d0*(psic(dffts%nl(1:npw))+conjg(psic(dffts%nlm(1:npw)))) 549 psi_g(1:npw,iv+1)=(0.d0,-0.5d0)*(psic(dffts%nl(1:npw))-conjg(psic(dffts%nlm(1:npw)))) 550 if(gstart==2) psi_g(1,iv)=dble(psi_g(1,iv)) 551 if(gstart==2) psi_g(1,iv+1)=dble(psi_g(1,iv+1)) 552 else 553 psi_g(1:npw,iv)=psic(dffts%nl(1:npw)) 554 if(gstart==2) psi_g(1,iv)=dble(psi_g(1,iv)) 555 endif 556 557!!project on conduction manifold 558 call start_clock('opsi_pc') 559 if(iv/=numv) then 560 call pc_operator(psi_g(:,iv),1,.false.) 561 call pc_operator(psi_g(:,iv+1),1,.false.) 562 else 563 call pc_operator(psi_g(:,iv),1,.false.) 564 endif 565 call stop_clock('opsi_pc') 566 enddo 567 568 write(stdout,*) 'DEBUG1' 569 FLUSH(stdout) 570!call (H-e)^-1 solver 571 if(.true.) then 572 psi_g2(1:npw,1:numv)=psi_g(1:npw,1:numv) 573 else 574 psi_g2(1:npw,1:numv)=tmp_psi(1:npw,1:numv) 575 endif 576 write(stdout,*) 'DEBUG1.5' 577 FLUSH(stdout) 578 call cgsolve_all_gamma (hpsi_pw4gww,cg_psi_pw4gww,et(1,1),psi_g,psi_g2, & 579 h_diag,npw,npw,ethr,1,kter,lconv_root,anorm,numv,1) 580 581 tmp_psi(1:npw,1:numv)=psi_g2(1:npw,1:numv) 582 write(stdout,*) 'DEBUG2',kter,lconv_root,anorm 583 FLUSH(stdout) 584 585 586 do iv=1,numv,2 587 !!project on conduction manifold 588 call start_clock('opsi_pc') 589 if(iv/=numv) then 590 call pc_operator(psi_g2(:,iv),1,.false.) 591 call pc_operator(psi_g2(:,iv+1),1,.false.) 592 else 593 call pc_operator(psi_g2(:,iv),1,.false.) 594 endif 595 call stop_clock('opsi_pc') 596 597 598 !!fourier transform to R space 599 psic(:)=(0.d0,0.d0) 600 if(iv/=numv) then 601 psic(dffts%nl(1:npw)) = psi_g2(1:npw,iv)+(0.d0,1.d0)*psi_g2(1:npw,iv+1) 602 psic(dffts%nlm(1:npw)) = CONJG( psi_g2(1:npw,iv) )+(0.d0,1.d0)*conjg(psi_g2(1:npw,iv+1)) 603 else 604 psic(dffts%nl(1:npw)) = psi_g2(1:npw,iv) 605 psic(dffts%nlm(1:npw)) = CONJG( psi_g2(1:npw,iv) ) 606 endif 607 call start_clock('opsi_fft') 608 CALL invfft ('Wave', psic, dffts) 609 call stop_clock('opsi_fft') 610 if(iv/=numv) then 611 psi_r(:,1)= DBLE(psic(:)) 612 psi_r(:,2)= dimag(psic(:)) 613 else 614 psi_r(:,1)= DBLE(psic(:)) 615 endif 616!!product with psi_v 617 if(iv/=numv) then 618 psi_r(1:dffts%nnr,1)=psi_r(1:dffts%nnr,1)*v_states(1:dffts%nnr,iv) 619 psi_r(1:dffts%nnr,2)=psi_r(1:dffts%nnr,2)*v_states(1:dffts%nnr,iv+1) 620 else 621 psi_r(1:dffts%nnr,1)=psi_r(1:dffts%nnr,1)*v_states(1:dffts%nnr,iv) 622 endif 623!!fourier transform in G space!! sum up results 624!!TAKE CARE OF SIGN 625 if(iv/=numv) then 626 psic(:)=cmplx(psi_r(:,1),psi_r(:,2)) 627 else 628 psic(:)=cmplx(psi_r(:,1),0.d0) 629 endif 630 CALL fwfft ('Wave', psic, dffts) 631 if(iv/=numv) then 632 opsi(1:npw)=opsi(1:npw)-0.5d0*(psic(dffts%nl(1:npw))+conjg(psic(dffts%nlm(1:npw)))) 633 opsi(1:npw)=opsi(1:npw)-(0.d0,-0.5d0)*(psic(dffts%nl(1:npw))-conjg(psic(dffts%nlm(1:npw)))) 634 else 635 opsi(1:npw)=opsi(1:npw)-psic(dffts%nl(igk_k(1:npw,1))) 636 endif 637 enddo 638 deallocate(h_diag,s_diag) 639 deallocate(psi_g2) 640 call deallocate_bec_type(becp) 641 call stop_clock('opsi_total') 642 ! call print_clock('opsi_total') 643 ! call print_clock('opsi_fft') 644 ! call print_clock('opsi_pc') 645 646 else if(pmat_type==1) then 647!loop on v 648 call start_clock('opsi_total') 649 do iv=1,numv,2 650!!product with psi_v 651 if(iv/=numv) then 652 psi_r(1:dffts%nnr,1)=psi_v(1:dffts%nnr)*v_states(1:dffts%nnr, iv) 653 psi_r(1:dffts%nnr,2)=psi_v(1:dffts%nnr)*v_states(1:dffts%nnr, iv+1) 654 else 655 psi_r(1:dffts%nnr,1)=psi_v(1:dffts%nnr)*v_states(1:dffts%nnr, iv) 656 endif 657!!fourier transfrm to G 658 if(iv/=numv) then 659 psic(1:dffts%nnr)=cmplx(psi_r(1:dffts%nnr,1),psi_r(1:dffts%nnr,2)) 660 else 661 psic(1:dffts%nnr)=cmplx(psi_r(1:dffts%nnr,1),0.d0) 662 endif 663 call start_clock('opsi_fft') 664 CALL fwfft ('Wave', psic, dffts) 665 call stop_clock('opsi_fft') 666 if(iv/=numv) then 667 psi_g(1:npw,1)=0.5d0*(psic(dffts%nl(1:npw))+conjg(psic(dffts%nlm(1:npw)))) 668 psi_g(1:npw,2)=(0.d0,-0.5d0)*(psic(dffts%nl(1:npw))-conjg(psic(dffts%nlm(1:npw)))) 669 if(gstart==2) psi_g(1,1)=dble(psi_g(1,1)) 670 if(gstart==2) psi_g(1,2)=dble(psi_g(1,2)) 671 else 672 psi_g(1:npw,1)=psic(dffts%nl(1:npw)) 673 if(gstart==2) psi_g(1,1)=dble(psi_g(1,1)) 674 endif 675 676!!project on conduction manifold 677 call start_clock('opsi_pc') 678 if(iv/=numv) then 679 call pc_operator(psi_g(:,1),1,.false.) 680 call pc_operator(psi_g(:,2),1,.false.) 681 else 682 call pc_operator(psi_g(:,1),1,.false.) 683 endif 684 !call pc_operator_test(psi_g) 685 if(l_freq) then 686 687 ! call cgsolve_all_gamma (h_psi, cg_psi, e, d0psi, dpsi, h_diag, & 688 ! ndmx, ndim, ethr, ik, kter, conv_root, anorm, nbnd, npol) 689 690 691 endif 692 call stop_clock('opsi_pc') 693!!fourier transform to R space 694 psic(:)=(0.d0,0.d0) 695 if(iv/=numv) then 696 psic(dffts%nl(1:npw)) = psi_g(1:npw,1)+(0.d0,1.d0)*psi_g(1:npw,2) 697 psic(dffts%nlm(1:npw)) = CONJG( psi_g(1:npw,1) )+(0.d0,1.d0)*conjg(psi_g(1:npw,2)) 698 else 699 psic(dffts%nl(1:npw)) = psi_g(1:npw,1) 700 psic(dffts%nlm(1:npw)) = CONJG( psi_g(1:npw,1) ) 701 endif 702 call start_clock('opsi_fft') 703 CALL invfft ('Wave', psic, dffts) 704 call stop_clock('opsi_fft') 705 if(iv/=numv) then 706 psi_r(:,1)= DBLE(psic(:)) 707 psi_r(:,2)= dimag(psic(:)) 708 else 709 psi_r(:,1)= DBLE(psic(:)) 710 endif 711!!product with psi_v 712 713 if(iv/=numv) then 714 psi_r(1:dffts%nnr,1)=psi_r(1:dffts%nnr,1)*v_states(1:dffts%nnr,iv) 715 psi_r(1:dffts%nnr,2)=psi_r(1:dffts%nnr,2)*v_states(1:dffts%nnr,iv+1) 716 else 717 psi_r(1:dffts%nnr,1)=psi_r(1:dffts%nnr,1)*v_states(1:dffts%nnr,iv) 718 endif 719!!fourier transform in G space 720!! sum up results 721!!TAKE CARE OF SIGN 722 723 if(iv/=numv) then 724 psic(:)=cmplx(psi_r(:,1),psi_r(:,2)) 725 else 726 psic(:)=cmplx(psi_r(:,1),0.d0) 727 endif 728 CALL fwfft ('Wave', psic, dffts) 729 if(iv/=numv) then 730 opsi(1:npw)=opsi(1:npw)-0.5d0*(psic(dffts%nl(1:npw))+conjg(psic(dffts%nlm(1:npw)))) 731 opsi(1:npw)=opsi(1:npw)-(0.d0,-0.5d0)*(psic(dffts%nl(1:npw))-conjg(psic(dffts%nlm(1:npw)))) 732 else 733 opsi(1:npw)=opsi(1:npw)-psic(dffts%nl(igk_k(1:npw,1))) 734 endif 735 736 enddo 737 call stop_clock('opsi_total') 738 ! call print_clock('opsi_total') 739 ! call print_clock('opsi_fft') 740 ! call print_clock('opsi_pc') 741 else if(pmat_type==2) then 742 psi_r(:,1)=0.d0 743 do iv=1,numv 744 psi_r(1:dffts%nnr,1)=psi_r(1:dffts%nnr,1)+psi_v(1:dffts%nnr)*v_states(1:dffts%nnr, iv) 745 enddo 746 747 psic(:)=cmplx(psi_r(:,1),0.d0) 748 CALL fwfft ('Wave', psic, dffts) 749 psi_g(1:npw,1)=psic(dffts%nl(igk_k(1:npw,1))) 750 if(gstart==2) psi_g(1,1)=dble(psi_g(1,1)) 751 752 call pc_operator(psi_g(:,1),1,.false.) 753 psi_g(1:npw,1)=psi_g(1:npw,1)*hdiag(1:npw) 754 call pc_operator(psi_g(:,1),1,.false.) 755 756 psic(dffts%nl(igk_k(1:npw,1))) = psi_g(1:npw,1) 757 psic(dffts%nlm(igk_k(1:npw,1))) = CONJG( psi_g(1:npw,1) ) 758 CALL invfft ('Wave', psic, dffts) 759 psi_r(:,1)=dble(psic(:)) 760 761 psi_v(:)= psi_r(:,1) 762 psi_r(:,1)=0.d0 763 do iv=1,numv 764 psi_r(1:dffts%nnr,1)=psi_r(1:dffts%nnr,1)+psi_v(1:dffts%nnr)*v_states(1:dffts%nnr, iv) 765 enddo 766 767 psic(:)=cmplx(psi_r(:,1),0.d0) 768 CALL fwfft ('Wave', psic, dffts) 769 opsi(1:npw)=-psic(dffts%nl(igk_k(1:npw,1))) 770 771 else!cases 3,4 772!form scalar products 773 allocate(p_terms(fcw_number),s_terms(fcw_number)) 774 call dgemm('T','N',fcw_number,1,2*npw,2.d0,fcw_state,2*npw,psi,2*npw,0.d0,p_terms,fcw_number) 775 if(gstart==2) then 776 do ii=1,fcw_number 777 p_terms(ii)=p_terms(ii)-dble(conjg(fcw_state(1,ii))*psi(1)) 778 enddo 779 endif 780 call mp_sum(p_terms(:),world_comm) 781 782!multiply to D matrix 783 l_blk= (fcw_number)/nproc 784 if(l_blk*nproc < (fcw_number)) l_blk = l_blk+1 785 nbegin=mpime*l_blk+1 786 nend=nbegin+l_blk-1 787 if(nend > fcw_number) nend=fcw_number 788 nsize=nend-nbegin+1 789 790 791 s_terms(:)=0.d0 792 if(nsize>0) then 793 call dgemm('T','N',nsize,1,fcw_number,1.d0,fcw_mat,fcw_number,p_terms,fcw_number,0.d0,& 794 &s_terms(nbegin:nend),nsize) 795 endif 796 797!collect from processors 798 call mp_sum(s_terms,world_comm) 799 800!multiply with gamma functions 801 call dgemm('N','N',2*npw,1,fcw_number,-1.d0,fcw_state,2*npw,s_terms,fcw_number,0.d0,opsi,2*npw) 802 803 804 deallocate(p_terms,s_terms) 805 endif 806 807 if(gstart==2) opsi(1)=(0.d0,0.d0) 808 ! 809 deallocate(psi_r, psi_g,psi_v) 810 deallocate(h_psi_g,s_psi_g) 811 ! 812 RETURN 813 ! 814END SUBROUTINE o_1psi_gamma 815 816 817SUBROUTINE evc_to_real(numv, v_states) 818!this subroutine fourier transform states from evc 819!to real space 820 821 USE io_global, ONLY : stdout, ionode, ionode_id 822 USE kinds, ONLY : DP 823 USE wvfct, ONLY : npwx, npw, nbnd 824 USE wavefunctions, ONLY : evc, psic 825 USE gvecs, ONLY : doublegrid 826 USE fft_base, ONLY : dfftp, dffts 827 USE fft_interfaces, ONLY : fwfft, invfft 828 USE klist, ONLY : igk_k 829 830 implicit none 831 832 INTEGER, INTENT(in) :: numv!number of states to be transformed 833 REAL(kind=DP), INTENT(out) :: v_states(dffts%nnr,numv)!target arrsy 834 835 INTEGER :: iv 836 837 do iv=1,numv,2 838 psic(:)=(0.d0,0.d0) 839 if(iv < numv) then 840 psic(dffts%nl(igk_k(1:npw,1))) = evc(1:npw,iv)+(0.d0,1.d0)*evc(1:npw,iv+1) 841 psic(dffts%nlm(igk_k(1:npw,1))) = CONJG( evc(1:npw,iv) )+(0.d0,1.d0)*CONJG(evc(1:npw,iv+1)) 842 else 843 psic(dffts%nl(igk_k(1:npw,1))) = evc(1:npw,iv) 844 psic(dffts%nlm(igk_k(1:npw,1))) = CONJG( evc(1:npw,iv) ) 845 endif 846 CALL invfft ('Wave', psic, dffts) 847 if(iv<numv) then 848 v_states(1:dffts%nnr,iv)=dble(psic(1:dffts%nnr)) 849 v_states(1:dffts%nnr,iv+1)=dimag(psic(1:dffts%nnr)) 850 else 851 v_states(1:dffts%nnr,iv)=dble(psic(1:dffts%nnr)) 852 endif 853 854 enddo 855return 856 857END SUBROUTINE evc_to_real 858 859 860SUBROUTINE o_basis_init(numpw,o_basis,numv,v_states,cutoff, ptype,fcw_number,fcw_state,fcw_mat,ethr) 861!this subroutine initializes the polarization basis to 862!random values and diagonalize with respect to the O matrix 863 864 865 USE gvect, ONLY : g 866 USE io_global, ONLY : stdout, ionode, ionode_id 867 USE kinds, ONLY : DP 868 USE wvfct, ONLY : g2kin, npwx, npw, nbnd 869 USE wavefunctions, ONLY : evc, psic 870 USE gvecs, ONLY : doublegrid 871 USE constants, ONLY : tpi 872 USE random_numbers, ONLY : randy 873 USE klist, ONLY : xk,igk_k 874 USE cell_base, ONLY : tpiba2 875 USE fft_base, ONLY : dfftp, dffts 876 USE fft_interfaces, ONLY : fwfft, invfft 877 878 879 implicit none 880 881 INTEGER, INTENT(in) :: numv!number of states to be transformed 882 REAL(kind=DP), INTENT(in) :: v_states(dffts%nnr,numv)!valence states in real space 883 INTEGER, INTENT(in) :: numpw!dimesion of the polarization basis 884 COMPLEX(kind=DP), INTENT(out) :: o_basis(npw,numpw)!polarization basis 885 REAL(kind=DP), INTENT(in) :: cutoff!cutoff for planewaves 886 INTEGER, INTENT(in) :: ptype!type of approximation for O operator 887 INTEGER, INTENT(in) :: fcw_number!number of "fake conduction" states for O matrix method 888 COMPLEX(kind=DP) :: fcw_state(npw,fcw_number)! "fake conduction" states for O matrix method 889 REAL(kind=DP) :: fcw_mat(fcw_number,fcw_number)! "fake conduction" matrix 890 REAL(kind=DP), INTENT(in) :: ethr!threshold o_1psi_gamma 891 892 893 INTEGER :: iw,ig 894 REAL(kind=DP) :: rr, arg 895 REAL(kind=DP), ALLOCATABLE :: e(:) 896 REAL(kind=DP), ALLOCATABLE :: hdiag(:) 897 898 899 900 allocate(hdiag(npw)) 901 902 g2kin(1:npw) = ( (g(1,igk_k(1:npw,1)) )**2 + & 903 ( g(2,igk_k(1:npw,1)) )**2 + & 904 ( g(3,igk_k(1:npw,1)) )**2 ) * tpiba2 905 906 907 do ig=1,npw 908 !if(g2kin(ig) <= 3.d0) then 909 if(g2kin(ig) <= cutoff) then 910 hdiag(ig)=1.d0!/(1.d0+g2kin(ig)) 911 else 912 hdiag(ig)=0.d0 913 endif 914 enddo 915 hdiag(:)=1.d0 916 917 918 allocate( e(numpw)) 919 920 DO iw = 1, numpw 921 ! 922 DO ig = 1, npw 923 ! 924 rr = randy() 925 926 arg = tpi * randy() 927 ! 928 o_basis(ig,iw) = CMPLX( rr*COS( arg ), rr*SIN( arg ) ) / & 929 ( g(1,igk_k(ig,1))**2 + g(2,igk_k(ig,1))**2 + g(3,igk_k(ig,1))**2 + 1.D0) 930 ! 931 END DO 932 ! 933 END DO 934 935 !CALL rinitcgg( npwx, npw, n_starting_wfc, & 936 ! nbnd, wfcatom, wfcatom, etatom ) 937 938 939 940 call o_rinitcgg( npwx, npw, numpw, numpw, o_basis, o_basis, e, numv, v_states,hdiag,ptype,fcw_number,fcw_state,fcw_mat,ethr) 941 do iw=1,numpw 942 write(stdout,*) 'E', iw, e(iw) 943 enddo 944 945 deallocate(e) 946 deallocate(hdiag) 947 return 948 949 END SUBROUTINE o_basis_init 950 951 952 SUBROUTINE o_basis_write(numpw, o_basis,lcutoff,ecutoff, l_pbc) 953 954!this subroutines writes the basis of the polarization on file 955 USE io_global, ONLY : stdout, ionode, ionode_id 956 USE kinds, ONLY : DP 957 USE wvfct, ONLY : npwx, npw, nbnd 958 USE wavefunctions, ONLY : evc, psic 959 USE gvecs, ONLY : doublegrid 960 USE wavefunctions, ONLY : psic 961 USE io_files, ONLY : diropn, prefix 962 USE gvect, ONLY : ngm, gg,gstart 963 USE cell_base, ONLY: tpiba2 964 USE wannier_gw, ONLY : max_ngm 965 USE fft_base, ONLY : dfftp, dffts 966 USE fft_interfaces, ONLY : fwfft, invfft 967 USE klist, ONLY : igk_k 968 969 implicit none 970 971 INTEGER, EXTERNAL :: find_free_unit 972 973 INTEGER, INTENT(inout) :: numpw!dimension of the polarization basis 974 COMPLEX(kind=DP), INTENT(in) :: o_basis(npw,numpw) 975 REAL(kind=DP), INTENT(in) :: ecutoff!cutoff in Rydberg for g sum 976 LOGICAL, INTENT(in) :: lcutoff !if true uses cutoff on G defined by ecutoff 977 LOGICAL, INTENT(in) :: l_pbc!if true PBC are assumed and element G=0 is added as the last one 978 979 INTEGER :: iw, iungprod, ig,ngm_max 980 LOGICAL :: exst 981 COMPLEX(kind=DP), ALLOCATABLE :: tmp_g(:) 982 983 allocate(tmp_g(max_ngm)) 984 985 if(lcutoff) then 986 ngm_max=0 987 do ig=1,ngm 988 if(gg(ig)*tpiba2 >= ecutoff) exit 989 ngm_max=ngm_max+1 990 enddo 991 else 992 ngm_max=ngm 993 endif 994 995 write(stdout,*) 'NGM MAX:', ngm_max, ngm 996 997 iungprod = find_free_unit() 998 CALL diropn( iungprod, 'wiwjwfc_red', max_ngm*2, exst ) 999 1000 1001 1002 do iw=1,numpw 1003 1004 psic(:)=(0.d0,0.d0) 1005 psic(dffts%nl(igk_k(1:npw,1))) = o_basis(1:npw,iw) 1006 psic(dffts%nlm(igk_k(1:npw,1))) = CONJG( o_basis(1:npw,iw)) 1007 tmp_g(1:max_ngm)=psic(dfftp%nl(1:max_ngm)) 1008 if(gstart==2) tmp_g(1)=(0.d0,0.d0) 1009 CALL davcio(tmp_g, max_ngm*2,iungprod,iw,1) 1010 enddo 1011 1012 if(l_pbc) then 1013 numpw=numpw+1 1014 tmp_g(1:max_ngm)=(0.d0,0.d0) 1015 if(gstart==2) tmp_g(1)=(1.d0,0.d0) 1016 CALL davcio(tmp_g, max_ngm*2,iungprod,numpw,1) 1017 endif 1018 close(iungprod) 1019 deallocate(tmp_g) 1020 1021 return 1022 1023 END SUBROUTINE o_basis_write 1024 1025 1026!---------------------------------------------------------------------------- 1027SUBROUTINE o_1psi_gamma_real( numv, v_states, psi, opsi) 1028 !---------------------------------------------------------------------------- 1029 ! 1030 ! 1031 !this subroutines applies the O oprator to a state psi 1032 !IT WORKS ONLY FOR NORMCONSERVING PSEUDOPOTENTIALS 1033 !the valence states in G space must be in evc 1034 ! Gamma point version in real space 1035 1036 USE io_global, ONLY : stdout, ionode, ionode_id 1037 USE kinds, ONLY : DP 1038 USE wannier_gw 1039 USE gvect 1040 USE constants, ONLY : e2, pi, tpi, fpi 1041 USE cell_base, ONLY: at, alat, tpiba, omega, tpiba2 1042 USE wvfct, ONLY : npwx, npw, nbnd 1043 USE wavefunctions, ONLY : evc, psic 1044 USE mp, ONLY : mp_sum, mp_barrier, mp_bcast 1045 USE mp_world, ONLY : mpime, world_comm 1046 USE gvecs, ONLY : doublegrid 1047 USE fft_base, ONLY : dfftp, dffts 1048 USE fft_interfaces, ONLY : fwfft, invfft 1049 USE klist, ONLY : igk_k 1050 1051 USE kinds, ONLY : DP 1052 ! 1053 IMPLICIT NONE 1054 1055 INTEGER, INTENT(in) :: numv!number of valence states 1056 REAL(kind=DP), INTENT(in) :: v_states(dffts%nnr,numv)!valence states in real space 1057 COMPLEX(kind=DP), INTENT(in) :: psi(npw)!input wavefunction 1058 COMPLEX(kind=DP), INTENT(out) :: opsi(npw)!O|\psi> 1059 1060 REAL(kind=DP), ALLOCATABLE :: psi_r(:), psi_v(:), psi_w(:) 1061 COMPLEX(kind=DP), ALLOCATABLE :: psi_g(:) 1062 REAL(kind=DP), ALLOCATABLE :: prod(:) 1063 INTEGER :: iv 1064 1065 allocate(psi_r(dffts%nnr),psi_g(npw),psi_v(dffts%nnr)) 1066 allocate(prod(numv),psi_w(dffts%nnr)) 1067 1068!fourier transform psi to R space 1069 1070 opsi(1:npw)=(0.d0,0.d0) 1071 1072 psic(:)=(0.d0,0.d0) 1073 psic(dffts%nl(igk_k(1:npw,1))) = psi(1:npw) 1074 psic(dffts%nlm(igk_k(1:npw,1))) = CONJG( psi(1:npw) ) 1075 CALL invfft ('Wave', psic, dffts) 1076 psi_v(1:dffts%nnr)= DBLE(psic(1:dffts%nnr)) 1077 psi_w(:)=0.d0 1078!loop on v 1079 do iv=1,numv 1080 1081!!product with psi_v 1082 psi_r(1:dffts%nnr)=psi_v(1:dffts%nnr)*v_states(1:dffts%nnr, iv) 1083!!project on conduction manifold 1084 call dgemm('T','N',numv,1,dffts%nnr,1.d0,v_states,dffts%nnr,psi_r,dffts%nnr,0.d0,prod,numv) 1085 call mp_sum(prod(1:numv), world_comm) 1086 prod(:)=prod(:)/dble(dffts%nr1*dffts%nr2*dffts%nr3) 1087 call dgemm('N','N',dffts%nnr,1,numv,-1.d0,v_states,dffts%nnr,prod,numv,1.d0, psi_r,dffts%nnr) 1088! 1089 1090 1091 psi_r(1:dffts%nnr)=psi_r(1:dffts%nnr)*v_states(1:dffts%nnr,iv) 1092 1093!add up result (with sign) 1094 psi_w(:)=psi_w(:)-psi_r(:) 1095 1096 1097 1098 enddo 1099!!fourier transform in G space 1100 1101 1102 psic(:)=cmplx(psi_w(:),0.d0) 1103 CALL fwfft ('Wave', psic, dffts) 1104 opsi(1:npw)=psic(dffts%nl(igk_k(1:npw,1))) 1105 1106 1107 if(gstart==2) opsi(1)=dble(opsi(1)) 1108 ! 1109 deallocate(psi_r, psi_g,psi_v) 1110 deallocate(prod,psi_w) 1111 ! 1112 RETURN 1113 ! 1114END SUBROUTINE o_1psi_gamma_real 1115 1116 1117 1118 SUBROUTINE o_basis_test(numv,v_states,numpw, lcutoff,ecutoff) 1119 1120!this subroutines writes the basis of the polarization on file 1121 USE io_global, ONLY : stdout, ionode, ionode_id 1122 USE kinds, ONLY : DP 1123 USE wvfct, ONLY : npwx, npw, nbnd 1124 USE wavefunctions, ONLY : evc, psic 1125 USE gvecs, ONLY : doublegrid 1126 USE wavefunctions, ONLY : psic 1127 USE io_files, ONLY : prefix, tmp_dir, diropn 1128 USE gvect, ONLY : ngm, gg,gstart 1129 USE cell_base, ONLY: tpiba2 1130 USE wannier_gw, ONLY : max_ngm 1131 USE mp, ONLY : mp_sum 1132 USE mp_world, ONLY : world_comm 1133 USE fft_base, ONLY : dfftp, dffts 1134 USE fft_interfaces, ONLY : fwfft, invfft 1135 USE klist, ONLY : igk_k 1136 1137 1138 1139 implicit none 1140 1141 INTEGER, EXTERNAL :: find_free_unit 1142 INTEGER, INTENT(in) :: numv!number of valence states 1143 REAL(kind=DP), INTENT(in) :: v_states(dffts%nnr,numv)!valence states in real space 1144 INTEGER, INTENT(in) :: numpw!dimension of the polarization basis 1145 REAL(kind=DP), INTENT(in) :: ecutoff!cutoff in Rydberg for g sum 1146 LOGICAL, INTENT(in) :: lcutoff !if true uses cutoff on G defined by ecutoff 1147 1148 INTEGER :: iw, iungprod, ig,ngm_max 1149 LOGICAL :: exst 1150 COMPLEX(kind=DP), ALLOCATABLE :: tmp_g(:), psi1(:),psi2(:) 1151 REAL(kind=DP) :: sca 1152 1153 allocate(tmp_g(max_ngm),psi1(npw),psi2(npw)) 1154 1155 if(lcutoff) then 1156 ngm_max=0 1157 do ig=1,ngm 1158 if(gg(ig)*tpiba2 >= ecutoff) exit 1159 ngm_max=ngm_max+1 1160 enddo 1161 else 1162 ngm_max=ngm 1163 endif 1164 1165 write(stdout,*) 'NGM MAX:', ngm_max, ngm 1166 1167 iungprod = find_free_unit() 1168 CALL diropn( iungprod, 'wiwjwfc_red', max_ngm*2, exst ) 1169 1170 do iw=1,numpw 1171 CALL davcio(tmp_g, max_ngm*2,iungprod,iw,-1) 1172 psic(:)=(0.d0,0.d0) 1173 do ig=1,max_ngm 1174 psic(dffts%nl(ig))=tmp_g(ig) 1175 psic(dffts%nlm(ig))=conjg(tmp_g(ig)) 1176 enddo 1177 do ig=1,npw 1178 psi1(ig)=psic(dffts%nl(igk_k(ig,1))) 1179 enddo 1180 1181 call o_1psi_gamma_real( numv, v_states, psi1, psi2) 1182 sca=0.d0 1183 do ig=1,npw 1184 sca=sca+2.d0*dble(conjg(psi2(ig))*psi2(ig)) 1185 enddo 1186 if(gstart==2) sca=sca-dble(conjg(psi2(1))*psi2(1)) 1187 call mp_sum(sca,world_comm) 1188 sca=dsqrt(sca) 1189 psi2(:)=psi2(:)/sca 1190 sca=0.d0 1191 do ig=1,npw 1192 sca=sca+2.d0*dble(conjg(psi1(ig))*psi2(ig)) 1193 enddo 1194 if(gstart==2) sca=sca-dble(conjg(psi1(1))*psi2(1)) 1195 call mp_sum(sca,world_comm) 1196 write(stdout,*) 'o basis test:',iw,sca 1197 1198 enddo 1199 1200 1201 close(iungprod) 1202 deallocate(tmp_g,psi1,psi2) 1203 1204 return 1205 1206 END SUBROUTINE o_basis_test 1207 1208