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