OOPSE/libmdtools/neighborLists.F90

!! Module neighborLists
!! Impliments 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.7 2003-10-30 14:11:28 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(natoms, error)
    integer, intent(out) :: error
    integer :: natoms
    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(getNRow(plan_row)),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 * getNCol(plan_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(natoms),stat=alloc_error)
       if (alloc_error /= 0) then
          write(default_error,*) &
               "ExpandNeighborLists: Error in allocating point"
          error = -1
          return
       end if
       allocate(list(listMultiplier * natoms),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 * getNCol(plan_col) + oldSize
#else
       newSize = listMultiplier * natoms + 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(natoms, q, listSkin, update_nlist)
    integer :: natoms
    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
    integer :: nlocal


    nlocal = natoms
    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) /= nlocal) then
       update_nlist = .true.
       return
    end if


#ifdef MPI 
    
    dispmx_tmp = 0.0E0_DP
    do i = 1, nlocal
       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, nlocal
       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(natoms, q)

    integer :: natoms
    real(kind = dp ), dimension(3,natoms), intent(in)  :: q
    integer :: list_size


    !! get size of list
    list_size = natoms
    
    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:	844
Committed:	Thu Oct 30 14:11:28 2003 UTC (20 years, 8 months ago) by gezelter
File size:	6937 byte(s)
Log Message:	Fixed bug that size(q0) was being queried before q0 was allocated.
#	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	gezelter	844	!! @version $Id: neighborLists.F90,v 1.7 2003-10-30 14:11:28 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			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	gezelter	747	#ifdef IS_MPI !! // MPI
63	mmeineke	377	if (.not. associated(point) .and. &
64			.not. associated(list) ) then
65	gezelter	747	allocate(point(getNRow(plan_row)),stat=alloc_error)
66	mmeineke	377	if (alloc_error /= 0) then
67	chuckv	480	write(default_error,*) &
68	gezelter	747	"ExpandNeighborLists: Error in allocating MPI point"
69			error = -1
70	mmeineke	377	return
71			end if
72	gezelter	747	allocate(list(listMultiplier * getNCol(plan_col)),stat=alloc_error)
73	mmeineke	377	if (alloc_error /= 0) then
74	chuckv	480	write(default_error,*) &
75	gezelter	747	"ExpandNeighborLists: Error in allocating MPI list"
76	mmeineke	377	error = -1
77			return
78			end if
79			listSize = size(list)
80			nAllocations = nAllocations + 1
81			return
82			end if
83	gezelter	747	#else !! // NONMPI
84	mmeineke	377	if (.not. associated(point) .and. &
85			.not. associated(list) ) then
86	gezelter	747	allocate(point(natoms),stat=alloc_error)
87	mmeineke	377	if (alloc_error /= 0) then
88	chuckv	480	write(default_error,*) &
89	gezelter	747	"ExpandNeighborLists: Error in allocating point"
90			error = -1
91	mmeineke	377	return
92			end if
93	gezelter	747	allocate(list(listMultiplier * natoms),stat=alloc_error)
94	mmeineke	377	if (alloc_error /= 0) then
95	chuckv	480	write(default_error,*) &
96	gezelter	747	"ExpandNeighborLists: Error in allocating list"
97	mmeineke	377	error = -1
98			return
99			end if
100			listSize = size(list)
101			nAllocations = nAllocations + 1
102			return
103			end if
104			#endif !! //MPI
105
106			! Expand the neighbor list
107
108			! Check to see if we have exceeded the maximum number of allocations.
109			if (nAllocations > maxAllocations) then
110	chuckv	480	write(default_error,*) &
111			"ExpandNeighborList: exceeded maximum number of re-allocations"
112	mmeineke	377	error = -1
113			return
114			else !! Expand the list.
115			oldSize = size(list)
116
117	chuckv	480
118	mmeineke	377	#ifdef IS_MPI !! MPI
119			newSize = listMultiplier * getNCol(plan_col) + oldSize
120			#else
121			newSize = listMultiplier * natoms + oldSize
122			#endif !! MPI
123
124
125
126			allocate(newList(newSize), stat=alloc_error)
127			if (alloc_error /= 0) then
128			write(,) "Error allocating new neighborlist"
129			error = -1
130			return
131			end if
132
133			!! Copy old list to new list
134			do i = 1, oldSize
135			newList(i) = list(i)
136			end do
137			!! Free old list
138			if(associated(list)) deallocate(list,stat=alloc_error)
139			if (alloc_error /= 0) then
140			error = -1
141			return
142			end if
143
144			!! Point list at new list
145
146			list => newList
147			end if
148
149			nAllocations = nAllocations + 1
150			listSize = size(list)
151			end subroutine expandNeighborList
152
153			!! checks to see if any long range particle has moved
154			!! through the neighbor list skin thickness.
155	mmeineke	626	subroutine checkNeighborList(natoms, q, listSkin, update_nlist)
156	mmeineke	377	integer :: natoms
157	mmeineke	626	real(kind = dp), intent(in) :: listSkin
158	mmeineke	377	real( kind = dp ), dimension(:,:) :: q
159			integer :: i
160			real( kind = DP ) :: dispmx
161			logical, intent(out) :: update_nlist
162			real( kind = DP ) :: dispmx_tmp
163			integer :: nlocal
164	chuckv	480
165
166	mmeineke	377	nlocal = natoms
167			dispmx = 0.0E0_DP
168			!! calculate the largest displacement of any atom in any direction
169
170			!! If we have changed the particle idents, then we need to update
171	gezelter	844	if (.not. allocated(q0)) then
172	mmeineke	377	update_nlist = .true.
173			return
174			end if
175	chuckv	480
176	gezelter	844	if (size(q0,2) /= nlocal) then
177			update_nlist = .true.
178			return
179			end if
180	chuckv	480
181	gezelter	844
182	mmeineke	377	#ifdef MPI
183
184			dispmx_tmp = 0.0E0_DP
185			do i = 1, nlocal
186			dispmx_tmp = max( abs ( q(1,i) - q0(1,i) ), dispmx_tmp )
187			dispmx_tmp = max( abs ( q(2,i) - q0(2,i) ), dispmx_tmp )
188			dispmx_tmp = max( abs ( q(3,i) - q0(3,i) ), dispmx_tmp )
189			end do
190			call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_double_precision, &
191			mpi_max,mpi_comm_world,mpi_err)
192
193			#else
194
195			dispmx = 0.0_DP
196			do i = 1, nlocal
197			dispmx = max( abs ( q(1,i) - q0(1,i) ), dispmx )
198			dispmx = max( abs ( q(2,i) - q0(2,i) ), dispmx )
199			dispmx = max( abs ( q(3,i) - q0(3,i) ), dispmx )
200			end do
201
202			#endif
203
204	chuckv	673
205	mmeineke	377	!! a conservative test of list skin crossings
206			dispmx = 2.0E0_DP * sqrt (3.0E0_DP * dispmx * dispmx)
207
208	mmeineke	626	update_nlist = (dispmx.gt.listSkin)
209	mmeineke	377
210	chuckv	673	end subroutine checkNeighborList
211	mmeineke	377
212
213			!! Saves neighbor list for comparison in check.
214			!! Save_neighborList will work even if the number of
215			!! local atoms has changed.
216	mmeineke	459	subroutine saveNeighborList(natoms, q)
217
218			integer :: natoms
219			real(kind = dp ), dimension(3,natoms), intent(in) :: q
220	mmeineke	377	integer :: list_size
221	mmeineke	459
222	chuckv	673
223	mmeineke	377
224			!! get size of list
225	mmeineke	459	list_size = natoms
226	mmeineke	377
227			if (.not. allocated(q0)) then
228			allocate(q0(3,list_size))
229	chuckv	673	else if( list_size > size(q0,2)) then
230	mmeineke	377	deallocate(q0)
231			allocate(q0(3,list_size))
232			endif
233			q0 = q
234			end subroutine saveNeighborList
235
236
237			function getNeighborListSize() result(returnListSize)
238			integer :: returnListSize
239			returnListSize = listSize
240			end function getNeighborListSize
241
242			end module neighborLists