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

# User Rev Content
1 mmeineke 377 !! 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 mmeineke 459 !! @version $Id: neighborLists.F90,v 1.2 2003-04-04 19:57:01 mmeineke Exp $,
10 mmeineke 377
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 mmeineke 459 real( kind = dp ), dimension(:,:), allocatable, save :: q0
36 mmeineke 377 !! 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 mmeineke 459 subroutine saveNeighborList(natoms, q)
200    
201     integer :: natoms
202     real(kind = dp ), dimension(3,natoms), intent(in) :: q
203 mmeineke 377 integer :: list_size
204 mmeineke 459
205 mmeineke 377
206     !! get size of list
207 mmeineke 459 list_size = natoms
208 mmeineke 377
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