ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/doForces.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2262 by chuckv, Sun Jul 3 20:53:43 2005 UTC vs.
Revision 2269 by chuckv, Tue Aug 9 19:40:56 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.23 2005-07-03 20:53:43 chuckv Exp $, $Date: 2005-07-03 20:53:43 $, $Name: not supported by cvs2svn $, $Revision: 1.23 $
48 > !! @version $Id: doForces.F90,v 1.27 2005-08-09 19:40:56 chuckv Exp $, $Date: 2005-08-09 19:40:56 $, $Name: not supported by cvs2svn $, $Revision: 1.27 $
49  
50  
51   module doForces
# Line 116 | Line 116 | module doForces
116    logical, save :: SIM_uses_PBC
117    logical, save :: SIM_uses_molecular_cutoffs
118  
119 –  !!!GO AWAY---------
120 –  !!!!!real(kind=dp), save :: rlist, rlistsq
119  
120    public :: init_FF
121    public :: do_force_loop
# Line 126 | Line 124 | module doForces
124    !public :: setInteractionHash
125    !public :: getInteractionHash
126    public :: createInteractionMap
127 <  public :: createRcuts
127 >  public :: createGroupCutoffs
128  
129   #ifdef PROFILE
130    public :: getforcetime
# Line 134 | Line 132 | module doForces
132    real :: forceTimeInitial, forceTimeFinal
133    integer :: nLoops
134   #endif
137 –
138 –  type, public :: Interaction
139 –     integer :: InteractionHash
140 –     real(kind=dp) :: rList = 0.0_dp
141 –     real(kind=dp) :: rListSq = 0.0_dp
142 –  end type Interaction
135    
136 <  type(Interaction), dimension(:,:),allocatable :: InteractionMap
137 <  
136 > !! Variables for cutoff mapping and interaction mapping
137 > ! Bit hash to determine pair-pair interactions.
138 >  integer, dimension(:,:),allocatable :: InteractionHash
139 > !! Cuttoffs in OOPSE are handled on a Group-Group pair basis.
140 > ! Largest cutoff for atypes for all potentials
141 >  real(kind=dp), dimension(:), allocatable :: atypeMaxCuttoff
142 > ! Largest cutoff for groups
143 >  real(kind=dp), dimension(:), allocatable :: groupMaxCutoff
144 > ! Group to Gtype transformation Map
145 >  integer,dimension(:), allocatable :: groupToGtype
146 > ! Group Type Max Cutoff
147 >  real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff
148 > ! GroupType definition
149 >  type ::gtype
150 >     real(kind=dp) :: rcut ! Group Cutoff
151 >     real(kind=dp) :: rcutsq ! Group Cutoff Squared
152 >     real(kind=dp) :: rlistsq ! List cutoff Squared    
153 >  end type gtype
154  
155 +  type(gtype), dimension(:,:), allocatable :: gtypeCutoffMap
156    
157   contains
158  
159  
160    subroutine createInteractionMap(status)
161      integer :: nAtypes
162 <    integer :: status
162 >    integer, intent(out) :: status
163      integer :: i
164      integer :: j
165      integer :: ihash
166      real(kind=dp) :: myRcut
167 < ! Test Types
167 >    !! Test Types
168      logical :: i_is_LJ
169      logical :: i_is_Elect
170      logical :: i_is_Sticky
# Line 171 | Line 180 | contains
180      logical :: j_is_EAM
181      logical :: j_is_Shape
182      
183 <    
183 >    status = 0  
184 >
185      if (.not. associated(atypes)) then
186 <       call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
186 >       call handleError("atype", "atypes was not present before call of createDefaultInteractionHash!")
187         status = -1
188         return
189      endif
# Line 185 | Line 195 | contains
195         return
196      end if
197  
198 <    if (.not. allocated(InteractionMap)) then
199 <       allocate(InteractionMap(nAtypes,nAtypes))
198 >    if (.not. allocated(InteractionHash)) then
199 >       allocate(InteractionHash(nAtypes,nAtypes))
200      endif
201          
202      do i = 1, nAtypes
# Line 240 | Line 250 | contains
250            if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
251  
252  
253 <          InteractionMap(i,j)%InteractionHash = iHash
254 <          InteractionMap(j,i)%InteractionHash = iHash
253 >          InteractionHash(i,j) = iHash
254 >          InteractionHash(j,i) = iHash
255  
256         end do
257  
258      end do
259 +
260 +    haveInteractionMap = .true.
261    end subroutine createInteractionMap
262  
263 < ! Query each potential and return the cutoff for that potential. We build the neighbor list based on the largest cutoff value for that atype. Each potential can decide whether to calculate the force for that atype based upon it's own cutoff.
252 <  subroutine createRcuts(defaultRList)
263 >  subroutine createGroupCutoffs(skinThickness,defaultrList,stat)
264      real(kind=dp), intent(in), optional :: defaultRList
265 +    real(kind-dp), intent(in), :: skinThickenss
266 +  ! Query each potential and return the cutoff for that potential. We
267 +  ! build the neighbor list based on the largest cutoff value for that
268 +  ! atype. Each potential can decide whether to calculate the force for
269 +  ! that atype based upon it's own cutoff.
270 +  
271 +
272 +    real(kind=dp), intent(in), optional :: defaultRCut, defaultSkinThickness
273 +
274      integer :: iMap
275      integer :: map_i,map_j
276      real(kind=dp) :: thisRCut = 0.0_dp
277      real(kind=dp) :: actualCutoff = 0.0_dp
278 +    integer, intent(out) :: stat
279      integer :: nAtypes
280 +    integer :: myStatus
281  
282 <    if(.not. allocated(InteractionMap)) return
282 >    stat = 0
283 >    if (.not. haveInteractionMap) then
284  
285 <    nAtypes = getSize(atypes)
286 < ! If we pass a default rcut, set all atypes to that cutoff distance
285 >       call createInteractionMap(myStatus)
286 >
287 >       if (myStatus .ne. 0) then
288 >          write(default_error, *) 'createInteractionMap failed in doForces!'
289 >          stat = -1
290 >          return
291 >       endif
292 >    endif
293 >
294 >    nAtypes = getSize(atypes)
295 >    !! If we pass a default rcut, set all atypes to that cutoff distance
296      if(present(defaultRList)) then
297 <       InteractionMap(:,:)%rList = defaultRList
298 <       InteractionMap(:,:)%rListSq = defaultRList*defaultRList
297 >       InteractionMap(:,:)%rCut = defaultRCut
298 >       InteractionMap(:,:)%rCutSq = defaultRCut*defaultRCut
299 >       InteractionMap(:,:)%rListSq = (defaultRCut+defaultSkinThickness)**2
300         haveRlist = .true.
301         return
302      end if
# Line 273 | Line 306 | contains
306            iMap = InteractionMap(map_i, map_j)%InteractionHash
307            
308            if ( iand(iMap, LJ_PAIR).ne.0 ) then
309 < !            thisRCut = getLJCutOff(map_i,map_j)
309 >             ! thisRCut = getLJCutOff(map_i,map_j)
310               if (thisRcut > actualCutoff) actualCutoff = thisRcut
311            endif
312            
313            if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
314 < !            thisRCut = getElectrostaticCutOff(map_i,map_j)
314 >             ! thisRCut = getElectrostaticCutOff(map_i,map_j)
315               if (thisRcut > actualCutoff) actualCutoff = thisRcut
316            endif
317            
318            if ( iand(iMap, STICKY_PAIR).ne.0 ) then
319 < !             thisRCut = getStickyCutOff(map_i,map_j)
319 >             ! thisRCut = getStickyCutOff(map_i,map_j)
320                if (thisRcut > actualCutoff) actualCutoff = thisRcut
321             endif
322            
323             if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
324 < !              thisRCut = getStickyPowerCutOff(map_i,map_j)
324 >              ! thisRCut = getStickyPowerCutOff(map_i,map_j)
325                if (thisRcut > actualCutoff) actualCutoff = thisRcut
326             endif
327            
328             if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
329 < !              thisRCut = getGayberneCutOff(map_i,map_j)
329 >              ! thisRCut = getGayberneCutOff(map_i,map_j)
330                if (thisRcut > actualCutoff) actualCutoff = thisRcut
331             endif
332            
# Line 316 | Line 349 | contains
349   !              thisRCut = getShapeLJCutOff(map_i,map_j)
350                if (thisRcut > actualCutoff) actualCutoff = thisRcut
351             endif
352 <           InteractionMap(map_i, map_j)%rList = actualCutoff
353 <           InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
352 >           InteractionMap(map_i, map_j)%rCut = actualCutoff
353 >           InteractionMap(map_i, map_j)%rCutSq = actualCutoff * actualCutoff
354 >           InteractionMap(map_i, map_j)%rListSq = (actualCutoff + skinThickness)**2
355 >
356 >           InteractionMap(map_j, map_i)%rCut = InteractionMap(map_i, map_j)%rCut
357 >           InteractionMap(map_j, map_i)%rCutSq = InteractionMap(map_i, map_j)%rCutSq
358 >           InteractionMap(map_j, map_i)%rListSq = InteractionMap(map_i, map_j)%rListSq
359          end do
360       end do
361 <          haveRlist = .true.
324 <  end subroutine createRcuts
361 >     ! now the groups
362  
363  
327 – !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
328 – !!$  subroutine setRlistDF( this_rlist )
329 – !!$
330 – !!$   real(kind=dp) :: this_rlist
331 – !!$
332 – !!$    rlist = this_rlist
333 – !!$    rlistsq = rlist * rlist
334 – !!$
335 – !!$    haveRlist = .true.
336 – !!$
337 – !!$  end subroutine setRlistDF
364  
365 +     haveRlist = .true.
366 +   end subroutine createGroupCutoffs
367  
368    subroutine setSimVariables()
369      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
# Line 367 | Line 395 | contains
395      error = 0
396  
397      if (.not. haveInteractionMap) then
398 <
399 <       myStatus = 0
372 <
398 >      
399 >       myStatus = 0      
400         call createInteractionMap(myStatus)
401 <
401 >      
402         if (myStatus .ne. 0) then
403            write(default_error, *) 'createInteractionMap failed in doForces!'
404            error = -1
# Line 722 | Line 749 | contains
749         iend = nGroups - 1
750   #endif
751         outer: do i = istart, iend
752 +
753 + #ifdef IS_MPI
754 +             me_i = atid_row(i)
755 + #else
756 +             me_i = atid(i)
757 + #endif
758  
759            if (update_nlist) point(i) = nlist + 1
760  
# Line 750 | Line 783 | contains
783               endif
784  
785   #ifdef IS_MPI
786 +             me_j = atid_col(j)
787               call get_interatomic_vector(q_group_Row(:,i), &
788                    q_group_Col(:,j), d_grp, rgrpsq)
789   #else
790 +             me_j = atid(j)
791               call get_interatomic_vector(q_group(:,i), &
792                    q_group(:,j), d_grp, rgrpsq)
793   #endif
# Line 975 | Line 1010 | contains
1010   #else
1011               me_i = atid(i)
1012   #endif
1013 <             iMap = InteractionMap(me_i, me_j)%InteractionHash
1013 >             iMap = InteractionHash(me_i,me_j)
1014              
1015               if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1016  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines