ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/neighborLists.F90
Revision: 459
Committed: Fri Apr 4 19:57:01 2003 UTC (21 years, 3 months ago) by mmeineke
File size: 6438 byte(s)
Log Message:
fixed a memory read bug in neighborlist

File Contents

# Content
1 !! Module neighborLists
2 !! Impliments verlet neighbor lists for force modules.
3 !! Automagically expands neighbor list if size too small
4 !! up to maxAllocations times. If after maxAllocations we try to
5 !! expand the neighbor list, we get an error message and quit.
6 !! @author Charles F. Vardeman II
7 !! @author Matthew Meineke
8 !! @author J. Daniel Gezelter
9 !! @version $Id: neighborLists.F90,v 1.2 2003-04-04 19:57:01 mmeineke Exp $,
10
11 module neighborLists
12
13 use definitions
14 #ifdef IS_MPI
15 use mpiSimulation
16 #endif
17
18 implicit none
19 PRIVATE
20
21 !--------------MODULE VARIABLES---------------------->
22 !! Parameter for size > # of long range particles neighbor list
23 !! should be.
24 integer, parameter :: listMultiplier = 80
25 !! Maximum number of times we should reallocate neighbor list.
26 integer, parameter :: maxAllocations = 5
27 !! Number of times we have allocated the neighbor list.
28 integer, save :: nAllocations = 0
29 !! Pointer array to location in list for atom i.
30 integer, dimension(:),public, pointer :: point => null()
31 !! Neighbor list for atom i.
32 integer, dimension(:),public, pointer :: list => null()
33 !! Position array of previous positions for check. Allocated first time
34 !! into saveNeighborList.
35 real( kind = dp ), dimension(:,:), allocatable, save :: q0
36 !! Current list size
37 integer, save :: listSize
38 !--------------MODULE ACCESS-------------------------->
39 public :: expandNeighborList
40 public :: checkNeighborList
41 public :: saveNeighborList
42 public :: getNeighborListSize
43
44 contains
45
46
47 subroutine expandNeighborList(natoms, error)
48 integer, intent(out) :: error
49 integer :: natoms
50 integer :: alloc_error
51 integer :: oldSize = 0
52 integer :: newSize = 0
53 integer :: i
54 integer, dimension(:), pointer :: newList => null()
55 error = 0
56
57
58 !! First time through we should allocate point and list.
59 !! If one is associated and one is not, something is wrong
60 !! and return a error.
61
62 #ifndef IS_MPI !!/Non MPI
63 if (.not. associated(point) .and. &
64 .not. associated(list) ) then
65 allocate(point(natoms),stat=alloc_error)
66 if (alloc_error /= 0) then
67 error = -1
68 return
69 end if
70 allocate(list(listMultiplier * natoms),stat=alloc_error)
71 if (alloc_error /= 0) then
72 error = -1
73 return
74 end if
75 listSize = size(list)
76 nAllocations = nAllocations + 1
77 return
78 end if
79 #else !!// MPI
80 if (.not. associated(point) .and. &
81 .not. associated(list) ) then
82 allocate(point(getNRow(plan_row)),stat=alloc_error)
83 if (alloc_error /= 0) then
84 error = -1
85 return
86 end if
87 allocate(list(listMultiplier * getNCol(plan_col)),stat=alloc_error)
88 if (alloc_error /= 0) then
89 error = -1
90 return
91 end if
92 listSize = size(list)
93 nAllocations = nAllocations + 1
94 return
95 end if
96 #endif !! //MPI
97
98 ! Expand the neighbor list
99
100 ! Check to see if we have exceeded the maximum number of allocations.
101 if (nAllocations > maxAllocations) then
102 error = -1
103 return
104 else !! Expand the list.
105 oldSize = size(list)
106
107 #ifdef IS_MPI !! MPI
108 newSize = listMultiplier * getNCol(plan_col) + oldSize
109 #else
110 newSize = listMultiplier * natoms + oldSize
111 #endif !! MPI
112
113
114
115 allocate(newList(newSize), stat=alloc_error)
116 if (alloc_error /= 0) then
117 write(*,*) "Error allocating new neighborlist"
118 error = -1
119 return
120 end if
121
122 !! Copy old list to new list
123 do i = 1, oldSize
124 newList(i) = list(i)
125 end do
126 !! Free old list
127 if(associated(list)) deallocate(list,stat=alloc_error)
128 if (alloc_error /= 0) then
129 error = -1
130 return
131 end if
132
133 !! Point list at new list
134
135 list => newList
136 end if
137
138 nAllocations = nAllocations + 1
139 listSize = size(list)
140 end subroutine expandNeighborList
141
142 !! checks to see if any long range particle has moved
143 !! through the neighbor list skin thickness.
144 subroutine checkNeighborList(natoms, q, rcut, rlist, update_nlist)
145 integer :: natoms
146 real(kind = dp), intent(in) :: rcut, rlist
147 real( kind = dp ), dimension(:,:) :: q
148 integer :: i
149 real( kind = DP ) :: dispmx
150 logical, intent(out) :: update_nlist
151 real( kind = DP ) :: dispmx_tmp
152 real( kind = dp ) :: skin_thickness
153 integer :: nlocal
154
155 nlocal = natoms
156 skin_thickness = rlist - rcut
157 dispmx = 0.0E0_DP
158 !! calculate the largest displacement of any atom in any direction
159
160 !! If we have changed the particle idents, then we need to update
161 if (.not. allocated(q0) .or. size(q0) /= nlocal) then
162 update_nlist = .true.
163 return
164 end if
165
166 #ifdef MPI
167
168 dispmx_tmp = 0.0E0_DP
169 do i = 1, nlocal
170 dispmx_tmp = max( abs ( q(1,i) - q0(1,i) ), dispmx_tmp )
171 dispmx_tmp = max( abs ( q(2,i) - q0(2,i) ), dispmx_tmp )
172 dispmx_tmp = max( abs ( q(3,i) - q0(3,i) ), dispmx_tmp )
173 end do
174 call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_double_precision, &
175 mpi_max,mpi_comm_world,mpi_err)
176
177 #else
178
179 dispmx = 0.0_DP
180 do i = 1, nlocal
181 dispmx = max( abs ( q(1,i) - q0(1,i) ), dispmx )
182 dispmx = max( abs ( q(2,i) - q0(2,i) ), dispmx )
183 dispmx = max( abs ( q(3,i) - q0(3,i) ), dispmx )
184 end do
185
186 #endif
187
188 !! a conservative test of list skin crossings
189 dispmx = 2.0E0_DP * sqrt (3.0E0_DP * dispmx * dispmx)
190
191 update_nlist = (dispmx.gt.(skin_thickness))
192
193 end subroutine checkNeighborList
194
195
196 !! Saves neighbor list for comparison in check.
197 !! Save_neighborList will work even if the number of
198 !! local atoms has changed.
199 subroutine saveNeighborList(natoms, q)
200
201 integer :: natoms
202 real(kind = dp ), dimension(3,natoms), intent(in) :: q
203 integer :: list_size
204
205
206 !! get size of list
207 list_size = natoms
208
209 if (.not. allocated(q0)) then
210 allocate(q0(3,list_size))
211 else if( list_size > size(q0)) then
212 deallocate(q0)
213 allocate(q0(3,list_size))
214 endif
215 q0 = q
216 end subroutine saveNeighborList
217
218
219 function getNeighborListSize() result(returnListSize)
220 integer :: returnListSize
221 returnListSize = listSize
222 end function getNeighborListSize
223
224 end module neighborLists