-
Notifications
You must be signed in to change notification settings - Fork 105
/
Copy pathlink_list_module.f90
343 lines (257 loc) · 13.7 KB
/
link_list_module.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
! link_list_module.f90
! Link list handling routines for MC or MD simulation
MODULE link_list_module
!------------------------------------------------------------------------------------------------!
! This software was written in 2016/17 !
! by Michael P. Allen <[email protected]>/<[email protected]> !
! and Dominic J. Tildesley <[email protected]> ("the authors"), !
! to accompany the book "Computer Simulation of Liquids", second edition, 2017 ("the text"), !
! published by Oxford University Press ("the publishers"). !
! !
! LICENCE !
! Creative Commons CC0 Public Domain Dedication. !
! To the extent possible under law, the authors have dedicated all copyright and related !
! and neighboring rights to this software to the PUBLIC domain worldwide. !
! This software is distributed without any warranty. !
! You should have received a copy of the CC0 Public Domain Dedication along with this software. !
! If not, see <http://creativecommons.org/publicdomain/zero/1.0/>. !
! !
! DISCLAIMER !
! The authors and publishers make no warranties about the software, and disclaim liability !
! for all uses of the software, to the fullest extent permitted by applicable law. !
! The authors and publishers do not recommend use of this software for any purpose. !
! It is made freely available, solely to clarify points made in the text. When using or citing !
! the software, you should not imply endorsement by the authors or publishers. !
!------------------------------------------------------------------------------------------------!
USE, INTRINSIC :: iso_fortran_env, ONLY : output_unit, error_unit
IMPLICIT NONE
PRIVATE
! Public routines
PUBLIC :: initialize_list, extend_list, finalize_list, make_list, check_list
PUBLIC :: move_in_list, create_in_list, destroy_in_list, c_index, neighbours
! Public (protected) data
INTEGER, PROTECTED, SAVE, PUBLIC :: sc ! dimensions of head array
INTEGER, DIMENSION(:,:,:), ALLOCATABLE, PROTECTED, SAVE, PUBLIC :: head ! head(0:sc-1,0:sc-1,0:sc-1)
INTEGER, DIMENSION(:), ALLOCATABLE, PROTECTED, SAVE, PUBLIC :: list ! list(n)
INTEGER, DIMENSION(:,:), ALLOCATABLE, PROTECTED, SAVE, PUBLIC :: c ! c(3,n) 3D cell index of each atom
! Private data
INTEGER, DIMENSION(:), ALLOCATABLE :: tmp_list ! for dynamic reallocation of list(n)
INTEGER, DIMENSION(:,:), ALLOCATABLE :: tmp_c ! for dynamic reallocation of c(3,n)
CONTAINS
SUBROUTINE initialize_list ( n, r_cut_box ) ! Routine to allocate list arrays
IMPLICIT NONE
INTEGER, INTENT(in) :: n ! Number of particles
REAL, INTENT(in) :: r_cut_box ! rcut/box
WRITE ( unit=output_unit, fmt='(a,t40,f15.6)') 'Link cells based on r_cut/box', r_cut_box
sc = FLOOR ( 1.0 / r_cut_box ) ! Number of cells in each dimension
WRITE ( unit=output_unit, fmt='(a,t40,i15)') 'Number of cells in box', sc
IF ( sc < 3 ) THEN
WRITE ( unit=error_unit, fmt='(a,i15)') 'System is too small to use link cells', sc
STOP 'Error in initialize_list'
END IF
ALLOCATE ( list(n), c(3,n) )
ALLOCATE ( head(0:sc-1,0:sc-1,0:sc-1) )
END SUBROUTINE initialize_list
SUBROUTINE extend_list ( n_new )
IMPLICIT NONE
INTEGER, INTENT(in) :: n_new ! New size of arrays
INTEGER :: n_old
n_old = SIZE ( list )
IF ( n_new <= n_old ) THEN
WRITE ( unit=error_unit, fmt='(a,2i15)') 'Unexpected array size ', n_old, n_new
STOP 'Error in extend_list'
END IF
! Reallocate list array, preserving existing data
ALLOCATE ( tmp_list(n_new) )
tmp_list(1:n_old) = list
tmp_list(n_old+1:) = 0
CALL MOVE_ALLOC ( tmp_list, list )
! Reallocate c array, preserving existing data
ALLOCATE ( tmp_c(3,n_new) )
tmp_c(:,1:n_old) = c
tmp_c(:,n_old+1:) = 0
CALL MOVE_ALLOC ( tmp_c, c )
END SUBROUTINE extend_list
SUBROUTINE finalize_list ! Routine to deallocate list arrays
IMPLICIT NONE
DEALLOCATE ( list, c )
DEALLOCATE ( head )
END SUBROUTINE finalize_list
SUBROUTINE make_list ( n, r_cut_box, r ) ! Routine to make list
IMPLICIT NONE
INTEGER, INTENT(in) :: n ! Number of atoms
REAL, INTENT(in) :: r_cut_box ! rcut/box
REAL, DIMENSION(3,n), INTENT(in) :: r ! Atom coordinates
INTEGER :: i
! Allow for possibility of box size change
sc = FLOOR ( 1.0 / r_cut_box )
IF ( sc < 3 ) THEN
WRITE ( unit=error_unit, fmt='(a,i15)') 'System is too small to use link cells', sc
STOP 'Error in make_list'
END IF
IF ( sc > SIZE ( head, dim=1 ) ) THEN
DEALLOCATE ( head )
ALLOCATE ( head(0:sc-1,0:sc-1,0:sc-1) )
END IF
head(:,:,:) = 0
DO i = 1, n ! Loop over all atoms
c(:,i) = c_index ( r(:,i) ) ! Index function allocating atom i to cell
CALL create_in_list ( i, c(:,i) ) ! This does the work of adding atom i to list
END DO ! End loop over all atoms
END SUBROUTINE make_list
FUNCTION c_index ( ri ) RESULT ( ci )
IMPLICIT NONE
INTEGER, DIMENSION(3) :: ci ! Returns 3D cell index in range 0..sc-1, calculated from
REAL, DIMENSION(3), INTENT(in) :: ri ! position in box = 1 units
! We do not want to do any periodic imaging here, so as to cope with
! Lees-Edwards boundaries as well as normal ones
! But we must check that ri is within bounds
IF ( ANY ( ABS(ri) > 0.5 ) ) THEN ! Should never happen
WRITE ( unit=error_unit, fmt='(a,3f15.6)') 'Atom not in main box', ri
STOP 'Error in c_index'
END IF
ci(:) = FLOOR ( ( ri(:) + 0.5 ) * REAL(sc) ) ! The index formula
! Guard against small chance of roundoff error
WHERE ( ci(:) < 0 ) ci(:) = 0
WHERE ( ci(:) > sc-1 ) ci(:) = sc-1
END FUNCTION c_index
SUBROUTINE create_in_list ( i, ci ) ! Routine to create atom i in cell ci
IMPLICIT NONE
INTEGER, INTENT(in) :: i ! Index of atom
INTEGER, DIMENSION(3), INTENT(in) :: ci ! 3D index of cell in which i lies
list(i) = head(ci(1),ci(2),ci(3)) ! Transfer old head to list
head(ci(1),ci(2),ci(3)) = i ! Atom i becomes new head for this list
c(:,i) = ci(:) ! Store 3D index in array
END SUBROUTINE create_in_list
SUBROUTINE destroy_in_list ( i, ci ) ! Routine to destroy atom i in cell ci
IMPLICIT NONE
INTEGER, INTENT(in) :: i ! Index of atom
INTEGER, DIMENSION(3), INTENT(in) :: ci ! 3D index of cell in which i lies
INTEGER :: this, next
this = head(ci(1),ci(2),ci(3)) ! Locate head of list corresponding to cell
IF ( this == i ) THEN ! Atom i is the head atom in that cell
head(ci(1),ci(2),ci(3)) = list(i) ! Simply point head at next atom, we're done
ELSE ! Atom i lies further down the list
DO ! Loop traversing link-list
next = list(this) ! Look ahead to the next entry
IF ( next == i ) THEN ! Found our atom, just link over it
list(this) = list(i)
EXIT ! Leave the loop
ELSE IF ( next == 0 ) THEN ! This should never happen
WRITE ( unit=error_unit, fmt='(a,4i15)') 'Could not find particle in its cell', i, ci
STOP 'Error in destroy_in_list'
ELSE ! Move on to the next
this = next ! Keep this index for next iteration
END IF
END DO ! End loop traversing link-list
END IF
END SUBROUTINE destroy_in_list
SUBROUTINE move_in_list ( i, ci ) ! Routine to move atom i from current cell to ci
IMPLICIT NONE
INTEGER, INTENT(in) :: i
INTEGER, DIMENSION(3), INTENT(in) :: ci
IF ( ALL ( ci(:) == c(:,i) ) ) RETURN ! No need to do anything
CALL destroy_in_list ( i, c(:,i) ) ! Remove atom i from old cell
CALL create_in_list ( i, ci(:) ) ! Add atom i to new cell
END SUBROUTINE move_in_list
SUBROUTINE check_list ( n, r ) ! Routine to check consistency of cell lists
IMPLICIT NONE
INTEGER, INTENT(in) :: n ! Number of atoms
REAL, DIMENSION(3,n), INTENT(in) :: r ! Atom positions
INTEGER, DIMENSION(3) :: ci
INTEGER :: ci1, ci2, ci3, i
DO i = 1, n ! Loop to check that each atom's cell is correct
ci = c_index ( r(:,i) ) ! Index function
IF ( ANY ( ci(:) /= c(:,i) ) ) THEN
WRITE ( unit=error_unit, fmt='(a,7i10)') 'Inconsistency 1 found:', i, ci, c(:,i)
STOP 'Error in check_list'
END IF
END DO ! End loop to check that each atom's cell is correct
! Triple loop over cells
DO ci1 = 0, sc-1
DO ci2 = 0, sc-1
DO ci3 = 0, sc-1
ci = [ ci1, ci2, ci3 ] ! Store 3D cell index in array
i = head(ci1,ci2,ci3) ! Locate head of list corresponding to cell
DO ! Loop traversing link list
IF ( i == 0 ) EXIT ! Reached end of list
IF ( ANY ( ci(:) /= c(:,i) ) ) THEN
WRITE ( unit=error_unit, fmt='(a,7i10)') 'Inconsistency 2 found:', i, ci, c(:,i)
STOP 'Error in check_list'
END IF
i = list(i) ! Move on to next list entry
END DO ! End loop traversing link list
END DO
END DO
END DO
! End triple loop over cells
END SUBROUTINE check_list
FUNCTION neighbours ( n, i, ci, half ) RESULT ( j_list )
IMPLICIT NONE
INTEGER, INTENT(in) :: n ! Number of atoms
INTEGER, INTENT(in) :: i ! Atom whose neighbours are required
INTEGER, DIMENSION(3), INTENT(in) :: ci ! Cell of atom of interest
LOGICAL, INTENT(in) :: half ! Determining the range of neighbours searched
INTEGER, DIMENSION(n) :: j_list ! Resulting list of indices
! This routine uses the link-list cell structure to fill out the array j_list
! with possible neighbours of atom i, padding with zeroes
! If half==.false., cell ci and all 26 surrounding cells are searched.
! If half==.true., cell ci, and just 13 of the neighbour cells, are searched
! and moreover, in ci, we only look down-list making use of list(i)
! There is a subtlety: using list(i) assumes that our interest is in the cells that
! are neighbours of c(:,i), i.e. that ci(:)==c(:,i), and we check for this explicitly.
! In other words, we assume that atom i has not moved since list was constructed.
! If half==.false., particle i might be in a very different position, and ci might be
! very different to c(:,i) but in that case we make no use of list(i), in normal use
! We have a cubic cell lattice
! Set up vectors to each cell in the 3x3x3 neighbourhood of a given cell
! To work properly, these are listed with inversion symmetry about (0,0,0)
INTEGER, PARAMETER :: nk = 13
INTEGER, DIMENSION(3,-nk:nk), PARAMETER :: d = RESHAPE( [ &
& -1,-1,-1, 0,-1,-1, 1,-1,-1, &
& -1, 1,-1, 0, 1,-1, 1, 1,-1, &
& -1, 0,-1, 1, 0,-1, 0, 0,-1, &
& 0,-1, 0, 1,-1, 0, -1,-1, 0, &
& -1, 0, 0, 0, 0, 0, 1, 0, 0, &
& 1, 1, 0, -1, 1, 0, 0, 1, 0, &
& 0, 0, 1, -1, 0, 1, 1, 0, 1, &
& -1,-1, 1, 0,-1, 1, 1,-1, 1, &
& -1, 1, 1, 0, 1, 1, 1, 1, 1 ], [ 3, 2*nk+1 ] )
INTEGER :: k1, k2, k, j, nj
INTEGER, DIMENSION(3) :: cj
IF ( half ) THEN ! Check half neighbour cells and j downlist from i in current cell
k1 = 0
k2 = nk
IF ( ANY ( ci(:) /= c(:,i) ) ) THEN ! should never happen
WRITE ( unit=error_unit, fmt='(a,6i15)' ) 'Cell mismatch ', ci(:), c(:,i)
STOP 'Error in get_neighbours'
END IF
ELSE ! Check every atom other than i in all cells
k1 = -nk
k2 = nk
END IF
j_list = 0 ! Initialize with zero values everywhere
nj = 0 ! Next position in list to be filled
DO k = k1, k2 ! Begin loop over neighbouring cells
cj(:) = ci(:) + d(:,k) ! Neighbour cell index
cj(:) = MODULO ( cj(:), sc ) ! Periodic boundary correction
IF ( k == 0 .AND. half ) THEN
j = list(i) ! Check down-list from i in i-cell
ELSE
j = head(cj(1),cj(2),cj(3)) ! Check entire j-cell
END IF
DO ! Begin loop over j atoms in list
IF ( j == 0 ) EXIT ! Exhausted list
IF ( j /= i ) THEN ! Skip self
nj = nj + 1 ! Increment count of j atoms
IF ( nj >= n ) THEN ! Check more than n-1 neighbours (should never happen)
WRITE ( unit=error_unit, fmt='(a,2i15)' ) 'Neighbour error for j_list', nj, n
STOP 'Impossible error in get_neighbours'
END IF ! End check more than n-1 neighbours
j_list(nj) = j ! Store new j atom
END IF
j = list(j) ! Next atom in j cell
END DO ! End loop over j atoms in list
END DO ! End loop over neighbouring cells
END FUNCTION neighbours
END MODULE link_list_module