1! 2! Copyright (C) 2002-2007 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!----------------------------------------------------------------------- 9SUBROUTINE vofrho_x( nfi, rhor, drhor, rhog, drhog, rhos, rhoc, tfirst, & 10 tlast, ei1, ei2, ei3, irb, eigrb, sfac, tau0, fion ) 11!----------------------------------------------------------------------- 12! computes: the one-particle potential v in real space, 13! the total energy etot, 14! the forces fion acting on the ions, 15! the derivative of total energy to cell parameters h 16! rhor input : electronic charge on dense real space grid 17! (plus core charge if present) 18! rhog input : electronic charge in g space (up to density cutoff) 19! rhos input : electronic charge on smooth real space grid 20! rhor output: total potential on dense real space grid 21! rhos output: total potential on smooth real space grid 22! 23 USE kinds, ONLY: dp 24 USE control_flags, ONLY: iprint, iverbosity, thdyn, tpre, tfor, & 25 tprnfor, iesr, textfor 26 USE io_global, ONLY: stdout 27 USE ions_base, ONLY: nsp, na, nat, rcmax, compute_eextfor 28 USE cell_base, ONLY: omega, r_to_s 29 USE cell_base, ONLY: alat, at, tpiba2, h, ainv 30 USE cell_base, ONLY: ibrav, isotropic !True if volume option is chosen for cell_dofree 31 USE cp_main_variables, ONLY: iprint_stdout !print control 32 USE gvect, ONLY: gstart, gg, g 33 USE electrons_base, ONLY: nspin 34 USE constants, ONLY: pi, fpi, au_gpa, e2 35 USE energies, ONLY: etot, eself, enl, ekin, epseu, esr, eht, & 36 exc, eextfor, exx 37 USE local_pseudo, ONLY: vps, dvps, rhops 38 USE uspp, ONLY: nlcc_any 39 USE smallbox_gvec 40 USE dener, ONLY: detot, dekin, dps, dh, dsr, dxc, denl, denlc,& 41 detot6, dekin6, dps6, dh6, dsr6, dxc6, denl6 42 USE mp, ONLY: mp_sum 43 USE mp_global, ONLY: intra_bgrp_comm 44 USE funct, ONLY: dft_is_meta, dft_is_nonlocc, nlc, get_inlc,& 45 dft_is_hybrid, exx_is_active 46 USE vdW_DF, ONLY: vdW_DF_stress 47 use rVV10, ONLY: rVV10_stress 48 USE pres_ai_mod, ONLY: abivol, abisur, v_vol, P_ext, volclu, & 49 Surf_t, surfclu 50 USE fft_interfaces, ONLY: fwfft, invfft 51 USE sic_module, ONLY: self_interaction, sic_alpha 52 USE energies, ONLY: self_exc, self_ehte 53 USE cp_interfaces, ONLY: pseudo_stress, compute_gagb, stress_hartree, & 54 add_drhoph, stress_local, force_loc, self_vofhar 55 USE fft_base, ONLY: dfftp, dffts 56 USE ldaU_cp, ONLY: e_hubbard 57 USE control_flags, ONLY: ts_vdw 58 USE tsvdw_module, ONLY: tsvdw_calculate 59 USE tsvdw_module, ONLY: EtsvdW,UtsvdW,FtsvdW,HtsvdW 60 USE mp_global, ONLY: me_image 61 USE exx_module, ONLY: dexx_dh, exxalfa ! exx_wf related 62 USE fft_rho 63 USE fft_helper_subroutines 64 65 USE plugin_variables, ONLY: plugin_etot 66 67 IMPLICIT NONE 68! 69 LOGICAL :: tlast, tfirst 70 INTEGER :: nfi 71 REAL(DP) :: rhor(:,:), drhor(:,:,:,:), rhos(:,:), fion(:,:) 72 REAL(DP) :: rhoc(:), tau0(:,:) 73 ! COMPLEX(DP) ei1(-nr1:nr1,nat), ei2(-nr2:nr2,nat), ei3(-nr3:nr3,nat) 74 COMPLEX(DP) :: ei1(:,:), ei2(:,:), ei3(:,:) 75 COMPLEX(DP) :: eigrb(:,:) 76 COMPLEX(DP) :: rhog(:,:), drhog(:,:,:,:) 77 COMPLEX(DP) :: sfac(:,:) 78 INTEGER :: irb(:,:) 79 ! 80 INTEGER iss, isup, isdw, ig, ir, i, j, k, ij, is, ia, inlc 81 REAL(DP) :: vtxc, vave, ebac, wz, eh, ehpre, enlc 82 COMPLEX(DP) fp, fm, ci, drhop, zpseu, zh 83 COMPLEX(DP), ALLOCATABLE :: rhotmp(:), vtemp(:) 84 COMPLEX(DP), ALLOCATABLE :: drhot(:,:) 85 REAL(DP), ALLOCATABLE :: gagb(:,:), rhosave(:,:), newrhosave(:,:), rhocsave(:) 86 ! 87 REAL(DP), ALLOCATABLE :: fion1( :, : ) 88 REAL(DP), ALLOCATABLE :: stmp( :, : ) 89 ! 90 COMPLEX(DP), ALLOCATABLE :: self_vloc(:) 91 COMPLEX(DP) :: self_rhoeg 92 REAL(DP) :: self_ehtet, fpibg 93 LOGICAL :: ttsic 94 REAL(DP) :: detmp( 3, 3 ), desr( 6 ), deps( 6 ) 95 REAL(DP) :: detmp2( 3, 3 ) 96 REAL(DP) :: ht( 3, 3 ) 97 REAL(DP) :: deht( 6 ) 98 COMPLEX(DP) :: screen_coul( 1 ) 99 REAL(DP) :: dexx(3,3) ! stress tensor from exact exchange exx_wf related 100! 101 INTEGER, DIMENSION(6), PARAMETER :: alpha = (/ 1,2,3,2,3,3 /) 102 INTEGER, DIMENSION(6), PARAMETER :: beta = (/ 1,1,1,2,2,3 /) 103 104 ! ... dalbe(:) = delta( alpha(:), beta(:) ) 105 REAL(DP), DIMENSION(6), PARAMETER :: dalbe = & 106 (/ 1.0_DP, 0.0_DP, 0.0_DP, 1.0_DP, 0.0_DP, 1.0_DP /) 107 108 109 CALL start_clock( 'vofrho' ) 110 111 ! 112 ! TS-vdW calculation (RAD) 113 ! 114 IF (ts_vdw) THEN 115 ! 116 CALL start_clock( 'ts_vdw' ) 117 ALLOCATE (stmp(3,nat), rhocsave(dfftp%nnr) ) 118 stmp(:,:) = tau0(:,:) 119 ! 120 IF ( nspin==2 ) THEN 121 rhocsave(:) = rhor(:,1) + rhor(:,2) 122 ELSE 123 rhocsave(:) = rhor(:,1) 124 ENDIF 125 ! 126 CALL tsvdw_calculate(stmp,rhocsave) 127 ! 128 DEALLOCATE (rhocsave,stmp) 129 CALL stop_clock( 'ts_vdw' ) 130 ! 131 END IF 132 ! 133 ci = ( 0.0d0, 1.0d0 ) 134 ! 135 ! wz = factor for g.neq.0 because of c*(g)=c(-g) 136 ! 137 wz = 2.0d0 138 ! 139 ht = TRANSPOSE( h ) 140 ! 141 ALLOCATE( vtemp( dfftp%ngm ) ) 142 ALLOCATE( rhotmp( dfftp%ngm ) ) 143 ! 144 IF ( tpre ) THEN 145 ALLOCATE( drhot( dfftp%ngm, 6 ) ) 146 ALLOCATE( gagb( 6, dfftp%ngm ) ) 147 CALL compute_gagb( gagb, g, dfftp%ngm, tpiba2 ) 148 END IF 149! 150! ab-initio pressure and surface tension contributions to the potential 151! 152 if (abivol.or.abisur) call vol_clu(rhor,rhog,nfi) 153 ! 154 ! compute plugin contributions to the potential, add it later 155 ! 156 CALL plugin_get_potential(rhor,nfi) 157 ! 158 ! compute plugin contributions to the energy 159 ! 160 plugin_etot = 0.0_dp 161 ! 162 CALL plugin_energy(rhor,plugin_etot) 163 ! 164 ttsic = ( ABS( self_interaction ) /= 0 ) 165 ! 166 IF( ttsic ) ALLOCATE( self_vloc( dfftp%ngm ) ) 167 ! 168 ! first routine in which fion is calculated: annihilation 169 ! 170 fion = 0.d0 171 ! 172 ! forces on ions, ionic term in real space 173 ! 174 IF( tprnfor .OR. tfor .OR. tfirst .OR. tpre ) THEN 175 ! 176 ALLOCATE( stmp( 3, nat ) ) 177 ! 178 CALL r_to_s( tau0, stmp, nat, ainv ) 179 ! 180 CALL vofesr( iesr, esr, dsr6, fion, stmp, tpre, h ) 181 ! 182 call mp_sum( fion, intra_bgrp_comm ) 183 ! 184 DEALLOCATE( stmp ) 185 ! 186 END IF 187! 188!$omp parallel default(shared), private(ig,is,ij,i,j,k) 189!$omp workshare 190 rhotmp( 1:dfftp%ngm ) = rhog( 1:dfftp%ngm, 1 ) 191!$omp end workshare 192 IF( nspin == 2 ) THEN 193!$omp workshare 194 rhotmp( 1:dfftp%ngm ) = rhotmp( 1:dfftp%ngm ) + rhog( 1:dfftp%ngm, 2 ) 195!$omp end workshare 196 END IF 197 ! 198 IF( tpre ) THEN 199!$omp do 200 DO ij = 1, 6 201 i = alpha( ij ) 202 j = beta( ij ) 203 drhot( :, ij ) = 0.0d0 204 DO k = 1, 3 205 drhot( :, ij ) = drhot( :, ij ) + drhog( :, 1, i, k ) * ht( k, j ) 206 END DO 207 END DO 208!$omp end do 209 IF( nspin == 2 ) THEN 210!$omp do 211 DO ij = 1, 6 212 i = alpha( ij ) 213 j = beta( ij ) 214 DO k = 1, 3 215 drhot( :, ij ) = drhot( :, ij ) + drhog( :, 2, i, k ) * ht( k, j ) 216 END DO 217 END DO 218!$omp end do 219 ENDIF 220 END IF 221 ! 222 ! calculation local potential energy 223 ! 224!$omp master 225 zpseu = 0.0d0 226!$omp end master 227 ! 228!$omp do 229 DO ig = 1, SIZE(vtemp) 230 vtemp(ig)=(0.d0,0.d0) 231 END DO 232 DO is=1,nsp 233!$omp do 234 DO ig=1,dffts%ngm 235 vtemp(ig)=vtemp(ig)+CONJG(rhotmp(ig))*sfac(ig,is)*vps(ig,is) 236 END DO 237 END DO 238!$omp do reduction(+:zpseu) 239 DO ig=1,dffts%ngm 240 zpseu = zpseu + vtemp(ig) 241 END DO 242!$omp end parallel 243 244 epseu = wz * DBLE(zpseu) 245 ! 246 IF (gstart == 2) epseu = epseu - DBLE( vtemp(1) ) 247 ! 248 CALL mp_sum( epseu, intra_bgrp_comm ) 249 epseu = epseu * omega 250 ! 251 IF( tpre ) THEN 252 ! 253 CALL stress_local( dps6, epseu, gagb, sfac, rhotmp, drhot, omega ) 254 ! 255 END IF 256 ! 257 ! 258 ! calculation hartree energy 259 ! 260 ! 261 self_ehtet = 0.d0 262 ! 263 IF( ttsic ) self_vloc = 0.d0 264 265 zh = 0.0d0 266 267!$omp parallel default(shared), private(ig,is) 268 269 DO is=1,nsp 270!$omp do 271 DO ig=1,dffts%ngm 272 rhotmp(ig)=rhotmp(ig)+sfac(ig,is)*rhops(ig,is) 273 END DO 274 END DO 275 ! 276!$omp do 277 DO ig = gstart, dfftp%ngm 278 vtemp(ig) = CONJG( rhotmp( ig ) ) * rhotmp( ig ) / gg( ig ) 279 END DO 280 281!$omp do reduction(+:zh) 282 DO ig = gstart, dfftp%ngm 283 zh = zh + vtemp(ig) 284 END DO 285 286!$omp end parallel 287 288 eh = DBLE( zh ) * wz * 0.5d0 * fpi / tpiba2 289! 290 CALL mp_sum( eh, intra_bgrp_comm ) 291 ! 292 IF ( ttsic ) THEN 293 ! 294 CALL self_vofhar( .false., self_ehte, self_vloc, rhog, omega, h ) 295 ! 296 eh = eh - self_ehte / omega 297 ! 298 END IF 299 ! 300 IF(tpre) THEN 301 ! 302 CALL add_drhoph( drhot, sfac, gagb ) 303 CALL stress_hartree(dh6, eh*omega, sfac, rhotmp, drhot, gagb, omega ) 304 ! 305 DEALLOCATE( gagb ) 306 DEALLOCATE( drhot ) 307 ! 308 END IF 309 ! 310 ! forces on ions, ionic term in reciprocal space 311 ! 312 ALLOCATE( fion1( 3, nat ) ) 313 ! 314 fion1 = 0.d0 315 ! 316 IF( tprnfor .OR. tfor .OR. tpre) THEN 317 vtemp( 1:dfftp%ngm ) = rhog( 1:dfftp%ngm, 1 ) 318 IF( nspin == 2 ) THEN 319 vtemp( 1:dfftp%ngm ) = vtemp(1:dfftp%ngm) + rhog( 1:dfftp%ngm, 2 ) 320 END IF 321 CALL force_loc( .false., vtemp, fion1, rhops, vps, ei1, ei2, ei3, sfac, omega, screen_coul ) 322 END IF 323 ! 324 ! calculation hartree + local pseudo potential 325 ! 326 ! 327 IF (gstart == 2) vtemp(1)=(0.d0,0.d0) 328 329!$omp parallel default(shared), private(ig,is) 330!$omp do 331 DO ig=gstart,dfftp%ngm 332 vtemp(ig)=rhotmp(ig)*fpi/(tpiba2*gg(ig)) 333 END DO 334 ! 335 DO is=1,nsp 336!$omp do 337 DO ig=1,dffts%ngm 338 vtemp(ig)=vtemp(ig)+sfac(ig,is)*vps(ig,is) 339 END DO 340 END DO 341!$omp end parallel 342 343 DEALLOCATE (rhotmp) 344! 345! vtemp = v_loc(g) + v_h(g) 346! 347! =================================================================== 348! calculation exchange and correlation energy and potential 349! ------------------------------------------------------------------- 350 ! 351 ! ... UGLY HACK WARNING: rhor must be saved before exch_corr_h 352 ! ... overwrites it and before add_cc adds to it the core charge 353 ! ... We also need an allocated rhoc array even in absence of core charge 354 ! 355 IF ( dft_is_nonlocc() ) THEN 356 ALLOCATE ( rhosave(dfftp%nnr,nspin), rhocsave(dfftp%nnr) ) 357 ALLOCATE ( newrhosave(dfftp%nnr,nspin) ) 358 rhosave(:,:) = rhor(:,:) 359 IF ( SIZE(rhoc) == dfftp%nnr ) THEN 360 rhocsave(:)= rhoc(:) 361 ELSE 362 rhocsave(:)= 0.0_dp 363 ENDIF 364 END IF 365 ! 366 IF ( nlcc_any ) CALL add_cc( rhoc, rhog, rhor ) 367 CALL exch_corr_h( nspin, rhog, rhor, rhoc, sfac, exc, dxc, self_exc ) 368 ! 369 ! ... add non local corrections (if any) 370 ! 371 IF ( dft_is_nonlocc() ) THEN 372 ! 373 ! ... UGLY HACK WARNING: nlc adds nonlocal term IN RYDBERG A.U. 374 ! ... to input energy and potential (potential is stored in rhor) 375 ! 376 enlc = 0.0_dp 377 vtxc = 0.0_dp 378 rhor = rhor*e2 379 CALL nlc( rhosave, rhocsave, nspin, enlc, vtxc, rhor ) 380 rhor = rhor/e2 381 CALL mp_sum( enlc, intra_bgrp_comm ) 382 exc = exc + enlc / e2 383 ! 384 ! ... non-local XC contribution to stress brought to crystal axis, 385 ! ... transformed into energy derivative (Ha), added to dxc 386 ! 387 IF ( tpre ) THEN 388 CALL mp_sum( vtxc, intra_bgrp_comm ) 389 DO i=1,3 390 DO j=1,3 391 dxc( i, j ) = dxc( i, j ) + (enlc - vtxc)/e2*ainv( j, i ) 392 END DO 393 END DO 394 denlc(:,:) = 0.0_dp 395 ! 396 !^^ ... TEMPORARY FIX (newlsda) ... 397 IF ( nspin==2 ) THEN ! PH adjusted 05/2020 398 rhosave(:,1) = rhosave(:,1) + rhosave(:,2) 399 newrhosave(:,1) = rhosave(:,1) + rhosave(:,2) 400 newrhosave(:,2) = rhosave(:,1) - rhosave(:,2) 401 ! CALL errore('stres_vdW', 'LSDA+stress+vdW-DF not implemented',1) 402 END IF 403 !^^....................... 404 ! 405 inlc = get_inlc() 406 IF ( inlc > 0 .AND. inlc < 26 ) THEN 407 CALL vdW_DF_stress ( newrhosave, rhocsave, nspin, denlc ) 408 ELSEIF ( inlc == 26 ) then 409 CALL rVV10_stress ( rhosave(:,1), rhocsave, nspin, denlc ) 410 END IF 411 ! 412 dxc(:,:) = dxc(:,:) - omega/e2 * MATMUL(denlc,TRANSPOSE(ainv)) 413 END IF 414 DEALLOCATE ( rhocsave, rhosave, newrhosave ) 415 ELSE 416 denlc(:,:) = 0.0_dp 417 END IF 418 ! 419 ! Add TS-vdW wavefunction forces to rhor here... (RAD) 420 ! 421 IF (ts_vdw) THEN 422 ! 423 CALL fftx_add_field( rhor, UtsvdW, dfftp ) 424 ! 425 END IF 426 ! 427 ! add plugin contributions to potential here... 428 ! 429 CALL plugin_add_potential( rhor ) 430 ! 431! 432! rhor contains the xc potential in r-space 433! 434! =================================================================== 435! fourier transform of xc potential to g-space (dense grid) 436! ------------------------------------------------------------------- 437! 438 IF( abivol .or. abisur ) THEN 439 CALL rho_r2g ( dfftp, rhor, rhog, v_vol ) 440 ELSE 441 CALL rho_r2g ( dfftp, rhor, rhog ) 442 END IF 443 444 IF( nspin == 1 ) THEN 445 rhog( 1:dfftp%ngm, 1 ) = rhog( 1:dfftp%ngm, 1 ) + vtemp(1:dfftp%ngm) 446 ELSE 447 isup=1 448 isdw=2 449 rhog( 1:dfftp%ngm, isup ) = rhog( 1:dfftp%ngm, isup ) + vtemp(1:dfftp%ngm) 450 rhog( 1:dfftp%ngm, isdw ) = rhog( 1:dfftp%ngm, isdw ) + vtemp(1:dfftp%ngm) 451 IF( ttsic ) THEN 452 rhog( 1:dfftp%ngm, isup ) = rhog( 1:dfftp%ngm, isup ) - self_vloc(1:dfftp%ngm) 453 rhog( 1:dfftp%ngm, isdw ) = rhog( 1:dfftp%ngm, isdw ) - self_vloc(1:dfftp%ngm) 454 END IF 455 END IF 456 457 DEALLOCATE (vtemp) 458 IF( ttsic ) DEALLOCATE( self_vloc ) 459! 460! rhog contains now the total (local+Hartree+xc) potential in g-space 461! 462 IF( tprnfor .OR. tfor ) THEN 463 464 IF ( nlcc_any ) CALL force_cc( irb, eigrb, rhor, fion1 ) 465 466 CALL mp_sum( fion1, intra_bgrp_comm ) 467 ! 468 ! add g-space ionic and core correction contributions to fion 469 ! 470 fion = fion + fion1 471 ! 472 ! Add TS-vdW ion forces to fion here... (RAD) 473 ! 474 IF (ts_vdw) THEN 475 fion1(:,:) = FtsvdW(:,:) 476 fion = fion + fion1 477 !fion=fion+FtsvdW 478 END IF 479 ! 480 ! plugin patches on internal forces 481 ! 482 CALL plugin_int_forces(fion) 483 484 END IF 485 486 DEALLOCATE( fion1 ) 487! 488! 489! =================================================================== 490! fourier transform of total potential to r-space (dense grid) 491! ------------------------------------------------------------------- 492 493 CALL rho_g2r( dfftp, rhog, rhor ) 494 495 IF(nspin.EQ.1) THEN 496 vave=SUM(rhor(:,1))/DBLE( dfftp%nr1* dfftp%nr2* dfftp%nr3) 497 ELSE 498 isup=1 499 isdw=2 500 vave=(SUM(rhor(:,isup))+SUM(rhor(:,isdw))) / 2.0d0 / DBLE( dfftp%nr1 * dfftp%nr2 * dfftp%nr3 ) 501 END IF 502 503 CALL mp_sum( vave, intra_bgrp_comm ) 504 505 ! 506 ! fourier transform of total potential to r-space (smooth grid) 507 ! 508 CALL rho_g2r ( dffts, rhog, rhos ) 509 510 IF( dft_is_meta() ) CALL vofrho_meta( ) 511 512 ebac = 0.0d0 513 ! 514 eht = eh * omega + esr - eself 515 ! 516 eextfor = 0.0_DP 517 IF( textfor ) eextfor = compute_eextfor( tau0 ) 518 ! 519 ! etot is the total energy ; ekin, enl were calculated in rhoofr 520 ! 521 etot = ekin + eht + epseu + enl + exc + ebac +e_hubbard + eextfor 522 ! 523 ! Add EXX energy to etot here. exx_wf related 524 ! 525 IF(dft_is_hybrid().AND.exx_is_active()) THEN 526 ! 527 etot = etot - exxalfa*exx 528 ! 529 END IF 530 ! 531 ! Add TS-vdW energy to etot here... (RAD) 532 ! 533 IF (ts_vdw) etot=etot+EtsvdW 534 ! 535 if (abivol) etot = etot + P_ext*volclu 536 if (abisur) etot = etot + Surf_t*surfclu 537 ! 538 ! Add possible external contribution from plugins to the energy 539 ! 540 etot = etot + plugin_etot 541 ! 542 IF( tpre ) THEN 543 ! 544 detot6 = dekin6 + dh6 + dps6 + dsr6 545 ! 546 call mp_sum( detot6, intra_bgrp_comm ) 547 ! 548 DO k = 1, 6 549 detmp( alpha(k), beta(k) ) = detot6(k) 550 detmp( beta(k), alpha(k) ) = detmp( alpha(k), beta(k) ) 551 END DO 552 ! 553 detot = MATMUL( detmp(:,:), TRANSPOSE( ainv(:,:) ) ) 554 ! 555 detot = detot + denl + dxc 556 ! 557 ! Add TS-vdW cell derivatives to detot here... (RAD) 558 ! 559 IF (ts_vdw) THEN 560 ! 561 detot = detot + HtsvdW 562 ! 563 ! BS / RAD start print ts_vdW pressure ------------- 564 ! 565 IF(MOD(nfi,iprint_stdout).EQ.0) THEN 566 detmp = HtsvdW 567 detmp = -1.d0 * (MATMUL( detmp(:,:), TRANSPOSE(h) )) / omega 568 detmp = detmp * au_gpa ! GPa 569 WRITE( stdout,9013) (detmp(1,1)+detmp(2,2)+detmp(3,3))/3.0_DP , nfi 570 9013 FORMAT (/,5X,'TS-vdW Pressure (GPa)',F15.5,I7) 571 END IF 572 ! 573 END IF 574 ! 575 ! BS / RAD / HK 576 ! Adding the stress tensor from exact exchange here. exx_wf related 577 ! 578 IF(dft_is_hybrid().AND.exx_is_active()) THEN 579 ! 580 IF (isotropic .and. (ibrav.eq.1)) THEN 581 ! 582 ! BS / RAD 583 ! This part is dE/dV; so works only for cubic cells and isotropic change 584 ! in simulation cell while doing variable cell calculation .. 585 ! dE/dV = -(1/3) * (-exx * 0.25_DP) / V 586 ! dexx(3,3) = dE/dh = (dE/dV) * (dV/dh) = (dE/dV) * V * (Transpose h)^-1 587 ! 588 DO k = 1, 6 589 IF(alpha(k) .EQ. beta(k)) THEN 590 detmp( alpha(k), beta(k) ) = - (1.0_DP/3.0_DP) * (-exx * exxalfa) 591 ELSE 592 detmp( alpha(k), beta(k) ) = 0.0_DP 593 END IF 594 detmp( beta(k), alpha(k) ) = detmp( alpha(k), beta(k) ) 595 END DO 596 ! 597 dexx = MATMUL( detmp(:,:), TRANSPOSE( ainv(:,:) ) ) 598 ! 599 ELSE 600 ! 601 ! HK: general case is computed in exx_gs.f90 602 ! (notice that the negative sign comes from the definition of 603 ! exchange energy being as positive value) 604 ! 605 dexx = -dexx_dh*exxalfa 606 ! 607 END IF 608 ! 609 detot = detot + dexx 610 ! 611 IF(MOD(nfi,iprint_stdout).EQ.0) WRITE( stdout,9014) (-1.0_DP/3.0_DP)*& 612 (dexx(1,1)+dexx(2,2)+dexx(3,3))*au_gpa , nfi 613 9014 FORMAT (5X,'EXX Pressure (GPa)',F15.5,I7) 614 ! 615 END IF 616 ! BS / RAD : stress tensor from exx ends here ... 617 ! 618 ! BS / RAD start print total electronic pressure ------------- 619 ! 620 IF(MOD(nfi,iprint_stdout).EQ.0) THEN 621 detmp = detot 622 detmp = -1.d0 * (MATMUL( detmp(:,:), TRANSPOSE(h) )) / omega 623 detmp = detmp * au_gpa ! GPa 624 WRITE( stdout,9015) (detmp(1,1)+detmp(2,2)+detmp(3,3))/3.0_DP , nfi 625 9015 FORMAT (5X,'Total Electronic Pressure (GPa)',F15.5,I7) 626 END IF 627 ! 628 END IF 629 ! 630 CALL stop_clock( 'vofrho' ) 631 ! 632 IF ( tpre ) THEN 633 ! 634 IF( ( iverbosity > 1 ) .AND. ( MOD( nfi - 1, iprint) == 0 ) ) THEN 635 ! 636 WRITE( stdout,*) 637 WRITE( stdout,*) "From vofrho:" 638 WRITE( stdout,*) "cell parameters h" 639 WRITE( stdout,5555) (at(i,1)*alat, at(i,2)*alat, at(i,3)*alat,i=1,3) 640 ! 641 WRITE( stdout,*) 642 WRITE( stdout,*) "derivative of e(tot)" 643 WRITE( stdout,5555) ((detot(i,j),j=1,3),i=1,3) 644 WRITE( stdout,*) "kbar" 645 detmp = -1.0d0 * MATMUL( detot, TRANSPOSE( h ) ) / omega * au_gpa * 10.0d0 646 WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3) 647 ! 648 WRITE( stdout,*) 649 WRITE( stdout,*) "derivative of e(kin)" 650 WRITE( stdout,5555) ((dekin(i,j),j=1,3),i=1,3) 651 WRITE( stdout,*) "kbar" 652 detmp = -1.0d0 * MATMUL( dekin, TRANSPOSE( h ) ) / omega * au_gpa * 10.0d0 653 WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3) 654 ! 655 WRITE( stdout,*) "derivative of e(h)" 656 WRITE( stdout,5555) ((dh(i,j),j=1,3),i=1,3) 657 WRITE( stdout,*) "kbar" 658 detmp = -1.0d0 * MATMUL( dh, TRANSPOSE( h ) ) / omega * au_gpa * 10.0d0 659 WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3) 660 ! 661 WRITE( stdout,*) "derivative of e(sr)" 662 WRITE( stdout,5555) ((dsr(i,j),j=1,3),i=1,3) 663 WRITE( stdout,*) "kbar" 664 detmp = -1.0d0 * MATMUL( dsr, TRANSPOSE( h ) ) / omega * au_gpa * 10.0d0 665 WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3) 666 ! 667 WRITE( stdout,*) "derivative of e(ps)" 668 WRITE( stdout,5555) ((dps(i,j),j=1,3),i=1,3) 669 WRITE( stdout,*) "kbar" 670 detmp = -1.0d0 * MATMUL( dps, TRANSPOSE( h ) ) / omega * au_gpa * 10.0d0 671 WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3) 672 ! 673 WRITE( stdout,*) "derivative of e(nl)" 674 WRITE( stdout,5555) ((denl(i,j),j=1,3),i=1,3) 675 WRITE( stdout,*) "kbar" 676 detmp = -1.0d0 * MATMUL( denl, TRANSPOSE( h ) ) / omega * au_gpa * 10.0d0 677 WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3) 678 ! 679 WRITE( stdout,*) "derivative of e(xc)" 680 WRITE( stdout,5555) ((dxc(i,j),j=1,3),i=1,3) 681 WRITE( stdout,*) "kbar" 682 detmp = -1.0d0 * MATMUL( dxc, TRANSPOSE( h ) ) / omega * au_gpa * 10.0d0 683 WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3) 684 ! 685 IF (ts_vdw) THEN 686 WRITE( stdout,*) "derivative of e(TS-vdW)" 687 WRITE( stdout,5555) ((HtsvdW(i,j),j=1,3),i=1,3) 688 WRITE( stdout,*) "kbar" 689 detmp = -1.0d0 * MATMUL( HtsvdW, TRANSPOSE( h ) ) / omega * au_gpa * 10.0d0 690 WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3) 691 END IF 692 ENDIF 693 ENDIF 694 695 RETURN 696 6975555 FORMAT(1x,f12.5,1x,f12.5,1x,f12.5/ & 698 & 1x,f12.5,1x,f12.5,1x,f12.5/ & 699 & 1x,f12.5,1x,f12.5,1x,f12.5//) 700! 701 702END SUBROUTINE vofrho_x 703