UseTheForce/DarkSide/neighborLists.F90

!!
!! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
!!
!! The University of Notre Dame grants you ("Licensee") a
!! non-exclusive, royalty free, license to use, modify and
!! redistribute this software in source and binary code form, provided
!! that the following conditions are met:
!!
!! 1. Acknowledgement of the program authors must be made in any
!!    publication of scientific results based in part on use of the
!!    program.  An acceptable form of acknowledgement is citation of
!!    the article in which the program was described (Matthew
!!    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
!!    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
!!    Parallel Simulation Engine for Molecular Dynamics,"
!!    J. Comput. Chem. 26, pp. 252-271 (2005))
!!
!! 2. Redistributions of source code must retain the above copyright
!!    notice, this list of conditions and the following disclaimer.
!!
!! 3. Redistributions in binary form must reproduce the above copyright
!!    notice, this list of conditions and the following disclaimer in the
!!    documentation and/or other materials provided with the
!!    distribution.
!!
!! This software is provided "AS IS," without a warranty of any
!! kind. All express or implied conditions, representations and
!! warranties, including any implied warranty of merchantability,
!! fitness for a particular purpose or non-infringement, are hereby
!! excluded.  The University of Notre Dame and its licensors shall not
!! be liable for any damages suffered by licensee as a result of
!! using, modifying or distributing the software or its
!! derivatives. In no event will the University of Notre Dame or its
!! licensors be liable for any lost revenue, profit or data, or for
!! direct, indirect, special, consequential, incidental or punitive
!! damages, however caused and regardless of the theory of liability,
!! arising out of the use of or inability to use software, even if the
!! University of Notre Dame has been advised of the possibility of
!! such damages.
!!


!! Module neighborLists
!! Implements 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 2006-05-17 19:54:26 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 = 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(nGroups, error)
    integer, intent(out) :: error
    integer :: nGroups
    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(getNgroupsInRow(plan_group_row)+1),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 * getNgroupsInCol(plan_group_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(nGroups),stat=alloc_error)
       if (alloc_error /= 0) then
          write(default_error,*) &
               "ExpandNeighborLists: Error in allocating point"
          error = -1
          return
       end if
       allocate(list(listMultiplier * nGroups),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 * getNgroupsInCol(plan_group_col) + oldSize
#else
       newSize = listMultiplier * nGroups + 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(nGroups, q, listSkin, update_nlist)
    integer :: nGroups
    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

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


#ifdef MPI 

    dispmx_tmp = 0.0E0_DP
    do i = 1, nGroups
       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
#ifdef SINGLE_PRECISION
    call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_real, & 
         mpi_max,mpi_comm_world,mpi_err)
#else
    call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_double_precision, & 
         mpi_max,mpi_comm_world,mpi_err)
#endif

#else

    dispmx = 0.0_DP
    do i = 1, nGroups
       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(nGroups, q)

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

    !! get size of list
    list_size = nGroups

    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:	2758
Committed:	Wed May 17 19:54:27 2006 UTC (18 years, 3 months ago) by gezelter
File size:	9078 byte(s)
Log Message:	Adding single precision capabilities to the fortran side
#	User	Rev	Content
1	gezelter	1930	!!
2			!! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			!!
4			!! The University of Notre Dame grants you ("Licensee") a
5			!! non-exclusive, royalty free, license to use, modify and
6			!! redistribute this software in source and binary code form, provided
7			!! that the following conditions are met:
8			!!
9			!! 1. Acknowledgement of the program authors must be made in any
10			!! publication of scientific results based in part on use of the
11			!! program. An acceptable form of acknowledgement is citation of
12			!! the article in which the program was described (Matthew
13			!! A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14			!! J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15			!! Parallel Simulation Engine for Molecular Dynamics,"
16			!! J. Comput. Chem. 26, pp. 252-271 (2005))
17			!!
18			!! 2. Redistributions of source code must retain the above copyright
19			!! notice, this list of conditions and the following disclaimer.
20			!!
21			!! 3. Redistributions in binary form must reproduce the above copyright
22			!! notice, this list of conditions and the following disclaimer in the
23			!! documentation and/or other materials provided with the
24			!! distribution.
25			!!
26			!! This software is provided "AS IS," without a warranty of any
27			!! kind. All express or implied conditions, representations and
28			!! warranties, including any implied warranty of merchantability,
29			!! fitness for a particular purpose or non-infringement, are hereby
30			!! excluded. The University of Notre Dame and its licensors shall not
31			!! be liable for any damages suffered by licensee as a result of
32			!! using, modifying or distributing the software or its
33			!! derivatives. In no event will the University of Notre Dame or its
34			!! licensors be liable for any lost revenue, profit or data, or for
35			!! direct, indirect, special, consequential, incidental or punitive
36			!! damages, however caused and regardless of the theory of liability,
37			!! arising out of the use of or inability to use software, even if the
38			!! University of Notre Dame has been advised of the possibility of
39			!! such damages.
40			!!
41
42
43	gezelter	1490	!! Module neighborLists
44			!! Implements verlet neighbor lists for force modules.
45			!! Automagically expands neighbor list if size too small
46			!! up to maxAllocations times. If after maxAllocations we try to
47			!! expand the neighbor list, we get an error message and quit.
48			!! @author Charles F. Vardeman II
49			!! @author Matthew Meineke
50			!! @author J. Daniel Gezelter
51	gezelter	2758	!! @version $Id: neighborLists.F90,v 1.5 2006-05-17 19:54:26 gezelter Exp $,
52	gezelter	1490
53			module neighborLists
54	gezelter	2204
55	gezelter	1490	use definitions
56			#ifdef IS_MPI
57			use mpiSimulation
58			#endif
59	gezelter	2204
60	gezelter	1490	implicit none
61			PRIVATE
62	gezelter	2204
63	gezelter	1490	!--------------MODULE VARIABLES---------------------->
64			!! Parameter for size > # of long range particles neighbor list
65			!! should be.
66			integer, parameter :: listMultiplier = 80
67			!! Maximum number of times we should reallocate neighbor list.
68	chrisfen	2548	integer, parameter :: maxAllocations = 10
69	gezelter	1490	!! Number of times we have allocated the neighbor list.
70			integer, save :: nAllocations = 0
71			!! Pointer array to location in list for atom i.
72			integer, dimension(:),public, pointer :: point => null()
73			!! Neighbor list for atom i.
74			integer, dimension(:),public, pointer :: list => null()
75			!! Position array of previous positions for check. Allocated first time
76			!! into saveNeighborList.
77			real( kind = dp ), dimension(:,:), allocatable, save :: q0
78			!! Current list size
79			integer, save :: listSize
80			!--------------MODULE ACCESS-------------------------->
81			public :: expandNeighborList
82			public :: checkNeighborList
83			public :: saveNeighborList
84			public :: getNeighborListSize
85	gezelter	2204
86	gezelter	1490	contains
87
88
89			subroutine expandNeighborList(nGroups, error)
90			integer, intent(out) :: error
91			integer :: nGroups
92			integer :: alloc_error
93			integer :: oldSize = 0
94			integer :: newSize = 0
95			integer :: i
96			integer, dimension(:), pointer :: newList => null()
97			error = 0
98
99			!! First time through we should allocate point and list.
100			!! If one is associated and one is not, something is wrong
101			!! and return a error.
102
103			#ifdef IS_MPI !! // MPI
104			if (.not. associated(point) .and. &
105			.not. associated(list) ) then
106			allocate(point(getNgroupsInRow(plan_group_row)+1),stat=alloc_error)
107			if (alloc_error /= 0) then
108			write(default_error,*) &
109			"ExpandNeighborLists: Error in allocating MPI point"
110	gezelter	2204	error = -1
111	gezelter	1490	return
112			end if
113			allocate(list(listMultiplier * getNgroupsInCol(plan_group_col)),stat=alloc_error)
114			if (alloc_error /= 0) then
115			write(default_error,*) &
116			"ExpandNeighborLists: Error in allocating MPI list"
117			error = -1
118			return
119			end if
120			listSize = size(list)
121			nAllocations = nAllocations + 1
122			return
123			end if
124			#else !! // NONMPI
125			if (.not. associated(point) .and. &
126			.not. associated(list) ) then
127			allocate(point(nGroups),stat=alloc_error)
128			if (alloc_error /= 0) then
129			write(default_error,*) &
130			"ExpandNeighborLists: Error in allocating point"
131			error = -1
132			return
133			end if
134			allocate(list(listMultiplier * nGroups),stat=alloc_error)
135			if (alloc_error /= 0) then
136			write(default_error,*) &
137			"ExpandNeighborLists: Error in allocating list"
138			error = -1
139			return
140			end if
141			listSize = size(list)
142			nAllocations = nAllocations + 1
143			return
144			end if
145			#endif !! //MPI
146	gezelter	2204
147	gezelter	1490	! Expand the neighbor list
148	gezelter	2204
149	gezelter	1490	! Check to see if we have exceeded the maximum number of allocations.
150			if (nAllocations > maxAllocations) then
151			write(default_error,*) &
152			"ExpandNeighborList: exceeded maximum number of re-allocations"
153			error = -1
154			return
155			else !! Expand the list.
156			oldSize = size(list)
157
158
159			#ifdef IS_MPI !! MPI
160			newSize = listMultiplier * getNgroupsInCol(plan_group_col) + oldSize
161			#else
162			newSize = listMultiplier * nGroups + oldSize
163			#endif !! MPI
164
165
166	gezelter	2204
167	gezelter	1490	allocate(newList(newSize), stat=alloc_error)
168			if (alloc_error /= 0) then
169			write(,) "Error allocating new neighborlist"
170			error = -1
171			return
172			end if
173
174			!! Copy old list to new list
175			do i = 1, oldSize
176			newList(i) = list(i)
177			end do
178			!! Free old list
179			if(associated(list)) deallocate(list,stat=alloc_error)
180			if (alloc_error /= 0) then
181			error = -1
182			return
183			end if
184	gezelter	2204
185	gezelter	1490	!! Point list at new list
186	gezelter	2204
187	gezelter	1490	list => newList
188			end if
189
190			nAllocations = nAllocations + 1
191			listSize = size(list)
192			end subroutine expandNeighborList
193	gezelter	2204
194	gezelter	1490	!! checks to see if any long range particle has moved
195			!! through the neighbor list skin thickness.
196			subroutine checkNeighborList(nGroups, q, listSkin, update_nlist)
197			integer :: nGroups
198			real(kind = dp), intent(in) :: listSkin
199			real( kind = dp ), dimension(:,:) :: q
200			integer :: i
201			real( kind = DP ) :: dispmx
202			logical, intent(out) :: update_nlist
203			real( kind = DP ) :: dispmx_tmp
204
205			dispmx = 0.0E0_DP
206			!! calculate the largest displacement of any atom in any direction
207	gezelter	2204
208	gezelter	1490	!! If we have changed the particle idents, then we need to update
209			if (.not. allocated(q0)) then
210			update_nlist = .true.
211			return
212			end if
213
214			if (size(q0,2) /= nGroups) then
215			update_nlist = .true.
216			return
217			end if
218
219
220			#ifdef MPI
221	gezelter	2204
222	gezelter	1490	dispmx_tmp = 0.0E0_DP
223			do i = 1, nGroups
224			dispmx_tmp = max( abs ( q(1,i) - q0(1,i) ), dispmx_tmp )
225			dispmx_tmp = max( abs ( q(2,i) - q0(2,i) ), dispmx_tmp )
226			dispmx_tmp = max( abs ( q(3,i) - q0(3,i) ), dispmx_tmp )
227			end do
228	gezelter	2758	#ifdef SINGLE_PRECISION
229			call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_real, &
230			mpi_max,mpi_comm_world,mpi_err)
231			#else
232	gezelter	1490	call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_double_precision, &
233			mpi_max,mpi_comm_world,mpi_err)
234	gezelter	2758	#endif
235	gezelter	1490
236			#else
237	gezelter	2204
238	gezelter	1490	dispmx = 0.0_DP
239			do i = 1, nGroups
240			dispmx = max( abs ( q(1,i) - q0(1,i) ), dispmx )
241			dispmx = max( abs ( q(2,i) - q0(2,i) ), dispmx )
242			dispmx = max( abs ( q(3,i) - q0(3,i) ), dispmx )
243			end do
244
245			#endif
246
247			!! a conservative test of list skin crossings
248			dispmx = 2.0E0_DP * sqrt (3.0E0_DP * dispmx * dispmx)
249
250			update_nlist = (dispmx.gt.listSkin)
251	gezelter	2204
252			end subroutine checkNeighborList
253
254
255	gezelter	1490	!! Saves neighbor list for comparison in check.
256			!! Save_neighborList will work even if the number of
257			!! local atoms has changed.
258			subroutine saveNeighborList(nGroups, q)
259
260			integer :: nGroups
261			real(kind = dp ), dimension(3,nGroups), intent(in) :: q
262			integer :: list_size
263	gezelter	2204
264	gezelter	1490	!! get size of list
265			list_size = nGroups
266	gezelter	2204
267	gezelter	1490	if (.not. allocated(q0)) then
268			allocate(q0(3,list_size))
269			else if( list_size > size(q0,2)) then
270			deallocate(q0)
271			allocate(q0(3,list_size))
272			endif
273			q0 = q
274			end subroutine saveNeighborList
275	gezelter	2204
276
277	gezelter	1490	function getNeighborListSize() result(returnListSize)
278			integer :: returnListSize
279			returnListSize = listSize
280			end function getNeighborListSize
281	gezelter	2204
282	gezelter	1490	end module neighborLists