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