1! 2! Copyright (C) 2002-2008 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! ... Gradient Correction & exchange and correlation 10!=----------------------------------------------------------------------------=! 11 12 subroutine exch_corr_h( nspin, rhog, rhor, rhoc, sfac, exc, dxc, self_exc ) 13! 14! calculate exch-corr potential, energy, and derivatives dxc(i,j) 15! of e(xc) with respect to to cell parameter h(i,j) 16! 17 use funct, only : dft_is_gradient, dft_is_meta 18 use fft_base, only : dfftp, dffts 19 use cell_base, only : ainv, omega, h 20 use ions_base, only : nsp 21 use control_flags, only : tpre, iverbosity 22 use core, only : drhocg 23 use gvect, only : g 24 use uspp, only : nlcc_any 25 use mp, only : mp_sum 26 use metagga_cp, ONLY : kedtaur 27 USE io_global, ONLY : stdout 28 USE mp_global, ONLY : intra_bgrp_comm 29 use kinds, ONLY : DP 30 use constants, ONLY : au_gpa 31 USE sic_module, ONLY : self_interaction, sic_alpha 32 USE cp_interfaces, ONLY : denlcc 33 use cp_main_variables, only : drhor 34! 35 implicit none 36 37 ! input 38 ! 39 integer nspin 40 ! 41 ! rhog contains the charge density in G space 42 ! rhor contains the charge density in R space 43 ! 44 complex(DP) :: rhog( dfftp%ngm, nspin ) 45 complex(DP) :: sfac( dffts%ngm, nsp ) 46 ! 47 ! output 48 ! rhor contains the exchange-correlation potential 49 ! 50 real(DP) :: rhor( dfftp%nnr, nspin ), rhoc( dfftp%nnr ) 51 real(DP) :: dxc( 3, 3 ), exc 52 real(DP) :: dcc( 3, 3 ), drc( 3, 3 ) 53 ! 54 ! local 55 ! 56 integer :: i, j, ir, iss 57 real(DP) :: dexc(3,3) 58 real(DP), allocatable :: gradr(:,:,:) 59 ! 60 !sic 61 REAL(DP) :: self_exc 62 REAL(DP), ALLOCATABLE :: self_rho( :,: ), self_gradr(:,:,:) 63 complex(DP), ALLOCATABLE :: self_rhog( :,: ) 64 LOGICAL :: ttsic 65 real(DP) :: detmp(3,3) 66 ! 67 ! filling of gradr with the gradient of rho using fft's 68 ! 69 if ( dft_is_gradient() ) then 70 ! 71 allocate( gradr( 3, dfftp%nnr, nspin ) ) 72 do iss = 1, nspin 73 CALL fft_gradient_g2r ( dfftp, rhog(1,iss), g, gradr(1,1,iss) ) 74 end do 75 ! 76 else 77 ! 78 allocate( gradr( 3, 1, 2 ) ) 79 ! 80 end if 81 82 83 ttsic = (self_interaction /= 0 ) 84 ! 85 IF ( ttsic ) THEN 86 ! 87 IF ( dft_is_meta() ) CALL errore ('exch_corr_h', & 88 'SIC and meta-GGA not together', 1) 89 IF ( tpre ) CALL errore( 'exch_corr_h', 'SIC and stress not implemented', 1) 90 91 ! allocate the sic_arrays 92 ! 93 ALLOCATE( self_rho( dfftp%nnr, nspin ) ) 94 ALLOCATE( self_rhog(dfftp%ngm, nspin ) ) 95 96 self_rho(:, :) = rhor( :, :) 97 IF( dft_is_gradient() ) THEN 98 ALLOCATE( self_gradr( 3, dfftp%nnr, nspin ) ) 99 self_gradr(:, :, :) = gradr(:, :, :) 100 ENDIF 101 self_rhog(:, :) = rhog( :, :) 102 ! 103 END IF 104! 105 self_exc = 0.d0 106! 107 if( dft_is_meta() ) then 108 ! 109 call tpssmeta( dfftp%nnr, nspin, gradr, rhor, kedtaur, exc ) 110 ! 111 else 112 ! 113 CALL exch_corr_cp(dfftp%nnr, nspin, gradr, rhor, exc) 114 ! 115 IF ( ttsic ) THEN 116 CALL exch_corr_cp(dfftp%nnr, nspin, self_gradr, self_rho, self_exc) 117 self_exc = sic_alpha * (exc - self_exc) 118 exc = exc - self_exc 119 END IF 120! 121 end if 122 123 call mp_sum( exc, intra_bgrp_comm ) 124 IF ( ttsic ) call mp_sum( self_exc, intra_bgrp_comm ) 125 126 exc = exc * omega / DBLE( dfftp%nr1 * dfftp%nr2 * dfftp%nr3 ) 127 IF ( ttsic ) self_exc = self_exc * omega/DBLE(dfftp%nr1 * dfftp%nr2 *dfftp%nr3 ) 128 129 ! WRITE(*,*) 'Debug: calcolo exc', exc, 'eself', self_exc 130 ! 131 ! exchange-correlation contribution to pressure 132 ! 133 dxc = 0.0d0 134 ! 135 if ( tpre ) then 136 ! 137 ! Add term: Vxc( r ) * Drhovan( r )_ij - Vxc( r ) * rho( r ) * ((H^-1)^t)_ij 138 ! 139 do iss = 1, nspin 140 do j=1,3 141 do i=1,3 142 do ir=1,dfftp%nnr 143 dxc(i,j) = dxc(i,j) + rhor( ir, iss ) * drhor( ir, iss, i, j ) 144 end do 145 end do 146 end do 147 end do 148 ! 149 dxc = dxc * omega / DBLE( dfftp%nr1*dfftp%nr2*dfftp%nr3 ) 150 ! 151 call mp_sum ( dxc, intra_bgrp_comm ) 152 ! 153 do j = 1, 3 154 do i = 1, 3 155 dxc( i, j ) = dxc( i, j ) + exc * ainv( j, i ) 156 end do 157 end do 158 ! 159 ! DEBUG 160 ! 161 ! write (stdout,*) "derivative of e(xc)" 162 ! write (stdout,5555) ((dxc(i,j),j=1,3),i=1,3) 163 ! 164 IF( iverbosity > 1 ) THEN 165 DO i=1,3 166 DO j=1,3 167 detmp(i,j)=exc*ainv(j,i) 168 END DO 169 END DO 170 WRITE( stdout,*) "derivative of e(xc) - diag - kbar" 171 detmp = -1.0d0 * MATMUL( detmp, TRANSPOSE( h ) ) / omega * au_gpa * 10.0d0 172 WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3) 173 END IF 174 ! 175 end if 176 ! 177 if (dft_is_gradient()) then 178 ! 179 ! Add second part of the xc-potential to rhor 180 ! Compute contribution to the stress dexc 181 ! 182 call gradh( nspin, gradr, rhog, rhor, dexc) 183 ! 184 if (tpre) then 185 ! 186 call mp_sum ( dexc, intra_bgrp_comm ) 187 ! 188 dxc = dxc + dexc 189 ! 190 end if 191 ! 192 end if 193 ! 194 195 IF( ttsic ) THEN 196! 197 IF (dft_is_gradient()) then 198 199 call gradh( nspin, self_gradr, self_rhog, self_rho, dexc) 200 201 gradr(:,:, 1) = (1.d0 - sic_alpha ) * gradr(:,:, 1) 202 gradr(:,:, 2) = (1.d0 - sic_alpha ) * gradr(:,:, 2) + & 203 & sic_alpha * ( self_gradr(:,:,1) + self_gradr(:,:,2) ) 204 ENDIF 205 206 rhor(:, 1) = (1.d0 - sic_alpha ) * rhor(:, 1) 207 rhor(:, 2) = (1.d0 - sic_alpha ) * rhor(:, 2) + & 208 & sic_alpha * ( self_rho(:,1) + self_rho(:,2) ) 209 210 IF(ALLOCATED(self_gradr)) DEALLOCATE(self_gradr) 211 IF(ALLOCATED(self_rhog)) DEALLOCATE(self_rhog) 212 IF(ALLOCATED(self_rho)) DEALLOCATE(self_rho) 213! 214 ENDIF 215 216 IF( tpre ) THEN 217 ! 218 dcc = 0.0d0 219 ! 220 IF( nlcc_any ) CALL denlcc( dfftp%nnr, nspin, rhor, sfac, drhocg, dcc ) 221 ! 222 ! DEBUG 223 ! 224 ! write (stdout,*) "derivative of e(xc) - nlcc part" 225 ! write (stdout,5555) ((dcc(i,j),j=1,3),i=1,3) 226 ! 227 dxc = dxc + dcc 228 ! 229 do iss = 1, nspin 230 drc = 0.0d0 231 IF( nlcc_any ) THEN 232 do j=1,3 233 do i=1,3 234 do ir=1,dfftp%nnr 235 drc(i,j) = drc(i,j) + rhor( ir, iss ) * rhoc( ir ) * ainv(j,i) 236 end do 237 end do 238 end do 239 call mp_sum ( drc, intra_bgrp_comm ) 240 END IF 241 dxc = dxc - drc * ( 1.0d0 / nspin ) * omega / ( dfftp%nr1*dfftp%nr2*dfftp%nr3 ) 242 end do 243 ! 244 END IF 245 ! 246 247 IF( ALLOCATED( gradr ) ) DEALLOCATE( gradr ) 248 2495555 format(1x,f12.5,1x,f12.5,1x,f12.5/ & 250 & 1x,f12.5,1x,f12.5,1x,f12.5/ & 251 & 1x,f12.5,1x,f12.5,1x,f12.5//) 252 ! 253 return 254 end subroutine exch_corr_h 255 256 257!=----------------------------------------------------------------------------=! 258 259 subroutine gradh( nspin, gradr, rhog, rhor, dexc ) 260! _________________________________________________________________ 261! 262! calculate the second part of gradient corrected xc potential 263! plus the gradient-correction contribution to pressure 264! 265 USE kinds, ONLY: DP 266 use control_flags, only: iprint, tpre 267 use gvect, only: g 268 use cell_base, only: ainv, tpiba, omega 269 use cp_main_variables, only: drhog 270 USE fft_interfaces, ONLY: fwfft, invfft 271 USE fft_base, ONLY: dfftp 272 USE fft_helper_subroutines, ONLY: fftx_threed2oned, fftx_oned2threed 273! 274 implicit none 275! input 276 integer nspin 277 real(DP) :: gradr( 3, dfftp%nnr, nspin ), rhor( dfftp%nnr, nspin ), dexc( 3, 3 ) 278 complex(DP) :: rhog( dfftp%ngm, nspin ) 279! 280 complex(DP), allocatable:: v(:), vp(:), vm(:) 281 complex(DP), allocatable:: x(:), vtemp(:) 282 complex(DP) :: ci, fp, fm 283 integer :: iss, ig, ir, i,j 284! 285 allocate(v(dfftp%nnr)) 286 allocate(x(dfftp%ngm)) 287 allocate(vp(dfftp%ngm)) 288 allocate(vm(dfftp%ngm)) 289 allocate(vtemp(dfftp%ngm)) 290 ! 291 ci=(0.0d0,1.0d0) 292 ! 293 dexc = 0.0d0 294 ! 295 do iss=1, nspin 296! _________________________________________________________________ 297! second part xc-potential: 3 forward ffts 298! 299 do ir=1,dfftp%nnr 300 v(ir)=CMPLX(gradr(1,ir,iss),0.d0,kind=DP) 301 end do 302 call fwfft('Rho',v, dfftp ) 303 CALL fftx_threed2oned( dfftp, v, vp ) 304 do ig=1,dfftp%ngm 305 x(ig)=ci*tpiba*g(1,ig)*vp(ig) 306 end do 307! 308 if(tpre) then 309 do i=1,3 310 do j=1,3 311 do ig=1,dfftp%ngm 312 vtemp(ig) = omega*ci*CONJG(vp(ig))* & 313 & tpiba*(-rhog(ig,iss)*g(i,ig)*ainv(j,1)+ & 314 & g(1,ig)*drhog(ig,iss,i,j)) 315 end do 316 dexc(i,j) = dexc(i,j) + DBLE(SUM(vtemp))*2.0d0 317 end do 318 end do 319 endif 320! 321 do ir=1,dfftp%nnr 322 v(ir)=CMPLX(gradr(2,ir,iss),gradr(3,ir,iss),kind=DP) 323 end do 324 call fwfft('Rho',v, dfftp ) 325 CALL fftx_threed2oned( dfftp, v, vp, vm ) 326! 327 do ig=1,dfftp%ngm 328 x(ig) = x(ig) + ci*tpiba*g(2,ig)*vp(ig) 329 x(ig) = x(ig) + ci*tpiba*g(3,ig)*vm(ig) 330 end do 331! 332 if(tpre) then 333 do i=1,3 334 do j=1,3 335 do ig=1,dfftp%ngm 336 vtemp(ig) = omega*ci*( & 337 & CONJG(vp(ig))*tpiba*(-rhog(ig,iss)*g(i,ig)*ainv(j,2)+g(2,ig)*drhog(ig,iss,i,j))+ & 338 & CONJG(vm(ig))*tpiba*(-rhog(ig,iss)*g(i,ig)*ainv(j,3)+g(3,ig)*drhog(ig,iss,i,j)) ) 339 end do 340 dexc(i,j) = dexc(i,j) + 2.0d0*DBLE(SUM(vtemp)) 341 end do 342 end do 343 endif 344! _________________________________________________________________ 345! second part xc-potential: 1 inverse fft 346! 347 CALL fftx_oned2threed( dfftp, v, x ) 348 call invfft('Rho',v, dfftp ) 349 do ir=1,dfftp%nnr 350 rhor(ir,iss)=rhor(ir,iss)-DBLE(v(ir)) 351 end do 352 end do 353! 354 deallocate(vtemp) 355 deallocate(x) 356 deallocate(v) 357 deallocate(vp) 358 deallocate(vm) 359! 360 return 361 end subroutine gradh 362 363!=----------------------------------------------------------------------------=! 364! 365! For CP we need a further small interface subroutine 366! 367!=----------------------------------------------------------------------------=! 368 369subroutine exch_corr_cp(nnr,nspin,grhor,rhor,etxc) 370 use kinds, only: DP 371 use funct, only: dft_is_gradient, get_igcc 372 use xc_lda_lsda, only: xc 373 use xc_gga, only: xc_gcx, change_threshold_gga 374 implicit none 375 integer, intent(in) :: nnr 376 integer, intent(in) :: nspin 377 real(DP) :: grhor( 3, nnr, nspin ) 378 real(DP) :: rhor( nnr, nspin ) 379 real(DP) :: etxc 380 381 real(DP), parameter :: epsr = 1.0d-10 382 real(DP), parameter :: e2=1.0_dp 383 integer :: ir, is, k, ipol, neg(3) 384 real(DP) :: grup, grdw 385 real(DP), allocatable :: v(:,:) 386 real(DP), allocatable :: h(:,:,:) 387 real(DP), allocatable :: rhox (:,:)!^ 388 real(DP), allocatable :: ex(:), ec(:) 389 real(DP), allocatable :: vx(:,:), vc(:,:) 390 REAL(DP), allocatable :: sx(:), sc(:) 391 REAL(DP), allocatable :: v1x(:,:), v2x(:,:), v1c(:,:), v2c(:,:) 392 real(dp), allocatable :: v2c_ud(:) 393 real(dp) :: zetas 394 ! 395 logical :: debug_xc = .false. 396 logical :: igcc_is_lyp 397 ! 398 allocate( v( nnr, nspin ) ) 399 if( dft_is_gradient() ) then 400 allocate( h( nnr, nspin, nspin ) ) 401 else 402 allocate( h( 1, 1, 1 ) ) 403 endif 404 ! 405 igcc_is_lyp = (get_igcc() == 3) 406 ! 407 etxc = 0.0d0 408 ! 409 allocate ( ex(nnr), ec(nnr), vx(nnr,nspin), vc(nnr,nspin) ) 410 IF ( nspin == 1 ) THEN 411 ! 412 ! spin-unpolarized case 413 ! 414 CALL xc( nnr, nspin, nspin, rhor, ex, ec, vx, vc ) 415 ! 416 v(:,nspin) = e2 * (vx(:,1) + vc(:,1) ) 417 etxc = e2 * SUM( (ex + ec)*rhor(:,nspin) ) 418 ! 419 ELSE 420 ! 421 ! spin-polarized case 422 ! 423 neg(1) = 0 424 neg(2) = 0 425 neg(3) = 0 426 ! 427 allocate ( rhox(nnr,2) ) !^ 428 rhox(:,1) = rhor(:,1) + rhor(:,2) 429 rhox(:,2) = rhor(:,1) - rhor(:,2) 430 ! 431 CALL xc( nnr, 2, 2, rhox, ex, ec, vx, vc ) 432 ! 433 DO ir = 1, nnr 434 ! 435 DO is = 1, nspin 436 v(ir,is) = e2 * (vx(ir,is) + vc(ir,is)) 437 ENDDO 438 etxc = etxc + e2 * (ex(ir) + ec(ir)) * rhox(ir,1) 439 ! 440 zetas = rhox(ir,2) / rhox(ir,1) 441 IF (rhor(ir,1) < 0.d0) neg(1) = neg(1) + 1 442 IF (rhor(ir,2) < 0.d0) neg(2) = neg(2) + 1 443 IF (ABS(zetas) > 1.d0) neg(3) = neg(3) + 1 444 ! 445 ENDDO 446 ! 447 deallocate ( rhox ) !^ 448 ! 449 ENDIF 450 deallocate ( vc, vx, ec, ex ) 451 ! 452 if( debug_xc ) then 453 open(unit=17,form='unformatted') 454 write(17) nnr, nspin 455 write(17) rhor 456 write(17) grhor 457 close(17) 458 debug_xc = .false. 459 end if 460 ! 461 ! gradient corrections 462 ! 463 if ( dft_is_gradient() ) then 464 ! 465 call change_threshold_gga( epsr ) 466 ! 467 allocate ( sx(nnr), sc(nnr), v1x(nnr,nspin), v1c(nnr,nspin), & 468 v2x(nnr,nspin), v2c(nnr,nspin) ) 469 if (nspin == 1) then 470 ! 471 ! ... This is the spin-unpolarised case 472 ! 473 call xc_gcx( nnr, nspin, rhor, grhor, sx, sc, v1x, v2x, v1c, v2c ) 474 ! 475 do k = 1, nnr 476 ! first term of the gradient correction: D(rho*Exc)/D(rho) 477 v(k,1) = v(k,1) + e2 * (v1x(k,1) + v1c(k,1)) 478 ! HERE h contains D(rho*Exc)/D(|grad rho|) / |grad rho| 479 h(k, 1, 1) = e2 * (v2x(k,1) + v2c(k,1)) 480 etxc = etxc + e2 * (sx(k) + sc(k)) 481 enddo 482 ! 483 else 484 ! 485 ! ... Spin-polarised case 486 ! 487 allocate (v2c_ud(nnr)) 488 call xc_gcx( nnr, 2, rhor, grhor, sx, sc, v1x, v2x, v1c, v2c, v2c_ud ) 489 ! 490 ! first term of the gradient correction : D(rho*Exc)/D(rho) 491 ! 492 v = v + e2*( v1x + v1c ) 493 ! 494 ! HERE h contains D(rho*Exc)/D(|grad rho|) / |grad rho| 495 ! 496 h(:,1,1) = e2 * (v2x(:,1) + v2c(:,1)) ! Spin UP-UP 497 h(:,1,2) = e2 * v2c_ud(:) ! Spin UP-DW 498 h(:,2,1) = e2 * v2c_ud(:) ! Spin DW-UP 499 h(:,2,2) = e2 * (v2x(:,2) + v2c(:,2)) ! Spin DW-DW 500 ! 501 etxc = etxc + e2 * SUM( sx(:)+sc(:) ) 502 ! 503 deallocate (v2c_ud) 504 endif 505 ! 506 deallocate ( v2c, v2x, v1c, v1x, sc, sx ) 507 end if 508 509 if( dft_is_gradient() ) then 510 ! 511 if( nspin == 1 ) then 512 ! 513 ! h contains D(rho*Exc)/D(|grad rho|) * (grad rho) / |grad rho| 514 ! 515!$omp parallel default(none), shared(nnr,grhor,h), private(ipol,k) 516 do ipol = 1, 3 517!$omp do 518 do k = 1, nnr 519 grhor (ipol, k, 1) = h (k, 1, 1) * grhor (ipol, k, 1) 520 enddo 521!$omp end do 522 end do 523!$omp end parallel 524 ! 525 ! 526 else 527 ! 528!$omp parallel default(none), shared(nnr,grhor,h), private(ipol,k,grup,grdw) 529 do ipol = 1, 3 530!$omp do 531 do k = 1, nnr 532 grup = grhor (ipol, k, 1) 533 grdw = grhor (ipol, k, 2) 534 grhor (ipol, k, 1) = h (k, 1, 1) * grup + h (k, 1, 2) * grdw 535 grhor (ipol, k, 2) = h (k, 2, 2) * grdw + h (k, 2, 1) * grup 536 enddo 537!$omp end do 538 enddo 539!$omp end parallel 540 ! 541 end if 542 ! 543 end if 544 545 rhor = v 546 547 deallocate( v ) 548 deallocate( h ) 549 550 return 551end subroutine exch_corr_cp 552