mdtools/libmdCode/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.3 2003-03-12 19:10:54 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 save_neighborList.
  real( kind = dp ), dimension(:,:), allocatable , save :: q0
  !! Current list size
  integer, save :: listSize
  !--------------MODULE ACCESS-------------------------->
  public :: expandNeighborList
  public :: checkNeighborList
  public :: save_neighborList
  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.
#ifndef IS_MPI !!/Non MPI
    if (.not. associated(point) .and. &
         .not. associated(list) ) then
       allocate(point(natoms),stat=alloc_error)
       if (alloc_error /= 0) then
          error = -1
          return
       end if
       allocate(list(listMultiplier * natoms),stat=alloc_error)
       if (alloc_error /= 0) then
          error = -1
          return
       end if
       nAllocations = nAllocations + 1
    else
       error = -1
       return
    end if
#else !!// MPI
    if (.not. associated(point) .and. &
         .not. associated(list) ) then
       allocate(point(getNRow(plan_row)),stat=alloc_error)
       if (alloc_error /= 0) then
          error = -1
          return
       end if
       allocate(list(listMultiplier * getNCol(plan_col)),stat=alloc_error)
       if (alloc_error /= 0) then
          error = -1
          return
       end if
       nAllocations = nAllocations + 1
       return
    else
       error = -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
       error = -1
       return
    else !! Expand the list.
       oldSize = size(list)
       
#ifndef IS_MPI !!Not MPI
       newSize = listMultiplier * natoms + oldSize
       allocate(newList(newSize), stat=alloc_error)
       if (alloc_error /= 0) then
          error = -1
          return
       end if
#else !! IS_MPI
       newSize = listMultiplier * getNCol(plan_col) + oldSize
       allocate(newList(newSize), stat = alloc_error)
       if (alloc_error /= 0) then
          error = -1
          return
       end if
#endif !! IS_MPI
       !! Copy old list to new list
       do i = 1, oldSize
          newList(i) = list(i)
       end do
       !! Free old 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
    
    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, rcut, rlist, update_nlist)
    integer :: natoms
    real(kind = dp), intent(in) :: rcut, rlist
    real( kind = dp ), dimension(:,:)  :: q 
    integer :: i
    real( kind = DP ) :: dispmx 
    logical, intent(out) :: update_nlist
    real( kind = DP ) :: dispmx_tmp
    real( kind = dp ) :: skin_thickness
    integer :: nlocal
    
    nlocal = natoms
    skin_thickness = rlist - rcut
    dispmx = 0.0E0_DP
    !! calculate the largest displacement of any atom in any direction
    
    
#ifdef MPI 
    
    !! If we have changed the particle idents, then we need to update    
    if (.not. allocated(q0) .or. &
         size(q0) /= nlocal) then
       update_nlist = .true.
       return
    end if
    
    dispmx_tmp = 0.0E0_DP
    do i = 1, nlocal
       dispmx_tmp = max( abs ( q(1,i) - q0(1,i) ), dispmx )
       dispmx_tmp = max( abs ( q(2,i) - q0(2,i) ), dispmx )
       dispmx_tmp = max( abs ( q(3,i) - q0(3,i) ), dispmx )
    end do
    call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_double_precision, & 
         mpi_max,mpi_comm_world,mpi_err)
#else
    
    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.(skin_thickness))
    
  end subroutine checkNeighborList
  
  
  !! Saves neighbor list for comparison in check.
  !! Save_neighborList will work even if the number of
  !! local atoms has changed.
  subroutine save_neighborList(q)
    real(kind = dp ), dimension(:,:), intent(in)  :: q
    integer :: list_size
    
    !! get size of list
    list_size = size(q)
    
    if (.not. allocated(q0)) then
       allocate(q0(3,list_size))
    else if( list_size > size(q0)) then
       deallocate(q0)
       allocate(q0(3,list_size))
    endif
    q0 = q
  end subroutine save_neighborList
  
  
  function getNeighborListSize() result(returnListSize)
    integer :: returnListSize
    returnListSize = listSize
  end function getNeighborListSize
  

end module neighborLists
Revision:	325
Committed:	Wed Mar 12 19:10:54 2003 UTC (21 years, 5 months ago) by gezelter
File size:	6434 byte(s)
Log Message:	Massive rewrite
#	User	Rev	Content
1	mmeineke	285	!! 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	325	!! @version $Id: neighborLists.F90,v 1.3 2003-03-12 19:10:54 gezelter Exp $,
10	mmeineke	285
11			module neighborLists
12	gezelter	325
13			use definitions
14	mmeineke	285	#ifdef IS_MPI
15			use mpiSimulation
16			#endif
17	gezelter	325
18	mmeineke	285	implicit none
19			PRIVATE
20	gezelter	325
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	mmeineke	285	integer, save :: nAllocations = 0
29	gezelter	325	!! Pointer array to location in list for atom i.
30	mmeineke	285	integer, dimension(:),public, pointer :: point => null()
31	gezelter	325	!! Neighbor list for atom i.
32	mmeineke	285	integer, dimension(:),public, pointer :: list => null()
33	gezelter	325	!! Position array of previous positions for check. Allocated first time
34			!! into save_neighborList.
35			real( kind = dp ), dimension(:,:), allocatable , save :: q0
36			!! Current list size
37	mmeineke	285	integer, save :: listSize
38	gezelter	325	!--------------MODULE ACCESS-------------------------->
39			public :: expandNeighborList
40			public :: checkNeighborList
41			public :: save_neighborList
42			public :: getNeighborListSize
43
44	mmeineke	285	contains
45
46
47	gezelter	325	subroutine expandNeighborList(natoms, error)
48	mmeineke	285	integer, intent(out) :: error
49	gezelter	325	integer :: natoms
50	mmeineke	285	integer :: alloc_error
51			integer :: oldSize = 0
52			integer :: newSize = 0
53	gezelter	325	integer :: i
54			integer, dimension(:), pointer :: newList => null()
55	mmeineke	285	error = 0
56	gezelter	325
57	mmeineke	285
58	gezelter	325	!! 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	mmeineke	285	#ifndef IS_MPI !!/Non MPI
62			if (.not. associated(point) .and. &
63			.not. associated(list) ) then
64	gezelter	325	allocate(point(natoms),stat=alloc_error)
65	mmeineke	285	if (alloc_error /= 0) then
66			error = -1
67			return
68			end if
69	gezelter	325	allocate(list(listMultiplier * natoms),stat=alloc_error)
70	mmeineke	285	if (alloc_error /= 0) then
71			error = -1
72			return
73			end if
74			nAllocations = nAllocations + 1
75			else
76			error = -1
77			return
78			end if
79			#else !!// MPI
80			if (.not. associated(point) .and. &
81			.not. associated(list) ) then
82	gezelter	325	allocate(point(getNRow(plan_row)),stat=alloc_error)
83	mmeineke	285	if (alloc_error /= 0) then
84			error = -1
85			return
86			end if
87	gezelter	325	allocate(list(listMultiplier * getNCol(plan_col)),stat=alloc_error)
88	mmeineke	285	if (alloc_error /= 0) then
89			error = -1
90			return
91			end if
92			nAllocations = nAllocations + 1
93			return
94			else
95			error = -1
96			return
97			end if
98			#endif !! //MPI
99	gezelter	325
100			! Expand the neighbor list
101
102			! Check to see if we have exceeded the maximum number of allocations.
103	mmeineke	285	if (nAllocations > maxAllocations) then
104			error = -1
105			return
106			else !! Expand the list.
107			oldSize = size(list)
108	gezelter	325
109	mmeineke	285	#ifndef IS_MPI !!Not MPI
110	gezelter	325	newSize = listMultiplier * natoms + oldSize
111	mmeineke	285	allocate(newList(newSize), stat=alloc_error)
112			if (alloc_error /= 0) then
113			error = -1
114			return
115			end if
116			#else !! IS_MPI
117	gezelter	325	newSize = listMultiplier * getNCol(plan_col) + oldSize
118	mmeineke	285	allocate(newList(newSize), stat = alloc_error)
119			if (alloc_error /= 0) then
120			error = -1
121			return
122			end if
123			#endif !! IS_MPI
124	gezelter	325	!! Copy old list to new list
125	mmeineke	285	do i = 1, oldSize
126			newList(i) = list(i)
127			end do
128	gezelter	325	!! Free old list
129	mmeineke	285	deallocate(list,stat=alloc_error)
130			if (alloc_error /= 0) then
131			error = -1
132			return
133			end if
134	gezelter	325
135			!! Point list at new list
136	mmeineke	285	list => newList
137			end if
138	gezelter	325
139	mmeineke	285	listSize = size(list)
140	gezelter	325	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	mmeineke	285	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	gezelter	325
155			nlocal = natoms
156			skin_thickness = rlist - rcut
157	mmeineke	285	dispmx = 0.0E0_DP
158			!! calculate the largest displacement of any atom in any direction
159	gezelter	325
160
161
162	mmeineke	285	#ifdef MPI
163	gezelter	325
164			!! If we have changed the particle idents, then we need to update
165	mmeineke	285	if (.not. allocated(q0) .or. &
166			size(q0) /= nlocal) then
167			update_nlist = .true.
168			return
169			end if
170	gezelter	325
171	mmeineke	285	dispmx_tmp = 0.0E0_DP
172			do i = 1, nlocal
173			dispmx_tmp = max( abs ( q(1,i) - q0(1,i) ), dispmx )
174			dispmx_tmp = max( abs ( q(2,i) - q0(2,i) ), dispmx )
175			dispmx_tmp = max( abs ( q(3,i) - q0(3,i) ), dispmx )
176			end do
177			call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_double_precision, &
178	gezelter	325	mpi_max,mpi_comm_world,mpi_err)
179	mmeineke	285	#else
180	gezelter	325
181	mmeineke	285	do i = 1, nlocal
182			dispmx = max( abs ( q(1,i) - q0(1,i) ), dispmx )
183			dispmx = max( abs ( q(2,i) - q0(2,i) ), dispmx )
184			dispmx = max( abs ( q(3,i) - q0(3,i) ), dispmx )
185			end do
186			#endif
187	gezelter	325
188	mmeineke	285	!! 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	gezelter	325	end subroutine checkNeighborList
194	mmeineke	285
195	gezelter	325
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	285	subroutine save_neighborList(q)
200			real(kind = dp ), dimension(:,:), intent(in) :: q
201			integer :: list_size
202	gezelter	325
203			!! get size of list
204	mmeineke	285	list_size = size(q)
205	gezelter	325
206	mmeineke	285	if (.not. allocated(q0)) then
207			allocate(q0(3,list_size))
208			else if( list_size > size(q0)) then
209			deallocate(q0)
210			allocate(q0(3,list_size))
211			endif
212			q0 = q
213			end subroutine save_neighborList
214
215	gezelter	325
216	chuckv	288	function getNeighborListSize() result(returnListSize)
217			integer :: returnListSize
218			returnListSize = listSize
219			end function getNeighborListSize
220	gezelter	325
221	chuckv	288
222	mmeineke	285	end module neighborLists