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