OOPSE/libmdtools/neighborLists.F90

!! Module neighborLists
!! Implements verlet neighbor lists for force modules.
!! Automagically expands neighbor list if size too small
!! up to maxAllocations times. If after maxAllocations we try to
!! expand the neighbor list, we get an error message and quit.
!! @author Charles F. Vardeman II
!! @author Matthew Meineke
!! @author J. Daniel Gezelter
!! @version $Id: neighborLists.F90,v 1.12 2004-06-01 21:45:22 gezelter Exp $, 

module neighborLists
  
  use definitions
#ifdef IS_MPI
  use mpiSimulation
#endif
  
  implicit none
  PRIVATE
  
  !--------------MODULE VARIABLES---------------------->
  !! Parameter for size > # of long range particles neighbor list
  !! should be.
  integer, parameter :: listMultiplier = 80
  !! Maximum number of times we should reallocate neighbor list.
  integer, parameter :: maxAllocations = 5
  !! Number of times we have allocated the neighbor list.
  integer, save       :: nAllocations = 0
  !! Pointer array to location in list for atom i.
  integer, dimension(:),public, pointer :: point => null()
  !! Neighbor list for atom i.
  integer, dimension(:),public, pointer :: list  => null()
  !! Position array of previous positions for check. Allocated first time
  !! into saveNeighborList.
  real( kind = dp ), dimension(:,:), allocatable, save  :: q0
  !! Current list size
  integer, save :: listSize
  !--------------MODULE ACCESS-------------------------->
  public :: expandNeighborList
  public :: checkNeighborList
  public :: saveNeighborList
  public :: getNeighborListSize
  
contains


  subroutine expandNeighborList(nGroups, error)
    integer, intent(out) :: error
    integer :: nGroups
    integer :: alloc_error
    integer :: oldSize = 0
    integer :: newSize = 0
    integer :: i
    integer, dimension(:), pointer :: newList => null()
    error = 0   

    !! First time through we should allocate point and list.
    !! If one is associated and one is not, something is wrong
    !! and return a error.

#ifdef IS_MPI !! // MPI
    if (.not. associated(point) .and. &
         .not. associated(list) ) then
       allocate(point(getNgroupsInRow(plan_group_row)+1),stat=alloc_error)
       if (alloc_error /= 0) then
          write(default_error,*) &
               "ExpandNeighborLists: Error in allocating MPI point"
           error = -1
          return
       end if
       allocate(list(listMultiplier * getNgroupsInCol(plan_group_col)),stat=alloc_error)
       if (alloc_error /= 0) then
          write(default_error,*) &
               "ExpandNeighborLists: Error in allocating MPI list"
          error = -1
          return
       end if
       listSize = size(list)
       nAllocations = nAllocations + 1
       return
    end if
#else !! // NONMPI
    if (.not. associated(point) .and. &
         .not. associated(list) ) then
       allocate(point(nGroups),stat=alloc_error)
       if (alloc_error /= 0) then
          write(default_error,*) &
               "ExpandNeighborLists: Error in allocating point"
          error = -1
          return
       end if
       allocate(list(listMultiplier * nGroups),stat=alloc_error)
       if (alloc_error /= 0) then
          write(default_error,*) &
               "ExpandNeighborLists: Error in allocating list"
          error = -1
          return
       end if
       listSize = size(list)
       nAllocations = nAllocations + 1
       return
    end if
#endif !! //MPI
    
    ! Expand the neighbor list
    
    ! Check to see if we have exceeded the maximum number of allocations.
    if (nAllocations > maxAllocations) then
       write(default_error,*) & 
            "ExpandNeighborList: exceeded maximum number of re-allocations"
       error = -1
       return
    else !! Expand the list.
       oldSize = size(list)


#ifdef IS_MPI !! MPI
       newSize = listMultiplier * getNgroupsInCol(plan_group_col) + oldSize
#else
       newSize = listMultiplier * nGroups + oldSize
#endif !! MPI

       
       allocate(newList(newSize), stat=alloc_error)
       if (alloc_error /= 0) then
          write(*,*) "Error allocating new neighborlist"
          error = -1
          return
       end if

       !! Copy old list to new list
       do i = 1, oldSize
          newList(i) = list(i)
       end do
       !! Free old list
       if(associated(list)) deallocate(list,stat=alloc_error)
       if (alloc_error /= 0) then
          error = -1
          return
       end if
       
       !! Point list at new list
       
       list => newList
    end if

    nAllocations = nAllocations + 1
    listSize = size(list)
  end subroutine expandNeighborList
  
  !! checks to see if any long range particle has moved
  !! through the neighbor list skin thickness.
  subroutine checkNeighborList(nGroups, q, listSkin, update_nlist)
    integer :: nGroups
    real(kind = dp), intent(in) :: listSkin
    real( kind = dp ), dimension(:,:)  :: q 
    integer :: i
    real( kind = DP ) :: dispmx 
    logical, intent(out) :: update_nlist
    real( kind = DP ) :: dispmx_tmp

    dispmx = 0.0E0_DP
    !! calculate the largest displacement of any atom in any direction
    
    !! If we have changed the particle idents, then we need to update    
    if (.not. allocated(q0)) then       
       update_nlist = .true.
       return
    end if

    if (size(q0,2) /= nGroups) then
       update_nlist = .true.
       return
    end if


#ifdef MPI 
    
    dispmx_tmp = 0.0E0_DP
    do i = 1, nGroups
       dispmx_tmp = max( abs ( q(1,i) - q0(1,i) ), dispmx_tmp )
       dispmx_tmp = max( abs ( q(2,i) - q0(2,i) ), dispmx_tmp )
       dispmx_tmp = max( abs ( q(3,i) - q0(3,i) ), dispmx_tmp )
    end do
    call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_double_precision, & 
         mpi_max,mpi_comm_world,mpi_err)

#else
    
    dispmx = 0.0_DP
    do i = 1, nGroups
       dispmx = max( abs ( q(1,i) - q0(1,i) ), dispmx )
       dispmx = max( abs ( q(2,i) - q0(2,i) ), dispmx )
       dispmx = max( abs ( q(3,i) - q0(3,i) ), dispmx )
    end do

#endif   

    !! a conservative test of list skin crossings
    dispmx = 2.0E0_DP * sqrt (3.0E0_DP * dispmx * dispmx)   

    update_nlist = (dispmx.gt.listSkin)
 
   end subroutine checkNeighborList
  
  
  !! Saves neighbor list for comparison in check.
  !! Save_neighborList will work even if the number of
  !! local atoms has changed.
  subroutine saveNeighborList(nGroups, q)

    integer :: nGroups
    real(kind = dp ), dimension(3,nGroups), intent(in)  :: q
    integer :: list_size
    
    !! get size of list
    list_size = nGroups
    
    if (.not. allocated(q0)) then
       allocate(q0(3,list_size))
    else if( list_size > size(q0,2)) then
       deallocate(q0)
       allocate(q0(3,list_size))
    endif
    q0 = q
  end subroutine saveNeighborList
  
  
  function getNeighborListSize() result(returnListSize)
    integer :: returnListSize
    returnListSize = listSize
  end function getNeighborListSize
    
end module neighborLists
Revision:	1217
Committed:	Tue Jun 1 21:45:22 2004 UTC (21 years, 5 months ago) by gezelter
File size:	6943 byte(s)
Log Message:	Bug fix (fixes of skipList and neighbor list under MPI)
#	User	Rev	Content
1	mmeineke	377	!! Module neighborLists
2	gezelter	1150	!! Implements verlet neighbor lists for force modules.
3	mmeineke	377	!! 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	gezelter	1217	!! @version $Id: neighborLists.F90,v 1.12 2004-06-01 21:45:22 gezelter 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	tim	1198	subroutine expandNeighborList(nGroups, error)
48	mmeineke	377	integer, intent(out) :: error
49	tim	1198	integer :: nGroups
50	mmeineke	377	integer :: alloc_error
51			integer :: oldSize = 0
52			integer :: newSize = 0
53			integer :: i
54			integer, dimension(:), pointer :: newList => null()
55	gezelter	845	error = 0
56	mmeineke	377
57			!! First time through we should allocate point and list.
58			!! If one is associated and one is not, something is wrong
59			!! and return a error.
60
61	gezelter	747	#ifdef IS_MPI !! // MPI
62	mmeineke	377	if (.not. associated(point) .and. &
63			.not. associated(list) ) then
64	gezelter	1217	allocate(point(getNgroupsInRow(plan_group_row)+1),stat=alloc_error)
65	mmeineke	377	if (alloc_error /= 0) then
66	chuckv	480	write(default_error,*) &
67	gezelter	747	"ExpandNeighborLists: Error in allocating MPI point"
68			error = -1
69	mmeineke	377	return
70			end if
71	tim	1198	allocate(list(listMultiplier * getNgroupsInCol(plan_group_col)),stat=alloc_error)
72	mmeineke	377	if (alloc_error /= 0) then
73	chuckv	480	write(default_error,*) &
74	gezelter	747	"ExpandNeighborLists: Error in allocating MPI list"
75	mmeineke	377	error = -1
76			return
77			end if
78			listSize = size(list)
79			nAllocations = nAllocations + 1
80			return
81			end if
82	gezelter	747	#else !! // NONMPI
83	mmeineke	377	if (.not. associated(point) .and. &
84			.not. associated(list) ) then
85	tim	1198	allocate(point(nGroups),stat=alloc_error)
86	mmeineke	377	if (alloc_error /= 0) then
87	chuckv	480	write(default_error,*) &
88	gezelter	747	"ExpandNeighborLists: Error in allocating point"
89			error = -1
90	mmeineke	377	return
91			end if
92	tim	1198	allocate(list(listMultiplier * nGroups),stat=alloc_error)
93	mmeineke	377	if (alloc_error /= 0) then
94	chuckv	480	write(default_error,*) &
95	gezelter	747	"ExpandNeighborLists: Error in allocating list"
96	mmeineke	377	error = -1
97			return
98			end if
99			listSize = size(list)
100			nAllocations = nAllocations + 1
101			return
102			end if
103			#endif !! //MPI
104
105			! Expand the neighbor list
106
107			! Check to see if we have exceeded the maximum number of allocations.
108			if (nAllocations > maxAllocations) then
109	chuckv	480	write(default_error,*) &
110			"ExpandNeighborList: exceeded maximum number of re-allocations"
111	mmeineke	377	error = -1
112			return
113			else !! Expand the list.
114			oldSize = size(list)
115
116	chuckv	480
117	mmeineke	377	#ifdef IS_MPI !! MPI
118	tim	1198	newSize = listMultiplier * getNgroupsInCol(plan_group_col) + oldSize
119	mmeineke	377	#else
120	tim	1198	newSize = listMultiplier * nGroups + oldSize
121	mmeineke	377	#endif !! MPI
122
123
124
125			allocate(newList(newSize), stat=alloc_error)
126			if (alloc_error /= 0) then
127			write(,) "Error allocating new neighborlist"
128			error = -1
129			return
130			end if
131
132			!! Copy old list to new list
133			do i = 1, oldSize
134			newList(i) = list(i)
135			end do
136			!! Free old list
137			if(associated(list)) deallocate(list,stat=alloc_error)
138			if (alloc_error /= 0) then
139			error = -1
140			return
141			end if
142
143			!! Point list at new list
144
145			list => newList
146			end if
147
148			nAllocations = nAllocations + 1
149			listSize = size(list)
150			end subroutine expandNeighborList
151
152			!! checks to see if any long range particle has moved
153			!! through the neighbor list skin thickness.
154	tim	1198	subroutine checkNeighborList(nGroups, q, listSkin, update_nlist)
155			integer :: nGroups
156	mmeineke	626	real(kind = dp), intent(in) :: listSkin
157	mmeineke	377	real( kind = dp ), dimension(:,:) :: q
158			integer :: i
159			real( kind = DP ) :: dispmx
160			logical, intent(out) :: update_nlist
161			real( kind = DP ) :: dispmx_tmp
162	chuckv	480
163	mmeineke	377	dispmx = 0.0E0_DP
164			!! calculate the largest displacement of any atom in any direction
165
166			!! If we have changed the particle idents, then we need to update
167	gezelter	844	if (.not. allocated(q0)) then
168	mmeineke	377	update_nlist = .true.
169			return
170			end if
171	chuckv	480
172	tim	1198	if (size(q0,2) /= nGroups) then
173	gezelter	844	update_nlist = .true.
174			return
175			end if
176	chuckv	480
177	gezelter	844
178	mmeineke	377	#ifdef MPI
179
180			dispmx_tmp = 0.0E0_DP
181	tim	1198	do i = 1, nGroups
182	mmeineke	377	dispmx_tmp = max( abs ( q(1,i) - q0(1,i) ), dispmx_tmp )
183			dispmx_tmp = max( abs ( q(2,i) - q0(2,i) ), dispmx_tmp )
184			dispmx_tmp = max( abs ( q(3,i) - q0(3,i) ), dispmx_tmp )
185			end do
186			call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_double_precision, &
187			mpi_max,mpi_comm_world,mpi_err)
188
189			#else
190
191			dispmx = 0.0_DP
192	tim	1198	do i = 1, nGroups
193	mmeineke	377	dispmx = max( abs ( q(1,i) - q0(1,i) ), dispmx )
194			dispmx = max( abs ( q(2,i) - q0(2,i) ), dispmx )
195			dispmx = max( abs ( q(3,i) - q0(3,i) ), dispmx )
196			end do
197
198	gezelter	845	#endif
199	chuckv	673
200	mmeineke	377	!! a conservative test of list skin crossings
201	gezelter	845	dispmx = 2.0E0_DP * sqrt (3.0E0_DP * dispmx * dispmx)
202
203	mmeineke	626	update_nlist = (dispmx.gt.listSkin)
204	chrisfen	872
205	chuckv	673	end subroutine checkNeighborList
206	mmeineke	377
207
208			!! Saves neighbor list for comparison in check.
209			!! Save_neighborList will work even if the number of
210			!! local atoms has changed.
211	tim	1198	subroutine saveNeighborList(nGroups, q)
212	mmeineke	459
213	tim	1198	integer :: nGroups
214			real(kind = dp ), dimension(3,nGroups), intent(in) :: q
215	mmeineke	377	integer :: list_size
216
217			!! get size of list
218	tim	1198	list_size = nGroups
219	mmeineke	377
220			if (.not. allocated(q0)) then
221			allocate(q0(3,list_size))
222	chuckv	673	else if( list_size > size(q0,2)) then
223	mmeineke	377	deallocate(q0)
224			allocate(q0(3,list_size))
225			endif
226			q0 = q
227			end subroutine saveNeighborList
228
229
230			function getNeighborListSize() result(returnListSize)
231			integer :: returnListSize
232			returnListSize = listSize
233			end function getNeighborListSize
234
235			end module neighborLists