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

Comparing trunk/OOPSE-3.0/src/UseTheForce/doForces.F90 (file contents):
Revision 2267 by tim, Fri Jul 29 17:34:06 2005 UTC vs.
Revision 2350 by chuckv, Mon Oct 10 21:20:46 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.25 2005-07-29 17:34:06 tim Exp $, $Date: 2005-07-29 17:34:06 $, $Name: not supported by cvs2svn $, $Revision: 1.25 $
48 > !! @version $Id: doForces.F90,v 1.51 2005-10-10 21:20:46 chuckv Exp $, $Date: 2005-10-10 21:20:46 $, $Name: not supported by cvs2svn $, $Revision: 1.51 $
49  
50  
51   module doForces
# Line 58 | Line 58 | module doForces
58    use lj
59    use sticky
60    use electrostatic_module
61 <  use reaction_field
61 >  use reaction_field_module
62    use gb_pair
63    use shapes
64    use vector_class
# Line 73 | Line 73 | module doForces
73  
74   #define __FORTRAN90
75   #include "UseTheForce/fSwitchingFunction.h"
76 + #include "UseTheForce/fCutoffPolicy.h"
77   #include "UseTheForce/DarkSide/fInteractionMap.h"
78 + #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
79  
80 +
81    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
82    INTEGER, PARAMETER:: PAIR_LOOP    = 2
83  
81  logical, save :: haveRlist = .false.
84    logical, save :: haveNeighborList = .false.
85    logical, save :: haveSIMvariables = .false.
86    logical, save :: haveSaneForceField = .false.
87 <  logical, save :: haveInteractionMap = .false.
87 >  logical, save :: haveInteractionHash = .false.
88 >  logical, save :: haveGtypeCutoffMap = .false.
89 >  logical, save :: haveDefaultCutoffs = .false.
90 >  logical, save :: haveRlist = .false.
91  
92    logical, save :: FF_uses_DirectionalAtoms
88  logical, save :: FF_uses_LennardJones
89  logical, save :: FF_uses_Electrostatics
90  logical, save :: FF_uses_Charges
93    logical, save :: FF_uses_Dipoles
92  logical, save :: FF_uses_Quadrupoles
93  logical, save :: FF_uses_Sticky
94  logical, save :: FF_uses_StickyPower
94    logical, save :: FF_uses_GayBerne
95    logical, save :: FF_uses_EAM
97  logical, save :: FF_uses_Shapes
98  logical, save :: FF_uses_FLARB
99  logical, save :: FF_uses_RF
96  
97    logical, save :: SIM_uses_DirectionalAtoms
102  logical, save :: SIM_uses_LennardJones
103  logical, save :: SIM_uses_Electrostatics
104  logical, save :: SIM_uses_Charges
105  logical, save :: SIM_uses_Dipoles
106  logical, save :: SIM_uses_Quadrupoles
107  logical, save :: SIM_uses_Sticky
108  logical, save :: SIM_uses_StickyPower
109  logical, save :: SIM_uses_GayBerne
98    logical, save :: SIM_uses_EAM
111  logical, save :: SIM_uses_Shapes
112  logical, save :: SIM_uses_FLARB
113  logical, save :: SIM_uses_RF
99    logical, save :: SIM_requires_postpair_calc
100    logical, save :: SIM_requires_prepair_calc
101    logical, save :: SIM_uses_PBC
117  logical, save :: SIM_uses_molecular_cutoffs
102  
103 <  !!!GO AWAY---------
120 <  !!!!!real(kind=dp), save :: rlist, rlistsq
103 >  integer, save :: electrostaticSummationMethod
104  
105    public :: init_FF
106 +  public :: setDefaultCutoffs
107    public :: do_force_loop
108 < !  public :: setRlistDF
109 <  !public :: addInteraction
110 <  !public :: setInteractionHash
111 <  !public :: getInteractionHash
112 <  public :: createInteractionMap
113 <  public :: createRcuts
108 >  public :: createInteractionHash
109 >  public :: createGtypeCutoffMap
110 >  public :: getStickyCut
111 >  public :: getStickyPowerCut
112 >  public :: getGayBerneCut
113 >  public :: getEAMCut
114 >  public :: getShapeCut
115  
116   #ifdef PROFILE
117    public :: getforcetime
# Line 134 | Line 119 | module doForces
119    real :: forceTimeInitial, forceTimeFinal
120    integer :: nLoops
121   #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
122    
123 <  type(Interaction), dimension(:,:),allocatable :: InteractionMap
124 <  
123 >  !! Variables for cutoff mapping and interaction mapping
124 >  ! Bit hash to determine pair-pair interactions.
125 >  integer, dimension(:,:), allocatable :: InteractionHash
126 >  real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
127 >  real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
128 >  real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
129  
130 +  integer, dimension(:), allocatable, target :: groupToGtypeRow
131 +  integer, dimension(:), pointer :: groupToGtypeCol => null()
132 +
133 +  real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
134 +  real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
135 +  type ::gtypeCutoffs
136 +     real(kind=dp) :: rcut
137 +     real(kind=dp) :: rcutsq
138 +     real(kind=dp) :: rlistsq
139 +  end type gtypeCutoffs
140 +  type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
141 +
142 +  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
143 +  real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
144 +  real(kind=dp),save :: listSkin
145    
146   contains
147  
148 <
151 <  subroutine createInteractionMap(status)
148 >  subroutine createInteractionHash(status)
149      integer :: nAtypes
150      integer, intent(out) :: status
151      integer :: i
152      integer :: j
153 <    integer :: ihash
157 <    real(kind=dp) :: myRcut
153 >    integer :: iHash
154      !! Test Types
155      logical :: i_is_LJ
156      logical :: i_is_Elect
# Line 170 | Line 166 | contains
166      logical :: j_is_GB
167      logical :: j_is_EAM
168      logical :: j_is_Shape
169 <    
170 <    status = 0
171 <    
169 >    real(kind=dp) :: myRcut
170 >
171 >    status = 0  
172 >
173      if (.not. associated(atypes)) then
174 <       call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
174 >       call handleError("atype", "atypes was not present before call of createInteractionHash!")
175         status = -1
176         return
177      endif
# Line 186 | Line 183 | contains
183         return
184      end if
185  
186 <    if (.not. allocated(InteractionMap)) then
187 <       allocate(InteractionMap(nAtypes,nAtypes))
186 >    if (.not. allocated(InteractionHash)) then
187 >       allocate(InteractionHash(nAtypes,nAtypes))
188      endif
189 +
190 +    if (.not. allocated(atypeMaxCutoff)) then
191 +       allocate(atypeMaxCutoff(nAtypes))
192 +    endif
193          
194      do i = 1, nAtypes
195         call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
# Line 241 | Line 242 | contains
242            if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
243  
244  
245 <          InteractionMap(i,j)%InteractionHash = iHash
246 <          InteractionMap(j,i)%InteractionHash = iHash
245 >          InteractionHash(i,j) = iHash
246 >          InteractionHash(j,i) = iHash
247  
248         end do
249  
250      end do
251  
252 <    haveInteractionMap = .true.
253 <  end subroutine createInteractionMap
252 >    haveInteractionHash = .true.
253 >  end subroutine createInteractionHash
254  
255 < ! 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.
255 <  subroutine createRcuts(defaultRList,stat)
256 <    real(kind=dp), intent(in), optional :: defaultRList
257 <    integer :: iMap
258 <    integer :: map_i,map_j
259 <    real(kind=dp) :: thisRCut = 0.0_dp
260 <    real(kind=dp) :: actualCutoff = 0.0_dp
261 <    integer, intent(out) :: stat
262 <    integer :: nAtypes
263 <    integer :: myStatus
255 >  subroutine createGtypeCutoffMap(stat)
256  
257 <    stat = 0
258 <    if (.not. haveInteractionMap) then
257 >    integer, intent(out), optional :: stat
258 >    logical :: i_is_LJ
259 >    logical :: i_is_Elect
260 >    logical :: i_is_Sticky
261 >    logical :: i_is_StickyP
262 >    logical :: i_is_GB
263 >    logical :: i_is_EAM
264 >    logical :: i_is_Shape
265 >    logical :: GtypeFound
266  
267 <       call createInteractionMap(myStatus)
267 >    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
268 >    integer :: n_in_i, me_i, ia, g, atom1
269 >    integer :: nGroupsInRow
270 >    integer :: nGroupsInCol
271 >    integer :: nGroupTypesRow,nGroupTypesCol
272 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin
273 >    real(kind=dp) :: biggestAtypeCutoff
274  
275 +    stat = 0
276 +    if (.not. haveInteractionHash) then
277 +       call createInteractionHash(myStatus)      
278         if (myStatus .ne. 0) then
279 <          write(default_error, *) 'createInteractionMap failed in doForces!'
279 >          write(default_error, *) 'createInteractionHash failed in doForces!'
280            stat = -1
281            return
282         endif
283      endif
284 <
285 <
284 > #ifdef IS_MPI
285 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
286 >    nGroupsInCol = getNgroupsInCol(plan_group_col)
287 > #endif
288      nAtypes = getSize(atypes)
289 <    !! If we pass a default rcut, set all atypes to that cutoff distance
290 <    if(present(defaultRList)) then
291 <       InteractionMap(:,:)%rList = defaultRList
292 <       InteractionMap(:,:)%rListSq = defaultRList*defaultRList
293 <       haveRlist = .true.
294 <       return
295 <    end if
296 <
297 <    do map_i = 1,nAtypes
298 <       do map_j = map_i,nAtypes
299 <          iMap = InteractionMap(map_i, map_j)%InteractionHash
289 > ! Set all of the initial cutoffs to zero.
290 >    atypeMaxCutoff = 0.0_dp
291 >    do i = 1, nAtypes
292 >       if (SimHasAtype(i)) then    
293 >          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
294 >          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
295 >          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
296 >          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
297 >          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
298 >          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
299 >          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
300            
301 <          if ( iand(iMap, LJ_PAIR).ne.0 ) then
302 <             ! thisRCut = getLJCutOff(map_i,map_j)
303 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
301 >
302 >          if (haveDefaultCutoffs) then
303 >             atypeMaxCutoff(i) = defaultRcut
304 >          else
305 >             if (i_is_LJ) then          
306 >                thisRcut = getSigma(i) * 2.5_dp
307 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
308 >             endif
309 >             if (i_is_Elect) then
310 >                thisRcut = defaultRcut
311 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
312 >             endif
313 >             if (i_is_Sticky) then
314 >                thisRcut = getStickyCut(i)
315 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
316 >             endif
317 >             if (i_is_StickyP) then
318 >                thisRcut = getStickyPowerCut(i)
319 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
320 >             endif
321 >             if (i_is_GB) then
322 >                thisRcut = getGayBerneCut(i)
323 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
324 >             endif
325 >             if (i_is_EAM) then
326 >                thisRcut = getEAMCut(i)
327 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
328 >             endif
329 >             if (i_is_Shape) then
330 >                thisRcut = getShapeCut(i)
331 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
332 >             endif
333            endif
334            
296          if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
297             ! thisRCut = getElectrostaticCutOff(map_i,map_j)
298             if (thisRcut > actualCutoff) actualCutoff = thisRcut
299          endif
335            
336 <          if ( iand(iMap, STICKY_PAIR).ne.0 ) then
337 <             ! thisRCut = getStickyCutOff(map_i,map_j)
338 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
304 <           endif
305 <          
306 <           if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
307 <              ! thisRCut = getStickyPowerCutOff(map_i,map_j)
308 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
309 <           endif
310 <          
311 <           if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
312 <              ! thisRCut = getGayberneCutOff(map_i,map_j)
313 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
314 <           endif
315 <          
316 <           if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
317 < !              thisRCut = getGaybrneLJCutOff(map_i,map_j)
318 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
319 <           endif
320 <          
321 <           if ( iand(iMap, EAM_PAIR).ne.0 ) then      
322 < !              thisRCut = getEAMCutOff(map_i,map_j)
323 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
324 <           endif
325 <          
326 <           if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
327 < !              thisRCut = getShapeCutOff(map_i,map_j)
328 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
329 <           endif
330 <          
331 <           if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
332 < !              thisRCut = getShapeLJCutOff(map_i,map_j)
333 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
334 <           endif
335 <           InteractionMap(map_i, map_j)%rList = actualCutoff
336 <           InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
337 <        end do
338 <     end do
339 <     haveRlist = .true.
340 <  end subroutine createRcuts
336 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
337 >             biggestAtypeCutoff = atypeMaxCutoff(i)
338 >          endif
339  
340 +       endif
341 +    enddo
342 +  
343  
344 < !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
345 < !!$  subroutine setRlistDF( this_rlist )
346 < !!$
347 < !!$   real(kind=dp) :: this_rlist
348 < !!$
349 < !!$    rlist = this_rlist
350 < !!$    rlistsq = rlist * rlist
351 < !!$
352 < !!$    haveRlist = .true.
353 < !!$
354 < !!$  end subroutine setRlistDF
344 >    
345 >    istart = 1
346 >    jstart = 1
347 > #ifdef IS_MPI
348 >    iend = nGroupsInRow
349 >    jend = nGroupsInCol
350 > #else
351 >    iend = nGroups
352 >    jend = nGroups
353 > #endif
354 >    
355 >    !! allocate the groupToGtype and gtypeMaxCutoff here.
356 >    if(.not.allocated(groupToGtypeRow)) then
357 >     !  allocate(groupToGtype(iend))
358 >       allocate(groupToGtypeRow(iend))
359 >    else
360 >       deallocate(groupToGtypeRow)
361 >       allocate(groupToGtypeRow(iend))
362 >    endif
363 >    if(.not.allocated(groupMaxCutoffRow)) then
364 >       allocate(groupMaxCutoffRow(iend))
365 >    else
366 >       deallocate(groupMaxCutoffRow)
367 >       allocate(groupMaxCutoffRow(iend))
368 >    end if
369 >
370 >    if(.not.allocated(gtypeMaxCutoffRow)) then
371 >       allocate(gtypeMaxCutoffRow(iend))
372 >    else
373 >       deallocate(gtypeMaxCutoffRow)
374 >       allocate(gtypeMaxCutoffRow(iend))
375 >    endif
376 >
377 >
378 > #ifdef IS_MPI
379 >       ! We only allocate new storage if we are in MPI because Ncol /= Nrow
380 >    if(.not.allocated(groupToGtypeCol)) then
381 >       allocate(groupToGtypeCol(jend))
382 >    else
383 >       deallocate(groupToGtypeCol)
384 >       allocate(groupToGtypeCol(jend))
385 >    end if
386 >
387 >    if(.not.allocated(groupToGtypeCol)) then
388 >       allocate(groupToGtypeCol(jend))
389 >    else
390 >       deallocate(groupToGtypeCol)
391 >       allocate(groupToGtypeCol(jend))
392 >    end if
393 >    if(.not.allocated(gtypeMaxCutoffCol)) then
394 >       allocate(gtypeMaxCutoffCol(jend))
395 >    else
396 >       deallocate(gtypeMaxCutoffCol)      
397 >       allocate(gtypeMaxCutoffCol(jend))
398 >    end if
399 >
400 >       groupMaxCutoffCol = 0.0_dp
401 >       gtypeMaxCutoffCol = 0.0_dp
402 >
403 > #endif
404 >       groupMaxCutoffRow = 0.0_dp
405 >       gtypeMaxCutoffRow = 0.0_dp
406 >
407 >
408 >    !! first we do a single loop over the cutoff groups to find the
409 >    !! largest cutoff for any atypes present in this group.  We also
410 >    !! create gtypes at this point.
411 >    
412 >    tol = 1.0d-6
413 >    nGroupTypesRow = 0
414 >
415 >    do i = istart, iend      
416 >       n_in_i = groupStartRow(i+1) - groupStartRow(i)
417 >       groupMaxCutoffRow(i) = 0.0_dp
418 >       do ia = groupStartRow(i), groupStartRow(i+1)-1
419 >          atom1 = groupListRow(ia)
420 > #ifdef IS_MPI
421 >          me_i = atid_row(atom1)
422 > #else
423 >          me_i = atid(atom1)
424 > #endif          
425 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
426 >             groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
427 >          endif          
428 >       enddo
429  
430 +       if (nGroupTypesRow.eq.0) then
431 +          nGroupTypesRow = nGroupTypesRow + 1
432 +          gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
433 +          groupToGtypeRow(i) = nGroupTypesRow
434 +       else
435 +          GtypeFound = .false.
436 +          do g = 1, nGroupTypesRow
437 +             if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
438 +                groupToGtypeRow(i) = g
439 +                GtypeFound = .true.
440 +             endif
441 +          enddo
442 +          if (.not.GtypeFound) then            
443 +             nGroupTypesRow = nGroupTypesRow + 1
444 +             gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
445 +             groupToGtypeRow(i) = nGroupTypesRow
446 +          endif
447 +       endif
448 +    enddo    
449  
450 + #ifdef IS_MPI
451 +    do j = jstart, jend      
452 +       n_in_j = groupStartCol(j+1) - groupStartCol(j)
453 +       groupMaxCutoffCol(j) = 0.0_dp
454 +       do ja = groupStartCol(j), groupStartCol(j+1)-1
455 +          atom1 = groupListCol(ja)
456 +
457 +          me_j = atid_col(atom1)
458 +
459 +          if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
460 +             groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
461 +          endif          
462 +       enddo
463 +
464 +       if (nGroupTypesCol.eq.0) then
465 +          nGroupTypesCol = nGroupTypesCol + 1
466 +          gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
467 +          groupToGtypeCol(j) = nGroupTypesCol
468 +       else
469 +          GtypeFound = .false.
470 +          do g = 1, nGroupTypesCol
471 +             if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
472 +                groupToGtypeCol(j) = g
473 +                GtypeFound = .true.
474 +             endif
475 +          enddo
476 +          if (.not.GtypeFound) then            
477 +             nGroupTypesCol = nGroupTypesCol + 1
478 +             gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
479 +             groupToGtypeCol(j) = nGroupTypesCol
480 +          endif
481 +       endif
482 +    enddo    
483 +
484 + #else
485 + ! Set pointers to information we just found
486 +    nGroupTypesCol = nGroupTypesRow
487 +    groupToGtypeCol => groupToGtypeRow
488 +    gtypeMaxCutoffCol => gtypeMaxCutoffRow
489 +    groupMaxCutoffCol => groupMaxCutoffRow
490 + #endif
491 +
492 +
493 +
494 +
495 +
496 +    !! allocate the gtypeCutoffMap here.
497 +    allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
498 +    !! then we do a double loop over all the group TYPES to find the cutoff
499 +    !! map between groups of two types
500 +    tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
501 +
502 +    do i = 1, nGroupTypesRow
503 +       do j = 1, nGroupTypesCol
504 +      
505 +          select case(cutoffPolicy)
506 +          case(TRADITIONAL_CUTOFF_POLICY)
507 +             thisRcut = tradRcut
508 +          case(MIX_CUTOFF_POLICY)
509 +             thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
510 +          case(MAX_CUTOFF_POLICY)
511 +             thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
512 +          case default
513 +             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
514 +             return
515 +          end select
516 +          gtypeCutoffMap(i,j)%rcut = thisRcut
517 +          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
518 +          skin = defaultRlist - defaultRcut
519 +          listSkin = skin ! set neighbor list skin thickness
520 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
521 +
522 +          ! sanity check
523 +
524 +          if (haveDefaultCutoffs) then
525 +             if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
526 +                call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
527 +             endif
528 +          endif
529 +       enddo
530 +    enddo
531 +    if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
532 +    if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
533 +    if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
534 + #ifdef IS_MPI
535 +    if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
536 +    if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
537 + #endif
538 +    groupMaxCutoffCol => null()
539 +    gtypeMaxCutoffCol => null()
540 +    
541 +    haveGtypeCutoffMap = .true.
542 +   end subroutine createGtypeCutoffMap
543 +
544 +   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
545 +     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
546 +     integer, intent(in) :: cutPolicy
547 +
548 +     defaultRcut = defRcut
549 +     defaultRsw = defRsw
550 +     defaultRlist = defRlist
551 +     cutoffPolicy = cutPolicy
552 +
553 +     haveDefaultCutoffs = .true.
554 +   end subroutine setDefaultCutoffs
555 +
556 +   subroutine setCutoffPolicy(cutPolicy)
557 +
558 +     integer, intent(in) :: cutPolicy
559 +     cutoffPolicy = cutPolicy
560 +     call createGtypeCutoffMap()
561 +   end subroutine setCutoffPolicy
562 +    
563 +    
564    subroutine setSimVariables()
565      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
358    SIM_uses_LennardJones = SimUsesLennardJones()
359    SIM_uses_Electrostatics = SimUsesElectrostatics()
360    SIM_uses_Charges = SimUsesCharges()
361    SIM_uses_Dipoles = SimUsesDipoles()
362    SIM_uses_Sticky = SimUsesSticky()
363    SIM_uses_StickyPower = SimUsesStickyPower()
364    SIM_uses_GayBerne = SimUsesGayBerne()
566      SIM_uses_EAM = SimUsesEAM()
366    SIM_uses_Shapes = SimUsesShapes()
367    SIM_uses_FLARB = SimUsesFLARB()
368    SIM_uses_RF = SimUsesRF()
567      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
568      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
569      SIM_uses_PBC = SimUsesPBC()
# Line 382 | Line 580 | contains
580  
581      error = 0
582  
583 <    if (.not. haveInteractionMap) then
583 >    if (.not. haveInteractionHash) then      
584 >       myStatus = 0      
585 >       call createInteractionHash(myStatus)      
586 >       if (myStatus .ne. 0) then
587 >          write(default_error, *) 'createInteractionHash failed in doForces!'
588 >          error = -1
589 >          return
590 >       endif
591 >    endif
592  
593 <       myStatus = 0
594 <
595 <       call createInteractionMap(myStatus)
390 <
593 >    if (.not. haveGtypeCutoffMap) then        
594 >       myStatus = 0      
595 >       call createGtypeCutoffMap(myStatus)      
596         if (myStatus .ne. 0) then
597 <          write(default_error, *) 'createInteractionMap failed in doForces!'
597 >          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
598            error = -1
599            return
600         endif
# Line 399 | Line 604 | contains
604         call setSimVariables()
605      endif
606  
607 <    if (.not. haveRlist) then
608 <       write(default_error, *) 'rList has not been set in doForces!'
609 <       error = -1
610 <       return
611 <    endif
607 >  !  if (.not. haveRlist) then
608 >  !     write(default_error, *) 'rList has not been set in doForces!'
609 >  !     error = -1
610 >  !     return
611 >  !  endif
612  
613      if (.not. haveNeighborList) then
614         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 428 | Line 633 | contains
633    end subroutine doReadyCheck
634  
635  
636 <  subroutine init_FF(use_RF_c, thisStat)
636 >  subroutine init_FF(thisESM, thisStat)
637  
638 <    logical, intent(in) :: use_RF_c
434 <
638 >    integer, intent(in) :: thisESM
639      integer, intent(out) :: thisStat  
640      integer :: my_status, nMatches
641      integer, pointer :: MatchList(:) => null()
# Line 440 | Line 644 | contains
644      !! assume things are copacetic, unless they aren't
645      thisStat = 0
646  
647 <    !! Fortran's version of a cast:
444 <    FF_uses_RF = use_RF_c
647 >    electrostaticSummationMethod = thisESM
648  
649      !! init_FF is called *after* all of the atom types have been
650      !! defined in atype_module using the new_atype subroutine.
# Line 450 | Line 653 | contains
653      !! interactions are used by the force field.    
654  
655      FF_uses_DirectionalAtoms = .false.
453    FF_uses_LennardJones = .false.
454    FF_uses_Electrostatics = .false.
455    FF_uses_Charges = .false.    
656      FF_uses_Dipoles = .false.
457    FF_uses_Sticky = .false.
458    FF_uses_StickyPower = .false.
657      FF_uses_GayBerne = .false.
658      FF_uses_EAM = .false.
461    FF_uses_Shapes = .false.
462    FF_uses_FLARB = .false.
659  
660      call getMatchingElementList(atypes, "is_Directional", .true., &
661           nMatches, MatchList)
662      if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
663  
468    call getMatchingElementList(atypes, "is_LennardJones", .true., &
469         nMatches, MatchList)
470    if (nMatches .gt. 0) FF_uses_LennardJones = .true.
471
472    call getMatchingElementList(atypes, "is_Electrostatic", .true., &
473         nMatches, MatchList)
474    if (nMatches .gt. 0) then
475       FF_uses_Electrostatics = .true.
476    endif
477
478    call getMatchingElementList(atypes, "is_Charge", .true., &
479         nMatches, MatchList)
480    if (nMatches .gt. 0) then
481       FF_uses_Charges = .true.  
482       FF_uses_Electrostatics = .true.
483    endif
484
664      call getMatchingElementList(atypes, "is_Dipole", .true., &
665           nMatches, MatchList)
666 <    if (nMatches .gt. 0) then
488 <       FF_uses_Dipoles = .true.
489 <       FF_uses_Electrostatics = .true.
490 <       FF_uses_DirectionalAtoms = .true.
491 <    endif
492 <
493 <    call getMatchingElementList(atypes, "is_Quadrupole", .true., &
494 <         nMatches, MatchList)
495 <    if (nMatches .gt. 0) then
496 <       FF_uses_Quadrupoles = .true.
497 <       FF_uses_Electrostatics = .true.
498 <       FF_uses_DirectionalAtoms = .true.
499 <    endif
500 <
501 <    call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
502 <         MatchList)
503 <    if (nMatches .gt. 0) then
504 <       FF_uses_Sticky = .true.
505 <       FF_uses_DirectionalAtoms = .true.
506 <    endif
507 <
508 <    call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
509 <         MatchList)
510 <    if (nMatches .gt. 0) then
511 <       FF_uses_StickyPower = .true.
512 <       FF_uses_DirectionalAtoms = .true.
513 <    endif
666 >    if (nMatches .gt. 0) FF_uses_Dipoles = .true.
667      
668      call getMatchingElementList(atypes, "is_GayBerne", .true., &
669           nMatches, MatchList)
670 <    if (nMatches .gt. 0) then
518 <       FF_uses_GayBerne = .true.
519 <       FF_uses_DirectionalAtoms = .true.
520 <    endif
670 >    if (nMatches .gt. 0) FF_uses_GayBerne = .true.
671  
672      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
673      if (nMatches .gt. 0) FF_uses_EAM = .true.
674  
525    call getMatchingElementList(atypes, "is_Shape", .true., &
526         nMatches, MatchList)
527    if (nMatches .gt. 0) then
528       FF_uses_Shapes = .true.
529       FF_uses_DirectionalAtoms = .true.
530    endif
675  
532    call getMatchingElementList(atypes, "is_FLARB", .true., &
533         nMatches, MatchList)
534    if (nMatches .gt. 0) FF_uses_FLARB = .true.
535
536    !! Assume sanity (for the sake of argument)
676      haveSaneForceField = .true.
677  
678 <    !! check to make sure the FF_uses_RF setting makes sense
678 >    !! check to make sure the reaction field setting makes sense
679  
680 <    if (FF_uses_dipoles) then
681 <       if (FF_uses_RF) then
680 >    if (FF_uses_Dipoles) then
681 >       if (electrostaticSummationMethod == REACTION_FIELD) then
682            dielect = getDielect()
683            call initialize_rf(dielect)
684         endif
685      else
686 <       if (FF_uses_RF) then          
686 >       if (electrostaticSummationMethod == REACTION_FIELD) then
687            write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
688            thisStat = -1
689            haveSaneForceField = .false.
# Line 552 | Line 691 | contains
691         endif
692      endif
693  
555    !sticky module does not contain check_sticky_FF anymore
556    !if (FF_uses_sticky) then
557    !   call check_sticky_FF(my_status)
558    !   if (my_status /= 0) then
559    !      thisStat = -1
560    !      haveSaneForceField = .false.
561    !      return
562    !   end if
563    !endif
564
694      if (FF_uses_EAM) then
695         call init_EAM_FF(my_status)
696         if (my_status /= 0) then
# Line 581 | Line 710 | contains
710         endif
711      endif
712  
584    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
585    endif
586
713      if (.not. haveNeighborList) then
714         !! Create neighbor lists
715         call expandNeighborList(nLocal, my_status)
# Line 647 | Line 773 | contains
773      integer :: localError
774      integer :: propPack_i, propPack_j
775      integer :: loopStart, loopEnd, loop
776 <    integer :: iMap
777 <    real(kind=dp) :: listSkin = 1.0  
776 >    integer :: iHash
777 >  
778  
779      !! initialize local variables  
780  
# Line 738 | Line 864 | contains
864         iend = nGroups - 1
865   #endif
866         outer: do i = istart, iend
741
742 #ifdef IS_MPI
743             me_i = atid_row(i)
744 #else
745             me_i = atid(i)
746 #endif
867  
868            if (update_nlist) point(i) = nlist + 1
869  
# Line 779 | Line 899 | contains
899               me_j = atid(j)
900               call get_interatomic_vector(q_group(:,i), &
901                    q_group(:,j), d_grp, rgrpsq)
902 < #endif
902 > #endif      
903  
904 <             if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
904 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
905                  if (update_nlist) then
906                     nlist = nlist + 1
907  
# Line 902 | Line 1022 | contains
1022                  endif
1023               end if
1024            enddo
1025 +
1026         enddo outer
1027  
1028         if (update_nlist) then
# Line 981 | Line 1102 | contains
1102  
1103      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1104  
1105 <       if (FF_uses_RF .and. SIM_uses_RF) then
1105 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1106  
1107   #ifdef IS_MPI
1108            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 999 | Line 1120 | contains
1120   #else
1121               me_i = atid(i)
1122   #endif
1123 <             iMap = InteractionMap(me_i, me_j)%InteractionHash
1123 >             iHash = InteractionHash(me_i,me_j)
1124              
1125 <             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1125 >             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1126  
1127                  mu_i = getDipoleMoment(me_i)
1128  
# Line 1065 | Line 1186 | contains
1186      real ( kind = dp ), intent(inout) :: rijsq
1187      real ( kind = dp )                :: r
1188      real ( kind = dp ), intent(inout) :: d(3)
1068    real ( kind = dp ) :: ebalance
1189      integer :: me_i, me_j
1190  
1191 <    integer :: iMap
1191 >    integer :: iHash
1192  
1193      r = sqrt(rijsq)
1194      vpair = 0.0d0
# Line 1082 | Line 1202 | contains
1202      me_j = atid(j)
1203   #endif
1204  
1205 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1205 >    iHash = InteractionHash(me_i, me_j)
1206  
1207 <    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1207 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1208         call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1209      endif
1210  
1211 <    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1211 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1212         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1213              pot, eFrame, f, t, do_pot)
1214  
1215 <       if (FF_uses_RF .and. SIM_uses_RF) then
1215 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1216  
1217            ! CHECK ME (RF needs to know about all electrostatic types)
1218            call accumulate_rf(i, j, r, eFrame, sw)
# Line 1101 | Line 1221 | contains
1221  
1222      endif
1223  
1224 <    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1224 >    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1225         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1226              pot, A, f, t, do_pot)
1227      endif
1228  
1229 <    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1229 >    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1230         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1231              pot, A, f, t, do_pot)
1232      endif
1233  
1234 <    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1234 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1235         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1236              pot, A, f, t, do_pot)
1237      endif
1238      
1239 <    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1239 >    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1240   !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1241   !           pot, A, f, t, do_pot)
1242      endif
1243  
1244 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1244 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1245         call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1246              do_pot)
1247      endif
1248  
1249 <    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1249 >    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1250         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1251              pot, A, f, t, do_pot)
1252      endif
1253  
1254 <    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1254 >    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1255         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1256              pot, A, f, t, do_pot)
1257      endif
# Line 1153 | Line 1273 | contains
1273      real ( kind = dp )                :: r, rc
1274      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1275  
1276 <    integer :: me_i, me_j, iMap
1276 >    integer :: me_i, me_j, iHash
1277  
1278 +    r = sqrt(rijsq)
1279 +
1280   #ifdef IS_MPI  
1281      me_i = atid_row(i)
1282      me_j = atid_col(j)  
# Line 1163 | Line 1285 | contains
1285      me_j = atid(j)  
1286   #endif
1287  
1288 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1288 >    iHash = InteractionHash(me_i, me_j)
1289  
1290 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1290 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1291              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1292      endif
1293      
# Line 1362 | Line 1484 | contains
1484  
1485    function FF_UsesDirectionalAtoms() result(doesit)
1486      logical :: doesit
1487 <    doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1366 <         FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1367 <         FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1487 >    doesit = FF_uses_DirectionalAtoms
1488    end function FF_UsesDirectionalAtoms
1489  
1490    function FF_RequiresPrepairCalc() result(doesit)
# Line 1374 | Line 1494 | contains
1494  
1495    function FF_RequiresPostpairCalc() result(doesit)
1496      logical :: doesit
1497 <    doesit = FF_uses_RF
1497 >    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1498    end function FF_RequiresPostpairCalc
1499  
1500   #ifdef PROFILE

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines