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.9 2003-11-21 19:31:05 chrisfen 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:	872
Committed:	Fri Nov 21 19:31:05 2003 UTC (20 years, 7 months ago) by chrisfen
File size:	6926 byte(s)
Log Message:	Fixed a bug in SimInfo ordering of radii
#	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	chrisfen	872	!! @version $Id: neighborLists.F90,v 1.9 2003-11-21 19:31:05 chrisfen 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	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	747	allocate(point(getNRow(plan_row)),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	gezelter	747	allocate(list(listMultiplier * getNCol(plan_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	gezelter	747	allocate(point(natoms),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	gezelter	747	allocate(list(listMultiplier * natoms),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			newSize = listMultiplier * getNCol(plan_col) + oldSize
119			#else
120			newSize = listMultiplier * natoms + oldSize
121			#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	mmeineke	626	subroutine checkNeighborList(natoms, q, listSkin, update_nlist)
155	mmeineke	377	integer :: natoms
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			integer :: nlocal
163	chuckv	480
164	mmeineke	377	nlocal = natoms
165			dispmx = 0.0E0_DP
166			!! calculate the largest displacement of any atom in any direction
167
168			!! If we have changed the particle idents, then we need to update
169	gezelter	844	if (.not. allocated(q0)) then
170	mmeineke	377	update_nlist = .true.
171			return
172			end if
173	chuckv	480
174	gezelter	844	if (size(q0,2) /= nlocal) then
175			update_nlist = .true.
176			return
177			end if
178	chuckv	480
179	gezelter	844
180	mmeineke	377	#ifdef MPI
181
182			dispmx_tmp = 0.0E0_DP
183			do i = 1, nlocal
184			dispmx_tmp = max( abs ( q(1,i) - q0(1,i) ), dispmx_tmp )
185			dispmx_tmp = max( abs ( q(2,i) - q0(2,i) ), dispmx_tmp )
186			dispmx_tmp = max( abs ( q(3,i) - q0(3,i) ), dispmx_tmp )
187			end do
188			call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_double_precision, &
189			mpi_max,mpi_comm_world,mpi_err)
190
191			#else
192
193			dispmx = 0.0_DP
194			do i = 1, nlocal
195			dispmx = max( abs ( q(1,i) - q0(1,i) ), dispmx )
196			dispmx = max( abs ( q(2,i) - q0(2,i) ), dispmx )
197			dispmx = max( abs ( q(3,i) - q0(3,i) ), dispmx )
198			end do
199
200	gezelter	845	#endif
201	chuckv	673
202	mmeineke	377	!! a conservative test of list skin crossings
203	gezelter	845	dispmx = 2.0E0_DP * sqrt (3.0E0_DP * dispmx * dispmx)
204
205	mmeineke	626	update_nlist = (dispmx.gt.listSkin)
206	chrisfen	872
207	chuckv	673	end subroutine checkNeighborList
208	mmeineke	377
209
210			!! Saves neighbor list for comparison in check.
211			!! Save_neighborList will work even if the number of
212			!! local atoms has changed.
213	mmeineke	459	subroutine saveNeighborList(natoms, q)
214
215			integer :: natoms
216			real(kind = dp ), dimension(3,natoms), intent(in) :: q
217	mmeineke	377	integer :: list_size
218
219			!! get size of list
220	mmeineke	459	list_size = natoms
221	mmeineke	377
222			if (.not. allocated(q0)) then
223			allocate(q0(3,list_size))
224	chuckv	673	else if( list_size > size(q0,2)) then
225	mmeineke	377	deallocate(q0)
226			allocate(q0(3,list_size))
227			endif
228			q0 = q
229			end subroutine saveNeighborList
230
231
232			function getNeighborListSize() result(returnListSize)
233			integer :: returnListSize
234			returnListSize = listSize
235			end function getNeighborListSize
236
237			end module neighborLists