1! 2! Copyright (C) 2002-2005 FPMD-CPV groups 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 SUBROUTINE pstress_conv( de3x3, de6, ainv ) 11!------------------------------------------------------------------------------! 12 13 USE kinds, ONLY: DP 14 USE mp_global, ONLY: intra_bgrp_comm 15 USE mp, ONLY: mp_sum 16 USE stress_param, ONLY: alpha, beta 17 18 IMPLICIT NONE 19 20 REAL(DP) :: de3x3(3,3) 21 REAL(DP), INTENT(IN) :: de6(6) 22 REAL(DP), INTENT(IN) :: ainv(3,3) 23 REAL(DP) :: tmp(3,3) 24 25 INTEGER :: k 26 27 DO k = 1, 6 28 tmp( alpha(k), beta(k) ) = de6(k) 29 tmp( beta(k), alpha(k) ) = tmp(alpha(k),beta(k)) 30 END DO 31 32 de3x3 = MATMUL( tmp(:,:), TRANSPOSE( ainv(:,:) ) ) 33 34 CALL mp_sum( de3x3, intra_bgrp_comm ) 35 36 37 RETURN 38 END SUBROUTINE 39 40 41 42!------------------------------------------------------------------------------! 43 SUBROUTINE pseudo_stress_x( deps, epseu, gagb, sfac, dvps, rhoeg, omega ) 44!------------------------------------------------------------------------------! 45 ! 46 USE kinds, ONLY: DP 47 USE ions_base, ONLY: nsp 48 USE electrons_base, ONLY: nspin 49 USE stress_param, ONLY: dalbe 50 USE cp_interfaces, ONLY: stress_local 51 USE fft_base, ONLY: dffts 52 53 IMPLICIT NONE 54 55 REAL(DP), INTENT(IN) :: omega 56 REAL(DP), INTENT(OUT) :: deps(:) 57 REAL(DP), INTENT(IN) :: gagb(:,:) 58 COMPLEX(DP), INTENT(IN) :: rhoeg(:,:) 59 COMPLEX(DP), INTENT(IN) :: sfac(:,:) 60 REAL(DP), INTENT(IN) :: dvps(:,:) 61 REAL(DP), INTENT(IN) :: epseu 62 63 INTEGER :: k 64 COMPLEX(DP) :: rhets, depst(6) 65 COMPLEX(DP), ALLOCATABLE :: rhoe( : ) 66 COMPLEX(DP), ALLOCATABLE :: drhoe( :, : ) 67 ! 68 ALLOCATE( drhoe( dffts%ngm, 6 ), rhoe( dffts%ngm ) ) 69 70 rhoe( 1:dffts%ngm ) = rhoeg( 1:dffts%ngm, 1 ) 71 IF( nspin > 1 ) rhoe( 1:dffts%ngm ) = rhoe( 1:dffts%ngm ) + rhoeg( 1:dffts%ngm, 2 ) 72 73 DO k = 1, 6 74 drhoe( 1:dffts%ngm, k ) = - rhoe( 1:dffts%ngm ) * dalbe( k ) 75 END DO 76 77 CALL stress_local( deps, epseu, gagb, sfac, rhoe, drhoe, omega ) 78 79 DEALLOCATE( drhoe, rhoe ) 80 81 RETURN 82 END SUBROUTINE pseudo_stress_x 83 84 85 86!------------------------------------------------------------------------------! 87 SUBROUTINE stress_local_x( deps, epseu, gagb, sfac, rhoe, drhoe, omega ) 88!------------------------------------------------------------------------------! 89 ! 90 USE kinds, ONLY: DP 91 USE ions_base, ONLY: nsp 92 USE gvect, ONLY: gstart 93 USE electrons_base, ONLY: nspin 94 USE local_pseudo, ONLY: vps, dvps 95 USE fft_base, ONLY: dffts 96 97 IMPLICIT NONE 98 99 REAL(DP), INTENT(IN) :: omega 100 REAL(DP), INTENT(OUT) :: deps(:) 101 REAL(DP), INTENT(IN) :: gagb(:,:) 102 COMPLEX(DP), INTENT(IN) :: rhoe(:) 103 COMPLEX(DP), INTENT(IN) :: drhoe(:,:) 104 COMPLEX(DP), INTENT(IN) :: sfac(:,:) 105 REAL(DP), INTENT(IN) :: epseu 106 107 INTEGER :: ig,k,is, ispin 108 COMPLEX(DP) :: dsvp, svp, depst(6) 109 REAL(DP) :: wz 110 ! 111 depst = (0.d0,0.d0) 112 113 wz = 2.0d0 114 115 DO ig = gstart, dffts%ngm 116 svp = 0.0d0 117 DO is = 1, nsp 118 svp = svp + sfac( ig, is ) * vps( ig, is ) 119 END DO 120 depst = depst + wz * CONJG( drhoe( ig, : ) ) * svp 121 END DO 122 IF( gstart == 2 ) THEN 123 svp = 0.0d0 124 DO is = 1, nsp 125 svp = svp + sfac( 1, is ) * vps( 1, is ) 126 END DO 127 depst = depst + CONJG( drhoe( 1, : ) ) * svp 128 END IF 129 130 DO ig = gstart, dffts%ngm 131 dsvp = 0.0d0 132 DO is = 1, nsp 133 dsvp = dsvp + sfac( ig, is ) * dvps( ig, is ) 134 END DO 135 DO k = 1, 6 136 depst( k ) = depst( k ) - wz * 2.0d0 * CONJG( rhoe( ig ) ) * dsvp * gagb( k, ig ) 137 END DO 138 END DO 139 140 deps = omega * DBLE( depst ) 141 142 RETURN 143 END SUBROUTINE stress_local_x 144 145 146 147 148!------------------------------------------------------------------------------! 149 SUBROUTINE stress_kin_x( dekin, c0_bgrp, occ_bgrp ) 150!------------------------------------------------------------------------------! 151 152! this routine computes the kinetic energy contribution to the stress 153! tensor 154! 155! dekin(:) = - 2 (sum over i) f(i) * 156! ( (sum over g) gagb(:,g) CONJG( c0(g,i) ) c0(g,i) ) 157! 158 159 USE kinds, ONLY: DP 160 USE gvecw, ONLY: q2sigma, ecfixed, qcutz, ngw 161 USE constants, ONLY: pi 162 USE gvect, ONLY: gstart, gg, g 163 USE cell_base, ONLY: tpiba2 164 USE electrons_base, ONLY: nspin, iupdwn_bgrp, nupdwn_bgrp 165 USE stress_param, ONLY: alpha, beta 166 167 IMPLICIT NONE 168 169! ... declare subroutine arguments 170 REAL(DP), INTENT(OUT) :: dekin(:) 171 COMPLEX(DP), INTENT(IN) :: c0_bgrp(:,:) 172 REAL(DP), INTENT(IN) :: occ_bgrp(:) 173 174! ... declare other variables 175 REAL(DP) :: sk(6), scg, efac 176 REAL(DP), ALLOCATABLE :: arg(:) 177 INTEGER :: ib, ig, ispin, iwfc 178 179! ... end of declarations 180! ---------------------------------------------- 181 182 dekin = 0.0_DP 183 ALLOCATE( arg( ngw ) ) 184 185 efac = 2.0d0 * qcutz / q2sigma / SQRT(pi) 186 IF( efac > 0.0d0 ) THEN 187 DO ig = gstart, ngw 188 arg(ig) = 1.0d0 + efac * exp( -( ( tpiba2 *gg(ig) - ecfixed ) / q2sigma )**2 ) 189 END DO 190 ELSE 191 arg = 1.0d0 192 END IF 193 194 ! ... compute kinetic energy contribution 195 196 DO ispin = 1, nspin 197 DO ib = 1, nupdwn_bgrp( ispin ) 198 sk = 0.0_DP 199 iwfc = ib + iupdwn_bgrp( ispin ) - 1 200 DO ig = gstart, ngw 201 scg = arg(ig) * CONJG( c0_bgrp( ig, iwfc ) ) * c0_bgrp( ig, iwfc ) 202 sk(1) = sk(1) + scg * g( alpha( 1 ), ig ) * g( beta( 1 ), ig ) 203 sk(2) = sk(2) + scg * g( alpha( 2 ), ig ) * g( beta( 2 ), ig ) 204 sk(3) = sk(3) + scg * g( alpha( 3 ), ig ) * g( beta( 3 ), ig ) 205 sk(4) = sk(4) + scg * g( alpha( 4 ), ig ) * g( beta( 4 ), ig ) 206 sk(5) = sk(5) + scg * g( alpha( 5 ), ig ) * g( beta( 5 ), ig ) 207 sk(6) = sk(6) + scg * g( alpha( 6 ), ig ) * g( beta( 6 ), ig ) 208 END DO 209 dekin = dekin + occ_bgrp( iwfc ) * sk * tpiba2 210 END DO 211 END DO 212 dekin = - 2.0_DP * dekin 213 DEALLOCATE(arg) 214 RETURN 215 END SUBROUTINE stress_kin_x 216 217 218 219 220!------------------------------------------------------------------------------! 221 SUBROUTINE add_drhoph_x( drhot, sfac, gagb ) 222!------------------------------------------------------------------------------! 223 ! 224 USE kinds, ONLY: DP 225 USE ions_base, ONLY: nsp, rcmax 226 USE local_pseudo, ONLY: rhops 227 USE stress_param, ONLY: dalbe 228 USE fft_base, ONLY: dffts 229 ! 230 IMPLICIT NONE 231 ! 232 COMPLEX(DP), INTENT(INOUT) :: drhot( :, : ) 233 COMPLEX(DP), INTENT(IN) :: sfac( :, : ) 234 REAL(DP), INTENT(IN) :: gagb( :, : ) 235 ! 236 INTEGER :: ij, is, ig 237 COMPLEX(DP) :: drhop 238 ! 239 DO ij = 1, 6 240 IF( dalbe( ij ) > 0.0d0 ) THEN 241 DO is = 1, nsp 242 DO ig = 1, dffts%ngm 243 drhot(ig,ij) = drhot(ig,ij) - sfac(ig,is)*rhops(ig,is) 244 ENDDO 245 END DO 246 END IF 247 END DO 248 DO ig = 1, dffts%ngm 249 drhop = 0.0d0 250 DO is = 1, nsp 251 drhop = drhop - sfac( ig, is ) * rhops(ig,is) * rcmax(is)**2 * 0.5D0 252 END DO 253 DO ij = 1, 6 254 drhot(ig,ij) = drhot(ig,ij) - drhop * gagb( ij, ig ) 255 END DO 256 END DO 257 RETURN 258 END SUBROUTINE add_drhoph_x 259 260 261 262 263!------------------------------------------------------------------------------! 264 SUBROUTINE stress_har_x(deht, ehr, sfac, rhoeg, gagb, omega ) 265!------------------------------------------------------------------------------! 266 267 use kinds, only: DP 268 use ions_base, only: nsp, rcmax 269 use mp_global, ONLY: me_bgrp, root_bgrp 270 USE constants, ONLY: fpi 271 USE cell_base, ONLY: tpiba2 272 USE gvect, ONLY: gstart 273 USE local_pseudo, ONLY: rhops 274 USE electrons_base, ONLY: nspin 275 USE stress_param, ONLY: dalbe 276 USE cp_interfaces, ONLY: add_drhoph, stress_hartree 277 USE fft_base, ONLY: dffts, dfftp 278 279 IMPLICIT NONE 280 281 REAL(DP), INTENT(OUT) :: DEHT(:) 282 REAL(DP), INTENT(IN) :: omega, EHR, gagb(:,:) 283 COMPLEX(DP), INTENT(IN) :: RHOEG(:,:) 284 COMPLEX(DP), INTENT(IN) :: sfac(:,:) 285 286 COMPLEX(DP) DEHC(6) 287 COMPLEX(DP) RHOP,DRHOP 288 COMPLEX(DP) RHET,RHOG,RHETS,RHOGS 289 COMPLEX(DP) CFACT 290 COMPLEX(DP), ALLOCATABLE :: rhot(:), drhot(:,:) 291 REAL(DP) hgm1 292 293 INTEGER ig, is, k, ispin 294 295 296 ALLOCATE( rhot( dfftp%ngm ) ) 297 ALLOCATE( drhot( dfftp%ngm, 6 ) ) 298 299 ! sum up spin components 300 ! 301 rhot( gstart:dfftp%ngm ) = rhoeg( gstart:dfftp%ngm, 1 ) 302 IF( nspin > 1 ) rhot( gstart:dfftp%ngm ) = rhot( gstart:dfftp%ngm ) + rhoeg( gstart:dfftp%ngm, 2 ) 303 ! 304 ! add Ionic pseudo charges rho_I 305 ! 306 DO is = 1, nsp 307 rhot( gstart:dffts%ngm ) = rhot( gstart:dffts%ngm ) + sfac( gstart:dffts%ngm, is ) * rhops( gstart:dffts%ngm, is ) 308 END DO 309 310 ! add drho_e / dh 311 ! 312 DO k = 1, 6 313 IF( dalbe( k ) > 0.0d0 ) THEN 314 drhot( :, k ) = - rhoeg( :, 1 ) 315 IF( nspin > 1 ) drhot( :, k ) = drhot( :, k ) + rhoeg( :, 2 ) 316 ELSE 317 drhot( :, k ) = 0.0d0 318 END IF 319 END DO 320 321 ! add drho_I / dh 322 ! 323 CALL add_drhoph( drhot, sfac, gagb ) 324 325 CALL stress_hartree(deht, ehr, sfac, rhot, drhot, gagb, omega ) 326 327 DEALLOCATE( rhot, drhot ) 328 329 RETURN 330 END SUBROUTINE stress_har_x 331 332 333 334 335!------------------------------------------------------------------------------! 336 SUBROUTINE stress_hartree_x(deht, ehr, sfac, rhot, drhot, gagb, omega ) 337!------------------------------------------------------------------------------! 338 339 ! This subroutine computes: d E_hartree / dh = 340 ! E_hartree * h^t + 341 ! 4pi omega rho_t * CONJG( rho_t ) / G^2 / G^2 * G_alpha * G_beta + 342 ! 4pi omega Re{ CONJG( rho_t ) * drho_t / G^2 } 343 ! where: 344 ! rho_t = rho_e + rho_I 345 ! drho_t = d rho_t / dh = -rho_e + d rho_hard / dh + d rho_I / dh 346 347 use kinds, only: DP 348 use ions_base, only: nsp, rcmax 349 use mp_global, ONLY: me_bgrp, root_bgrp 350 USE constants, ONLY: fpi 351 USE cell_base, ONLY: tpiba2 352 USE gvect, ONLY: gstart, gg 353 USE local_pseudo, ONLY: rhops 354 USE electrons_base, ONLY: nspin 355 USE stress_param, ONLY: dalbe 356 USE fft_base, ONLY: dfftp 357 358 IMPLICIT NONE 359 360 REAL(DP), INTENT(OUT) :: DEHT(:) 361 REAL(DP), INTENT(IN) :: omega, EHR, gagb(:,:) 362 COMPLEX(DP) :: rhot(:) ! total charge: Sum_spin ( rho_e + rho_I ) 363 COMPLEX(DP) :: drhot(:,:) 364 COMPLEX(DP), INTENT(IN) :: sfac(:,:) 365 366 COMPLEX(DP) DEHC(6) 367 COMPLEX(DP) CFACT 368 REAL(DP), ALLOCATABLE :: hgm1( : ) 369 REAL(DP) :: wz 370 371 INTEGER ig, is, k, iss 372 373 DEHC = (0.D0,0.D0) 374 DEHT = 0.D0 375 376 wz = 2.0d0 377 378 ALLOCATE( hgm1( dfftp%ngm ) ) 379 380 hgm1( 1 ) = 0.0d0 381 DO ig = gstart, dfftp%ngm 382 hgm1( ig ) = 1.D0 / gg(ig) / tpiba2 383 END DO 384 385 ! Add term rho_t * CONJG( rho_t ) / G^2 * G_alpha * G_beta / G^2 386 387 DO ig = gstart, dfftp%ngm 388 cfact = rhot( ig ) * CONJG( rhot( ig ) ) * hgm1( ig ) ** 2 389 dehc = dehc + cfact * gagb(:,ig) 390 END DO 391 392 ! Add term 2 * Re{ CONJG( rho_t ) * drho_t / G^2 } 393 394 DO ig = gstart, dfftp%ngm 395 DO k = 1, 6 396 dehc( k ) = dehc( k ) + rhot( ig ) * CONJG( drhot( ig, k ) ) * hgm1( ig ) 397 END DO 398 END DO 399 400 ! term: E_h * h^t 401 402 if ( me_bgrp == root_bgrp ) then 403 deht = wz * fpi * omega * DBLE(dehc) + ehr * dalbe 404 else 405 deht = wz * fpi * omega * DBLE(dehc) 406 end if 407 408 DEALLOCATE( hgm1 ) 409 410 RETURN 411 END SUBROUTINE stress_hartree_x 412 413 414 415!------------------------------------------------------------------------------! 416 SUBROUTINE stress_debug(dekin, deht, dexc, desr, deps, denl, htm1) 417!------------------------------------------------------------------------------! 418 419 USE kinds, ONLY: DP 420 USE io_global, ONLY: stdout 421 USE mp_global, ONLY: intra_bgrp_comm 422 USE mp, ONLY: mp_sum 423 USE stress_param, ONLY: alpha, beta 424 425 IMPLICIT NONE 426 427 REAL(DP) :: dekin(6), deht(6), dexc(6), desr(6), deps(6), denl(6) 428 REAL(DP) :: detot(6), htm1(3,3) 429 REAL(DP) :: detmp(3,3) 430 431 INTEGER :: k, i, j 432 433 detot = dekin + deht + dexc + desr + deps + denl 434 435 DO k=1,6 436 detmp(alpha(k),beta(k)) = detot(k) 437 detmp(beta(k),alpha(k)) = detmp(alpha(k),beta(k)) 438 END DO 439 CALL mp_sum( detmp, intra_bgrp_comm ) 440 detmp = MATMUL( detmp(:,:), htm1(:,:) ) 441 WRITE( stdout,*) "derivative of e(tot)" 442 WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3) 443 444 DO k=1,6 445 detmp(alpha(k),beta(k)) = dekin(k) 446 detmp(beta(k),alpha(k)) = detmp(alpha(k),beta(k)) 447 END DO 448 CALL mp_sum( detmp, intra_bgrp_comm ) 449 detmp = MATMUL( detmp(:,:), htm1(:,:) ) 450 WRITE( stdout,*) "derivative of e(kin)" 451 WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3) 452 453 DO k=1,6 454 detmp(alpha(k),beta(k)) = deht(k) + desr(k) 455 detmp(beta(k),alpha(k)) = detmp(alpha(k),beta(k)) 456 END DO 457 CALL mp_sum( detmp, intra_bgrp_comm ) 458 detmp = MATMUL( detmp(:,:), htm1(:,:) ) 459 WRITE( stdout,*) "derivative of e(electrostatic)" 460 WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3) 461 462 DO k=1,6 463 detmp(alpha(k),beta(k)) = deht(k) 464 detmp(beta(k),alpha(k)) = detmp(alpha(k),beta(k)) 465 END DO 466 CALL mp_sum( detmp, intra_bgrp_comm ) 467 detmp = MATMUL( detmp(:,:), htm1(:,:) ) 468 WRITE( stdout,*) "derivative of e(h)" 469 WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3) 470 471 DO k=1,6 472 detmp(alpha(k),beta(k)) = desr(k) 473 detmp(beta(k),alpha(k)) = detmp(alpha(k),beta(k)) 474 END DO 475 CALL mp_sum( detmp, intra_bgrp_comm ) 476 detmp = MATMUL( detmp(:,:), htm1(:,:) ) 477 WRITE( stdout,*) "derivative of e(sr)" 478 WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3) 479 480 DO k=1,6 481 detmp(alpha(k),beta(k)) = deps(k) 482 detmp(beta(k),alpha(k)) = detmp(alpha(k),beta(k)) 483 END DO 484 CALL mp_sum( detmp, intra_bgrp_comm ) 485 detmp = MATMUL( detmp(:,:), htm1(:,:) ) 486 WRITE( stdout,*) "derivative of e(ps)" 487 WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3) 488 489 DO k=1,6 490 detmp(alpha(k),beta(k)) = denl(k) 491 detmp(beta(k),alpha(k)) = detmp(alpha(k),beta(k)) 492 END DO 493 CALL mp_sum( detmp, intra_bgrp_comm ) 494 detmp = MATMUL( detmp(:,:), htm1(:,:) ) 495 WRITE( stdout,*) "derivative of e(nl)" 496 WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3) 497 498 DO k=1,6 499 detmp(alpha(k),beta(k)) = dexc(k) 500 detmp(beta(k),alpha(k)) = detmp(alpha(k),beta(k)) 501 END DO 502 CALL mp_sum( detmp, intra_bgrp_comm ) 503 detmp = MATMUL( detmp(:,:), htm1(:,:) ) 504 WRITE( stdout,*) "derivative of e(xc)" 505 WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3) 506 5075555 format(1x,f12.5,1x,f12.5,1x,f12.5/ & 508 & 1x,f12.5,1x,f12.5,1x,f12.5/ & 509 & 1x,f12.5,1x,f12.5,1x,f12.5//) 510 511 RETURN 512 END SUBROUTINE stress_debug 513 514 515 516 517!------------------------------------------------------------------------------! 518 SUBROUTINE compute_gagb_x( gagb, g, ngm, tpiba2 ) 519!------------------------------------------------------------------------------! 520 521 ! ... compute G_alpha * G_beta 522 523 USE kinds, ONLY: DP 524 USE stress_param, ONLY: alpha, beta 525 526 IMPLICIT NONE 527 528 INTEGER, INTENT(IN) :: ngm 529 REAL(DP), INTENT(IN) :: g(:,:) 530 REAL(DP), INTENT(OUT) :: gagb(:,:) 531 REAL(DP), INTENT(IN) :: tpiba2 532 533 INTEGER :: k, ig 534 535!$omp parallel do default(shared), private(k) 536 DO ig = 1, ngm 537 DO k = 1, 6 538 gagb( k, ig ) = g( alpha( k ), ig ) * g( beta( k ), ig ) * tpiba2 539 END DO 540 END DO 541 542 END SUBROUTINE compute_gagb_x 543