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 KBLDNEWH 23 24 USE CONSTANTS_MOD 25 USE SETUPARRAY 26 USE NEBLISTARRAY 27 USE XBOARRAY 28 USE NONOARRAY 29 USE UNIVARRAY 30 USE KSPACEARRAY 31 USE MYPRECISION 32 33 IMPLICIT NONE 34 35 INTEGER :: I, J, NEWJ, K, L, II, JJ, KK, MM, MP, NN, SUBI 36 INTEGER :: IBRA, IKET, LBRA, LKET, MBRA, MKET 37 INTEGER :: INDEX, INDI, INDJ 38 INTEGER :: SWITCH, PREVJ 39 INTEGER :: PBCI, PBCJ, PBCK 40 INTEGER :: BASISI(5), BASISJ(5), LBRAINC, LKETINC 41 INTEGER :: NORBI, KCOUNT 42 INTEGER :: KX, KY, KZ 43 REAL(LATTEPREC) :: KX0, KY0, KZ0 44 REAL(LATTEPREC) :: ALPHA, BETA, COSBETA, PHI, TMP, PERM 45 REAL(LATTEPREC) :: RIJ(3), MAGR2, MAGR, MAGRP, RCUTTB 46 REAL(LATTEPREC) :: MAXARRAY(20), MAXRCUT, MAXRCUT2 47 REAL(LATTEPREC) :: ANGFACTOR 48 REAL(LATTEPREC) :: KPOINT(3), KDOTL 49 REAL(LATTEPREC) :: AMMBRA, WIGLBRAMBRA 50 REAL(LATTEPREC) :: B1(3), B2(3), B3(3), A1A2XA3, K0(3) 51 COMPLEX(LATTEPREC) :: BLOCH, KHTMP, KSTMP 52 REAL(LATTEPREC), EXTERNAL :: UNIVSCALE, WIGNERD, SLMMP, TLMMP, AM, BM 53 IF (EXISTERROR) RETURN 54 55 HK = CMPLX(ZERO) 56 57 ! Computing the reciprocal lattice vectors 58 59 B1(1) = BOX(2,2)*BOX(3,3) - BOX(3,2)*BOX(2,3) 60 B1(2) = BOX(3,1)*BOX(2,3) - BOX(2,1)*BOX(3,3) 61 B1(3) = BOX(2,1)*BOX(3,2) - BOX(3,1)*BOX(2,2) 62 63 A1A2XA3 = BOX(1,1)*B1(1) + BOX(1,2)*B1(2) + BOX(1,3)*B1(3) 64 65 ! B1 = (A2 X A3)/(A1.(A2 X A3)) 66 67 B1 = B1/A1A2XA3 68 69 ! B2 = (A3 x A1)/(A1(A2 X A3)) 70 71 B2(1) = (BOX(3,2)*BOX(1,3) - BOX(1,2)*BOX(3,3))/A1A2XA3 72 B2(2) = (BOX(1,1)*BOX(3,3) - BOX(3,1)*BOX(1,3))/A1A2XA3 73 B2(3) = (BOX(3,1)*BOX(1,2) - BOX(1,1)*BOX(3,2))/A1A2XA3 74 75 ! B3 = (A1 x A2)/(A1(A2 X A3)) 76 77 B3(1) = (BOX(1,2)*BOX(2,3) - BOX(2,2)*BOX(1,3))/A1A2XA3 78 B3(2) = (BOX(2,1)*BOX(1,3) - BOX(1,1)*BOX(2,3))/A1A2XA3 79 B3(3) = (BOX(1,1)*BOX(2,2) - BOX(2,1)*BOX(1,2))/A1A2XA3 80 81 INDEX = 0 82 83 ! Build diagonal elements 84 DO I = 1, NATS 85 86 SELECT CASE(BASIS(ELEMPOINTER(I))) 87 88 CASE("s") 89 90 INDEX = INDEX + 1 91 HK(INDEX, INDEX, 1) = CMPLX(HES(ELEMPOINTER(I))) 92 93 CASE("p") 94 95 DO SUBI = 1, 3 96 INDEX = INDEX + 1 97 HK(INDEX,INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I))) 98 ENDDO 99 100 CASE("d") 101 102 DO SUBI = 1, 5 103 INDEX = INDEX + 1 104 HK(INDEX,INDEX, 1) = CMPLX(HED(ELEMPOINTER(I))) 105 ENDDO 106 107 CASE("f") 108 109 DO SUBI = 1, 7 110 INDEX = INDEX + 1 111 HK(INDEX,INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I))) 112 ENDDO 113 114 CASE("sp") 115 116 DO SUBI = 1, 4 117 118 INDEX = INDEX + 1 119 IF (SUBI .EQ. 1) THEN 120 HK(INDEX,INDEX, 1) = CMPLX(HES(ELEMPOINTER(I))) 121 ELSE 122 HK(INDEX,INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I))) 123 ENDIF 124 125 ENDDO 126 127 CASE("sd") 128 129 DO SUBI = 1, 6 130 131 INDEX = INDEX + 1 132 IF (SUBI .EQ. 1) THEN 133 HK(INDEX,INDEX, 1) = CMPLX(HES(ELEMPOINTER(I))) 134 ELSE 135 HK(INDEX,INDEX, 1) = CMPLX(HED(ELEMPOINTER(I))) 136 ENDIF 137 138 ENDDO 139 140 CASE("sf") 141 142 DO SUBI = 1, 8 143 144 INDEX = INDEX + 1 145 IF (SUBI .EQ. 1) THEN 146 HK(INDEX,INDEX, 1) = CMPLX(HES(ELEMPOINTER(I))) 147 ELSE 148 HK(INDEX,INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I))) 149 ENDIF 150 151 ENDDO 152 153 CASE("pd") 154 155 DO SUBI = 1, 8 156 157 INDEX = INDEX + 1 158 IF (SUBI .LE. 3) THEN 159 HK(INDEX,INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I))) 160 ELSE 161 HK(INDEX,INDEX, 1) = CMPLX(HED(ELEMPOINTER(I))) 162 ENDIF 163 164 ENDDO 165 166 CASE("pf") 167 168 DO SUBI = 1, 10 169 170 INDEX = INDEX + 1 171 IF (SUBI .LE. 3) THEN 172 HK(INDEX,INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I))) 173 ELSE 174 HK(INDEX,INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I))) 175 ENDIF 176 177 ENDDO 178 179 CASE("df") 180 181 DO SUBI = 1, 12 182 183 INDEX = INDEX + 1 184 IF (SUBI .LE. 5) THEN 185 HK(INDEX,INDEX, 1) = CMPLX(HED(ELEMPOINTER(I))) 186 ELSE 187 HK(INDEX,INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I))) 188 ENDIF 189 190 ENDDO 191 192 CASE("spd") 193 194 DO SUBI = 1, 9 195 196 INDEX = INDEX + 1 197 IF (SUBI .EQ. 1) THEN 198 HK(INDEX,INDEX, 1) = CMPLX(HES(ELEMPOINTER(I))) 199 ELSEIF (SUBI .GT. 1 .AND. SUBI .LE. 4) THEN 200 HK(INDEX,INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I))) 201 ELSE 202 HK(INDEX,INDEX, 1) = CMPLX(HED(ELEMPOINTER(I))) 203 ENDIF 204 205 ENDDO 206 207 CASE("spf") 208 209 DO SUBI = 1, 11 210 211 INDEX = INDEX + 1 212 IF (SUBI .EQ. 1) THEN 213 HK(INDEX,INDEX, 1) = CMPLX(HES(ELEMPOINTER(I))) 214 ELSEIF (SUBI .GT. 1 .AND. SUBI .LE. 4) THEN 215 HK(INDEX,INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I))) 216 ELSE 217 HK(INDEX,INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I))) 218 ENDIF 219 220 ENDDO 221 222 CASE("sdf") 223 224 DO SUBI = 1, 13 225 226 INDEX = INDEX + 1 227 IF (SUBI .EQ. 1) THEN 228 HK(INDEX,INDEX, 1) = CMPLX(HES(ELEMPOINTER(I))) 229 ELSEIF (SUBI .GT. 1 .AND. SUBI .LE. 6) THEN 230 HK(INDEX,INDEX, 1) = CMPLX(HED(ELEMPOINTER(I))) 231 ELSE 232 HK(INDEX,INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I))) 233 ENDIF 234 235 ENDDO 236 237 238 CASE("pdf") 239 240 DO SUBI = 1, 15 241 242 INDEX = INDEX + 1 243 IF (SUBI .LE. 3) THEN 244 HK(INDEX,INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I))) 245 ELSEIF (SUBI .GT. 3 .AND. SUBI .LE. 8) THEN 246 HK(INDEX,INDEX, 1) = CMPLX(HED(ELEMPOINTER(I))) 247 ELSE 248 HK(INDEX,INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I))) 249 ENDIF 250 251 ENDDO 252 253 CASE("spdf") 254 255 DO SUBI = 1, 16 256 257 INDEX = INDEX + 1 258 IF (SUBI .EQ. 1) THEN 259 HK(INDEX, INDEX, 1) = CMPLX(HES(ELEMPOINTER(I))) 260 ELSEIF (SUBI .GT. 1 .AND. SUBI .LE. 4) THEN 261 HK(INDEX, INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I))) 262 ELSEIF (SUBI .GT. 1 .AND. SUBI .LE. 4) THEN 263 HK(INDEX, INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I))) 264 ELSEIF (SUBI .GT. 4 .AND. SUBI .LE. 9) THEN 265 HK(INDEX, INDEX, 1) = CMPLX(HED(ELEMPOINTER(I))) 266 ELSE 267 HK(INDEX, INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I))) 268 ENDIF 269 270 ENDDO 271 272 END SELECT 273 274 ENDDO 275 276 DO I = 2, NKTOT 277 DO J = 1, HDIM 278 HK(J,J,I) = HK(J,J,1) 279 ENDDO 280 ENDDO 281 282 ! We assign the diagonal elements in ADDQDEP 283 284 285 IF (BASISTYPE .EQ. "NONORTHO") THEN 286 287 SK = CMPLX(ZERO) 288 289 DO I = 1, NKTOT 290 DO J = 1, HDIM 291 SK(J,J,I) = CMPLX(ONE,ZERO) 292 ENDDO 293 ENDDO 294 295 ENDIF 296 297 K0 = PI*(ONE - REAL(NKX))/(REAL(NKX))*B1 + & 298 PI*(ONE - REAL(NKY))/(REAL(NKY))*B2 + & 299 PI*(ONE - REAL(NKZ))/(REAL(NKZ))*B3 - PI*KSHIFT 300 301!$OMP PARALLEL DO DEFAULT (NONE) & 302!$OMP SHARED(NATS, BASIS, ELEMPOINTER, TOTNEBTB, NEBTB) & 303!$OMP SHARED(CR, BOX, B1, B2, B3, NOINT, ATELE, ELE1, ELE2, HK, SK) & 304!$OMP SHARED(HCUT, SCUT, MATINDLIST, BASISTYPE, K0, NKX, NKY, NKZ) & 305!$OMP PRIVATE(I, J, K, NEWJ, BASISI, BASISJ, INDI, INDJ, PBCI, PBCJ, PBCK) & 306!$OMP PRIVATE(RIJ, MAGR2, MAGR, MAGRP, PHI, ALPHA, BETA, COSBETA) & 307!$OMP PRIVATE(LBRAINC, LBRA, MBRA, L, LKETINC, LKET, MKET) & 308!$OMP PRIVATE(BLOCH, KDOTL, KPOINT, KCOUNT, KHTMP, KSTMP)& 309!$OMP PRIVATE(RCUTTB, IBRA, IKET, AMMBRA, WIGLBRAMBRA, ANGFACTOR, MP) 310!!$OMP REDUCTION(+:HK) 311 312 DO I = 1, NATS 313 314 ! Build the lists of orbitals on each atom 315 316 SELECT CASE(BASIS(ELEMPOINTER(I))) 317 318 CASE("s") 319 BASISI(1) = 0 320 BASISI(2) = -1 321 CASE("p") 322 BASISI(1) = 1 323 BASISI(2) = -1 324 CASE("d") 325 BASISI(1) = 2 326 BASISI(2) = -1 327 CASE("f") 328 BASISI(1) = 3 329 BASISI(2) = -1 330 CASE("sp") 331 BASISI(1) = 0 332 BASISI(2) = 1 333 BASISI(3) = -1 334 CASE("sd") 335 BASISI(1) = 0 336 BASISI(2) = 2 337 BASISI(3) = -1 338 CASE("sf") 339 BASISI(1) = 0 340 BASISI(2) = 3 341 BASISI(3) = -1 342 CASE("pd") 343 BASISI(1) = 1 344 BASISI(2) = 2 345 BASISI(3) = -1 346 CASE("pf") 347 BASISI(1) = 1 348 BASISI(2) = 3 349 BASISI(3) = -1 350 CASE("df") 351 BASISI(1) = 2 352 BASISI(2) = 3 353 BASISI(3) = -1 354 CASE("spd") 355 BASISI(1) = 0 356 BASISI(2) = 1 357 BASISI(3) = 2 358 BASISI(4) = -1 359 CASE("spf") 360 BASISI(1) = 0 361 BASISI(2) = 1 362 BASISI(3) = 3 363 BASISI(4) = -1 364 CASE("sdf") 365 BASISI(1) = 0 366 BASISI(2) = 2 367 BASISI(3) = 3 368 BASISI(4) = -1 369 CASE("pdf") 370 BASISI(1) = 1 371 BASISI(2) = 2 372 BASISI(3) = 3 373 BASISI(4) = -1 374 CASE("spdf") 375 BASISI(1) = 0 376 BASISI(2) = 1 377 BASISI(3) = 2 378 BASISI(4) = 3 379 BASISI(5) = -1 380 END SELECT 381 382 INDI = MATINDLIST(I) 383 384 ! open loop over neighbors J of atom I 385 DO NEWJ = 1, TOTNEBTB(I) 386 387 J = NEBTB(1, NEWJ, I) 388 389 PBCI = NEBTB(2, NEWJ, I) 390 PBCJ = NEBTB(3, NEWJ, I) 391 PBCK = NEBTB(4, NEWJ, I) 392 393 RIJ(1) = CR(1,J) + REAL(PBCI)*BOX(1,1) + REAL(PBCJ)*BOX(2,1) + & 394 REAL(PBCK)*BOX(3,1) - CR(1,I) 395 396 RIJ(2) = CR(2,J) + REAL(PBCI)*BOX(1,2) + REAL(PBCJ)*BOX(2,2) + & 397 REAL(PBCK)*BOX(3,2) - CR(2,I) 398 399 RIJ(3) = CR(3,J) + REAL(PBCI)*BOX(1,3) + REAL(PBCJ)*BOX(2,3) + & 400 REAL(PBCK)*BOX(3,3) - CR(3,I) 401 402 MAGR2 = RIJ(1)*RIJ(1) + RIJ(2)*RIJ(2) + RIJ(3)*RIJ(3) 403 404 RCUTTB = ZERO 405 406 DO K = 1, NOINT 407 408 IF ( (ATELE(I) .EQ. ELE1(K) .AND. & 409 ATELE(J) .EQ. ELE2(K)) .OR. & 410 (ATELE(J) .EQ. ELE1(K) .AND. & 411 ATELE(I) .EQ. ELE2(K) )) THEN 412 413 IF (HCUT(K) .GT. RCUTTB ) RCUTTB = HCUT(K) 414 415 IF (BASISTYPE .EQ. "NONORTHO") THEN 416 IF (SCUT(K) .GT. RCUTTB ) RCUTTB = SCUT(K) 417 ENDIF 418 419 ENDIF 420 421 ENDDO 422 423 IF (MAGR2 .LT. RCUTTB*RCUTTB) THEN 424 425 MAGR = SQRT(MAGR2) 426 427 SELECT CASE(BASIS(ELEMPOINTER(J))) 428 CASE("s") 429 BASISJ(1) = 0 430 BASISJ(2) = -1 431 CASE("p") 432 BASISJ(1) = 1 433 BASISJ(2) = -1 434 CASE("d") 435 BASISJ(1) = 2 436 BASISJ(2) = -1 437 CASE("f") 438 BASISJ(1) = 3 439 BASISJ(2) = -1 440 CASE("sp") 441 BASISJ(1) = 0 442 BASISJ(2) = 1 443 BASISJ(3) = -1 444 CASE("sd") 445 BASISJ(1) = 0 446 BASISJ(2) = 2 447 BASISJ(3) = -1 448 CASE("sf") 449 BASISJ(1) = 0 450 BASISJ(2) = 3 451 BASISJ(3) = -1 452 CASE("pd") 453 BASISJ(1) = 1 454 BASISJ(2) = 2 455 BASISJ(3) = -1 456 CASE("pf") 457 BASISJ(1) = 1 458 BASISJ(2) = 3 459 BASISJ(3) = -1 460 CASE("df") 461 BASISJ(1) = 2 462 BASISJ(2) = 3 463 BASISJ(3) = -1 464 CASE("spd") 465 BASISJ(1) = 0 466 BASISJ(2) = 1 467 BASISJ(3) = 2 468 BASISJ(4) = -1 469 CASE("spf") 470 BASISJ(1) = 0 471 BASISJ(2) = 1 472 BASISJ(3) = 3 473 BASISJ(4) = -1 474 CASE("sdf") 475 BASISJ(1) = 0 476 BASISJ(2) = 2 477 BASISJ(3) = 3 478 BASISJ(4) = -1 479 CASE("pdf") 480 BASISJ(1) = 1 481 BASISJ(2) = 2 482 BASISJ(3) = 3 483 BASISJ(4) = -1 484 CASE("spdf") 485 BASISJ(1) = 0 486 BASISJ(2) = 1 487 BASISJ(3) = 2 488 BASISJ(4) = 3 489 BASISJ(5) = -1 490 END SELECT 491 492 INDJ = MATINDLIST(J) 493 494 MAGRP = SQRT(RIJ(1) * RIJ(1) + RIJ(2) * RIJ(2)) 495 496 ! transform to system in which z-axis is aligned with RIJ, 497 IF (ABS(RIJ(1)) .GT. 1.0E-12) THEN 498 499 IF (RIJ(1) .GT. ZERO .AND. RIJ(2) .GE. ZERO) THEN 500 PHI = ZERO 501 ELSEIF (RIJ(1) .GT. ZERO .AND. RIJ(2) .LT. ZERO) THEN 502 PHI = TWO * PI 503 ELSE 504 PHI = PI 505 ENDIF 506 ALPHA = ATAN(RIJ(2) / RIJ(1)) + PHI 507 508 ELSEIF (ABS(RIJ(2)) .GT. 1.0E-12) THEN 509 510 IF (RIJ(2) .GT. 1.0E-12) THEN 511 ALPHA = PI / TWO 512 ELSE 513 ALPHA = THREE * PI / TWO 514 ENDIF 515 516 ELSE 517 518 ! pathological case: beta=0 and alpha undefined, but 519 ! this doesn't matter for matrix elements 520 521 ALPHA = ZERO 522 523 ENDIF 524 525 COSBETA = RIJ(3)/MAGR 526 BETA = ACOS( RIJ(3) / MAGR ) 527 528 ! Build matrix elements using eqns (1)-(9) in PRB 72 165107 529 530 ! The loops over LBRA and LKET need to take into account 531 ! the orbitals assigned to each atom, e.g., sd rather than 532 ! spd... 533 534 IBRA = INDI + 1 535 536 LBRAINC = 1 537 DO WHILE (BASISI(LBRAINC) .NE. -1) 538 539 LBRA = BASISI(LBRAINC) 540 LBRAINC = LBRAINC + 1 541 542 DO MBRA = -LBRA, LBRA 543 544 ! We can calculate these two outside the 545 ! MKET loop... 546 547 AMMBRA = AM(MBRA, ALPHA) 548 WIGLBRAMBRA = WIGNERD(LBRA, ABS(MBRA), 0, COSBETA) 549 550 IKET = INDJ + 1 551 552 LKETINC = 1 553 DO WHILE (BASISJ(LKETINC) .NE. -1) 554 555 LKET = BASISJ(LKETINC) 556 LKETINC = LKETINC + 1 557 558 DO MKET = -LKET, LKET 559 560 ! This is the sigma bonds (mp = 0) 561 562 ! Hamiltonian build 563 564 ! Pre-compute the angular part so we can use it 565 ! again later if we're building the S matrix too 566 567 ANGFACTOR = TWO * AMMBRA * & 568 AM(MKET, ALPHA) * & 569 WIGLBRAMBRA * & 570 WIGNERD(LKET, ABS(MKET), 0, COSBETA) 571 572 KHTMP = CMPLX(ANGFACTOR * & 573 UNIVSCALE(I, J, LBRA, LKET, & 574 0, MAGR, "H")) 575 576 IF (BASISTYPE .EQ. "NONORTHO") THEN 577 KSTMP = CMPLX(ANGFACTOR * & 578 UNIVSCALE(I, J, LBRA, LKET, & 579 0, MAGR, "S")) 580 ENDIF 581 582 KCOUNT = 0 583 584 !Loop over k-points 585 586 DO KX = 1, NKX 587 588 DO KY = 1, NKY 589 590 DO KZ = 1, NKZ 591 592 KPOINT = TWO*PI*(REAL(KX-1)*B1/REAL(NKX) + & 593 REAL(KY-1)*B2/REAL(NKY) + & 594 REAL(KZ-1)*B3/REAL(NKZ)) + K0 595 596 KCOUNT = KCOUNT+1 597 598 KDOTL = KPOINT(1)*RIJ(1) + KPOINT(2)*RIJ(2) + & 599 KPOINT(3)*RIJ(3) 600 601 BLOCH = EXP(CMPLX(ZERO,KDOTL)) 602 603 HK(IBRA, IKET, KCOUNT) = & 604 HK(IBRA, IKET, KCOUNT) + & 605 BLOCH*KHTMP 606 607 IF (BASISTYPE .EQ. "NONORTHO") THEN 608 SK(IBRA, IKET, KCOUNT) = & 609 SK(IBRA, IKET, KCOUNT) + & 610 BLOCH*KSTMP 611 ENDIF 612 613 ENDDO 614 ENDDO 615 ENDDO 616 617 618 ! everything else 619 620 DO MP = 1, MIN(LBRA, LKET) 621 622 ANGFACTOR = & 623 SLMMP(LBRA, MBRA, MP, ALPHA, COSBETA)* & 624 SLMMP(LKET, MKET, MP, ALPHA, COSBETA)+ & 625 TLMMP(LBRA, MBRA, MP, ALPHA, COSBETA)* & 626 TLMMP(LKET, MKET, MP, ALPHA, COSBETA) 627 628 KHTMP = CMPLX(ANGFACTOR * & 629 UNIVSCALE(I, J, LBRA, LKET, & 630 MP, MAGR, "H")) 631 632 IF (BASISTYPE .EQ. "NONORTHO") THEN 633 KSTMP = CMPLX(ANGFACTOR * & 634 UNIVSCALE(I, J, LBRA, LKET, & 635 MP, MAGR, "S")) 636 ENDIF 637 638 KCOUNT = 0 639 640 DO KX = 1, NKX 641 DO KY = 1, NKY 642 DO KZ = 1, NKZ 643 644 KPOINT = TWO*PI*(REAL(KX-1)*B1/REAL(NKX) + & 645 REAL(KY-1)*B2/REAL(NKY) + & 646 REAL(KZ-1)*B3/REAL(NKZ)) + K0 647 648 KCOUNT = KCOUNT+1 649 650 KDOTL = KPOINT(1)*RIJ(1) + & 651 KPOINT(2)*RIJ(2) + KPOINT(3)*RIJ(3) 652 653 BLOCH = EXP(CMPLX(ZERO,KDOTL)) 654 655 HK(IBRA, IKET, KCOUNT) = & 656 HK(IBRA, IKET, KCOUNT) + & 657 BLOCH*KHTMP 658 659 IF (BASISTYPE .EQ. "NONORTHO") THEN 660 SK(IBRA, IKET, KCOUNT) = & 661 SK(IBRA, IKET, KCOUNT) + & 662 BLOCH*KSTMP 663 ENDIF 664 665 ENDDO 666 ENDDO 667 ENDDO 668 669 ENDDO 670 671 IKET = IKET + 1 672 673 ENDDO 674 675 ENDDO 676 677 IBRA = IBRA + 1 678 679 ENDDO 680 ENDDO 681 ENDIF 682 ENDDO 683 ENDDO 684 685!$OMP END PARALLEL DO 686 687 ! Save the diagonal elements: it will help a lot when we add in the partial charges 688 689 IF (BASISTYPE .EQ. "ORTHO") THEN 690 691 DO I = 1, NKTOT 692 DO J = 1,HDIM 693 HKDIAG(J,I) = HK(J,J,I) 694 ENDDO 695 ENDDO 696 697 ELSEIF (BASISTYPE .EQ. "NONORTHO") THEN 698 699 ! keep a copy of the charge-independent Hamiltonian 700 701 HK0 = HK 702 703 CALL KGENX ! Get the factors for the Lowdin transform 704 705 ENDIF 706 707 708 IF (DEBUGON .EQ. 1 ) THEN 709 710 OPEN(UNIT=31, STATUS="UNKNOWN", FILE="myS0_k_R.dat") 711 712 OPEN(UNIT=32, STATUS="UNKNOWN", FILE="myS0_k_I.dat") 713 714 DO K = 1, NKTOT 715 WRITE(31,*) K 716 WRITE(32,*) K 717 DO I = 1, HDIM 718 WRITE(31,10) (REAL(SK(I,J,K)), J = 1, HDIM) 719 WRITE(32,10) (AIMAG(SK(I,J,K)), J = 1, HDIM) 720 ENDDO 721 ENDDO 722 723 CLOSE(31) 724 CLOSE(32) 725 72610 FORMAT(648F12.6) 727 728 ENDIF 729 730 RETURN 731 732END SUBROUTINE KBLDNEWH 733