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

Comparing trunk/OOPSE-2.0/src/UseTheForce/doForces.F90 (file contents):
Revision 2261 by gezelter, Tue Jun 28 13:58:45 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.22 2005-06-28 13:58:45 gezelter Exp $, $Date: 2005-06-28 13:58:45 $, $Name: not supported by cvs2svn $, $Revision: 1.22 $
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 82 | Line 82 | module doForces
82    logical, save :: haveNeighborList = .false.
83    logical, save :: haveSIMvariables = .false.
84    logical, save :: haveSaneForceField = .false.
85 +  logical, save :: haveInteractionMap = .false.
86  
87    logical, save :: FF_uses_DirectionalAtoms
88    logical, save :: FF_uses_LennardJones
# Line 115 | Line 116 | module doForces
116    logical, save :: SIM_uses_PBC
117    logical, save :: SIM_uses_molecular_cutoffs
118  
118 –  !!!GO AWAY---------
119 –  !!!!!real(kind=dp), save :: rlist, rlistsq
119  
120    public :: init_FF
121    public :: do_force_loop
122   !  public :: setRlistDF
123 +  !public :: addInteraction
124 +  !public :: setInteractionHash
125 +  !public :: getInteractionHash
126 +  public :: createInteractionMap
127 +  public :: createGroupCutoffs
128  
129   #ifdef PROFILE
130    public :: getforcetime
# Line 128 | Line 132 | module doForces
132    real :: forceTimeInitial, forceTimeFinal
133    integer :: nLoops
134   #endif
135 +  
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, public :: Interaction
133 <     integer :: InteractionHash
134 <     real(kind=dp) :: rCut
135 <  end type Interaction
155 >  type(gtype), dimension(:,:), allocatable :: gtypeCutoffMap
156    
137 –  type(Interaction), dimension(:,:),allocatable :: InteractionMap
138 –  
139 –  !public :: addInteraction
140 –  !public :: setInteractionHash
141 –  !public :: getInteractionHash
142 –  public :: createInteractionMap
143 –  
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 167 | 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 181 | 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 236 | 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 +  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 < !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
275 < !!$  subroutine setRlistDF( this_rlist )
276 < !!$
277 < !!$   real(kind=dp) :: this_rlist
278 < !!$
279 < !!$    rlist = this_rlist
280 < !!$    rlistsq = rlist * rlist
256 < !!$
257 < !!$    haveRlist = .true.
258 < !!$
259 < !!$  end subroutine setRlistDF
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 +    stat = 0
283 +    if (.not. haveInteractionMap) then
284  
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(:,:)%rCut = defaultRCut
298 +       InteractionMap(:,:)%rCutSq = defaultRCut*defaultRCut
299 +       InteractionMap(:,:)%rListSq = (defaultRCut+defaultSkinThickness)**2
300 +       haveRlist = .true.
301 +       return
302 +    end if
303 +
304 +    do map_i = 1,nAtypes
305 +       do map_j = map_i,nAtypes
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)
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)
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)
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)
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)
330 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
331 +           endif
332 +          
333 +           if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
334 + !              thisRCut = getGaybrneLJCutOff(map_i,map_j)
335 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
336 +           endif
337 +          
338 +           if ( iand(iMap, EAM_PAIR).ne.0 ) then      
339 + !              thisRCut = getEAMCutOff(map_i,map_j)
340 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
341 +           endif
342 +          
343 +           if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
344 + !              thisRCut = getShapeCutOff(map_i,map_j)
345 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
346 +           endif
347 +          
348 +           if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
349 + !              thisRCut = getShapeLJCutOff(map_i,map_j)
350 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
351 +           endif
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 +     ! now the groups
362 +
363 +
364 +
365 +     haveRlist = .true.
366 +   end subroutine createGroupCutoffs
367 +
368    subroutine setSimVariables()
369      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
370      SIM_uses_LennardJones = SimUsesLennardJones()
# Line 289 | Line 395 | contains
395      error = 0
396  
397      if (.not. haveInteractionMap) then
398 <
399 <       myStatus = 0
294 <
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 553 | Line 658 | contains
658      integer :: localError
659      integer :: propPack_i, propPack_j
660      integer :: loopStart, loopEnd, loop
661 <
661 >    integer :: iMap
662      real(kind=dp) :: listSkin = 1.0  
663  
664      !! initialize local variables  
# Line 644 | 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 672 | 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
794  
795 <             if (rgrpsq < rlistsq) then
795 >             if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
796                  if (update_nlist) then
797                     nlist = nlist + 1
798  
# Line 897 | 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  
# Line 1015 | Line 1128 | contains
1128      endif
1129      
1130      if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1131 <       call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1132 <            pot, A, f, t, do_pot)
1131 > !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1132 > !           pot, A, f, t, do_pot)
1133      endif
1134  
1135      if ( iand(iMap, EAM_PAIR).ne.0 ) then      

Diff Legend

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