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.6 2003-03-20 19:50:42 chuckv 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 = 10
  !! 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
    
    write(*,*) "Inside expand nlist ", nAllocations
    !! 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)
       write(*,*) "Expanding list time, size = ", nAllocations,oldSize       
#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
       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
    
    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
    
    !! 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
    
#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.(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 saveNeighborList(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 saveNeighborList
  
  
  function getNeighborListSize() result(returnListSize)
    integer :: returnListSize
    returnListSize = listSize
  end function getNeighborListSize
    
end module neighborLists
Revision:	374
Committed:	Thu Mar 20 19:50:42 2003 UTC (21 years, 5 months ago) by chuckv
File size:	6588 byte(s)
Log Message:	* empty log message *
#	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	chuckv	374	!! @version $Id: neighborLists.F90,v 1.6 2003-03-20 19:50:42 chuckv Exp $,
10	mmeineke	285
11			module neighborLists
12	gezelter	334
13	gezelter	325	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	chuckv	374	integer, parameter :: maxAllocations = 10
27	gezelter	325	!! 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	gezelter	334	!! into saveNeighborList.
35	gezelter	325	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	gezelter	334	public :: saveNeighborList
42	gezelter	325	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	chuckv	374	write(,) "Inside expand nlist ", nAllocations
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	chuckv	374
62	mmeineke	285	#ifndef IS_MPI !!/Non MPI
63			if (.not. associated(point) .and. &
64			.not. associated(list) ) then
65	gezelter	325	allocate(point(natoms),stat=alloc_error)
66	mmeineke	285	if (alloc_error /= 0) then
67			error = -1
68			return
69			end if
70	gezelter	325	allocate(list(listMultiplier * natoms),stat=alloc_error)
71	mmeineke	285	if (alloc_error /= 0) then
72			error = -1
73			return
74			end if
75			nAllocations = nAllocations + 1
76			else
77			error = -1
78			return
79			end if
80			#else !!// MPI
81			if (.not. associated(point) .and. &
82			.not. associated(list) ) then
83	gezelter	325	allocate(point(getNRow(plan_row)),stat=alloc_error)
84	mmeineke	285	if (alloc_error /= 0) then
85			error = -1
86			return
87			end if
88	gezelter	325	allocate(list(listMultiplier * getNCol(plan_col)),stat=alloc_error)
89	mmeineke	285	if (alloc_error /= 0) then
90			error = -1
91			return
92			end if
93			nAllocations = nAllocations + 1
94			return
95			else
96			error = -1
97			return
98			end if
99			#endif !! //MPI
100	gezelter	325
101			! Expand the neighbor list
102
103			! Check to see if we have exceeded the maximum number of allocations.
104	mmeineke	285	if (nAllocations > maxAllocations) then
105			error = -1
106			return
107			else !! Expand the list.
108			oldSize = size(list)
109	chuckv	374	write(,) "Expanding list time, size = ", nAllocations,oldSize
110	mmeineke	285	#ifndef IS_MPI !!Not MPI
111	gezelter	325	newSize = listMultiplier * natoms + oldSize
112	mmeineke	285	allocate(newList(newSize), stat=alloc_error)
113			if (alloc_error /= 0) then
114			error = -1
115			return
116			end if
117			#else !! IS_MPI
118	gezelter	325	newSize = listMultiplier * getNCol(plan_col) + oldSize
119	mmeineke	285	allocate(newList(newSize), stat = alloc_error)
120			if (alloc_error /= 0) then
121			error = -1
122			return
123			end if
124			#endif !! IS_MPI
125	gezelter	325	!! Copy old list to new list
126	mmeineke	285	do i = 1, oldSize
127			newList(i) = list(i)
128			end do
129	gezelter	325	!! Free old list
130	chuckv	368	if(associated(list)) deallocate(list,stat=alloc_error)
131	mmeineke	285	if (alloc_error /= 0) then
132			error = -1
133			return
134			end if
135	gezelter	325
136			!! Point list at new list
137	mmeineke	285	list => newList
138			end if
139	gezelter	325
140	mmeineke	285	listSize = size(list)
141	gezelter	325	end subroutine expandNeighborList
142
143			!! checks to see if any long range particle has moved
144			!! through the neighbor list skin thickness.
145			subroutine checkNeighborList(natoms, q, rcut, rlist, update_nlist)
146			integer :: natoms
147			real(kind = dp), intent(in) :: rcut, rlist
148	mmeineke	285	real( kind = dp ), dimension(:,:) :: q
149			integer :: i
150			real( kind = DP ) :: dispmx
151			logical, intent(out) :: update_nlist
152			real( kind = DP ) :: dispmx_tmp
153			real( kind = dp ) :: skin_thickness
154			integer :: nlocal
155	gezelter	325
156			nlocal = natoms
157			skin_thickness = rlist - rcut
158	mmeineke	285	dispmx = 0.0E0_DP
159			!! calculate the largest displacement of any atom in any direction
160	gezelter	325
161			!! If we have changed the particle idents, then we need to update
162	chuckv	368	if (.not. allocated(q0) .or. size(q0) /= nlocal) then
163	mmeineke	285	update_nlist = .true.
164			return
165			end if
166	gezelter	325
167	chuckv	368	#ifdef MPI
168
169	mmeineke	285	dispmx_tmp = 0.0E0_DP
170			do i = 1, nlocal
171	chuckv	368	dispmx_tmp = max( abs ( q(1,i) - q0(1,i) ), dispmx_tmp )
172			dispmx_tmp = max( abs ( q(2,i) - q0(2,i) ), dispmx_tmp )
173			dispmx_tmp = max( abs ( q(3,i) - q0(3,i) ), dispmx_tmp )
174	mmeineke	285	end do
175			call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_double_precision, &
176	gezelter	325	mpi_max,mpi_comm_world,mpi_err)
177	chuckv	368
178	mmeineke	285	#else
179	gezelter	325
180	chuckv	368	dispmx = 0.0_DP
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	chuckv	368
187	mmeineke	285	#endif
188	gezelter	325
189	mmeineke	285	!! a conservative test of list skin crossings
190			dispmx = 2.0E0_DP * sqrt (3.0E0_DP * dispmx * dispmx)
191
192			update_nlist = (dispmx.gt.(skin_thickness))
193
194	gezelter	325	end subroutine checkNeighborList
195	mmeineke	285
196	gezelter	325
197			!! Saves neighbor list for comparison in check.
198			!! Save_neighborList will work even if the number of
199			!! local atoms has changed.
200	gezelter	334	subroutine saveNeighborList(q)
201	mmeineke	285	real(kind = dp ), dimension(:,:), intent(in) :: q
202			integer :: list_size
203	gezelter	325
204			!! get size of list
205	mmeineke	285	list_size = size(q)
206	gezelter	325
207	mmeineke	285	if (.not. allocated(q0)) then
208			allocate(q0(3,list_size))
209			else if( list_size > size(q0)) then
210			deallocate(q0)
211			allocate(q0(3,list_size))
212			endif
213			q0 = q
214	gezelter	334	end subroutine saveNeighborList
215	mmeineke	285
216	gezelter	325
217	chuckv	288	function getNeighborListSize() result(returnListSize)
218			integer :: returnListSize
219			returnListSize = listSize
220			end function getNeighborListSize
221	gezelter	334
222	mmeineke	285	end module neighborLists