1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2! Copyright 2010. Los Alamos National Security, LLC. This material was ! 3! produced under U.S. Government contract DE-AC52-06NA25396 for Los Alamos ! 4! National Laboratory (LANL), which is operated by Los Alamos National ! 5! Security, LLC for the U.S. Department of Energy. The U.S. Government has ! 6! rights to use, reproduce, and distribute this software. NEITHER THE ! 7! GOVERNMENT NOR LOS ALAMOS NATIONAL SECURITY, LLC MAKES ANY WARRANTY, ! 8! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS ! 9! SOFTWARE. If software is modified to produce derivative works, such ! 10! modified software should be clearly marked, so as not to confuse it ! 11! with the version available from LANL. ! 12! ! 13! Additionally, this program is free software; you can redistribute it ! 14! and/or modify it under the terms of the GNU General Public License as ! 15! published by the Free Software Foundation; version 2.0 of the License. ! 16! Accordingly, this program is distributed in the hope that it will be ! 17! useful, but WITHOUT ANY WARRANTY; without even the implied warranty of ! 18! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General ! 19! Public License for more details. ! 20!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 21 22SUBROUTINE PULAY 23 24 USE CONSTANTS_MOD 25 USE SETUPARRAY 26 USE UNIVARRAY 27 USE NONOARRAY 28 USE NEBLISTARRAY 29 USE SPINARRAY 30 USE VIRIALARRAY 31 USE MYPRECISION 32 33 IMPLICIT NONE 34 35 INTEGER :: I, J, K, L, M, N, KK, INDI, INDJ 36 INTEGER :: LBRA, MBRA, LKET, MKET 37 INTEGER :: PREVJ, NEWJ 38 INTEGER :: PBCI, PBCJ, PBCK 39 INTEGER :: BASISI(5), BASISJ(5), LBRAINC, LKETINC 40 INTEGER :: SPININDI, SPININDJ 41 REAL(LATTEPREC) :: ALPHA, BETA, PHI, RHO, RHODIFF, COSBETA 42 REAL(LATTEPREC) :: RIJ(3), DC(3) 43 REAL(LATTEPREC) :: MAGR, MAGR2, MAGRP, MAGRP2 44 REAL(LATTEPREC) :: FTMP_PULAY(3), FTMP_COUL(3), FTMP_SPIN(3) 45 REAL(LATTEPREC) :: MAXRCUT, MAXRCUT2 46 REAL(LATTEPREC) :: WSPINI, WSPINJ 47 REAL(LATTEPREC) :: MYDFDA, MYDFDB, MYDFDR, RCUTTB 48 49 REAL(LATTEPREC), EXTERNAL :: DFDA, DFDB, DFDR 50 LOGICAL PATH 51 IF (EXISTERROR) RETURN 52 53 FPUL = ZERO 54 VIRPUL = ZERO 55 56 ! We first have to make the matrix S^-1 H rho = X^2 H rho 57 58 IF (SPINON .EQ. 0) THEN 59 60 IF (KBT .GT. 0.0000001) THEN 61 62#ifdef DOUBLEPREC 63 64#ifdef GPUON 65 66 ! XMAT * XMAT = S^-1 67 68 CALL mmlatte(HDIM, 0, 0, ONE, ZERO, XMAT, XMAT, X2HRHO) 69 70 ! S^-1 * H 71 72 CALL mmlatte(HDIM, 0, 0, ONE, ZERO, X2HRHO, H, NONOTMP) 73 74 ! (S^-1 * H)*RHO 75 76 CALL mmlatte(HDIM, 0, 0, ONE, ZERO, NONOTMP, BO, X2HRHO) 77 78#elif defined(GPUOFF) 79 80 ! XMAT * XMAT = S^-1 81 82 CALL DGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 83 XMAT, HDIM, XMAT, HDIM, ZERO, X2HRHO, HDIM) 84 85 ! S^-1 * H 86 87 CALL DGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 88 X2HRHO, HDIM, H, HDIM, ZERO, NONOTMP, HDIM) 89 90 ! (S^-1 * H)*RHO 91 92 CALL DGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 93 NONOTMP, HDIM, BO, HDIM, ZERO, X2HRHO, HDIM) 94 95#endif 96 97#elif defined(SINGLEPREC) 98 99 ! XMAT * XMAT = S^-1 100 101 CALL SGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 102 XMAT, HDIM, XMAT, HDIM, ZERO, X2HRHO, HDIM) 103 104 ! S^-1 * H 105 106 CALL SGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 107 X2HRHO, HDIM, H, HDIM, ZERO, NONOTMP, HDIM) 108 109 ! (S^-1 * H)*RHO 110 111 CALL SGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 112 NONOTMP, HDIM, BO, HDIM, ZERO, X2HRHO, HDIM) 113 114#endif 115 116 ELSE 117 118 ! Te = 0 : Fp = 2Tr[rho H rho dS/dR] 119 120 ! Be careful - we're working with bo = 2rho, so we need 121 ! the factor of 1/2... 122 123#ifdef DOUBLEPREC 124 125#ifdef GPUON 126 127 CALL mmlatte(HDIM, 0, 0, ONE, ZERO, BO, H, NONOTMP) 128 129 CALL mmlatte(HDIM, 0, 0, HALF, ZERO, NONOTMP, BO, X2HRHO) 130 131#elif defined(GPUOFF) 132 133 134 CALL DGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 135 BO, HDIM, H, HDIM, ZERO, NONOTMP, HDIM) 136 137 CALL DGEMM('N', 'N', HDIM, HDIM, HDIM, HALF, & 138 NONOTMP, HDIM, BO, HDIM, ZERO, X2HRHO, HDIM) 139 140#endif 141 142#elif defined(SINGLEPREC) 143 144 CALL SGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 145 BO, HDIM, H, HDIM, ZERO, NONOTMP, HDIM) 146 147 CALL SGEMM('N', 'N', HDIM, HDIM, HDIM, HALF, & 148 NONOTMP, HDIM, BO, HDIM, ZERO, X2HRHO, HDIM) 149 150#endif 151 152 ENDIF 153 154 ELSE 155 156 IF (KBT .GT. 0.0000001) THEN 157 158#ifdef DOUBLEPREC 159 160#ifdef GPUON 161 162 ! XMAT * XMAT = S^-1 163 164 CALL mmlatte(HDIM, 0, 0, ONE, ZERO, XMAT, XMAT, SPINTMP) 165 166 ! Hup * rhoup 167 168 CALL mmlatte(HDIM, 0, 0, ONE, ZERO, HUP, RHOUP, NONOTMP) 169 170 ! (Hup * rhoup) + Hdown*rhodown 171 172 CALL mmlatte(HDIM, 0, 0, ONE, ONE, HDOWN, RHODOWN, NONOTMP) 173 174 CALL mmlatte(HDIM, 0, 0, ONE, ZERO, SPINTMP, NONOTMP, X2HRHO) 175 176#elif defined(GPUOFF) 177 178 179 ! XMAT * XMAT = S^-1 180 181 CALL DGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 182 XMAT, HDIM, XMAT, HDIM, ZERO, SPINTMP, HDIM) 183 184 ! Hup * rhoup 185 186 CALL DGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 187 HUP, HDIM, RHOUP, HDIM, ZERO, NONOTMP, HDIM) 188 189 ! (Hup * rhoup) + Hdown*rhodown 190 191 CALL DGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 192 HDOWN, HDIM, RHODOWN, HDIM, ONE, NONOTMP, HDIM) 193 194 195 CALL DGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 196 SPINTMP, HDIM, NONOTMP, HDIM, ZERO, X2HRHO, HDIM) 197 198#endif 199 200#elif defined(SINGLEPREC) 201 202 ! XMAT * XMAT = S^-1 203 204 CALL SGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 205 XMAT, HDIM, XMAT, HDIM, ZERO, SPINTMP, HDIM) 206 207 ! Hup * rhoup 208 209 CALL SGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 210 HUP, HDIM, RHOUP, HDIM, ZERO, NONOTMP, HDIM) 211 212 ! (Hup * rhoup) + Hdown*rhodown 213 214 CALL SGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 215 HDOWN, HDIM, RHODOWN, HDIM, ONE, NONOTMP, HDIM) 216 217 218 CALL SGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 219 SPINTMP, HDIM, NONOTMP, HDIM, ZERO, X2HRHO, HDIM) 220 221#endif 222 223 ELSE 224 225#ifdef DOUBLEPREC 226 227 ! At zero K: tr(rho H rho) 228 229#ifdef GPUON 230 231 CALL mmlatte(HDIM, 0, 0, ONE, ZERO, RHOUP, HUP, NONOTMP) 232 233 CALL mmlatte(HDIM, 0, 0, ONE, ZERO, NONOTMP, RHOUP, X2HRHO) 234 235 CALL mmlatte(HDIM, 0, 0, ONE, ZERO, RHODOWN, HDOWN, NONOTMP) 236 237 CALL mmlatte(HDIM, 0, 0, ONE, ONE, NONOTMP, RHODOWN, X2HRHO) 238 239#elif defined(GPUOFF) 240 241 CALL DGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 242 RHOUP, HDIM, HUP, HDIM, ZERO, NONOTMP, HDIM) 243 244 CALL DGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 245 NONOTMP, HDIM, RHOUP, HDIM, ZERO, X2HRHO, HDIM) 246 247 CALL DGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 248 RHODOWN, HDIM, HDOWN, HDIM, ZERO, NONOTMP, HDIM) 249 250 CALL DGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 251 NONOTMP, HDIM, RHODOWN, HDIM, ONE, X2HRHO, HDIM) 252 253#endif 254 255#elif defined(SINGLEPREC) 256 257 CALL SGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 258 RHOUP, HDIM, HUP, HDIM, ZERO, NONOTMP, HDIM) 259 260 CALL SGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 261 NONOTMP, HDIM, RHOUP, HDIM, ZERO, X2HRHO, HDIM) 262 263 CALL SGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 264 RHODOWN, HDIM, HDOWN, HDIM, ZERO, NONOTMP, HDIM) 265 266 CALL SGEMM('N', 'N', HDIM, HDIM, HDIM, ONE, & 267 NONOTMP, HDIM, RHODOWN, HDIM, ONE, X2HRHO, HDIM) 268 269#endif 270 271 ENDIF 272 273 ENDIF 274 275 276!$OMP PARALLEL DO DEFAULT (NONE) & 277!$OMP SHARED(NATS, BASIS, ELEMPOINTER, TOTNEBTB, NEBTB) & 278!$OMP SHARED(CR, BOX, X2HRHO, NOINT, ATELE, ELE1, ELE2) & 279!$OMP SHARED(BO, RHOUP, RHODOWN, SPINON, ELECTRO) & 280!$OMP SHARED(HCUT, SCUT, MATINDLIST, BASISTYPE) & 281!$OMP SHARED(DELTASPIN, WSS, WPP, WDD, WFF, SPININDLIST) & 282!$OMP SHARED(HUBBARDU, DELTAQ, COULOMBV) & 283!$OMP SHARED(LCNSHIFT) & 284!$OMP SHARED(FPUL) & 285!$OMP PRIVATE(I, J, K, NEWJ, BASISI, BASISJ, INDI, INDJ, PBCI, PBCJ, PBCK) & 286!$OMP PRIVATE(RIJ, MAGR2, MAGR, MAGRP2, MAGRP, PATH, PHI, ALPHA, BETA, COSBETA) & 287!$OMP PRIVATE(FTMP_PULAY, FTMP_COUL, FTMP_SPIN) & 288!$OMP PRIVATE(DC, LBRAINC, LBRA, MBRA, L, LKETINC, LKET, MKET, RHO, RHODIFF) & 289!$OMP PRIVATE(MYDFDA, MYDFDB, MYDFDR, RCUTTB) & 290!$OMP PRIVATE(WSPINI, WSPINJ, SPININDI, SPININDJ) & 291!$OMP REDUCTION(+:VIRPUL) 292 293 294 DO I = 1, NATS 295 296 ! Build list of orbitals on atom I 297 SELECT CASE(BASIS(ELEMPOINTER(I))) 298 299 CASE("s") 300 BASISI(1) = 0 301 BASISI(2) = -1 302 CASE("p") 303 BASISI(1) = 1 304 BASISI(2) = -1 305 CASE("d") 306 BASISI(1) = 2 307 BASISI(2) = -1 308 CASE("f") 309 BASISI(1) = 3 310 BASISI(2) = -1 311 CASE("sp") 312 BASISI(1) = 0 313 BASISI(2) = 1 314 BASISI(3) = -1 315 CASE("sd") 316 BASISI(1) = 0 317 BASISI(2) = 2 318 BASISI(3) = -1 319 CASE("sf") 320 BASISI(1) = 0 321 BASISI(2) = 3 322 BASISI(3) = -1 323 CASE("pd") 324 BASISI(1) = 1 325 BASISI(2) = 2 326 BASISI(3) = -1 327 CASE("pf") 328 BASISI(1) = 1 329 BASISI(2) = 3 330 BASISI(3) = -1 331 CASE("df") 332 BASISI(1) = 2 333 BASISI(2) = 3 334 BASISI(3) = -1 335 CASE("spd") 336 BASISI(1) = 0 337 BASISI(2) = 1 338 BASISI(3) = 2 339 BASISI(4) = -1 340 CASE("spf") 341 BASISI(1) = 0 342 BASISI(2) = 1 343 BASISI(3) = 3 344 BASISI(4) = -1 345 CASE("sdf") 346 BASISI(1) = 0 347 BASISI(2) = 2 348 BASISI(3) = 3 349 BASISI(4) = -1 350 CASE("pdf") 351 BASISI(1) = 1 352 BASISI(2) = 2 353 BASISI(3) = 3 354 BASISI(4) = -1 355 CASE("spdf") 356 BASISI(1) = 0 357 BASISI(2) = 1 358 BASISI(3) = 2 359 BASISI(4) = 3 360 BASISI(5) = -1 361 END SELECT 362 363 INDI = MATINDLIST(I) 364 IF (SPINON .EQ. 1) SPININDI = SPININDLIST(I) 365 366 DO NEWJ = 1, TOTNEBTB(I) 367 368 J = NEBTB(1, NEWJ, I) 369 PBCI = NEBTB(2, NEWJ, I) 370 PBCJ = NEBTB(3, NEWJ, I) 371 PBCK = NEBTB(4, NEWJ, I) 372 373 RIJ(1) = CR(1,J) + REAL(PBCI)*BOX(1,1) + REAL(PBCJ)*BOX(2,1) + & 374 REAL(PBCK)*BOX(3,1) - CR(1,I) 375 376 RIJ(2) = CR(2,J) + REAL(PBCI)*BOX(1,2) + REAL(PBCJ)*BOX(2,2) + & 377 REAL(PBCK)*BOX(3,2) - CR(2,I) 378 379 RIJ(3) = CR(3,J) + REAL(PBCI)*BOX(1,3) + REAL(PBCJ)*BOX(2,3) + & 380 REAL(PBCK)*BOX(3,3) - CR(3,I) 381 382 MAGR2 = RIJ(1)*RIJ(1) + RIJ(2)*RIJ(2) + RIJ(3)*RIJ(3) 383 384 ! Building the forces is expensive - use the cut-off 385 386 RCUTTB = ZERO 387 DO K = 1, NOINT 388 389 IF ( (ATELE(I) .EQ. ELE1(K) .AND. & 390 ATELE(J) .EQ. ELE2(K)) .OR. & 391 (ATELE(J) .EQ. ELE1(K) .AND. & 392 ATELE(I) .EQ. ELE2(K) )) THEN 393 394 IF (SCUT(K) .GT. RCUTTB ) RCUTTB = SCUT(K) 395 396 ENDIF 397 398 ENDDO 399 400 IF (MAGR2 .LT. RCUTTB*RCUTTB) THEN 401 402 MAGR = SQRT(MAGR2) 403 404 ! Build list of orbitals on atom J 405 406 SELECT CASE(BASIS(ELEMPOINTER(J))) 407 CASE("s") 408 BASISJ(1) = 0 409 BASISJ(2) = -1 410 CASE("p") 411 BASISJ(1) = 1 412 BASISJ(2) = -1 413 CASE("d") 414 BASISJ(1) = 2 415 BASISJ(2) = -1 416 CASE("f") 417 BASISJ(1) = 3 418 BASISJ(2) = -1 419 CASE("sp") 420 BASISJ(1) = 0 421 BASISJ(2) = 1 422 BASISJ(3) = -1 423 CASE("sd") 424 BASISJ(1) = 0 425 BASISJ(2) = 2 426 BASISJ(3) = -1 427 CASE("sf") 428 BASISJ(1) = 0 429 BASISJ(2) = 3 430 BASISJ(3) = -1 431 CASE("pd") 432 BASISJ(1) = 1 433 BASISJ(2) = 2 434 BASISJ(3) = -1 435 CASE("pf") 436 BASISJ(1) = 1 437 BASISJ(2) = 3 438 BASISJ(3) = -1 439 CASE("df") 440 BASISJ(1) = 2 441 BASISJ(2) = 3 442 BASISJ(3) = -1 443 CASE("spd") 444 BASISJ(1) = 0 445 BASISJ(2) = 1 446 BASISJ(3) = 2 447 BASISJ(4) = -1 448 CASE("spf") 449 BASISJ(1) = 0 450 BASISJ(2) = 1 451 BASISJ(3) = 3 452 BASISJ(4) = -1 453 CASE("sdf") 454 BASISJ(1) = 0 455 BASISJ(2) = 2 456 BASISJ(3) = 3 457 BASISJ(4) = -1 458 CASE("pdf") 459 BASISJ(1) = 1 460 BASISJ(2) = 2 461 BASISJ(3) = 3 462 BASISJ(4) = -1 463 CASE("spdf") 464 BASISJ(1) = 0 465 BASISJ(2) = 1 466 BASISJ(3) = 2 467 BASISJ(4) = 3 468 BASISJ(5) = -1 469 END SELECT 470 471 INDJ = MATINDLIST(J) 472 IF (SPINON .EQ. 1) SPININDJ = SPININDLIST(J) 473 474 475 MAGRP2 = RIJ(1)*RIJ(1) + RIJ(2)*RIJ(2) 476 MAGRP = SQRT(MAGRP2) 477 478 ! transform to system in which z-axis is aligned with RIJ 479 480 PATH = .FALSE. 481 IF (ABS(RIJ(1)) .GT. 1E-12) THEN 482 IF (RIJ(1) .GT. ZERO .AND. RIJ(2) .GE. ZERO) THEN 483 PHI = ZERO 484 ELSEIF (RIJ(1) .GT. ZERO .AND. RIJ(2) .LT. ZERO) THEN 485 PHI = TWO * PI 486 ELSE 487 PHI = PI 488 ENDIF 489 ALPHA = ATAN(RIJ(2) / RIJ(1)) + PHI 490 ELSEIF (ABS(RIJ(2)) .GT. 1E-12) THEN 491 IF (RIJ(2) .GT. 1E-12) THEN 492 ALPHA = PI / TWO 493 ELSE 494 ALPHA = THREE * PI / TWO 495 ENDIF 496 ELSE 497 ! pathological case: pole in alpha at beta=0 498 PATH = .TRUE. 499 ENDIF 500 501 COSBETA = RIJ(3)/MAGR 502 BETA = ACOS(RIJ(3) / MAGR) 503 504 DC = RIJ/MAGR 505 506 ! build forces using PRB 72 165107 eq. (12) - the sign of the 507 ! dfda contribution seems to be wrong, but gives the right 508 ! answer(?) 509 510 FTMP_PULAY = ZERO 511 FTMP_COUL = ZERO 512 FTMP_SPIN = ZERO 513 514 K = INDI 515 516 LBRAINC = 1 517 DO WHILE (BASISI(LBRAINC) .NE. -1) 518 519 LBRA = BASISI(LBRAINC) 520 LBRAINC = LBRAINC + 1 521 522 IF (SPINON .EQ. 1) THEN 523 524 SELECT CASE(LBRA) 525 CASE(0) 526 WSPINI = WSS(ELEMPOINTER(I)) 527 CASE(1) 528 WSPINI = WPP(ELEMPOINTER(I)) 529 CASE(2) 530 WSPINI = WDD(ELEMPOINTER(I)) 531 CASE(3) 532 WSPINI = WFF(ELEMPOINTER(I)) 533 END SELECT 534 535 WSPINI = WSPINI*DELTASPIN(SPININDI + LBRAINC - 1) 536 537 ENDIF 538 539 540 DO MBRA = -LBRA, LBRA 541 542 K = K + 1 543 L = INDJ 544 545 LKETINC = 1 546 DO WHILE (BASISJ(LKETINC) .NE. -1) 547 548 LKET = BASISJ(LKETINC) 549 LKETINC = LKETINC + 1 550 551 IF (SPINON .EQ. 1) THEN 552 553 SELECT CASE(LKET) 554 CASE(0) 555 WSPINJ = WSS(ELEMPOINTER(J)) 556 CASE(1) 557 WSPINJ = WPP(ELEMPOINTER(J)) 558 CASE(2) 559 WSPINJ = WDD(ELEMPOINTER(J)) 560 CASE(3) 561 WSPINJ = WFF(ELEMPOINTER(J)) 562 END SELECT 563 564 WSPINJ = WSPINJ*DELTASPIN(SPININDJ + LKETINC - 1) 565 566 ENDIF 567 568 DO MKET = -LKET, LKET 569 570 L = L + 1 571 572 SELECT CASE(SPINON) 573 CASE(0) 574 RHO = BO(L, K) 575 CASE(1) 576 RHO = RHOUP(L, K) + RHODOWN(L, K) 577 RHODIFF = (RHOUP(L, K) - RHODOWN(L, K))*(WSPINI + WSPINJ) 578 END SELECT 579 580! RHO = X2HRHO(L, K) 581 582 IF (.NOT. PATH) THEN 583 584 ! Unroll loops and pre-compute 585 586 MYDFDA = DFDA(I, J, LBRA, LKET, MBRA, & 587 MKET, MAGR, ALPHA, COSBETA, "S") 588 589 MYDFDB = DFDB(I, J, LBRA, LKET, MBRA, & 590 MKET, MAGR, ALPHA, COSBETA, "S") 591 592 MYDFDR = DFDR(I, J, LBRA, LKET, MBRA, & 593 MKET, MAGR, ALPHA, COSBETA, "S") 594 595 ! 596 ! d/d_alpha 597 ! 598 599 FTMP_PULAY(1) = FTMP_PULAY(1) + X2HRHO(L, K) * & 600 (-RIJ(2) / MAGRP2 * MYDFDA) 601 602 FTMP_PULAY(2) = FTMP_PULAY(2) + X2HRHO(L, K) * & 603 (RIJ(1)/ MAGRP2 * MYDFDA) 604 605 606 FTMP_COUL(1) = FTMP_COUL(1) + RHO * & 607 (-RIJ(2) / MAGRP2 * MYDFDA) 608 609 FTMP_COUL(2) = FTMP_COUL(2) + RHO * & 610 (RIJ(1)/ MAGRP2 * MYDFDA) 611 612 IF (SPINON .EQ. 1) THEN 613 614 FTMP_SPIN(1) = FTMP_SPIN(1) + RHODIFF * & 615 (-RIJ(2) / MAGRP2 * MYDFDA) 616 617 FTMP_SPIN(2) = FTMP_SPIN(2) + RHODIFF * & 618 (RIJ(1)/ MAGRP2 * MYDFDA) 619 620 ENDIF 621 622 623 ! 624 ! d/d_beta 625 ! 626 627 FTMP_PULAY(1) = FTMP_PULAY(1) + X2HRHO(L, K) * & 628 (((((RIJ(3) * RIJ(1)) / & 629 MAGR2)) / MAGRP) * MYDFDB) 630 631 632 FTMP_PULAY(2) = FTMP_PULAY(2) + X2HRHO(L, K) * & 633 (((((RIJ(3) * RIJ(2)) / & 634 MAGR2)) / MAGRP) * MYDFDB) 635 636 FTMP_PULAY(3) = FTMP_PULAY(3) - X2HRHO(L, K) * & 637 (((ONE - ((RIJ(3) * RIJ(3)) / & 638 MAGR2)) / MAGRP) * MYDFDB) 639 640 641 FTMP_COUL(1) = FTMP_COUL(1) + RHO * & 642 (((((RIJ(3) * RIJ(1)) / & 643 MAGR2)) / MAGRP) * MYDFDB) 644 645 FTMP_COUL(2) = FTMP_COUL(2) + RHO * & 646 (((((RIJ(3) * RIJ(2)) / & 647 MAGR2)) / MAGRP) * MYDFDB) 648 649 FTMP_COUL(3) = FTMP_COUL(3) - RHO * & 650 (((ONE - ((RIJ(3) * RIJ(3)) / & 651 MAGR2)) / MAGRP) * MYDFDB) 652 653 IF (SPINON .EQ. 1) THEN 654 655 FTMP_SPIN(1) = FTMP_SPIN(1) + RHODIFF * & 656 (((((RIJ(3) * RIJ(1)) / & 657 MAGR2)) / MAGRP) * MYDFDB) 658 659 FTMP_SPIN(2) = FTMP_SPIN(2) + RHODIFF * & 660 (((((RIJ(3) * RIJ(2)) / & 661 MAGR2)) / MAGRP) * MYDFDB) 662 663 FTMP_SPIN(3) = FTMP_SPIN(3) - RHODIFF * & 664 (((ONE - ((RIJ(3) * RIJ(3)) / & 665 MAGR2)) / MAGRP) * MYDFDB) 666 667 ENDIF 668 669 670 ! 671 ! d/dR 672 ! 673 674 FTMP_PULAY(1) = FTMP_PULAY(1) - X2HRHO(L, K) * DC(1) * & 675 MYDFDR 676 677 FTMP_PULAY(2) = FTMP_PULAY(2) - X2HRHO(L, K) * DC(2) * & 678 MYDFDR 679 680 FTMP_PULAY(3) = FTMP_PULAY(3) - X2HRHO(L, K) * DC(3) * & 681 MYDFDR 682 683 684 FTMP_COUL(1) = FTMP_COUL(1) - RHO * DC(1) * & 685 MYDFDR 686 687 FTMP_COUL(2) = FTMP_COUL(2) - RHO * DC(2) * & 688 MYDFDR 689 690 FTMP_COUL(3) = FTMP_COUL(3) - RHO * DC(3) * & 691 MYDFDR 692 693 IF (SPINON .EQ. 1) THEN 694 695 FTMP_SPIN(1) = FTMP_SPIN(1) - RHODIFF * DC(1) * & 696 MYDFDR 697 698 FTMP_SPIN(2) = FTMP_SPIN(2) - RHODIFF * DC(2) * & 699 MYDFDR 700 701 FTMP_SPIN(3) = FTMP_SPIN(3) - RHODIFF * DC(3) * & 702 MYDFDR 703 704 ENDIF 705 706 707 ELSE 708 709 ! pathological configuration in which beta=0 710 ! or pi => alpha undefined 711 712 ! fixed: MJC 12/17/13 713 714 MYDFDB = DFDB(I, J, LBRA, LKET, & 715 MBRA, MKET, MAGR, ZERO, COSBETA, "S") / MAGR 716 717 FTMP_PULAY(1) = FTMP_PULAY(1) - X2HRHO(L, K) * (COSBETA * MYDFDB) 718 719 FTMP_COUL(1) = FTMP_COUL(1) - RHO * (COSBETA * MYDFDB) 720 721 IF (SPINON .EQ. 1) THEN 722 FTMP_SPIN(1) = FTMP_SPIN(1) - RHODIFF * (COSBETA * MYDFDB) 723 ENDIF 724 725 MYDFDB = DFDB(I, J, LBRA, LKET, & 726 MBRA, MKET, MAGR, PI/TWO, COSBETA, "S") / MAGR 727 728 FTMP_PULAY(2) = FTMP_PULAY(2) - X2HRHO(L, K) * (COSBETA * MYDFDB) 729 730 FTMP_COUL(2) = FTMP_COUL(2) - RHO * (COSBETA * MYDFDB) 731 732 IF (SPINON .EQ. 1) THEN 733 FTMP_SPIN(2) = FTMP_SPIN(2) - RHODIFF * (COSBETA * MYDFDB) 734 ENDIF 735 736 MYDFDR = DFDR(I, J, LBRA, LKET, MBRA, & 737 MKET, MAGR, ZERO, COSBETA, "S") 738 739 FTMP_PULAY(3) = FTMP_PULAY(3) - X2HRHO(L, K) * COSBETA * MYDFDR 740 741 FTMP_COUL(3) = FTMP_COUL(3) - RHO * COSBETA * MYDFDR 742 743 IF (SPINON .EQ. 1) THEN 744 FTMP_SPIN(3) = FTMP_SPIN(3) - RHODIFF * COSBETA * MYDFDR 745 ENDIF 746 747 ENDIF 748 749 ENDDO 750 ENDDO 751 ENDDO 752 ENDDO 753 754 FPUL(1,I) = FPUL(1,I) - TWO * FTMP_PULAY(1) 755 FPUL(2,I) = FPUL(2,I) - TWO * FTMP_PULAY(2) 756 FPUL(3,I) = FPUL(3,I) - TWO * FTMP_PULAY(3) 757 758 VIRPUL(1) = VIRPUL(1) - RIJ(1) * FTMP_PULAY(1) 759 VIRPUL(2) = VIRPUL(2) - RIJ(2) * FTMP_PULAY(2) 760 VIRPUL(3) = VIRPUL(3) - RIJ(3) * FTMP_PULAY(3) 761 VIRPUL(4) = VIRPUL(4) - RIJ(1) * FTMP_PULAY(2) 762 VIRPUL(5) = VIRPUL(5) - RIJ(2) * FTMP_PULAY(3) 763 VIRPUL(6) = VIRPUL(6) - RIJ(3) * FTMP_PULAY(1) 764 765 IF (ELECTRO .EQ. 1) THEN 766 767 FTMP_COUL = FTMP_COUL * ( HUBBARDU(ELEMPOINTER(J))*DELTAQ(J) + COULOMBV(J) & 768 +HUBBARDU(ELEMPOINTER(I))*DELTAQ(I) + COULOMBV(I)) 769 770 FPUL(1,I) = FPUL(1,I) + FTMP_COUL(1) 771 FPUL(2,I) = FPUL(2,I) + FTMP_COUL(2) 772 FPUL(3,I) = FPUL(3,I) + FTMP_COUL(3) 773 774 ! with the factor of 2... \ 775 776 777 VIRPUL(1) = VIRPUL(1) + RIJ(1)*FTMP_COUL(1)/TWO 778 VIRPUL(2) = VIRPUL(2) + RIJ(2)*FTMP_COUL(2)/TWO 779 VIRPUL(3) = VIRPUL(3) + RIJ(3)*FTMP_COUL(3)/TWO 780 VIRPUL(4) = VIRPUL(4) + RIJ(1)*FTMP_COUL(2)/TWO 781 VIRPUL(5) = VIRPUL(5) + RIJ(2)*FTMP_COUL(3)/TWO 782 VIRPUL(6) = VIRPUL(6) + RIJ(3)*FTMP_COUL(1)/TWO 783 784 ELSE 785 786 FTMP_COUL = FTMP_COUL * (LCNSHIFT(I) + LCNSHIFT(J)) 787 788 FPUL(1,I) = FPUL(1,I) + FTMP_COUL(1) 789 FPUL(2,I) = FPUL(2,I) + FTMP_COUL(2) 790 FPUL(3,I) = FPUL(3,I) + FTMP_COUL(3) 791 792 ! with the factor of 2... \ 793 794 795 VIRPUL(1) = VIRPUL(1) + RIJ(1)*FTMP_COUL(1)/TWO 796 VIRPUL(2) = VIRPUL(2) + RIJ(2)*FTMP_COUL(2)/TWO 797 VIRPUL(3) = VIRPUL(3) + RIJ(3)*FTMP_COUL(3)/TWO 798 VIRPUL(4) = VIRPUL(4) + RIJ(1)*FTMP_COUL(2)/TWO 799 VIRPUL(5) = VIRPUL(5) + RIJ(2)*FTMP_COUL(3)/TWO 800 VIRPUL(6) = VIRPUL(6) + RIJ(3)*FTMP_COUL(1)/TWO 801 802 ENDIF 803 804 IF (SPINON .EQ. 1) THEN 805 806 FPUL(1,I) = FPUL(1,I) + FTMP_SPIN(1) 807 FPUL(2,I) = FPUL(2,I) + FTMP_SPIN(2) 808 FPUL(3,I) = FPUL(3,I) + FTMP_SPIN(3) 809 810 ! with the factor of 2... \ 811 812 VIRPUL(1) = VIRPUL(1) + RIJ(1)*FTMP_SPIN(1)/TWO 813 VIRPUL(2) = VIRPUL(2) + RIJ(2)*FTMP_SPIN(2)/TWO 814 VIRPUL(3) = VIRPUL(3) + RIJ(3)*FTMP_SPIN(3)/TWO 815 VIRPUL(4) = VIRPUL(4) + RIJ(1)*FTMP_SPIN(2)/TWO 816 VIRPUL(5) = VIRPUL(5) + RIJ(2)*FTMP_SPIN(3)/TWO 817 VIRPUL(6) = VIRPUL(6) + RIJ(3)*FTMP_SPIN(1)/TWO 818 819 820 ENDIF 821 822 ENDIF 823 824 ENDDO 825 826 ENDDO 827 828!$OMP END PARALLEL DO 829 830 ! DO I= 1, NATS 831 ! WRITE(6,10) I, FPUL(1,I), FPUL(2,I), FPUL(3,I) 832 ! ENDDO 833 834 !10 FORMAT(I4, 3F12.6) 835 836 RETURN 837 838END SUBROUTINE PULAY 839