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 KPULAY 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 KSPACEARRAY 32 USE MYPRECISION 33 34 IMPLICIT NONE 35 36 INTEGER :: I, J, K, L, M, N, KK, INDI, INDJ 37 INTEGER :: LBRA, MBRA, LKET, MKET 38 INTEGER :: PREVJ, NEWJ 39 INTEGER :: PBCI, PBCJ, PBCK 40 INTEGER :: BASISI(5), BASISJ(5), LBRAINC, LKETINC 41 INTEGER :: KX, KY, KZ, KCOUNT 42 REAL(LATTEPREC) :: ALPHA, BETA, PHI, COSBETA 43 REAL(LATTEPREC) :: RIJ(3), DC(3) 44 REAL(LATTEPREC) :: MAGR, MAGR2, MAGRP, MAGRP2 45 REAL(LATTEPREC) :: MAXRCUT, MAXRCUT2 46 REAL(LATTEPREC) :: MYDFDA, MYDFDB, MYDFDR, RCUTTB 47 REAL(LATTEPREC) :: KPOINT(3), KX0, KY0, KZ0, KDOTL 48 REAL(LATTEPREC) :: B1(3), B2(3), B3(3), MAG1, MAG2, MAG3, A1A2XA3, K0(3) 49 REAL(LATTEPREC), EXTERNAL :: DFDA, DFDB, DFDR 50 COMPLEX(LATTEPREC) :: FTMP_PULAY(3), FTMP_COUL(3) 51 COMPLEX(LATTEPREC) :: RHO_PULAY, RHO_COUL, CONJGBLOCH 52 COMPLEX(LATTEPREC), ALLOCATABLE :: KX2HRHO(:,:,:), KTMP(:,:) 53 COMPLEX(LATTEPREC), PARAMETER :: ZONE=CMPLX(ONE), ZZERO=CMPLX(ZERO), ZHALF=CMPLX(HALF) 54 55 LOGICAL PATH 56 IF (EXISTERROR) RETURN 57 58 ! These were allocated elsewhere. We'll use them to accumulate the complex forces 59 60 IF (SPINON .EQ. 1) THEN 61 CALL ERRORS("kpulay","Non-ortho k-space and spin polarization not yet implemented") 62 ENDIF 63 64 KF = CMPLX(ZERO) 65 VIRBONDK = CMPLX(ZERO) 66 67 ALLOCATE(KX2HRHO(HDIM, HDIM, NKTOT), KTMP(HDIM, HDIM)) 68 69 ! Computing the reciprocal lattice vectors 70 71 B1(1) = BOX(2,2)*BOX(3,3) - BOX(3,2)*BOX(2,3) 72 B1(2) = BOX(3,1)*BOX(2,3) - BOX(2,1)*BOX(3,3) 73 B1(3) = BOX(2,1)*BOX(3,2) - BOX(3,1)*BOX(2,2) 74 75 A1A2XA3 = BOX(1,1)*B1(1) + BOX(1,2)*B1(2) + BOX(1,3)*B1(3) 76 77 ! B1 = 2*PI*(A2 X A3)/(A1.(A2 X A3)) 78 79 B1 = B1/A1A2XA3 80 81 ! B2 = 2*PI*(A3 x A1)/(A1(A2 X A3)) 82 83 B2(1) = (BOX(3,2)*BOX(1,3) - BOX(1,2)*BOX(3,3))/A1A2XA3 84 B2(2) = (BOX(1,1)*BOX(3,3) - BOX(3,1)*BOX(1,3))/A1A2XA3 85 B2(3) = (BOX(3,1)*BOX(1,2) - BOX(1,1)*BOX(3,2))/A1A2XA3 86 87 ! B3 = 2*PI*(A1 x A2)/(A1(A2 X A3)) 88 89 B3(1) = (BOX(1,2)*BOX(2,3) - BOX(2,2)*BOX(1,3))/A1A2XA3 90 B3(2) = (BOX(2,1)*BOX(1,3) - BOX(1,1)*BOX(2,3))/A1A2XA3 91 B3(3) = (BOX(1,1)*BOX(2,2) - BOX(2,1)*BOX(1,2))/A1A2XA3 92 93 K0 = PI*(ONE - REAL(NKX))/(REAL(NKX))*B1 + & 94 PI*(ONE - REAL(NKY))/(REAL(NKY))*B2 + & 95 PI*(ONE - REAL(NKZ))/(REAL(NKZ))*B3 - PI*KSHIFT 96 97 ! We first have to make the matrix S^-1 H rho = X^2 H rho 98 99 IF (KBT .GT. 0.000001) THEN 100 101 ! Finite temperature 102 103 DO K = 1, NKTOT 104 105 CALL ZGEMM('N', 'N', HDIM, HDIM, HDIM, ZONE, & 106 KXMAT(:,:,K), HDIM, KXMAT(:,:,K), HDIM, ZZERO, KX2HRHO(:,:,K), HDIM) 107 108 CALL ZGEMM('N', 'N', HDIM, HDIM, HDIM, ZONE, & 109 KX2HRHO(:,:,K), HDIM, HK(:,:,K), HDIM, ZZERO, KTMP, HDIM) 110 111 ! (S^-1 * H)*RHO 112 113 CALL ZGEMM('N', 'N', HDIM, HDIM, HDIM, ZONE, & 114 KTMP, HDIM, KBO(:,:,K), HDIM, ZZERO, KX2HRHO(:,:,K), HDIM) 115 116 ENDDO 117 118 ELSE 119 120 ! Te = 0 : Fp = 2Tr[rho H rho dS/dR] 121 122 ! Be careful - we're working with bo = 2rho, so we need 123 ! the factor of 1/2... 124 125 DO K = 1, NKTOT 126 127 CALL ZGEMM('N', 'N', HDIM, HDIM, HDIM, ZONE, & 128 KBO(:,:,K), HDIM, HK(:,:,K), HDIM, ZZERO, KTMP, HDIM) 129 130 CALL ZGEMM('N', 'N', HDIM, HDIM, HDIM, ZHALF, & 131 KTMP, HDIM, KBO(:,:,K), HDIM, ZZERO, KX2HRHO(:,:,K), HDIM) 132 133 ENDDO 134 135 ENDIF 136 137!$OMP PARALLEL DO DEFAULT (NONE) & 138!$OMP SHARED(NATS, BASIS, ELEMPOINTER, TOTNEBTB, NEBTB) & 139!$OMP SHARED(CR, BOX, KX2HRHO, KBO, NOINT, ATELE, ELE1, ELE2) & 140!$OMP SHARED(HCUT, SCUT, MATINDLIST, BASISTYPE, ELECTRO) & 141!$OMP SHARED(K0, B1, B2, B3, NKX, NKY, NKZ, KF) & 142!$OMP SHARED(LCNSHIFT, HUBBARDU, DELTAQ, COULOMBV) & 143!$OMP PRIVATE(I, J, K, NEWJ, BASISI, BASISJ, INDI, INDJ, PBCI, PBCJ, PBCK) & 144!$OMP PRIVATE(RIJ, MAGR2, MAGR, MAGRP2, MAGRP, PATH, PHI, ALPHA, BETA, COSBETA, FTMP_PULAY, FTMP_COUL) & 145!$OMP PRIVATE(DC, LBRAINC, LBRA, MBRA, L, LKETINC, LKET, MKET, RHO_PULAY, RHO_COUL) & 146!$OMP PRIVATE(MYDFDA, MYDFDB, MYDFDR, RCUTTB, CONJGBLOCH, KDOTL) & 147!$OMP PRIVATE(KPOINT, KCOUNT) & 148!$OMP REDUCTION(+: VIRBONDK) 149 150 151 DO I = 1, NATS 152 153 ! Build list of orbitals on atom I 154 SELECT CASE(BASIS(ELEMPOINTER(I))) 155 156 CASE("s") 157 BASISI(1) = 0 158 BASISI(2) = -1 159 CASE("p") 160 BASISI(1) = 1 161 BASISI(2) = -1 162 CASE("d") 163 BASISI(1) = 2 164 BASISI(2) = -1 165 CASE("f") 166 BASISI(1) = 3 167 BASISI(2) = -1 168 CASE("sp") 169 BASISI(1) = 0 170 BASISI(2) = 1 171 BASISI(3) = -1 172 CASE("sd") 173 BASISI(1) = 0 174 BASISI(2) = 2 175 BASISI(3) = -1 176 CASE("sf") 177 BASISI(1) = 0 178 BASISI(2) = 3 179 BASISI(3) = -1 180 CASE("pd") 181 BASISI(1) = 1 182 BASISI(2) = 2 183 BASISI(3) = -1 184 CASE("pf") 185 BASISI(1) = 1 186 BASISI(2) = 3 187 BASISI(3) = -1 188 CASE("df") 189 BASISI(1) = 2 190 BASISI(2) = 3 191 BASISI(3) = -1 192 CASE("spd") 193 BASISI(1) = 0 194 BASISI(2) = 1 195 BASISI(3) = 2 196 BASISI(4) = -1 197 CASE("spf") 198 BASISI(1) = 0 199 BASISI(2) = 1 200 BASISI(3) = 3 201 BASISI(4) = -1 202 CASE("sdf") 203 BASISI(1) = 0 204 BASISI(2) = 2 205 BASISI(3) = 3 206 BASISI(4) = -1 207 CASE("pdf") 208 BASISI(1) = 1 209 BASISI(2) = 2 210 BASISI(3) = 3 211 BASISI(4) = -1 212 CASE("spdf") 213 BASISI(1) = 0 214 BASISI(2) = 1 215 BASISI(3) = 2 216 BASISI(4) = 3 217 BASISI(5) = -1 218 END SELECT 219 220 INDI = MATINDLIST(I) 221 222 DO NEWJ = 1, TOTNEBTB(I) 223 224 J = NEBTB(1, NEWJ, I) 225 PBCI = NEBTB(2, NEWJ, I) 226 PBCJ = NEBTB(3, NEWJ, I) 227 PBCK = NEBTB(4, NEWJ, I) 228 229 RIJ(1) = CR(1,J) + REAL(PBCI)*BOX(1,1) + REAL(PBCJ)*BOX(2,1) + & 230 REAL(PBCK)*BOX(3,1) - CR(1,I) 231 232 RIJ(2) = CR(2,J) + REAL(PBCI)*BOX(1,2) + REAL(PBCJ)*BOX(2,2) + & 233 REAL(PBCK)*BOX(3,2) - CR(2,I) 234 235 RIJ(3) = CR(3,J) + REAL(PBCI)*BOX(1,3) + REAL(PBCJ)*BOX(2,3) + & 236 REAL(PBCK)*BOX(3,3) - CR(3,I) 237 238 MAGR2 = RIJ(1)*RIJ(1) + RIJ(2)*RIJ(2) + RIJ(3)*RIJ(3) 239 240 ! Building the forces is expensive - use the cut-off 241 242 RCUTTB = ZERO 243 DO K = 1, NOINT 244 245 IF ( (ATELE(I) .EQ. ELE1(K) .AND. & 246 ATELE(J) .EQ. ELE2(K)) .OR. & 247 (ATELE(J) .EQ. ELE1(K) .AND. & 248 ATELE(I) .EQ. ELE2(K) )) THEN 249 250 IF (SCUT(K) .GT. RCUTTB ) RCUTTB = SCUT(K) 251 252 ENDIF 253 254 ENDDO 255 256 IF (MAGR2 .LT. RCUTTB*RCUTTB) THEN 257 258 MAGR = SQRT(MAGR2) 259 260 ! Build list of orbitals on atom J 261 262 SELECT CASE(BASIS(ELEMPOINTER(J))) 263 CASE("s") 264 BASISJ(1) = 0 265 BASISJ(2) = -1 266 CASE("p") 267 BASISJ(1) = 1 268 BASISJ(2) = -1 269 CASE("d") 270 BASISJ(1) = 2 271 BASISJ(2) = -1 272 CASE("f") 273 BASISJ(1) = 3 274 BASISJ(2) = -1 275 CASE("sp") 276 BASISJ(1) = 0 277 BASISJ(2) = 1 278 BASISJ(3) = -1 279 CASE("sd") 280 BASISJ(1) = 0 281 BASISJ(2) = 2 282 BASISJ(3) = -1 283 CASE("sf") 284 BASISJ(1) = 0 285 BASISJ(2) = 3 286 BASISJ(3) = -1 287 CASE("pd") 288 BASISJ(1) = 1 289 BASISJ(2) = 2 290 BASISJ(3) = -1 291 CASE("pf") 292 BASISJ(1) = 1 293 BASISJ(2) = 3 294 BASISJ(3) = -1 295 CASE("df") 296 BASISJ(1) = 2 297 BASISJ(2) = 3 298 BASISJ(3) = -1 299 CASE("spd") 300 BASISJ(1) = 0 301 BASISJ(2) = 1 302 BASISJ(3) = 2 303 BASISJ(4) = -1 304 CASE("spf") 305 BASISJ(1) = 0 306 BASISJ(2) = 1 307 BASISJ(3) = 3 308 BASISJ(4) = -1 309 CASE("sdf") 310 BASISJ(1) = 0 311 BASISJ(2) = 2 312 BASISJ(3) = 3 313 BASISJ(4) = -1 314 CASE("pdf") 315 BASISJ(1) = 1 316 BASISJ(2) = 2 317 BASISJ(3) = 3 318 BASISJ(4) = -1 319 CASE("spdf") 320 BASISJ(1) = 0 321 BASISJ(2) = 1 322 BASISJ(3) = 2 323 BASISJ(4) = 3 324 BASISJ(5) = -1 325 END SELECT 326 327 INDJ = MATINDLIST(J) 328 329 MAGRP2 = RIJ(1)*RIJ(1) + RIJ(2)*RIJ(2) 330 MAGRP = SQRT(MAGRP2) 331 332 333 ! transform to system in which z-axis is aligned with RIJ 334 335 PATH = .FALSE. 336 IF (ABS(RIJ(1)) .GT. 1E-12) THEN 337 IF (RIJ(1) .GT. ZERO .AND. RIJ(2) .GE. ZERO) THEN 338 PHI = ZERO 339 ELSEIF (RIJ(1) .GT. ZERO .AND. RIJ(2) .LT. ZERO) THEN 340 PHI = TWO * PI 341 ELSE 342 PHI = PI 343 ENDIF 344 ALPHA = ATAN(RIJ(2) / RIJ(1)) + PHI 345 ELSEIF (ABS(RIJ(2)) .GT. 1E-12) THEN 346 IF (RIJ(2) .GT. 1E-12) THEN 347 ALPHA = PI / TWO 348 ELSE 349 ALPHA = THREE * PI / TWO 350 ENDIF 351 ELSE 352 ! pathological case: pole in alpha at beta=0 353 PATH = .TRUE. 354 ENDIF 355 356 COSBETA = RIJ(3)/MAGR 357 BETA = ACOS(RIJ(3) / MAGR) 358 359 DC = RIJ/MAGR 360 361 ! build forces using PRB 72 165107 eq. (12) - the sign of the 362 ! dfda contribution seems to be wrong, but gives the right 363 ! answer(?) 364 365 FTMP_PULAY = ZERO 366 FTMP_COUL = ZERO 367 368 K = INDI 369 370 LBRAINC = 1 371 DO WHILE (BASISI(LBRAINC) .NE. -1) 372 373 LBRA = BASISI(LBRAINC) 374 LBRAINC = LBRAINC + 1 375 376 DO MBRA = -LBRA, LBRA 377 378 K = K + 1 379 L = INDJ 380 381 LKETINC = 1 382 DO WHILE (BASISJ(LKETINC) .NE. -1) 383 384 LKET = BASISJ(LKETINC) 385 LKETINC = LKETINC + 1 386 387 DO MKET = -LKET, LKET 388 389 L = L + 1 390 391 !RHO = X2HRHO(L, K) 392 393 IF (.NOT. PATH) THEN 394 395 ! Unroll loops and pre-compute 396 397 MYDFDA = DFDA(I, J, LBRA, LKET, MBRA, & 398 MKET, MAGR, ALPHA, COSBETA, "S") 399 400 MYDFDB = DFDB(I, J, LBRA, LKET, MBRA, & 401 MKET, MAGR, ALPHA, COSBETA, "S") 402 403 MYDFDR = DFDR(I, J, LBRA, LKET, MBRA, & 404 MKET, MAGR, ALPHA, COSBETA, "S") 405 406 407 KCOUNT = 0 408 409 DO KX = 1, NKX 410 411 DO KY = 1, NKY 412 413 DO KZ = 1, NKZ 414 415 KPOINT = TWO*PI*(REAL(KX-1)*B1/REAL(NKX) + & 416 REAL(KY-1)*B2/REAL(NKY) + & 417 REAL(KZ-1)*B3/REAL(NKZ)) + K0 418 419 KCOUNT = KCOUNT+1 420 421 KDOTL = KPOINT(1)*RIJ(1) + KPOINT(2)*RIJ(2) + & 422 KPOINT(3)*RIJ(3) 423 424 CONJGBLOCH = EXP(CMPLX(ZERO,-KDOTL)) 425 426 RHO_PULAY = KX2HRHO(K,L,KCOUNT)*CONJGBLOCH 427 428 RHO_COUL = KBO(K,L,KCOUNT)*CONJGBLOCH 429 430 ! 431 ! d/d_alpha 432 ! 433 434 FTMP_PULAY(1) = FTMP_PULAY(1) + RHO_PULAY * & 435 (-RIJ(2) / MAGRP2 * MYDFDA) 436 437 FTMP_PULAY(2) = FTMP_PULAY(2) + RHO_PULAY * & 438 (RIJ(1)/ MAGRP2 * MYDFDA) 439 440 FTMP_COUL(1) = FTMP_COUL(1) + RHO_COUL * & 441 (-RIJ(2) / MAGRP2 * MYDFDA) 442 443 FTMP_COUL(2) = FTMP_COUL(2) + RHO_COUL * & 444 (RIJ(1)/ MAGRP2 * MYDFDA) 445 446 447 ! 448 ! d/d_beta 449 ! 450 451 FTMP_PULAY(1) = FTMP_PULAY(1) + RHO_PULAY * & 452 (((((RIJ(3) * RIJ(1)) / & 453 MAGR2)) / MAGRP) * MYDFDB) 454 455 FTMP_PULAY(2) = FTMP_PULAY(2) + RHO_PULAY * & 456 (((((RIJ(3) * RIJ(2)) / & 457 MAGR2)) / MAGRP) * MYDFDB) 458 459 FTMP_PULAY(3) = FTMP_PULAY(3) - RHO_PULAY * & 460 (((ONE - ((RIJ(3) * RIJ(3)) / & 461 MAGR2)) / MAGRP) * MYDFDB) 462 463 FTMP_COUL(1) = FTMP_COUL(1) + RHO_COUL * & 464 (((((RIJ(3) * RIJ(1)) / & 465 MAGR2)) / MAGRP) * MYDFDB) 466 467 FTMP_COUL(2) = FTMP_COUL(2) + RHO_COUL * & 468 (((((RIJ(3) * RIJ(2)) / & 469 MAGR2)) / MAGRP) * MYDFDB) 470 471 FTMP_COUL(3) = FTMP_COUL(3) - RHO_COUL * & 472 (((ONE - ((RIJ(3) * RIJ(3)) / & 473 MAGR2)) / MAGRP) * MYDFDB) 474 475 476 477 ! 478 ! d/dR 479 ! 480 481 FTMP_PULAY(1) = FTMP_PULAY(1) - RHO_PULAY * DC(1) * & 482 MYDFDR 483 484 FTMP_PULAY(2) = FTMP_PULAY(2) - RHO_PULAY * DC(2) * & 485 MYDFDR 486 487 FTMP_PULAY(3) = FTMP_PULAY(3) - RHO_PULAY * DC(3) * & 488 MYDFDR 489 490 FTMP_COUL(1) = FTMP_COUL(1) - RHO_COUL * DC(1) * & 491 MYDFDR 492 493 FTMP_COUL(2) = FTMP_COUL(2) - RHO_COUL * DC(2) * & 494 MYDFDR 495 496 FTMP_COUL(3) = FTMP_COUL(3) - RHO_COUL * DC(3) * & 497 MYDFDR 498 499 500 ENDDO 501 ENDDO 502 ENDDO 503 504 ELSE 505 506 ! pathological configuration in which beta=0 507 ! or pi => alpha undefined 508 509 ! fixed: MJC 12/17/13 510 511 MYDFDB = DFDB(I, J, LBRA, LKET, & 512 MBRA, MKET, MAGR, ZERO, COSBETA, "S") / MAGR 513 514 515 MYDFDB = DFDB(I, J, LBRA, LKET, & 516 MBRA, MKET, MAGR, PI/TWO, COSBETA, "S") / MAGR 517 518 MYDFDR = DFDR(I, J, LBRA, LKET, MBRA, & 519 MKET, MAGR, ZERO, COSBETA, "S") 520 521 KCOUNT = 0 522 523 DO KX = 1, NKX 524 525 DO KY = 1, NKY 526 527 DO KZ = 1, NKZ 528 529 KPOINT = TWO*PI*(REAL(KX-1)*B1/REAL(NKX) + & 530 REAL(KY-1)*B2/REAL(NKY) + & 531 REAL(KZ-1)*B3/REAL(NKZ)) + K0 532 533 KCOUNT = KCOUNT+1 534 535 KDOTL = KPOINT(1)*RIJ(1) + KPOINT(2)*RIJ(2) + & 536 KPOINT(3)*RIJ(3) 537 538 CONJGBLOCH = EXP(CMPLX(ZERO,-KDOTL)) 539 540 RHO_PULAY = KX2HRHO(K,L,KCOUNT)*CONJGBLOCH 541 542 RHO_COUL = KBO(K,L,KCOUNT)*CONJGBLOCH 543 544 FTMP_PULAY(1) = FTMP_PULAY(1) - RHO_PULAY * COSBETA * MYDFDB 545 FTMP_PULAY(2) = FTMP_PULAY(2) - RHO_PULAY * COSBETA * MYDFDB 546 FTMP_PULAY(3) = FTMP_PULAY(3) - RHO_PULAY * COSBETA * MYDFDR 547 548 FTMP_COUL(1) = FTMP_COUL(1) - RHO_COUL * COSBETA * MYDFDB 549 FTMP_COUL(2) = FTMP_COUL(2) - RHO_COUL * COSBETA * MYDFDB 550 FTMP_COUL(3) = FTMP_COUL(3) - RHO_COUL * COSBETA * MYDFDR 551 552 ENDDO 553 ENDDO 554 ENDDO 555 556 ENDIF 557 558 ENDDO 559 ENDDO 560 ENDDO 561 ENDDO 562 563 KF(1,I) = KF(1,I) - TWO * FTMP_PULAY(1) 564 KF(2,I) = KF(2,I) - TWO * FTMP_PULAY(2) 565 KF(3,I) = KF(3,I) - TWO * FTMP_PULAY(3) 566 567 VIRBONDK(1) = VIRBONDK(1) - RIJ(1) * FTMP_PULAY(1) 568 VIRBONDK(2) = VIRBONDK(2) - RIJ(2) * FTMP_PULAY(2) 569 VIRBONDK(3) = VIRBONDK(3) - RIJ(3) * FTMP_PULAY(3) 570 VIRBONDK(4) = VIRBONDK(4) - RIJ(1) * FTMP_PULAY(2) 571 VIRBONDK(5) = VIRBONDK(5) - RIJ(2) * FTMP_PULAY(3) 572 VIRBONDK(6) = VIRBONDK(6) - RIJ(3) * FTMP_PULAY(1) 573 574 575 IF (ELECTRO .EQ. 1) THEN 576 577 FTMP_COUL = FTMP_COUL * ( HUBBARDU(ELEMPOINTER(J))*DELTAQ(J) + COULOMBV(J) & 578 +HUBBARDU(ELEMPOINTER(I))*DELTAQ(I) + COULOMBV(I)) 579 580 ELSE 581 582 FTMP_COUL = FTMP_COUL * (LCNSHIFT(I) + LCNSHIFT(J)) 583 584 ENDIF 585 586 587 KF(1,I) = KF(1,I) + FTMP_COUL(1) 588 KF(2,I) = KF(2,I) + FTMP_COUL(2) 589 KF(3,I) = KF(3,I) + FTMP_COUL(3) 590 591 ! with the factor of 2... 592 593 VIRBONDK(1) = VIRBONDK(1) + RIJ(1)*FTMP_COUL(1)/TWO 594 VIRBONDK(2) = VIRBONDK(2) + RIJ(2)*FTMP_COUL(2)/TWO 595 VIRBONDK(3) = VIRBONDK(3) + RIJ(3)*FTMP_COUL(3)/TWO 596 VIRBONDK(4) = VIRBONDK(4) + RIJ(1)*FTMP_COUL(2)/TWO 597 VIRBONDK(5) = VIRBONDK(5) + RIJ(2)*FTMP_COUL(3)/TWO 598 VIRBONDK(6) = VIRBONDK(6) + RIJ(3)*FTMP_COUL(1)/TWO 599 600 601 ENDIF 602 603 ENDDO 604 605 ENDDO 606 607!$OMP END PARALLEL DO 608 609 DEALLOCATE(KX2HRHO, KTMP) 610 611 ! PRINT*, KF(1,1), KF(2,1), KF(3,1) 612 613 ! DO I= 1, NATS 614 ! WRITE(6,10) I, FPUL(1,I), FPUL(2,I), FPUL(3,I) 615 ! ENDDO 616 617 !10 FORMAT(I4, 3F12.6) 618 619 FPUL = REAL(KF)/REAL(NKTOT) 620 VIRPUL = REAL(VIRBONDK)/REAL(NKTOT) 621 622 ! print*, VIRPUL(1) 623 RETURN 624 625END SUBROUTINE KPULAY 626