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 2266 by chuckv, Thu Jul 28 22:12:45 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.24 2005-07-28 22:12:45 chuckv Exp $, $Date: 2005-07-28 22:12:45 $, $Name: not supported by cvs2svn $, $Revision: 1.24 $
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
154 <    real(kind=dp) :: myRcut
158 < ! Test Types
153 >    integer :: iHash
154 >    !! Test Types
155      logical :: i_is_LJ
156      logical :: i_is_Elect
157      logical :: i_is_Sticky
# 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
250 –  end subroutine createInteractionMap
251  
252 < ! 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.
253 <  subroutine createRcuts(defaultRList,stat)
254 <    real(kind=dp), intent(in), optional :: defaultRList
255 <    integer :: iMap
256 <    integer :: map_i,map_j
257 <    real(kind=dp) :: thisRCut = 0.0_dp
258 <    real(kind=dp) :: actualCutoff = 0.0_dp
259 <    integer, intent(out) :: stat
260 <    integer :: nAtypes
261 <    integer :: myStatus
252 >    haveInteractionHash = .true.
253 >  end subroutine createInteractionHash
254  
255 <    stat = 0
264 <    if (.not. haveInteractionMap) then
255 >  subroutine createGtypeCutoffMap(stat)
256  
257 <       call createInteractionMap(myStatus)
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 +    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 + #ifdef IS_MPI
285 +    nGroupsInRow = getNgroupsInRow(plan_group_row)
286 +    nGroupsInCol = getNgroupsInCol(plan_group_col)
287 + #endif
288 +    nAtypes = getSize(atypes)
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 +
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 +          
335 +          
336 +          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
337 +             biggestAtypeCutoff = atypeMaxCutoff(i)
338 +          endif
339  
340 +       endif
341 +    enddo
342 +  
343  
344 <    nAtypes = getSize(atypes)
345 < ! If we pass a default rcut, set all atypes to that cutoff distance
346 <    if(present(defaultRList)) then
347 <       InteractionMap(:,:)%rList = defaultRList
348 <       InteractionMap(:,:)%rListSq = defaultRList*defaultRList
349 <       haveRlist = .true.
350 <       return
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 <    do map_i = 1,nAtypes
371 <       do map_j = map_i,nAtypes
372 <          iMap = InteractionMap(map_i, map_j)%InteractionHash
373 <          
374 <          if ( iand(iMap, LJ_PAIR).ne.0 ) then
375 < !            thisRCut = getLJCutOff(map_i,map_j)
291 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
292 <          endif
293 <          
294 <          if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
295 < !            thisRCut = getElectrostaticCutOff(map_i,map_j)
296 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
297 <          endif
298 <          
299 <          if ( iand(iMap, STICKY_PAIR).ne.0 ) then
300 < !             thisRCut = getStickyCutOff(map_i,map_j)
301 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
302 <           endif
303 <          
304 <           if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
305 < !              thisRCut = getStickyPowerCutOff(map_i,map_j)
306 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
307 <           endif
308 <          
309 <           if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
310 < !              thisRCut = getGayberneCutOff(map_i,map_j)
311 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
312 <           endif
313 <          
314 <           if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
315 < !              thisRCut = getGaybrneLJCutOff(map_i,map_j)
316 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
317 <           endif
318 <          
319 <           if ( iand(iMap, EAM_PAIR).ne.0 ) then      
320 < !              thisRCut = getEAMCutOff(map_i,map_j)
321 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
322 <           endif
323 <          
324 <           if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
325 < !              thisRCut = getShapeCutOff(map_i,map_j)
326 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
327 <           endif
328 <          
329 <           if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
330 < !              thisRCut = getShapeLJCutOff(map_i,map_j)
331 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
332 <           endif
333 <           InteractionMap(map_i, map_j)%rList = actualCutoff
334 <           InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
335 <        end do
336 <     end do
337 <          haveRlist = .true.
338 <  end subroutine createRcuts
370 >    if(.not.allocated(gtypeMaxCutoffRow)) then
371 >       allocate(gtypeMaxCutoffRow(iend))
372 >    else
373 >       deallocate(gtypeMaxCutoffRow)
374 >       allocate(gtypeMaxCutoffRow(iend))
375 >    endif
376  
377  
378 < !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
379 < !!$  subroutine setRlistDF( this_rlist )
380 < !!$
381 < !!$   real(kind=dp) :: this_rlist
382 < !!$
383 < !!$    rlist = this_rlist
384 < !!$    rlistsq = rlist * rlist
385 < !!$
386 < !!$    haveRlist = .true.
387 < !!$
388 < !!$  end subroutine setRlistDF
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()
356 –    SIM_uses_LennardJones = SimUsesLennardJones()
357 –    SIM_uses_Electrostatics = SimUsesElectrostatics()
358 –    SIM_uses_Charges = SimUsesCharges()
359 –    SIM_uses_Dipoles = SimUsesDipoles()
360 –    SIM_uses_Sticky = SimUsesSticky()
361 –    SIM_uses_StickyPower = SimUsesStickyPower()
362 –    SIM_uses_GayBerne = SimUsesGayBerne()
566      SIM_uses_EAM = SimUsesEAM()
364 –    SIM_uses_Shapes = SimUsesShapes()
365 –    SIM_uses_FLARB = SimUsesFLARB()
366 –    SIM_uses_RF = SimUsesRF()
567      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
568      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
569      SIM_uses_PBC = SimUsesPBC()
# Line 380 | 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)
388 <
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 397 | 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 426 | 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
432 <
638 >    integer, intent(in) :: thisESM
639      integer, intent(out) :: thisStat  
640      integer :: my_status, nMatches
641      integer, pointer :: MatchList(:) => null()
# Line 438 | Line 644 | contains
644      !! assume things are copacetic, unless they aren't
645      thisStat = 0
646  
647 <    !! Fortran's version of a cast:
442 <    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 448 | Line 653 | contains
653      !! interactions are used by the force field.    
654  
655      FF_uses_DirectionalAtoms = .false.
451 –    FF_uses_LennardJones = .false.
452 –    FF_uses_Electrostatics = .false.
453 –    FF_uses_Charges = .false.    
656      FF_uses_Dipoles = .false.
455 –    FF_uses_Sticky = .false.
456 –    FF_uses_StickyPower = .false.
657      FF_uses_GayBerne = .false.
658      FF_uses_EAM = .false.
459 –    FF_uses_Shapes = .false.
460 –    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  
466 –    call getMatchingElementList(atypes, "is_LennardJones", .true., &
467 –         nMatches, MatchList)
468 –    if (nMatches .gt. 0) FF_uses_LennardJones = .true.
469 –
470 –    call getMatchingElementList(atypes, "is_Electrostatic", .true., &
471 –         nMatches, MatchList)
472 –    if (nMatches .gt. 0) then
473 –       FF_uses_Electrostatics = .true.
474 –    endif
475 –
476 –    call getMatchingElementList(atypes, "is_Charge", .true., &
477 –         nMatches, MatchList)
478 –    if (nMatches .gt. 0) then
479 –       FF_uses_Charges = .true.  
480 –       FF_uses_Electrostatics = .true.
481 –    endif
482 –
664      call getMatchingElementList(atypes, "is_Dipole", .true., &
665           nMatches, MatchList)
666 <    if (nMatches .gt. 0) then
486 <       FF_uses_Dipoles = .true.
487 <       FF_uses_Electrostatics = .true.
488 <       FF_uses_DirectionalAtoms = .true.
489 <    endif
490 <
491 <    call getMatchingElementList(atypes, "is_Quadrupole", .true., &
492 <         nMatches, MatchList)
493 <    if (nMatches .gt. 0) then
494 <       FF_uses_Quadrupoles = .true.
495 <       FF_uses_Electrostatics = .true.
496 <       FF_uses_DirectionalAtoms = .true.
497 <    endif
498 <
499 <    call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
500 <         MatchList)
501 <    if (nMatches .gt. 0) then
502 <       FF_uses_Sticky = .true.
503 <       FF_uses_DirectionalAtoms = .true.
504 <    endif
505 <
506 <    call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
507 <         MatchList)
508 <    if (nMatches .gt. 0) then
509 <       FF_uses_StickyPower = .true.
510 <       FF_uses_DirectionalAtoms = .true.
511 <    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
516 <       FF_uses_GayBerne = .true.
517 <       FF_uses_DirectionalAtoms = .true.
518 <    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  
523 –    call getMatchingElementList(atypes, "is_Shape", .true., &
524 –         nMatches, MatchList)
525 –    if (nMatches .gt. 0) then
526 –       FF_uses_Shapes = .true.
527 –       FF_uses_DirectionalAtoms = .true.
528 –    endif
675  
530 –    call getMatchingElementList(atypes, "is_FLARB", .true., &
531 –         nMatches, MatchList)
532 –    if (nMatches .gt. 0) FF_uses_FLARB = .true.
533 –
534 –    !! 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 550 | Line 691 | contains
691         endif
692      endif
693  
553 –    !sticky module does not contain check_sticky_FF anymore
554 –    !if (FF_uses_sticky) then
555 –    !   call check_sticky_FF(my_status)
556 –    !   if (my_status /= 0) then
557 –    !      thisStat = -1
558 –    !      haveSaneForceField = .false.
559 –    !      return
560 –    !   end if
561 –    !endif
562 –
694      if (FF_uses_EAM) then
695         call init_EAM_FF(my_status)
696         if (my_status /= 0) then
# Line 579 | Line 710 | contains
710         endif
711      endif
712  
582 –    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
583 –    endif
584 –
713      if (.not. haveNeighborList) then
714         !! Create neighbor lists
715         call expandNeighborList(nLocal, my_status)
# Line 645 | 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 743 | Line 871 | contains
871  
872            if (update_nlist) then
873   #ifdef IS_MPI
746 –             me_i = atid_row(i)
874               jstart = 1
875               jend = nGroupsInCol
876   #else
750 –             me_i = atid(i)
877               jstart = i+1
878               jend = nGroups
879   #endif
# Line 773 | 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 896 | Line 1022 | contains
1022                  endif
1023               end if
1024            enddo
1025 +
1026         enddo outer
1027  
1028         if (update_nlist) then
# Line 975 | 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 993 | 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 1059 | Line 1186 | contains
1186      real ( kind = dp ), intent(inout) :: rijsq
1187      real ( kind = dp )                :: r
1188      real ( kind = dp ), intent(inout) :: d(3)
1062 –    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 1076 | 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 1095 | 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 1147 | 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 1157 | 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 1356 | Line 1484 | contains
1484  
1485    function FF_UsesDirectionalAtoms() result(doesit)
1486      logical :: doesit
1487 <    doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1360 <         FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1361 <         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 1368 | 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