1! 2! Copyright (C) 2002-2008 Quantm-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! AB INITIO COSTANT PRESSURE MOLECULAR DYNAMICS 9! ---------------------------------------------- 10! Car-Parrinello Parallel Program 11 12 13 14 15 SUBROUTINE potential_print_info( iunit ) 16 17 USE control_flags, ONLY: iesr 18 19 INTEGER, INTENT(IN) :: iunit 20 21 WRITE(iunit,50) 22 WRITE(iunit,115) (2*iesr+1),(2*iesr+1),(2*iesr+1) 23 24 50 FORMAT(//,3X,'Potentials Parameters',/,3X,'---------------------') 25 115 FORMAT( 3X,'Ewald sum over ',I1,'*',I1,'*',I1,' cells') 26 27 RETURN 28 END SUBROUTINE potential_print_info 29 30!=----------------------------------------------------------------------------=! 31 32 SUBROUTINE cluster_bc( screen_coul, hg, omega, hmat ) 33 34 USE kinds, ONLY: DP 35 USE mp_global, ONLY: me_bgrp 36 USE fft_base, ONLY: dfftp 37 USE fft_interfaces, ONLY: fwfft 38 USE constants, ONLY: gsmall, pi 39 USE cell_base, ONLY: tpiba2, s_to_r, alat 40 USE fft_rho 41 42 IMPLICIT NONE 43 44 REAL(DP), INTENT(IN) :: hg( dfftp%ngm ) 45 REAL(DP), INTENT(IN) :: omega, hmat( 3, 3 ) 46 COMPLEX(DP) :: screen_coul( dfftp%ngm ) 47 48 REAL(DP), EXTERNAL :: qe_erf 49 50 ! ... Locals 51 ! 52 REAL(DP), ALLOCATABLE :: grr(:,:) 53 COMPLEX(DP), ALLOCATABLE :: grg(:,:) 54 REAL(DP) :: rc, r(3), s(3), rmod, g2, rc2, arg, fact 55 INTEGER :: ig, i, j, k, ir 56 INTEGER :: ir1, ir2, ir3, nr3l 57 58 ir1 = 1 59 ir2 = 1 60 ir3 = 1 61 DO k = 1, me_bgrp 62 ir3 = ir3 + dfftp%nr3p( k ) 63 END DO 64 nr3l = dfftp%my_nr3p 65 66 ALLOCATE( grr( dfftp%nnr, 1 ) ) 67 ALLOCATE( grg( dfftp%nnr, 1 ) ) 68 69 grr = 0.0d0 70 71 ! ... Martyna and Tuckerman convergence criterium 72 ! 73 rc = 7.0d0 / alat 74 rc2 = rc**2 75 fact = omega / ( dfftp%nr1 * dfftp%nr2 * dfftp%nr3 ) 76 IF( MOD(dfftp%nr1 * dfftp%nr2 * dfftp%nr3, 2) /= 0 ) fact = -fact 77 78 DO k = 1, nr3l 79 s(3) = DBLE ( (k-1) + (ir3 - 1) ) / dfftp%nr3 - 0.5d0 80 DO j = 1, dfftp%nr2 81 s(2) = DBLE ( (j-1) + (ir2 - 1) ) / dfftp%nr2 - 0.5d0 82 DO i = 1, dfftp%nr1 83 s(1) = DBLE ( (i-1) + (ir1 - 1) ) / dfftp%nr1 - 0.5d0 84 CALL S_TO_R( S, R, hmat ) 85 rmod = SQRT( r(1)**2 + r(2)**2 + r(3)**2 ) 86 ir = i + (j-1)*dfftp%nr1x + (k-1)*dfftp%nr1x*dfftp%nr2x 87 IF( rmod < gsmall ) THEN 88 grr( ir, 1 ) = fact * 2.0d0 * rc / SQRT( pi ) 89 ELSE 90 grr( ir, 1 ) = fact * qe_erf( rc * rmod ) / rmod 91 END IF 92 END DO 93 END DO 94 END DO 95 96 ! grg = FFT( grr ) 97 98 CALL rho_r2g( dfftp, grr, grg ) 99 100 DO ig = 1, SIZE( screen_coul ) 101 IF( hg(ig) < gsmall ) THEN 102 screen_coul(ig) = grg(1,1) - ( - pi / rc2 ) 103 ELSE 104 g2 = tpiba2 * hg(ig) 105 arg = - g2 / ( 4.0d0 * rc2 ) 106 screen_coul(ig) = grg(ig,1) - ( 4.0d0 * pi * EXP( arg ) / g2 ) 107 END IF 108 END DO 109 110 DEALLOCATE( grr, grg ) 111 112 RETURN 113 END SUBROUTINE cluster_bc 114 115 116!=----------------------------------------------------------------------------=! 117 118 119 SUBROUTINE vofps_x( eps, vloc, rhoeg, vps, sfac, omega ) 120 121 ! this routine computes: 122 ! omega = ht%deth 123 ! vloc_ps(ig) = (sum over is) sfac(is,ig) * vps(ig,is) 124 ! 125 ! Eps = Fact * omega * (sum over ig) cmplx ( rho_e(ig) ) * vloc_ps(ig) 126 ! if Gamma symmetry Fact = 2 else Fact = 1 127 ! 128 129 USE kinds, ONLY: DP 130 USE io_global, ONLY: stdout 131 USE ions_base, ONLY: nsp 132 USE gvect, ONLY: gstart 133 USE mp_global, ONLY: intra_bgrp_comm 134 USE mp, ONLY: mp_sum 135 USE fft_base, ONLY: dfftp 136 137 IMPLICIT NONE 138 139 ! ... Arguments 140 141 REAL(DP), INTENT(IN) :: vps(:,:) 142 REAL(DP), INTENT(IN) :: omega 143 COMPLEX(DP), INTENT(OUT) :: vloc(:) 144 COMPLEX(DP), INTENT(IN) :: rhoeg(:) 145 COMPLEX(DP), INTENT(IN) :: sfac(:,:) 146 COMPLEX(DP), INTENT(OUT) :: eps 147 148 ! ... Locals 149 150 INTEGER :: is, ig 151 COMPLEX(DP) :: vp 152 153 ! ... Subroutine body ... 154 ! 155 eps = (0.D0,0.D0) 156 ! 157 DO ig = gstart, dfftp%ngm 158 159 vp = (0.D0,0.D0) 160 DO is = 1, nsp 161 vp = vp + sfac( ig, is ) * vps( ig, is ) 162 END DO 163 164 vloc(ig) = vp 165 eps = eps + vp * CONJG( rhoeg( ig ) ) 166 167 END DO 168 ! ... 169 ! ... G = 0 element 170 ! 171 IF ( gstart == 2 ) THEN 172 vp = (0.D0,0.D0) 173 DO is = 1, nsp 174 vp = vp + sfac( 1, is) * vps(1, is) 175 END DO 176 vloc(1) = VP 177 eps = eps + vp * CONJG( rhoeg(1) ) * 0.5d0 178 END IF 179 ! 180 eps = 2.D0 * eps * omega 181 ! 182 CALL mp_sum( eps, intra_bgrp_comm ) 183 184 RETURN 185 END SUBROUTINE vofps_x 186 187 188!=----------------------------------------------------------------------------=! 189 190 SUBROUTINE vofloc_x( tscreen, ehte, ehti, eh, vloc, rhoeg, & 191 rhops, vps, sfac, omega, screen_coul ) 192 193 ! this routine computes: 194 ! omega = ht%deth 195 ! rho_e(ig) = (sum over iss) rhoeg(ig,iss) 196 ! rho_I(ig) = (sum over is) sfac(is,ig) * rhops(ig,is) 197 ! vloc_h(ig) = fpi / ( g(ig) * tpiba2 ) * { rho_e(ig) + rho_I(ig) } 198 ! 199 ! Eh = Fact * omega * (sum over ig) * fpi / ( g(ig) * tpiba2 ) * 200 ! { rho_e(ig) + rho_I(ig) } * conjugate { rho_e(ig) + rho_I(ig) } 201 ! if Gamma symmetry Fact = 1 else Fact = 1/2 202 ! 203 ! Hatree potential and local pseudopotential 204 ! vloc(ig) = vloc_h(ig) + vloc_ps(ig) 205 ! 206 207 USE kinds, ONLY: DP 208 USE constants, ONLY: fpi 209 USE cell_base, ONLY: tpiba2, tpiba 210 USE io_global, ONLY: stdout 211 USE gvect, ONLY: gstart, gg 212 USE ions_base, ONLY: nsp 213 USE mp_global, ONLY: intra_bgrp_comm 214 USE mp, ONLY: mp_sum 215 USE fft_base, ONLY: dfftp 216 217 IMPLICIT NONE 218 219 ! ... Arguments 220 221 LOGICAL, INTENT(IN) :: tscreen 222 REAL(DP), INTENT(IN) :: rhops(:,:), vps(:,:) 223 COMPLEX(DP), INTENT(INOUT) :: vloc(:) 224 COMPLEX(DP), INTENT(IN) :: rhoeg(:) 225 COMPLEX(DP), INTENT(IN) :: sfac(:,:) 226 REAL(DP), INTENT(OUT) :: ehte, ehti 227 REAL(DP), INTENT(IN) :: omega 228 COMPLEX(DP), INTENT(OUT) :: eh 229 COMPLEX(DP), INTENT(IN) :: screen_coul(:) 230 231 ! ... Locals 232 233 INTEGER :: is, ig 234 REAL(DP) :: fpibg, cost 235 COMPLEX(DP) :: rhet, rhog, rp, vscreen 236 237 ! ... Subroutine body ... 238 239 eh = 0.0d0 240 ehte = 0.0d0 241 ehti = 0.0d0 242 243!$omp parallel do default(shared), private(rp,is,rhet,rhog,fpibg), reduction(+:eh,ehte,ehti) 244 DO ig = gstart, dfftp%ngm 245 246 rp = (0.D0,0.D0) 247 DO is = 1, nsp 248 rp = rp + sfac( ig, is ) * rhops( ig, is ) 249 END DO 250 251 rhet = rhoeg( ig ) 252 rhog = rhet + rp 253 254 IF( tscreen ) THEN 255 fpibg = fpi / ( gg(ig) * tpiba2 ) + screen_coul(ig) 256 ELSE 257 fpibg = fpi / ( gg(ig) * tpiba2 ) 258 END IF 259 260 vloc(ig) = vloc(ig) + fpibg * rhog 261 eh = eh + fpibg * rhog * CONJG(rhog) 262 ehte = ehte + fpibg * DBLE(rhet * CONJG(rhet)) 263 ehti = ehti + fpibg * DBLE( rp * CONJG(rp)) 264 265 END DO 266 ! ... 267 ! ... G = 0 element 268 ! 269 IF ( gstart == 2 ) THEN 270 rp = (0.D0,0.D0) 271 IF( tscreen ) THEN 272 vscreen = screen_coul(1) 273 ELSE 274 vscreen = 0.0d0 275 END IF 276 DO IS = 1, nsp 277 rp = rp + sfac( 1, is) * rhops(1, is) 278 END DO 279 rhet = rhoeg(1) 280 rhog = rhet + rp 281 vloc(1) = vloc(1) + vscreen * rhog 282 eh = eh + vscreen * rhog * CONJG(rhog) 283 ehte = ehte + vscreen * DBLE(rhet * CONJG(rhet)) 284 ehti = ehti + vscreen * DBLE( rp * CONJG(rp)) 285 END IF 286 ! ... 287 eh = eh * omega 288 ehte = ehte * omega 289 ehti = ehti * omega 290 ! ... 291 CALL mp_sum(eh , intra_bgrp_comm) 292 CALL mp_sum(ehte, intra_bgrp_comm) 293 CALL mp_sum(ehti, intra_bgrp_comm) 294 ! 295 RETURN 296 END SUBROUTINE vofloc_x 297 298 299 SUBROUTINE force_loc_x( tscreen, rhoeg, fion, rhops, vps, ei1, ei2, ei3, & 300 sfac, omega, screen_coul ) 301 302 ! this routine computes: 303 ! 304 ! Local contribution to the forces on the ions 305 ! eigrx(ig,isa) = ei1( mill(1,ig), isa) 306 ! eigry(ig,isa) = ei2( mill(2,ig), isa) 307 ! eigrz(ig,isa) = ei3( mill(3,ig), isa) 308 ! fpibg = fpi / ( g(ig) * tpiba2 ) 309 ! tx_h(ig,is) = fpibg * rhops(ig, is) * CONJG( rho_e(ig) + rho_I(ig) ) 310 ! tx_ps(ig,is) = vps(ig,is) * CONJG( rho_e(ig) ) 311 ! gx(ig) = cmplx(0.D0, gx(1,ig),kind=DP) * tpiba 312 ! fion(x,isa) = fion(x,isa) + 313 ! Fact * omega * ( sum over ig, iss) (tx_h(ig,is) + tx_ps(ig,is)) * 314 ! gx(ig) * eigrx(ig,isa) * eigry(ig,isa) * eigrz(ig,isa) 315 ! if Gamma symmetry Fact = 2.0 else Fact = 1 316 ! 317 318 USE kinds, ONLY: DP 319 USE constants, ONLY: fpi 320 USE cell_base, ONLY: tpiba2, tpiba 321 USE io_global, ONLY: stdout 322 USE gvect, ONLY: mill, gstart, g, gg 323 USE ions_base, ONLY: nat, nsp, ityp 324 USE fft_base, ONLY: dfftp, dffts 325 326 IMPLICIT NONE 327 328 ! ... Arguments 329 330 LOGICAL :: tscreen 331 REAL(DP) :: fion(:,:) 332 REAL(DP) :: rhops(:,:), vps(:,:) 333 COMPLEX(DP) :: rhoeg(:) 334 COMPLEX(DP), INTENT(IN) :: sfac(:,:) 335 COMPLEX(DP) :: ei1(-dfftp%nr1:dfftp%nr1,nat) 336 COMPLEX(DP) :: ei2(-dfftp%nr2:dfftp%nr2,nat) 337 COMPLEX(DP) :: ei3(-dfftp%nr3:dfftp%nr3,nat) 338 REAL(DP) :: omega 339 COMPLEX(DP) :: screen_coul(:) 340 341 ! ... Locals 342 343 INTEGER :: is, ia, ig, ig1, ig2, ig3 344 REAL(DP) :: fpibg 345 COMPLEX(DP) :: cxc, rhet, rhog, vp, rp, gxc, gyc, gzc 346 COMPLEX(DP) :: teigr, cnvg, cvn, tx, ty, tz 347 COMPLEX(DP), ALLOCATABLE :: ftmp(:,:) 348 349 ! ... Subroutine body ... 350 351 ALLOCATE( ftmp( 3, SIZE( fion, 2 ) ) ) 352 353 ftmp = 0.0d0 354 355 DO ig = gstart, dffts%ngm 356 357 RP = (0.D0,0.D0) 358 DO IS = 1, nsp 359 RP = RP + sfac( ig, is ) * rhops( ig, is ) 360 END DO 361 362 RHET = RHOEG( ig ) 363 RHOG = RHET + RP 364 365 IF( tscreen ) THEN 366 FPIBG = fpi / ( gg(ig) * tpiba2 ) + screen_coul(ig) 367 ELSE 368 FPIBG = fpi / ( gg(ig) * tpiba2 ) 369 END IF 370 371 ig1 = mill(1,IG) 372 ig2 = mill(2,IG) 373 ig3 = mill(3,IG) 374 GXC = CMPLX(0.D0,g(1,IG),kind=DP) 375 GYC = CMPLX(0.D0,g(2,IG),kind=DP) 376 GZC = CMPLX(0.D0,g(3,IG),kind=DP) 377 DO ia = 1, nat 378 is = ityp(ia) 379 CNVG = RHOPS(IG,is) * FPIBG * CONJG(rhog) 380 CVN = VPS(ig, is) * CONJG(rhet) 381 TX = (CNVG+CVN) * GXC 382 TY = (CNVG+CVN) * GYC 383 TZ = (CNVG+CVN) * GZC 384 TEIGR = ei1(IG1,ia) * ei2(IG2,ia) * ei3(IG3,ia) 385 ftmp(1,ia) = ftmp(1,ia) + TEIGR*TX 386 ftmp(2,ia) = ftmp(2,ia) + TEIGR*TY 387 ftmp(3,ia) = ftmp(3,ia) + TEIGR*TZ 388 END DO 389 390 END DO 391 ! 392 fion = fion + DBLE(ftmp) * 2.D0 * omega * tpiba 393 394 DEALLOCATE( ftmp ) 395 396 RETURN 397 END SUBROUTINE force_loc_x 398 399 400! 401!=----------------------------------------------------------------------------=! 402 SUBROUTINE vofesr( iesr, esr, desr, fion, taus, tstress, hmat ) 403!=----------------------------------------------------------------------------=! 404 405 USE kinds, ONLY : DP 406 USE constants, ONLY : sqrtpm1 407 USE cell_base, ONLY : s_to_r, pbcs 408 USE mp_global, ONLY : nproc_bgrp, me_bgrp, intra_bgrp_comm 409 USE mp, ONLY : mp_sum 410 USE ions_base, ONLY : rcmax, zv, nsp, na, nat, ityp 411 412 IMPLICIT NONE 413 414! ... ARGUMENTS 415 416 INTEGER, INTENT(IN) :: iesr 417 REAL(DP), INTENT(IN) :: taus(3,nat) 418 REAL(DP) :: ESR 419 REAL(DP) :: DESR(6) 420 REAL(DP) :: FION(3,nat) 421 LOGICAL, INTENT(IN) :: TSTRESS 422 REAL(DP), INTENT(in) :: hmat( 3, 3 ) 423 424 REAL(DP), EXTERNAL :: qe_erfc 425 426 INTEGER, EXTERNAL :: ldim_block, gind_block 427 428 429! ... LOCALS 430 431 INTEGER :: na_loc, ia_s, ia_e, igis 432 INTEGER :: k, i, j, is, ia, ib, ix, iy, iz 433 LOGICAL :: split, tzero, tshift 434 REAL(DP), ALLOCATABLE :: zv2(:,:) 435 REAL(DP), ALLOCATABLE :: rc(:,:) 436 REAL(DP), ALLOCATABLE :: fionloc(:,:) 437 REAL(DP) :: rxlm(3), sxlm(3) 438 REAL(DP) :: xlm,ylm,zlm, xlm0,ylm0,zlm0, erre2, rlm, arg, esrtzero 439 REAL(DP) :: addesr, addpre, repand, fxx 440 REAL(DP) :: rckj_m1 441 REAL(DP) :: zvk, zvj, zv2_kj 442 REAL(DP) :: fact_pre 443 444 INTEGER, DIMENSION(6), PARAMETER :: ALPHA = (/ 1,2,3,2,3,3 /) 445 INTEGER, DIMENSION(6), PARAMETER :: BETA = (/ 1,1,1,2,2,3 /) 446 447! ... SUBROUTINE BODY 448 449 ALLOCATE( rc( nsp, nsp ) ) 450 ALLOCATE( zv2( nsp, nsp ) ) 451 ALLOCATE( fionloc( 3, nat ) ) 452 rc = 0.0_DP 453 zv2 = 0.0_DP 454 fionloc = 0.0_DP 455 456 ! Here pre-compute some factors 457 458 DO j = 1, nsp 459 DO k = 1, nsp 460 zv2( k, j ) = zv( k ) * zv( j ) 461 rc ( k, j ) = SQRT( rcmax(k)**2 + rcmax(j)**2 ) 462 END DO 463 END DO 464 465 xlm = 1.0_DP 466 ylm = 1.0_DP 467 zlm = 1.0_DP 468 ESR = 0.0_DP 469 DESR = 0.0_DP 470 471 ! Distribute the atoms pairs to processors 472 473 NA_LOC = ldim_block( nat, nproc_bgrp, me_bgrp) 474 IA_S = gind_block( 1, nat, nproc_bgrp, me_bgrp ) 475 IA_E = IA_S + NA_LOC - 1 476 477 DO ia = ia_s, ia_e 478 k = ityp(ia) 479 DO ib = ia, nat 480 j = ityp(ib) 481 482 zv2_kj = zv2(k,j) 483 rckj_m1 = 1.0_DP / rc(k,j) 484 fact_pre = (2.0_DP * zv2_kj * sqrtpm1) * rckj_m1 485 486 IF( ia.EQ.ib ) THEN 487 ! ... same atoms 488 xlm=0.0_DP; ylm=0.0_DP; zlm=0.0_DP; 489 tzero=.TRUE. 490 ELSE 491 ! ... different atoms 492 xlm0= taus(1,ia) - taus(1,ib) 493 ylm0= taus(2,ia) - taus(2,ib) 494 zlm0= taus(3,ia) - taus(3,ib) 495 CALL pbcs(xlm0,ylm0,zlm0,xlm,ylm,zlm,1) 496 TZERO=.FALSE. 497 END IF 498 499 DO IX=-IESR,IESR 500 sxlm(1) = XLM + DBLE(IX) 501 DO IY=-IESR,IESR 502 sxlm(2) = YLM + DBLE(IY) 503 DO IZ=-IESR,IESR 504 TSHIFT= IX.EQ.0 .AND. IY.EQ.0 .AND. IZ.EQ.0 505 IF( .NOT. ( TZERO .AND. TSHIFT ) ) THEN 506 sxlm(3) = ZLM + DBLE(IZ) 507 CALL S_TO_R( sxlm, rxlm, hmat ) 508 ERRE2 = rxlm(1)**2 + rxlm(2)**2 + rxlm(3)**2 509 RLM = SQRT(ERRE2) 510 ARG = RLM * rckj_m1 511 IF (TZERO) THEN 512 ESRTZERO=0.5D0 513 ELSE 514 ESRTZERO=1.D0 515 END IF 516 ADDESR = ZV2_KJ * qe_erfc(ARG) / RLM 517 ESR = ESR + ESRTZERO*ADDESR 518 ADDPRE = FACT_PRE * EXP(-ARG*ARG) 519 REPAND = ESRTZERO*(ADDESR + ADDPRE)/ERRE2 520 ! 521 DO i = 1, 3 522 fxx = repand * rxlm( i ) 523 fionloc( i, ia ) = fionloc( i, ia ) + fxx 524 fionloc( i, ib ) = fionloc( i, ib ) - fxx 525 END DO 526 ! 527 IF( tstress ) THEN 528 DO i = 1, 6 529 fxx = repand * rxlm( alpha( i ) ) * rxlm( beta( i ) ) 530 desr( i ) = desr( i ) - fxx 531 END DO 532 END IF 533 END IF 534 END DO ! IZ 535 END DO ! IY 536 END DO ! IX 537 END DO 538 END DO 539 540! 541! each processor add its own contribution to the array FION 542! 543 DO ia = 1, nat 544 FION(1,ia) = FION(1,ia)+FIONLOC(1,ia) 545 FION(2,ia) = FION(2,ia)+FIONLOC(2,ia) 546 FION(3,ia) = FION(3,ia)+FIONLOC(3,ia) 547 END DO 548 549 CALL mp_sum(esr, intra_bgrp_comm) 550 551 DEALLOCATE(rc) 552 DEALLOCATE(zv2) 553 DEALLOCATE(fionloc) 554 555 RETURN 556!=----------------------------------------------------------------------------=! 557 END SUBROUTINE vofesr 558!=----------------------------------------------------------------------------=! 559 560 561 562!=----------------------------------------------------------------------------=! 563 SUBROUTINE self_vofhar_x( tscreen, self_ehte, vloc, rhoeg, omega, hmat ) 564!=----------------------------------------------------------------------------=! 565 566 ! adds the hartree part of the self interaction 567 568 USE kinds, ONLY: DP 569 USE constants, ONLY: fpi 570 USE control_flags, ONLY: gamma_only 571 USE cell_base, ONLY: tpiba2 572 USE gvect, ONLY: gstart, gg 573 USE sic_module, ONLY: sic_epsilon, sic_alpha 574 USE mp_global, ONLY: intra_bgrp_comm 575 USE mp, ONLY: mp_sum 576 USE fft_base, ONLY: dfftp 577 578 IMPLICIT NONE 579 580 ! ... Arguments 581 LOGICAL :: tscreen 582 COMPLEX(DP) :: vloc(:) 583 COMPLEX(DP) :: rhoeg(:,:) 584 REAL(DP) :: self_ehte 585 REAL(DP), INTENT(IN) :: omega 586 REAL(DP), INTENT(IN) :: hmat( 3, 3 ) 587 588 ! ... Locals 589 590 INTEGER :: ig 591 REAL(DP) :: fpibg 592 COMPLEX(DP) :: rhog 593 COMPLEX(DP) :: ehte 594 COMPLEX(DP) :: vscreen 595 COMPLEX(DP), ALLOCATABLE :: screen_coul(:) 596 597 ! ... Subroutine body ... 598 599 IF( tscreen ) THEN 600 ALLOCATE( screen_coul( dfftp%ngm ) ) 601 CALL cluster_bc( screen_coul, gg, omega, hmat ) 602 END IF 603 604 !== HARTREE == 605 606 ehte = 0.D0 607 608 DO IG = gstart, dfftp%ngm 609 610 rhog = rhoeg(ig,1) - rhoeg(ig,2) 611 612 IF( tscreen ) THEN 613 FPIBG = fpi / ( gg(ig) * tpiba2 ) + screen_coul(ig) 614 ELSE 615 FPIBG = fpi / ( gg(ig) * tpiba2 ) 616 END IF 617 618 vloc(ig) = fpibg * rhog 619 ehte = ehte + fpibg * rhog * CONJG(rhog) 620 621 END DO 622 623 ! ... G = 0 element 624 ! 625 IF ( gstart == 2 ) THEN 626 rhog = rhoeg(1,1) - rhoeg(1,2) 627 IF( tscreen ) THEN 628 vscreen = screen_coul(1) 629 ELSE 630 vscreen = 0.0d0 631 END IF 632 vloc(1) = vscreen * rhog 633 ehte = ehte + vscreen * rhog * CONJG(rhog) 634 END IF 635 636 ! ... 637 638 IF( .NOT. gamma_only ) THEN 639 ehte = ehte * 0.5d0 640 END IF 641 ! 642 self_ehte = DBLE(ehte) * omega * sic_epsilon 643 vloc = vloc * sic_epsilon 644 645 CALL mp_sum( self_ehte, intra_bgrp_comm ) 646 647 IF( ALLOCATED( screen_coul ) ) DEALLOCATE( screen_coul ) 648 649 RETURN 650!=----------------------------------------------------------------------------=! 651 END SUBROUTINE self_vofhar_x 652!=----------------------------------------------------------------------------=! 653