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