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 energies_xc( lda, n, m, psi, e_xc, e_h,ispin ) 13 !---------------------------------------------------------------------------- 14 ! 15 ! computes the expectation values of the exchange and correlation potential 16 ! and of hartree potential 17 ! ... input: 18 ! ... lda leading dimension of arrays psi, spsi, hpsi 19 ! ... n true dimension of psi, spsi, hpsi 20 ! ... m number of states psi 21 ! ... psi 22 ! 23 ! ... output: 24 ! e_xc 25 ! e_h 26 USE kinds, ONLY : DP 27 USE uspp, ONLY : vkb, nkb 28 USE gvecs, ONLY : doublegrid 29 USE gvect, ONLY : ngm, gstart, g, gg, gcutm 30 USE cell_base, ONLY : alat, omega 31 USE lsda_mod, ONLY : nspin 32 USE ldaU, ONLY : lda_plus_u 33 USE lsda_mod, ONLY : current_spin 34 USE gvect, ONLY : gstart 35 USE io_global, ONLY :stdout 36 USE scf, ONLY : rho, vltot, vrs, rho_core,rhog_core, scf_type 37 USE constants, ONLY :rytoev 38 USE io_files, ONLY : diropn 39 USE mp, ONLY : mp_sum, mp_barrier 40 USE mp_world, ONLY : world_comm 41 USE control_flags, ONLY : gamma_only 42 USE funct, ONLY : dft_is_meta 43 USE fft_base, ONLY : dfftp, dffts 44 USE fft_interfaces, ONLY : fwfft, invfft, fft_interpolate 45 46 USE exx, ONLY : vexx !Suriano 47 USE funct, ONLY : exx_is_active,dft_is_hybrid 48 USE klist, ONLY : igk_k 49 50 ! 51 IMPLICIT NONE 52 INTEGER, EXTERNAL :: find_free_unit 53 ! 54 ! ... input/output arguments 55 ! 56 INTEGER :: lda, n, m 57 COMPLEX(DP) :: psi(lda,m) 58 REAL(kind=DP) :: e_xc(m), e_h(m) 59 INTEGER, INTENT(in) :: ispin !spin 1,2 60 61 REAL(kind=DP), ALLOCATABLE :: vr(:,:) 62 ! 63 ! 64 CALL start_clock( 'h_psi' ) 65 allocate(vr(dfftp%nnr,nspin)) 66 ! 67 IF ( gamma_only ) THEN 68 ! 69 CALL energies_xc_gamma() 70 ! 71 ELSE 72 ! 73 CALL energies_xc_k( ) 74 ! 75 END IF 76 77 ! 78 CALL stop_clock( 'h_psi' ) 79 deallocate(vr) 80 ! 81 RETURN 82 ! 83 CONTAINS 84 ! 85 !----------------------------------------------------------------------- 86 SUBROUTINE energies_xc_k( ) 87 !----------------------------------------------------------------------- 88 ! 89 ! ... k-points version 90 ! 91 USE wavefunctions, ONLY : psic 92 USE becmod, ONLY : becp 93 ! 94 IMPLICIT NONE 95 ! 96 97 INTEGER :: ibnd, j,is, ig 98 REAL(dp) :: etxc,vtxc 99 REAL(kind=DP) :: ehart, charge 100 ! counters 101 ! 102 ! 103 ! 104 ! ... Here we apply the kinetic energy (k+G)^2 psi 105 ! 106 ! 107 ! 108 ! ... Here we add the Hubbard potential times psi 109 ! 110 ! 111 ! ... the local potential V_Loc psi. First the psi in real space 112!set exchange and correlation potential 113 if(.not.allocated(psic)) write(stdout,*) 'psic not allocated' 114 ! 115 if (dft_is_meta()) then 116! call v_xc_meta( rho, rho_core, rhog_core, etxc, vtxc, v%of_r, v%kin_r ) 117 else 118 CALL v_xc( rho, rho_core, rhog_core, etxc, vtxc, vr ) 119 endif 120 ! 121 do is=1,nspin 122 vrs(:,is)=vr(:,is) 123 if(doublegrid) call fft_interpolate(dfftp, vrs(:,is),dffts,vrs(:,is)) ! interpolate from dense to smooth 124 enddo 125 ! 126 DO ibnd = 1, m 127 ! 128 CALL start_clock( 'firstfft' ) 129 ! 130 psic(1:dffts%nnr) = ( 0.D0, 0.D0 ) 131 ! 132 psic(dffts%nl(igk_k(1:n,1))) = psi(1:n,ibnd) 133 ! 134 CALL invfft ('Wave', psic, dffts) 135 136 ! 137 CALL stop_clock( 'firstfft' ) 138 ! 139 ! ... product with the potential vrs = (vltot+vr) on the smooth grid 140 ! 141 psic(1:dffts%nnr) = psic(1:dffts%nnr) * vrs(1:dffts%nnr,1) 142 ! 143 ! ... back to reciprocal space 144 ! 145 CALL start_clock( 'secondfft' ) 146 ! 147 CALL fwfft ('Wave', psic, dffts) 148 ! 149 ! ... addition to the total product 150 ! 151 e_xc(ibnd)=0.d0 152 do ig=1,n 153 e_xc(ibnd)=e_xc(ibnd)+real(conjg(psi(ig,ibnd))*psic(dffts%nl(igk_k(ig,1)))) 154 enddo 155 call mp_sum(e_xc(ibnd),world_comm) 156 write(stdout,*) 'energies_xc :', ibnd, e_xc(ibnd)*rytoev 157! 158 CALL stop_clock( 'secondfft' ) 159 ! 160 END DO 161 vr(:,:)=0.d0 162 call v_h(rho%of_g , ehart, charge, vr ) 163 do is=1,nspin 164 vrs(:,is)=vr(:,is) 165 if(doublegrid) call fft_interpolate(dfftp, vrs(:,is), dffts, vrs(:,is)) ! interpolate from dense to smooth 166 enddo 167 168 169 DO ibnd = 1, m 170 171 CALL start_clock( 'firstfft' ) 172 psic(1:dffts%nnr) = ( 0.D0, 0.D0 ) 173 psic(dffts%nl(igk_k(1:n,1))) = psi(1:n,ibnd) 174 175 CALL invfft ('Wave', psic, dffts) 176 177 178 CALL stop_clock( 'firstfft' ) 179 180 181 psic(1:dffts%nnr) = psic(1:dffts%nnr) * vrs(1:dffts%nnr,1) 182 183 CALL start_clock( 'secondfft' ) 184 CALL fwfft ('Wave', psic, dffts) 185 e_h(ibnd)=0.d0 186 do ig=1,n 187 e_h(ibnd)=e_h(ibnd)+real(conjg(psi(ig,ibnd))*psic(dffts%nl(igk_k(ig,1)))) 188 enddo 189 call mp_sum(e_h(ibnd),world_comm) 190 write(stdout,*) 'energies_h :', ibnd, e_h(ibnd)*rytoev 191 192 CALL stop_clock( 'secondfft' ) 193 194 enddo! 195 196 ! 197 ! 198 ! 199 RETURN 200 ! 201 END SUBROUTINE energies_xc_k 202 SUBROUTINE energies_xc_gamma 203 204 USE uspp, ONLY : okvan 205 USE wannier_gw, ONLY : becp_gw, restart_gww,l_whole_s,l_verbose,& 206 &l_scissor,scissor,num_nbndv,num_nbnds 207 ! USE realus, ONLY : adduspos_gamma_r 208 USE wvfct, ONLY : npwx,npw,nbnd, et,g2kin 209 USE wavefunctions, ONLY : evc 210 USE klist, ONLY : xk 211 USE mp, ONLY : mp_sum 212 USE mp_world, ONLY : world_comm 213 USE gvect, ONLY : gstart,g 214 USE constants, ONLY : rytoev 215 USE becmod, ONLY : becp, calbec,allocate_bec_type,deallocate_bec_type 216 USE cell_base, ONLY : tpiba2 217 USE io_global, ONLY : ionode 218 USE io_files, ONLY :prefix,tmp_dir 219 USE exx, ONLY : exxalfa 220 221 implicit none 222 223 INTEGER, EXTERNAL :: find_free_unit 224 225 INTEGER :: ibnd,jbnd,ir,ig 226 INTEGER :: iunwfcreal,iunu 227 REAL(kind=DP) :: etxc,vtxc,ehart, charge 228 REAL(kind=DP), ALLOCATABLE :: psi_r(:),psi_rs(:) 229 LOGICAL :: exst 230 REAL(kind=DP), ALLOCATABLE :: rho_fake_core(:) 231 COMPLEX(kind=DP), ALLOCATABLE :: hpsi(:,:) 232 REAL(kind=DP), ALLOCATABLE :: exact_x(:) 233 REAL(kind=DP), ALLOCATABLE :: e_hub(:)!Hubbard contribution to KS energies 234 REAL(kind=DP), ALLOCATABLE :: et_off(:,:)!off-diagonal energies 235 236 237 allocate(exact_x(nbnd)) 238 allocate(e_hub(nbnd)) 239 if(l_whole_s) then 240 allocate (et_off(nbnd,nbnd)) 241 endif 242!if required calculates also the KS energies 243 ! if(restart_gww==-1) then 244 if(l_verbose) write(stdout,*) 'ATTENZIONE1' 245 FLUSH(stdout) 246 !allocate( becp%r( nkb, nbnd ) ) 247 call allocate_bec_type ( nkb, nbnd, becp) 248 if(l_verbose) write(stdout,*) 'ATTENZIONE2' 249 FLUSH(stdout) 250 251 IF ( nkb > 0 ) CALL init_us_2( npw, igk_k(1,1), xk(1,1), vkb ) 252 !call ccalbec( nkb, npwx, npw, nbnd, becp%r, vkb, evc ) 253 !if(nkb> 0)CALL calbec ( npw, vkb, psi, becp, nbnd) 254 255 if(l_verbose)write(stdout,*) 'ATTENZIONE3' 256 FLUSH(stdout) 257 258 allocate(hpsi(npwx,nbnd)) 259 if(l_verbose)write(stdout,*) 'ATTENZIONE4' 260 FLUSH(stdout) 261 262 g2kin(1:npw) = ( ( xk(1,1) + g(1,igk_k(1:npw,1)) )**2 + & 263 ( xk(2,1) + g(2,igk_k(1:npw,1)) )**2 + & 264 ( xk(3,1) + g(3,igk_k(1:npw,1)) )**2 ) * tpiba2 265 266 267 if(l_verbose)write(stdout,*) 'ATTENZIONE5' 268 FLUSH(stdout) 269 270 271 ! exxalfa=0.d0!ATTENZIONE 272 call h_psi( npwx, npw, nbnd, psi, hpsi ) 273 et(:,ispin)=0.d0 274 if(l_verbose)write(stdout,*) 'ATTENZIONE6' 275 if(l_verbose)write(stdout,*) 'EXXALFA', exxalfa 276 FLUSH(stdout) 277 278 do ibnd=1,nbnd 279 280 ! call dgemm('T','N',1,1,2*npw,2.d0,evc(:,ibnd),2*npwx,hpsi(:,ibnd),2*npwx,& 281 ! &0.d0,et(ibnd,1),1) 282 do ig=1,npw 283 et(ibnd,ispin)=et(ibnd,ispin)+2.d0*dble(conjg(evc(ig,ibnd))*hpsi(ig,ibnd)) 284 enddo 285 if(gstart==2) then 286 et(ibnd,ispin)=et(ibnd,ispin)-dble(conjg(evc(1,ibnd))*hpsi(1,ibnd)) 287 endif 288 enddo 289 call mp_sum(et(:,ispin),world_comm) 290 if(l_scissor) then 291 et(1:num_nbndv(ispin),ispin)=et(1:num_nbndv(ispin),ispin)+scissor(1)/rytoev 292 et(num_nbndv(ispin)+1:num_nbnds,ispin)=et(num_nbndv(ispin)+1:num_nbnds,ispin)+scissor(2)/rytoev 293 endif 294 295 if(l_verbose)write(stdout,*) 'ATTENZIONE7' 296 FLUSH(stdout) 297!if required calculate Hubbard U contribution to eigen-energies 298 e_hub(:)=0.d0 299 if ( lda_plus_u ) then 300 hpsi(:,:)=(0.d0,0.d0) 301 CALL vhpsi( npwx, npw, nbnd, psi, hpsi ) 302 303 do ibnd=1,nbnd 304 305 do ig=1,npw 306 e_hub(ibnd)=e_hub(ibnd)+2.d0*dble(conjg(psi(ig,ibnd))*hpsi(ig,ibnd)) 307 enddo 308 if(gstart==2) then 309 e_hub(ibnd)=e_hub(ibnd)-dble(conjg(psi(1,ibnd))*hpsi(1,ibnd)) 310 endif 311 enddo 312 call mp_sum(e_hub(:),world_comm) 313 do ibnd=1,nbnd 314 write(stdout,*) 'Hubbard U energy:',ibnd,e_hub(ibnd)*rytoev 315 enddo 316 FLUSH(stdout) 317 318 endif 319 do ibnd=1,nbnd 320 write(stdout,*) 'KS energy:',ibnd,et(ibnd,ispin)*rytoev 321 enddo 322 FLUSH(stdout) 323 324!in case of hybrid functionals and HF we have to calculated also the exact exchange part 325 326 if(dft_is_hybrid()) then 327!NOT_TO_BE_INCLUDED_START 328 hpsi(:,:)=(0.d0,0.d0) 329 call vexx( npwx, npw, nbnd, psi, hpsi ) 330 do ibnd=1,nbnd 331 call dgemm('T','N',1,1,2*npw,2.d0,evc(:,ibnd),2*npwx,hpsi(:,ibnd),2*npwx,& 332 &0.d0,exact_x(ibnd),1) 333 if(gstart==2) then 334 exact_x(ibnd)=exact_x(ibnd)-dble(conjg(evc(1,ibnd))*hpsi(1,ibnd)) 335 endif 336 call mp_sum(exact_x(ibnd),world_comm) 337 write(stdout,*) 'Exact exchange :',ibnd, exact_x(ibnd) 338 enddo 339!NOT_TO_BE_INCLUDED_END 340 end if 341 342 ! deallocate(hpsi,becp%r) 343 call deallocate_bec_type ( becp) 344! endif 345 346 347 348 allocate(psi_r(dfftp%nnr),psi_rs(dfftp%nnr)) 349 350 iunwfcreal=find_free_unit() 351 CALL diropn( iunwfcreal, 'real_whole', dffts%nnr, exst ) 352 353 354!calculate xc potential on fine grid 355 356 357 if(.not.allocated(vr)) write(stdout,*) 'vr not allocated' 358 allocate(rho_fake_core(dfftp%nnr)) 359 rho_fake_core(:)=0.d0 360 ! 361 if (dft_is_meta()) then 362 ! call v_xc_meta( rho, rho_core, rhog_core, etxc, vtxc, v%of_r, v%kin_r ) 363 else 364 CALL v_xc( rho, rho_core, rhog_core, etxc, vtxc, vr ) 365 endif 366 ! 367 deallocate(rho_fake_core) 368 369 370 if(l_whole_s) then 371!NOT_TO_BE_INCLUDED_START 372 allocate(hpsi(npwx,nbnd)) 373 hpsi(:,:)=(0.d0,0.d0) 374 CALL vloc_psi_gamma ( npwx, npw, nbnd, evc, vr(1,ispin), hpsi ) 375 call dgemm('T','N',nbnd,nbnd,2*npw,2.d0,evc,2*npwx,hpsi,2*npwx,& 376 &0.d0,et_off,nbnd) 377 if(gstart==2) then 378 do ibnd=1,nbnd 379 do jbnd=1,nbnd 380 et_off(ibnd,jbnd)=et_off(ibnd,jbnd)-dble(conjg(evc(1,ibnd))*hpsi(1,jbnd)) 381 enddo 382 enddo 383 endif 384 deallocate(hpsi) 385 call mp_sum(et_off,world_comm) 386!write on file 387 if(ionode) then 388 iunu = find_free_unit() 389 if(ispin==1) then 390 open(unit=iunu,file=trim(tmp_dir)//trim(prefix)//'.exc_off',status='unknown',form='unformatted') 391 else 392 open(unit=iunu,file=trim(tmp_dir)//trim(prefix)//'.exc_off2',status='unknown',form='unformatted') 393 endif 394 write(iunu) nbnd 395 do ibnd=1,nbnd 396 write(iunu) et_off(1:nbnd,ibnd) 397 enddo 398 close(iunu) 399 endif 400!NOT_TO_BE_INCLUDED_END 401 endif 402 403 404 405 406 do ibnd=1,m!loop on states 407!read from disk wfc on coarse grid 408 CALL davcio( psi_rs,dffts%nnr,iunwfcreal,ibnd+(ispin-1)*nbnd,-1) 409 if(doublegrid) then 410 call fft_interpolate(dffts, psi_rs, dfftp, psi_r) ! interpolate from smooth to dense 411 else 412 psi_r(:)=psi_rs(:) 413 endif 414 415 do ir=1,dfftp%nnr 416 psi_r(ir)=psi_r(ir)**2.d0 417 enddo 418 419 !if(okvan) call adduspos_gamma_r(ibnd,ibnd,psi_r,1,becp_gw(:,ibnd),becp_gw(:,ibnd)) 420 421 e_xc(ibnd)=0.d0 422 do ir=1,dfftp%nnr 423 e_xc(ibnd)=e_xc(ibnd)+psi_r(ir)*vr(ir,ispin)!the 1 is for the spin NOT IMPLEMENTED YET 424 enddo 425 e_xc(ibnd)=e_xc(ibnd)/dble(dfftp%nr1*dfftp%nr2*dfftp%nr3) 426 427 call mp_sum(e_xc(ibnd),world_comm) 428 429!ifrequired add the contribution from exact exchange for hybrids and HF 430 if(dft_is_hybrid()) then 431!NOT_TO_BE_INCLUDED_START 432 e_xc(ibnd)=e_xc(ibnd)+exact_x(ibnd) 433!NOT_TO_BE_INCLUDED_END 434 endif 435 436 write(stdout,*) 'Routine energies_xc :', ibnd, e_xc(ibnd)*rytoev 437 438!now hartree term 439 440 441 enddo 442 443!if required add to e_xc Hubbard U terms 444 445 if(lda_plus_u) then 446 447 e_xc(1:nbnd)=e_xc(1:nbnd)+e_hub(1:nbnd) 448 endif 449 450 451 vr(:,:)=0.d0 452 453 call v_h(rho%of_g , ehart, charge, vr ) 454 455 456 do ibnd=1,m!loop on states 457!read from disk wfc on coarse grid 458 CALL davcio( psi_rs,dffts%nnr,iunwfcreal,ibnd+(ispin-1)*nbnd,-1) 459 if(doublegrid) then 460 call fft_interpolate(dffts, psi_rs, dfftp, psi_r) ! interpolate from smooth to dense 461 else 462 psi_r(:)=psi_rs(:) 463 endif 464 465 do ir=1,dfftp%nnr 466 psi_r(ir)=psi_r(ir)**2.d0 467 enddo 468 469 !if(okvan) call adduspos_gamma_r(ibnd,ibnd,psi_r,1,becp_gw(:,ibnd),becp_gw(:,ibnd)) 470 471 e_h(ibnd)=0.d0 472 do ir=1,dfftp%nnr 473 e_h(ibnd)=e_h(ibnd)+psi_r(ir)*vr(ir,ispin) 474 enddo 475 e_h(ibnd)=e_h(ibnd)/dble(dfftp%nr1*dfftp%nr2*dfftp%nr3) 476 477 call mp_sum(e_h(ibnd),world_comm) 478 write(stdout,*) 'Routine energies_h :', ibnd, e_h(ibnd)*rytoev 479 480!now hartree term 481 482 483 enddo 484 485 486 487 488 489 deallocate(psi_r,psi_rs) 490 deallocate(exact_x) 491 close(iunwfcreal) 492 deallocate(e_hub) 493 if(l_whole_s) then 494!NOT_TO_BE_INCLUDED_START 495 deallocate(et_off) 496!NOT_TO_BE_INCLUDED_END 497 endif 498 499 return 500 501 END SUBROUTINE energies_xc_gamma 502 ! 503END SUBROUTINE energies_xc 504 505 506SUBROUTINE write_energies_xc(e_xc) 507 508 USE kinds, ONLY : DP 509 USE wannier_gw, ONLY : num_nbnds, l_verbose 510 USE io_files, ONLY : prefix,tmp_dir 511 USE io_global, ONLY : ionode 512 USE wvfct, ONLY : nbnd 513 USE lsda_mod, ONLY : nspin 514 515 implicit none 516 517 INTEGER, EXTERNAL :: find_free_unit 518 519 REAL(kind=DP) :: e_xc(nbnd,nspin)!exchange and correlation energies 520 521 INTEGER :: iunu, iw 522 523 524 if(ionode) then 525 iunu = find_free_unit() 526 527 open(unit=iunu,file=trim(tmp_dir)//trim(prefix)//'.dft_xc',status='unknown',form='unformatted') 528 529 write(iunu) nbnd 530 531 532 do iw=1,nbnd 533 write(iunu) e_xc(iw,1) 534 if(l_verbose) WRITE(*,*) 'SCRITTO e_XC 1', e_xc(iw,1) 535 enddo 536 if(nspin==2) then 537 do iw=1,nbnd 538 write(iunu) e_xc(iw,2) 539 if(l_verbose) WRITE(*,*) 'SCRITTO e_XC 2', e_xc(iw,2) 540 enddo 541 endif 542 close(iunu) 543 endif 544 545 return 546 547END SUBROUTINE write_energies_xc 548