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