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