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 GRADH 23 24 USE CONSTANTS_MOD 25 USE SETUPARRAY 26 USE NEBLISTARRAY 27 USE UNIVARRAY 28 USE SPINARRAY 29 USE VIRIALARRAY 30 USE MYPRECISION 31 32 IMPLICIT NONE 33 34 INTEGER :: I, J, K, L, M, N, KK, INDI, INDJ 35 INTEGER :: LBRA, MBRA, LKET, MKET 36 INTEGER :: PREVJ, NEWJ 37 INTEGER :: PBCI, PBCJ, PBCK 38 INTEGER :: BASISI(5), BASISJ(5), LBRAINC, LKETINC 39 REAL(LATTEPREC) :: ALPHA, BETA, PHI, RHO, COSBETA 40 REAL(LATTEPREC) :: RIJ(3), DC(3) 41 REAL(LATTEPREC) :: MAGR, MAGR2, MAGRP, MAGRP2, FTMP(3) 42 REAL(LATTEPREC) :: MAXRCUT, MAXRCUT2 43 REAL(LATTEPREC) :: MYDFDA, MYDFDB, MYDFDR, RCUTTB 44 REAL(LATTEPREC), EXTERNAL :: DFDA, DFDB, DFDR 45 LOGICAL PATH 46 IF (EXISTERROR) RETURN 47 48 F = ZERO 49 VIRBOND = ZERO 50 51!$OMP PARALLEL DO DEFAULT (NONE) & 52!$OMP SHARED(NATS, BASIS, ELEMPOINTER, TOTNEBTB, NEBTB) & 53!$OMP SHARED(CR, BOX, BO, RHOUP, RHODOWN, SPINON, NOINT, ATELE, ELE1, ELE2) & 54!$OMP SHARED(HCUT, SCUT, MATINDLIST, BASISTYPE) & 55!$OMP SHARED(F) & 56!$OMP PRIVATE(I, J, K, NEWJ, BASISI, BASISJ, INDI, INDJ, PBCI, PBCJ, PBCK) & 57!$OMP PRIVATE(RIJ, MAGR2, MAGR, MAGRP2, MAGRP, PATH, PHI, ALPHA, BETA, COSBETA, FTMP) & 58!$OMP PRIVATE(DC, LBRAINC, LBRA, MBRA, L, LKETINC, LKET, MKET, RHO) & 59!$OMP PRIVATE(MYDFDA, MYDFDB, MYDFDR, RCUTTB) & 60!$OMP REDUCTION(+:VIRBOND) 61 62 DO I = 1, NATS 63 64 ! Build list of orbitals on atom I 65 SELECT CASE(BASIS(ELEMPOINTER(I))) 66 67 CASE("s") 68 BASISI(1) = 0 69 BASISI(2) = -1 70 CASE("p") 71 BASISI(1) = 1 72 BASISI(2) = -1 73 CASE("d") 74 BASISI(1) = 2 75 BASISI(2) = -1 76 CASE("f") 77 BASISI(1) = 3 78 BASISI(2) = -1 79 CASE("sp") 80 BASISI(1) = 0 81 BASISI(2) = 1 82 BASISI(3) = -1 83 CASE("sd") 84 BASISI(1) = 0 85 BASISI(2) = 2 86 BASISI(3) = -1 87 CASE("sf") 88 BASISI(1) = 0 89 BASISI(2) = 3 90 BASISI(3) = -1 91 CASE("pd") 92 BASISI(1) = 1 93 BASISI(2) = 2 94 BASISI(3) = -1 95 CASE("pf") 96 BASISI(1) = 1 97 BASISI(2) = 3 98 BASISI(3) = -1 99 CASE("df") 100 BASISI(1) = 2 101 BASISI(2) = 3 102 BASISI(3) = -1 103 CASE("spd") 104 BASISI(1) = 0 105 BASISI(2) = 1 106 BASISI(3) = 2 107 BASISI(4) = -1 108 CASE("spf") 109 BASISI(1) = 0 110 BASISI(2) = 1 111 BASISI(3) = 3 112 BASISI(4) = -1 113 CASE("sdf") 114 BASISI(1) = 0 115 BASISI(2) = 2 116 BASISI(3) = 3 117 BASISI(4) = -1 118 CASE("pdf") 119 BASISI(1) = 1 120 BASISI(2) = 2 121 BASISI(3) = 3 122 BASISI(4) = -1 123 CASE("spdf") 124 BASISI(1) = 0 125 BASISI(2) = 1 126 BASISI(3) = 2 127 BASISI(4) = 3 128 BASISI(5) = -1 129 END SELECT 130 131 INDI = MATINDLIST(I) 132 133 DO NEWJ = 1, TOTNEBTB(I) 134 135 J = NEBTB(1, NEWJ, I) 136 PBCI = NEBTB(2, NEWJ, I) 137 PBCJ = NEBTB(3, NEWJ, I) 138 PBCK = NEBTB(4, NEWJ, I) 139 140 RIJ(1) = CR(1,J) + REAL(PBCI)*BOX(1,1) + REAL(PBCJ)*BOX(2,1) + & 141 REAL(PBCK)*BOX(3,1) - CR(1,I) 142 143 RIJ(2) = CR(2,J) + REAL(PBCI)*BOX(1,2) + REAL(PBCJ)*BOX(2,2) + & 144 REAL(PBCK)*BOX(3,2) - CR(2,I) 145 146 RIJ(3) = CR(3,J) + REAL(PBCI)*BOX(1,3) + REAL(PBCJ)*BOX(2,3) + & 147 REAL(PBCK)*BOX(3,3) - CR(3,I) 148 149 MAGR2 = RIJ(1)*RIJ(1) + RIJ(2)*RIJ(2) + RIJ(3)*RIJ(3) 150 151 ! Building the forces is expensive - use the cut-off 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 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 ! PRINT*, ALPHA, BETA 275 276 DC = RIJ/MAGR 277 278 ! build forces using PRB 72 165107 eq. (12) - the sign of the 279 ! dfda contribution seems to be wrong, but gives the right 280 ! answer(?) 281 282 FTMP = ZERO 283 K = INDI 284 285 LBRAINC = 1 286 DO WHILE (BASISI(LBRAINC) .NE. -1) 287 288 LBRA = BASISI(LBRAINC) 289 LBRAINC = LBRAINC + 1 290 291 DO MBRA = -LBRA, LBRA 292 293 K = K + 1 294 L = INDJ 295 296 LKETINC = 1 297 DO WHILE (BASISJ(LKETINC) .NE. -1) 298 299 LKET = BASISJ(LKETINC) 300 LKETINC = LKETINC + 1 301 302 DO MKET = -LKET, LKET 303 304 L = L + 1 305 306 SELECT CASE(SPINON) 307 CASE(0) 308 RHO = BO(L, K) 309 CASE(1) 310 RHO = RHOUP(L, K) + RHODOWN(L, K) 311 END SELECT 312 313 IF (.NOT. PATH) THEN 314 315 ! Unroll loops and pre-compute 316 317 MYDFDA = DFDA(I, J, LBRA, LKET, MBRA, & 318 MKET, MAGR, ALPHA, COSBETA, "H") 319 320 MYDFDB = DFDB(I, J, LBRA, LKET, MBRA, & 321 MKET, MAGR, ALPHA, COSBETA, "H") 322 323 MYDFDR = DFDR(I, J, LBRA, LKET, MBRA, & 324 MKET, MAGR, ALPHA, COSBETA, "H") 325 326 ! 327 ! d/d_alpha 328 ! 329 330 FTMP(1) = FTMP(1) + RHO * & 331 (-RIJ(2) / MAGRP2 * MYDFDA) 332 333 FTMP(2) = FTMP(2) + RHO * & 334 (RIJ(1)/ MAGRP2 * MYDFDA) 335 336 ! 337 ! d/d_beta 338 ! 339 340 FTMP(1) = FTMP(1) + RHO * & 341 (((((RIJ(3) * RIJ(1)) / & 342 MAGR2)) / MAGRP) * MYDFDB) 343 344 FTMP(2) = FTMP(2) + RHO * & 345 (((((RIJ(3) * RIJ(2)) / & 346 MAGR2)) / MAGRP) * MYDFDB) 347 348 FTMP(3) = FTMP(3) - RHO * & 349 (((ONE - ((RIJ(3) * RIJ(3)) / & 350 MAGR2)) / MAGRP) * MYDFDB) 351 352 ! 353 ! d/dR 354 ! 355 356 FTMP(1) = FTMP(1) - RHO * DC(1) * & 357 MYDFDR 358 359 FTMP(2) = FTMP(2) - RHO * DC(2) * & 360 MYDFDR 361 362 FTMP(3) = FTMP(3) - RHO * DC(3) * & 363 MYDFDR 364 365 366 ELSE 367 368 ! pathological configuration in which beta=0 369 ! or pi => alpha undefined 370 371 ! fixed: MJC 12/17/13 372 373 MYDFDB = DFDB(I, J, LBRA, LKET, & 374 MBRA, MKET, MAGR, ZERO, COSBETA, "H") / MAGR 375 376 FTMP(1) = FTMP(1) - RHO * (COSBETA * MYDFDB) 377 378 MYDFDB = DFDB(I, J, LBRA, LKET, & 379 MBRA, MKET, MAGR, PI/TWO, COSBETA, "H") / MAGR 380 381 FTMP(2) = FTMP(2) - RHO * (COSBETA * MYDFDB) 382 383 MYDFDR = DFDR(I, J, LBRA, LKET, MBRA, & 384 MKET, MAGR, ZERO, COSBETA, "H") 385 386 FTMP(3) = FTMP(3) - RHO * COSBETA * MYDFDR 387 388 ENDIF 389 390 ENDDO 391 ENDDO 392 ENDDO 393 ENDDO 394 395 F(1,I) = F(1,I) + FTMP(1) 396 F(2,I) = F(2,I) + FTMP(2) 397 F(3,I) = F(3,I) + FTMP(3) 398 399 VIRBOND(1) = VIRBOND(1) + RIJ(1) * FTMP(1) 400 VIRBOND(2) = VIRBOND(2) + RIJ(2) * FTMP(2) 401 VIRBOND(3) = VIRBOND(3) + RIJ(3) * FTMP(3) 402 VIRBOND(4) = VIRBOND(4) + RIJ(1) * FTMP(2) 403 VIRBOND(5) = VIRBOND(5) + RIJ(2) * FTMP(3) 404 VIRBOND(6) = VIRBOND(6) + RIJ(3) * FTMP(1) 405 406 ENDIF 407 408 ENDDO 409 410 ! INDI = INDI + NORBI 411 412 ENDDO 413 414!$OMP END PARALLEL DO 415 416 RETURN 417 418END SUBROUTINE GRADH 419