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 NEBLISTS(AMIALLO) 23 24 USE CONSTANTS_MOD 25 USE SETUPARRAY 26 USE UNIVARRAY 27 USE PPOTARRAY 28 USE NEBLISTARRAY 29 USE COULOMBARRAY 30 USE MYPRECISION 31 32 IMPLICIT NONE 33 34 INTEGER :: I, J, K, L, M 35 INTEGER :: A, B, C, MYATOMI, MYATOMJ 36 INTEGER :: PBCX, PBCY, PBCZ 37 INTEGER :: BOXX, BOXY, BOXZ 38 INTEGER :: II, JJ, KK 39 INTEGER, INTENT(IN) :: AMIALLO 40 INTEGER :: XRANGE, YRANGE, ZRANGE 41 INTEGER :: MAXNEBTB, MAXNEBPP, MAXNEBCOUL 42 INTEGER :: NCELL(3), NUMCELL 43 INTEGER :: IPIV(3), INFO 44 INTEGER :: MYCELL, BOXID(3), COUNT 45! INTEGER, SAVE :: ALLOCEST 46 INTEGER, ALLOCATABLE :: TOTINCELL(:), CELLLIST(:,:) 47! INTEGER, ALLOCATABLE :: DIMTB(:), DIMPP(:), DIMCOUL(:) 48 REAL(LATTEPREC) :: RIJ(3), MAGR2 49 REAL(LATTEPREC) :: MAGA(3) 50 REAL(LATTEPREC) :: RCUTTB, RCUTCOUL, PPMAX, MAXCUT2 51 REAL(LATTEPREC) :: WORK(3), BOXINV(3,3), S(3) 52 REAL(LATTEPREC), PARAMETER :: MINR = 0.01 53 54 55 IF (PBCON .EQ. 1) CALL PBC 56 57 TOTNEBTB = 0 58 IF (PPOTON .GT. 0) TOTNEBPP = 0 59 IF (ELECTRO .EQ. 1) TOTNEBCOUL = 0 60 61 62 IF (AMIALLO .NE. 0) THEN 63 64 ! Reallocate the neighbor lists based on their size the last time 65 DEALLOCATE(NEBTB) 66 ALLOCATE(NEBTB( 4, MAXDIMTB, NATS )) 67 68 IF (PPOTON .NE. 0) THEN 69 DEALLOCATE(NEBPP) 70 ALLOCATE(NEBPP( 4, MAXDIMPP, NATS )) 71 ENDIF 72 73 IF (ELECTRO .NE. 1) THEN 74 DEALLOCATE(NEBCOUL) 75 ALLOCATE(NEBCOUL( 4, MAXDIMCOUL, NATS)) 76 ENDIF 77 78 ELSE 79 80 ! This bit is only done on the first neighbor list build 81 82 ! Let's get the cut-offs for our interactions 83 84 RCUTTB = ZERO 85 PPMAX = ZERO 86 87 ! Find the maximum cut off 88 89 DO K = 1, NOINT 90 91 IF (BOND(8,K) .GT. RCUTTB ) RCUTTB = BOND(8,K) 92 93 IF (BASISTYPE .EQ. "NONORTHO") THEN 94 IF (OVERL(8,K) .GT. RCUTTB ) RCUTTB = OVERL(8,K) 95 ENDIF 96 97 ENDDO 98 99 IF (PPOTON .GT. 0) THEN 100 101 DO K = 1, NOPPS 102 103 IF (PPOTON .EQ. 1 .AND. POTCOEF(10,K) .GT. PPMAX ) PPMAX = POTCOEF(10,K) 104 105 IF (PPOTON .EQ. 2 .AND. PPR(PPTABLENGTH(K), K) .GT. PPMAX) & 106 PPMAX = PPR(PPTABLENGTH(K), K) 107 108 ENDDO 109 110 ENDIF 111 112 RCUTTB = RCUTTB + SKIN 113 RCUTTB2 = RCUTTB*RCUTTB 114 115 IF (PPOTON .GT. 0) THEN 116 PPMAX = PPMAX + SKIN 117 PPMAX2 = PPMAX*PPMAX 118 ELSE 119 PPMAX = ZERO 120 PPMAX2 = ZERO 121 ENDIF 122 123 RCUTCOUL = COULCUT + SKIN 124 RCUTCOUL2 = RCUTCOUL * RCUTCOUL 125 126 IF (ELECTRO .EQ. 0) RCUTCOUL = ZERO 127 128 MAXCUT = MAX(RCUTTB, PPMAX, RCUTCOUL) 129 130 MAXCUT2 = MAXCUT*MAXCUT 131 132 ! Now let's estimate the size of the arrays we need for to 133 ! store the neighbor lists, plus some 134 135 IF (PBCON .EQ. 1) THEN 136 137 XRANGE = INT(MAXCUT/BOX(1,1)) + 1 138 YRANGE = INT(MAXCUT/BOX(2,2)) + 1 139 ZRANGE = INT(MAXCUT/BOX(3,3)) + 1 140! print*, maxcut, xrange, yrange, zrange 141 142 ! Here we're hoping atom 1 is in a typical environment 143 144 COUNT = 0 145 DO J = 1, NATS 146 DO II = -XRANGE, XRANGE 147 DO JJ = -YRANGE, YRANGE 148 DO KK = -ZRANGE, ZRANGE 149 150 RIJ(1) = CR(1,J) + REAL(II)*BOX(1,1) + & 151 REAL(JJ)*BOX(2,1) + REAL(KK)*BOX(3,1) - CR(1,1) 152 153 RIJ(2) = CR(2,J) + REAL(II)*BOX(1,2) + & 154 REAL(JJ)*BOX(2,2) + REAL(KK)*BOX(3,2) - CR(2,1) 155 156 RIJ(3) = CR(3,J) + REAL(II)*BOX(1,3) + & 157 REAL(JJ)*BOX(2,3) + REAL(KK)*BOX(3,3) - CR(3,1) 158 159 MAGR2 = RIJ(1)*RIJ(1) + RIJ(2)*RIJ(2) + RIJ(3)*RIJ(3) 160 161 IF (MAGR2 .LE. MAXCUT2) COUNT = COUNT + 1 162 163 ENDDO 164 ENDDO 165 ENDDO 166 ENDDO 167 168 MAXDIMTB = 2*COUNT 169 MAXDIMPP = 2*COUNT 170 MAXDIMCOUL = 2*COUNT 171 172 ELSEIF (PBCON .EQ. 0) THEN 173 174 DIMLIST = 0 175 DO I = 1, NATS 176 COUNT = 0 177 DO J = 1, NATS 178 179 RIJ(1) = CR(1,J) - CR(1,I) 180 RIJ(2) = CR(2,J) - CR(2,I) 181 RIJ(3) = CR(3,J) - CR(3,I) 182 183 MAGR2 = RIJ(1)*RIJ(1) + RIJ(2)*RIJ(2) + RIJ(3)*RIJ(3) 184 185 IF (MAGR2 .LE. MAXCUT2) COUNT = COUNT + 1 186 187 ENDDO 188 189 IF (COUNT .GT. DIMLIST) DIMLIST = COUNT 190 191 ENDDO 192 193 MAXDIMTB = DIMLIST 194 MAXDIMPP = DIMLIST 195 MAXDIMCOUL = DIMLIST 196 197 ENDIF 198 199 ALLOCATE ( NEBTB( 4, MAXDIMTB, NATS ) ) 200 IF (PPOTON .GT. 0) ALLOCATE(NEBPP( 4, MAXDIMPP, NATS ) ) 201 IF (ELECTRO .EQ. 1) ALLOCATE(NEBCOUL( 4, MAXDIMCOUL, NATS )) 202 203 ENDIF 204 205 ! Now build the neighbor list 206 207 ! With periodic boundaries first: 208 209 IF (PBCON .EQ. 1) THEN 210 211 MAGA(1) = SQRT(BOX(1,1)*BOX(1,1) + BOX(1,2)*BOX(1,2) + BOX(1,3)*BOX(1,3)) 212 MAGA(2) = SQRT(BOX(2,1)*BOX(2,1) + BOX(2,2)*BOX(2,2) + BOX(2,3)*BOX(2,3)) 213 MAGA(3) = SQRT(BOX(3,1)*BOX(3,1) + BOX(3,2)*BOX(3,2) + BOX(3,3)*BOX(3,3)) 214 215 XRANGE = INT(MAXCUT/MAGA(1)) + 1 216 YRANGE = INT(MAXCUT/MAGA(2)) + 1 217 ZRANGE = INT(MAXCUT/MAGA(3)) + 1 218 219 ! This gives the number of sub-cells along each lattice vector 220 221 NCELL(1) = MAX(INT(MAGA(1)/MAXCUT),1) 222 NCELL(2) = MAX(INT(MAGA(2)/MAXCUT),1) 223 NCELL(3) = MAX(INT(MAGA(3)/MAXCUT),1) 224 225 NUMCELL = NCELL(1)*NCELL(2)*NCELL(3) 226 227! PRINT*, NCELL(1), NCELL(2), NCELL(3), NUMCELL 228 229 IF (AMIALLO .EQ. 0) ALLOCEST = 2*NATS/NUMCELL 230 231 ALLOCATE(TOTINCELL(NUMCELL), CELLLIST(ALLOCEST, NUMCELL)) 232 233 TOTINCELL = 0 234 235 BOXINV = BOX 236 237 CALL DGETRF(3, 3, BOXINV, 3, IPIV, INFO) 238 239 CALL DGETRI(3, BOXINV, 3, IPIV, WORK, 3, INFO) 240 241 ! Put the atoms into the sub-cells 242 243 DO I = 1, NATS 244 245 CALL DGEMV('T', 3, 3, ONE, BOXINV, 3, CR(1,I), 1, ZERO, S, 1) 246 247 BOXID(1) = INT(S(1)*NCELL(1)) 248 BOXID(2) = INT(S(2)*NCELL(2)) 249 BOXID(3) = INT(S(3)*NCELL(3)) 250 251 MYCELL = BOXID(1) + NCELL(1)*BOXID(2) + NCELL(1)*NCELL(2)*BOXID(3) + 1 252 253 TOTINCELL(MYCELL) = TOTINCELL(MYCELL) + 1 254 CELLLIST(TOTINCELL(MYCELL), MYCELL) = I 255 256 ENDDO 257 258 ALLOCEST = 2*MAXVAL(TOTINCELL) 259 260 ! Loop over the subcells and build the lists 261 262 DO II = 1, NUMCELL 263 264 ! Indices of subcell II 265 266 BOXZ = (II - 1)/(NCELL(1)*NCELL(2)) 267 BOXY = (II - 1 - BOXZ*NCELL(1)*NCELL(2))/NCELL(1) 268 BOXX = (II - 1 - NCELL(1)*BOXY - NCELL(1)*NCELL(2)*BOXZ) 269! print*, boxx, boxy, boxz, totincell(ii), celllist(1,ii), NUMCELL 270 271 ! Loop over atoms in cell II 272 273 DO I = 1, TOTINCELL(II) 274 275 MYATOMI = CELLLIST(I,II) 276 277 ! Loop over the neighboring subcells 278 279 DO A = BOXZ - ZRANGE, BOXZ + ZRANGE 280 DO B = BOXY - YRANGE, BOXY + YRANGE 281 DO C = BOXX - XRANGE, BOXX + XRANGE 282 283 PBCX = 0 284 PBCY = 0 285 PBCZ = 0 286 287 BOXID(1) = C 288 BOXID(2) = B 289 BOXID(3) = A 290 291 IF (A .LT. 0) THEN 292 PBCZ = A 293 BOXID(3) = NCELL(3) - 1 294 ELSEIF (A .GE. NCELL(3)) THEN 295 PBCZ = A - NCELL(3) + 1 296 BOXID(3) = 0 297 ENDIF 298 299 IF (B .LT. 0) THEN 300 PBCY = B 301 BOXID(2) = NCELL(2) - 1 302 ELSEIF (B .GE. NCELL(2)) THEN 303 PBCY = B - NCELL(2) + 1 304 BOXID(2) = 0 305 ENDIF 306 307 IF (C .LT. 0) THEN 308 PBCX = C 309 BOXID(1) = NCELL(1) - 1 310 ELSEIF (C .GE. NCELL(1)) THEN 311 PBCX = C - NCELL(1) + 1 312 BOXID(1) = 0 313 ENDIF 314 315 MYCELL = BOXID(1) + NCELL(1)*BOXID(2) + NCELL(1)*NCELL(2)*BOXID(3) + 1 316 317 ! Loop over the atoms in the neighboring cell 318 319 DO J = 1, TOTINCELL(MYCELL) 320 321 MYATOMJ = CELLLIST(J, MYCELL) 322 323 RIJ(1) = CR(1,MYATOMJ) + REAL(PBCX)*BOX(1,1) + & 324 REAL(PBCY)*BOX(2,1) + REAL(PBCZ)*BOX(3,1) - CR(1,MYATOMI) 325 326 RIJ(2) = CR(2,MYATOMJ) + REAL(PBCX)*BOX(1,2) + & 327 REAL(PBCY)*BOX(2,2) + REAL(PBCZ)*BOX(3,2) - CR(2,MYATOMI) 328 329 RIJ(3) = CR(3,MYATOMJ) + REAL(PBCX)*BOX(1,3) + & 330 REAL(PBCY)*BOX(2,3) + REAL(PBCZ)*BOX(3,3) - CR(3,MYATOMI) 331 332 MAGR2 = RIJ(1)*RIJ(1) + RIJ(2)*RIJ(2) + RIJ(3)*RIJ(3) 333 334 IF (MAGR2 .GT. MINR .AND. MAGR2 .LT. RCUTTB2) THEN 335 336 TOTNEBTB(MYATOMI) = TOTNEBTB(MYATOMI) + 1 337 338 IF (TOTNEBTB(MYATOMI) .GT. MAXDIMTB) THEN 339 CALL ERRORS("neblists","NUMBER OF NEIGHBORS EXCEEDS ARRAY DIMENSION (TB)") 340 ENDIF 341 342 NEBTB( 1, TOTNEBTB(MYATOMI), MYATOMI ) = MYATOMJ 343 NEBTB( 2, TOTNEBTB(MYATOMI), MYATOMI ) = PBCX 344 NEBTB( 3, TOTNEBTB(MYATOMI), MYATOMI ) = PBCY 345 NEBTB( 4, TOTNEBTB(MYATOMI), MYATOMI ) = PBCZ 346 347 ENDIF 348 349 IF (PPOTON .NE. 0) THEN 350 351 IF (MAGR2 .GT. MINR .AND. MAGR2 .LT. PPMAX2) THEN 352 353 TOTNEBPP(MYATOMI) = TOTNEBPP(MYATOMI) + 1 354 355 IF (TOTNEBPP(MYATOMI) .GT. MAXDIMPP) THEN 356 CALL ERRORS("neblists","NUMBER OF NEIGHBORS EXCEEDS ARRAY DIMENSION (PP)") 357 ENDIF 358 359 NEBPP( 1, TOTNEBPP(MYATOMI), MYATOMI ) = MYATOMJ 360 NEBPP( 2, TOTNEBPP(MYATOMI), MYATOMI ) = PBCX 361 NEBPP( 3, TOTNEBPP(MYATOMI), MYATOMI ) = PBCY 362 NEBPP( 4, TOTNEBPP(MYATOMI), MYATOMI ) = PBCZ 363 364 ENDIF 365 366 ENDIF 367 368 IF (ELECTRO .NE. 0) THEN 369 370 IF (MAGR2 .GT. MINR .AND. MAGR2 .LT. RCUTCOUL2) THEN 371 372 TOTNEBCOUL(MYATOMI) = TOTNEBCOUL(MYATOMI) + 1 373 374 IF (TOTNEBCOUL(MYATOMI) .GT. MAXDIMCOUL) THEN 375 CALL ERRORS("neblists","NUMBER OF NEIGHBORS EXCEEDS ARRAY DIMENSION (COUL)") 376 ENDIF 377 378 NEBCOUL( 1, TOTNEBCOUL(MYATOMI), MYATOMI ) = MYATOMJ 379 NEBCOUL( 2, TOTNEBCOUL(MYATOMI), MYATOMI ) = PBCX 380 NEBCOUL( 3, TOTNEBCOUL(MYATOMI), MYATOMI ) = PBCY 381 NEBCOUL( 4, TOTNEBCOUL(MYATOMI), MYATOMI ) = PBCZ 382 383 ENDIF 384 385 ENDIF 386 387 388 ENDDO 389 ENDDO 390 ENDDO 391 392 ENDDO 393 394 ENDDO 395 396 ENDDO 397 398 DEALLOCATE(TOTINCELL, CELLLIST) 399 400 ELSEIF (PBCON .EQ. 0) THEN 401 402 ! Now we're doing building the neighbor lists for gas-phase systems 403 404 DO I = 1, NATS 405 DO J = 1, NATS 406 407 RIJ(1) = CR(1,J) - CR(1,I) 408 409 RIJ(2) = CR(2,J) - CR(2,I) 410 411 RIJ(3) = CR(3,J) - CR(3,I) 412 413 MAGR2 = RIJ(1)*RIJ(1) + RIJ(2)*RIJ(2) + RIJ(3)*RIJ(3) 414 415 IF (MAGR2 .GT. MINR .AND. MAGR2 .LT. RCUTTB2) THEN 416 417 TOTNEBTB(I) = TOTNEBTB(I) + 1 418 NEBTB( 1, TOTNEBTB(I), I ) = J 419 NEBTB( 2, TOTNEBTB(I), I ) = 0 420 NEBTB( 3, TOTNEBTB(I), I ) = 0 421 NEBTB( 4, TOTNEBTB(I), I ) = 0 422 423 ENDIF 424 425 IF (PPOTON .NE. 0) THEN 426 427 IF (MAGR2 .GT. MINR .AND. MAGR2 .LT. PPMAX2) THEN 428 429 TOTNEBPP(I) = TOTNEBPP(I) + 1 430 NEBPP( 1, TOTNEBPP(I), I ) = J 431 NEBPP( 2, TOTNEBPP(I), I ) = 0 432 NEBPP( 3, TOTNEBPP(I), I ) = 0 433 NEBPP( 4, TOTNEBPP(I), I ) = 0 434 435 ENDIF 436 437 ENDIF 438 439 IF (ELECTRO .NE. 0) THEN 440 441 IF (MAGR2 .GT. MINR .AND. MAGR2 .LT. RCUTCOUL2) THEN 442 443 TOTNEBCOUL(I) = TOTNEBCOUL(I) + 1 444 NEBCOUL( 1, TOTNEBCOUL(I), I ) = J 445 NEBCOUL( 2, TOTNEBCOUL(I), I ) = 0 446 NEBCOUL( 3, TOTNEBCOUL(I), I ) = 0 447 NEBCOUL( 4, TOTNEBCOUL(I), I ) = 0 448 449 ENDIF 450 451 ENDIF 452 453 ENDDO 454 ENDDO 455 456 ENDIF 457 458 ! Let's get the dimensions of the arrays about right for the next 459 ! loop through here 460 461 MAXDIMTB = 2*MAXVAL(TOTNEBTB) 462 IF (PPOTON .NE. 0) MAXDIMPP = 2*MAXVAL(TOTNEBPP) 463 IF (ELECTRO .NE. 0) MAXDIMCOUL = 2*MAXVAL(TOTNEBCOUL) 464 465 RETURN 466 467END SUBROUTINE NEBLISTS 468