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