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.5 2003-03-20 00:02:39 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 = 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.
#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
       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:	368
Committed:	Thu Mar 20 00:02:39 2003 UTC (21 years, 5 months ago) by chuckv
File size:	6465 byte(s)
Log Message:	Fixed bugs. Single version now runs w/o segfault. Still a conservation of energy bug. do_Forces.F90: Fixed pot not being passed to do_pair. neighborLists.F90: Fixed bugs in checkNeighborLists atype_module.F90: Fixed bug with allocating atypes on each new_atype call.Now checks to see if atypes is null, then calles initialize(16). vector_class.F90: Fixed some bugs with how MatchList was being allocated.
#	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	368	!! @version $Id: neighborLists.F90,v 1.5 2003-03-20 00:02:39 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			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	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	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	chuckv	368	if(associated(list)) deallocate(list,stat=alloc_error)
130	mmeineke	285	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			!! If we have changed the particle idents, then we need to update
161	chuckv	368	if (.not. allocated(q0) .or. size(q0) /= nlocal) then
162	mmeineke	285	update_nlist = .true.
163			return
164			end if
165	gezelter	325
166	chuckv	368	#ifdef MPI
167
168	mmeineke	285	dispmx_tmp = 0.0E0_DP
169			do i = 1, nlocal
170	chuckv	368	dispmx_tmp = max( abs ( q(1,i) - q0(1,i) ), dispmx_tmp )
171			dispmx_tmp = max( abs ( q(2,i) - q0(2,i) ), dispmx_tmp )
172			dispmx_tmp = max( abs ( q(3,i) - q0(3,i) ), dispmx_tmp )
173	mmeineke	285	end do
174			call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_double_precision, &
175	gezelter	325	mpi_max,mpi_comm_world,mpi_err)
176	chuckv	368
177	mmeineke	285	#else
178	gezelter	325
179	chuckv	368	dispmx = 0.0_DP
180	mmeineke	285	do i = 1, nlocal
181			dispmx = max( abs ( q(1,i) - q0(1,i) ), dispmx )
182			dispmx = max( abs ( q(2,i) - q0(2,i) ), dispmx )
183			dispmx = max( abs ( q(3,i) - q0(3,i) ), dispmx )
184			end do
185	chuckv	368
186	mmeineke	285	#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	gezelter	334	subroutine saveNeighborList(q)
200	mmeineke	285	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	gezelter	334	end subroutine saveNeighborList
214	mmeineke	285
215	gezelter	325
216	chuckv	288	function getNeighborListSize() result(returnListSize)
217			integer :: returnListSize
218			returnListSize = listSize
219			end function getNeighborListSize
220	gezelter	334
221	mmeineke	285	end module neighborLists