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 2267 by tim, Fri Jul 29 17:34:06 2005 UTC vs.
Revision 2357 by chuckv, Wed Oct 12 20:18:17 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.56 2005-10-12 20:18:01 chuckv Exp $, $Date: 2005-10-12 20:18:01 $, $Name: not supported by cvs2svn $, $Revision: 1.56 $
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 >    else
189 >       deallocate(InteractionHash)
190 >       allocate(InteractionHash(nAtypes,nAtypes))
191      endif
192 +
193 +    if (.not. allocated(atypeMaxCutoff)) then
194 +       allocate(atypeMaxCutoff(nAtypes))
195 +    else
196 +       deallocate(atypeMaxCutoff)
197 +       allocate(atypeMaxCutoff(nAtypes))
198 +    endif
199          
200      do i = 1, nAtypes
201         call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
# Line 241 | Line 248 | contains
248            if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
249  
250  
251 <          InteractionMap(i,j)%InteractionHash = iHash
252 <          InteractionMap(j,i)%InteractionHash = iHash
251 >          InteractionHash(i,j) = iHash
252 >          InteractionHash(j,i) = iHash
253  
254         end do
255  
256      end do
257  
258 <    haveInteractionMap = .true.
259 <  end subroutine createInteractionMap
258 >    haveInteractionHash = .true.
259 >  end subroutine createInteractionHash
260  
261 < ! 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
261 >  subroutine createGtypeCutoffMap(stat)
262  
263 <    stat = 0
264 <    if (.not. haveInteractionMap) then
263 >    integer, intent(out), optional :: stat
264 >    logical :: i_is_LJ
265 >    logical :: i_is_Elect
266 >    logical :: i_is_Sticky
267 >    logical :: i_is_StickyP
268 >    logical :: i_is_GB
269 >    logical :: i_is_EAM
270 >    logical :: i_is_Shape
271 >    logical :: GtypeFound
272  
273 <       call createInteractionMap(myStatus)
273 >    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
274 >    integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
275 >    integer :: nGroupsInRow
276 >    integer :: nGroupsInCol
277 >    integer :: nGroupTypesRow,nGroupTypesCol
278 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin
279 >    real(kind=dp) :: biggestAtypeCutoff
280  
281 +    stat = 0
282 +    if (.not. haveInteractionHash) then
283 +       call createInteractionHash(myStatus)      
284         if (myStatus .ne. 0) then
285 <          write(default_error, *) 'createInteractionMap failed in doForces!'
285 >          write(default_error, *) 'createInteractionHash failed in doForces!'
286            stat = -1
287            return
288         endif
289      endif
290 <
291 <
290 > #ifdef IS_MPI
291 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
292 >    nGroupsInCol = getNgroupsInCol(plan_group_col)
293 > #endif
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
299 <       haveRlist = .true.
300 <       return
301 <    end if
302 <
303 <    do map_i = 1,nAtypes
304 <       do map_j = map_i,nAtypes
305 <          iMap = InteractionMap(map_i, map_j)%InteractionHash
295 > ! Set all of the initial cutoffs to zero.
296 >    atypeMaxCutoff = 0.0_dp
297 >    do i = 1, nAtypes
298 >       if (SimHasAtype(i)) then    
299 >          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
300 >          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
301 >          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
302 >          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
303 >          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
304 >          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
305 >          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
306            
307 <          if ( iand(iMap, LJ_PAIR).ne.0 ) then
308 <             ! thisRCut = getLJCutOff(map_i,map_j)
309 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
310 <          endif
311 <          
312 <          if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
313 <             ! thisRCut = getElectrostaticCutOff(map_i,map_j)
314 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
307 >
308 >          if (haveDefaultCutoffs) then
309 >             atypeMaxCutoff(i) = defaultRcut
310 >          else
311 >             if (i_is_LJ) then          
312 >                thisRcut = getSigma(i) * 2.5_dp
313 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
314 >             endif
315 >             if (i_is_Elect) then
316 >                thisRcut = defaultRcut
317 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
318 >             endif
319 >             if (i_is_Sticky) then
320 >                thisRcut = getStickyCut(i)
321 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
322 >             endif
323 >             if (i_is_StickyP) then
324 >                thisRcut = getStickyPowerCut(i)
325 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
326 >             endif
327 >             if (i_is_GB) then
328 >                thisRcut = getGayBerneCut(i)
329 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
330 >             endif
331 >             if (i_is_EAM) then
332 >                thisRcut = getEAMCut(i)
333 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
334 >             endif
335 >             if (i_is_Shape) then
336 >                thisRcut = getShapeCut(i)
337 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
338 >             endif
339            endif
340            
341 <          if ( iand(iMap, STICKY_PAIR).ne.0 ) then
342 <             ! thisRCut = getStickyCutOff(map_i,map_j)
343 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
344 <           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
341 >          
342 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
343 >             biggestAtypeCutoff = atypeMaxCutoff(i)
344 >          endif
345  
346 +       endif
347 +    enddo
348 +  
349  
350 < !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
351 < !!$  subroutine setRlistDF( this_rlist )
352 < !!$
353 < !!$   real(kind=dp) :: this_rlist
354 < !!$
355 < !!$    rlist = this_rlist
356 < !!$    rlistsq = rlist * rlist
357 < !!$
358 < !!$    haveRlist = .true.
359 < !!$
360 < !!$  end subroutine setRlistDF
350 >    
351 >    istart = 1
352 >    jstart = 1
353 > #ifdef IS_MPI
354 >    iend = nGroupsInRow
355 >    jend = nGroupsInCol
356 > #else
357 >    iend = nGroups
358 >    jend = nGroups
359 > #endif
360 >    
361 >    !! allocate the groupToGtype and gtypeMaxCutoff here.
362 >    if(.not.allocated(groupToGtypeRow)) then
363 >     !  allocate(groupToGtype(iend))
364 >       allocate(groupToGtypeRow(iend))
365 >    else
366 >       deallocate(groupToGtypeRow)
367 >       allocate(groupToGtypeRow(iend))
368 >    endif
369 >    if(.not.allocated(groupMaxCutoffRow)) then
370 >       allocate(groupMaxCutoffRow(iend))
371 >    else
372 >       deallocate(groupMaxCutoffRow)
373 >       allocate(groupMaxCutoffRow(iend))
374 >    end if
375  
376 +    if(.not.allocated(gtypeMaxCutoffRow)) then
377 +       allocate(gtypeMaxCutoffRow(iend))
378 +    else
379 +       deallocate(gtypeMaxCutoffRow)
380 +       allocate(gtypeMaxCutoffRow(iend))
381 +    endif
382  
383 +
384 + #ifdef IS_MPI
385 +       ! We only allocate new storage if we are in MPI because Ncol /= Nrow
386 +    if(.not.associated(groupToGtypeCol)) then
387 +       allocate(groupToGtypeCol(jend))
388 +    else
389 +       deallocate(groupToGtypeCol)
390 +       allocate(groupToGtypeCol(jend))
391 +    end if
392 +
393 +    if(.not.associated(groupToGtypeCol)) then
394 +       allocate(groupToGtypeCol(jend))
395 +    else
396 +       deallocate(groupToGtypeCol)
397 +       allocate(groupToGtypeCol(jend))
398 +    end if
399 +    if(.not.associated(gtypeMaxCutoffCol)) then
400 +       allocate(gtypeMaxCutoffCol(jend))
401 +    else
402 +       deallocate(gtypeMaxCutoffCol)      
403 +       allocate(gtypeMaxCutoffCol(jend))
404 +    end if
405 +
406 +       groupMaxCutoffCol = 0.0_dp
407 +       gtypeMaxCutoffCol = 0.0_dp
408 +
409 + #endif
410 +       groupMaxCutoffRow = 0.0_dp
411 +       gtypeMaxCutoffRow = 0.0_dp
412 +
413 +
414 +    !! first we do a single loop over the cutoff groups to find the
415 +    !! largest cutoff for any atypes present in this group.  We also
416 +    !! create gtypes at this point.
417 +    
418 +    tol = 1.0d-6
419 +    nGroupTypesRow = 0
420 +
421 +    do i = istart, iend      
422 +       n_in_i = groupStartRow(i+1) - groupStartRow(i)
423 +       groupMaxCutoffRow(i) = 0.0_dp
424 +       do ia = groupStartRow(i), groupStartRow(i+1)-1
425 +          atom1 = groupListRow(ia)
426 + #ifdef IS_MPI
427 +          me_i = atid_row(atom1)
428 + #else
429 +          me_i = atid(atom1)
430 + #endif          
431 +          if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
432 +             groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
433 +          endif          
434 +       enddo
435 +
436 +       if (nGroupTypesRow.eq.0) then
437 +          nGroupTypesRow = nGroupTypesRow + 1
438 +          gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
439 +          groupToGtypeRow(i) = nGroupTypesRow
440 +       else
441 +          GtypeFound = .false.
442 +          do g = 1, nGroupTypesRow
443 +             if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
444 +                groupToGtypeRow(i) = g
445 +                GtypeFound = .true.
446 +             endif
447 +          enddo
448 +          if (.not.GtypeFound) then            
449 +             nGroupTypesRow = nGroupTypesRow + 1
450 +             gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
451 +             groupToGtypeRow(i) = nGroupTypesRow
452 +          endif
453 +       endif
454 +    enddo    
455 +
456 + #ifdef IS_MPI
457 +    do j = jstart, jend      
458 +       n_in_j = groupStartCol(j+1) - groupStartCol(j)
459 +       groupMaxCutoffCol(j) = 0.0_dp
460 +       do ja = groupStartCol(j), groupStartCol(j+1)-1
461 +          atom1 = groupListCol(ja)
462 +
463 +          me_j = atid_col(atom1)
464 +
465 +          if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
466 +             groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
467 +          endif          
468 +       enddo
469 +
470 +       if (nGroupTypesCol.eq.0) then
471 +          nGroupTypesCol = nGroupTypesCol + 1
472 +          gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
473 +          groupToGtypeCol(j) = nGroupTypesCol
474 +       else
475 +          GtypeFound = .false.
476 +          do g = 1, nGroupTypesCol
477 +             if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
478 +                groupToGtypeCol(j) = g
479 +                GtypeFound = .true.
480 +             endif
481 +          enddo
482 +          if (.not.GtypeFound) then            
483 +             nGroupTypesCol = nGroupTypesCol + 1
484 +             gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
485 +             groupToGtypeCol(j) = nGroupTypesCol
486 +          endif
487 +       endif
488 +    enddo    
489 +
490 + #else
491 + ! Set pointers to information we just found
492 +    nGroupTypesCol = nGroupTypesRow
493 +    groupToGtypeCol => groupToGtypeRow
494 +    gtypeMaxCutoffCol => gtypeMaxCutoffRow
495 +    groupMaxCutoffCol => groupMaxCutoffRow
496 + #endif
497 +
498 +
499 +
500 +
501 +
502 +    !! allocate the gtypeCutoffMap here.
503 +    allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
504 +    !! then we do a double loop over all the group TYPES to find the cutoff
505 +    !! map between groups of two types
506 +    tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
507 +
508 +    do i = 1, nGroupTypesRow
509 +       do j = 1, nGroupTypesCol
510 +      
511 +          select case(cutoffPolicy)
512 +          case(TRADITIONAL_CUTOFF_POLICY)
513 +             thisRcut = tradRcut
514 +          case(MIX_CUTOFF_POLICY)
515 +             thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
516 +          case(MAX_CUTOFF_POLICY)
517 +             thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
518 +          case default
519 +             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
520 +             return
521 +          end select
522 +          gtypeCutoffMap(i,j)%rcut = thisRcut
523 +          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
524 +          skin = defaultRlist - defaultRcut
525 +          listSkin = skin ! set neighbor list skin thickness
526 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
527 +
528 +          ! sanity check
529 +
530 +          if (haveDefaultCutoffs) then
531 +             if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
532 +                call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
533 +             endif
534 +          endif
535 +       enddo
536 +    enddo
537 +    if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
538 +    if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
539 +    if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
540 + #ifdef IS_MPI
541 +    if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
542 +    if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
543 + #endif
544 +    groupMaxCutoffCol => null()
545 +    gtypeMaxCutoffCol => null()
546 +    
547 +    haveGtypeCutoffMap = .true.
548 +   end subroutine createGtypeCutoffMap
549 +
550 +   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
551 +     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
552 +     integer, intent(in) :: cutPolicy
553 +
554 +     defaultRcut = defRcut
555 +     defaultRsw = defRsw
556 +     defaultRlist = defRlist
557 +     cutoffPolicy = cutPolicy
558 +
559 +     haveDefaultCutoffs = .true.
560 +   end subroutine setDefaultCutoffs
561 +
562 +   subroutine setCutoffPolicy(cutPolicy)
563 +
564 +     integer, intent(in) :: cutPolicy
565 +     cutoffPolicy = cutPolicy
566 +     call createGtypeCutoffMap()
567 +   end subroutine setCutoffPolicy
568 +    
569 +    
570    subroutine setSimVariables()
571      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()
572      SIM_uses_EAM = SimUsesEAM()
366 –    SIM_uses_Shapes = SimUsesShapes()
367 –    SIM_uses_FLARB = SimUsesFLARB()
368 –    SIM_uses_RF = SimUsesRF()
573      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
574      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
575      SIM_uses_PBC = SimUsesPBC()
# Line 382 | Line 586 | contains
586  
587      error = 0
588  
589 <    if (.not. haveInteractionMap) then
589 >    if (.not. haveInteractionHash) then      
590 >       myStatus = 0      
591 >       call createInteractionHash(myStatus)      
592 >       if (myStatus .ne. 0) then
593 >          write(default_error, *) 'createInteractionHash failed in doForces!'
594 >          error = -1
595 >          return
596 >       endif
597 >    endif
598  
599 <       myStatus = 0
600 <
601 <       call createInteractionMap(myStatus)
390 <
599 >    if (.not. haveGtypeCutoffMap) then        
600 >       myStatus = 0      
601 >       call createGtypeCutoffMap(myStatus)      
602         if (myStatus .ne. 0) then
603 <          write(default_error, *) 'createInteractionMap failed in doForces!'
603 >          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
604            error = -1
605            return
606         endif
# Line 399 | Line 610 | contains
610         call setSimVariables()
611      endif
612  
613 <    if (.not. haveRlist) then
614 <       write(default_error, *) 'rList has not been set in doForces!'
615 <       error = -1
616 <       return
617 <    endif
613 >  !  if (.not. haveRlist) then
614 >  !     write(default_error, *) 'rList has not been set in doForces!'
615 >  !     error = -1
616 >  !     return
617 >  !  endif
618  
619      if (.not. haveNeighborList) then
620         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 428 | Line 639 | contains
639    end subroutine doReadyCheck
640  
641  
642 <  subroutine init_FF(use_RF_c, thisStat)
642 >  subroutine init_FF(thisESM, thisStat)
643  
644 <    logical, intent(in) :: use_RF_c
434 <
644 >    integer, intent(in) :: thisESM
645      integer, intent(out) :: thisStat  
646      integer :: my_status, nMatches
647      integer, pointer :: MatchList(:) => null()
# Line 440 | Line 650 | contains
650      !! assume things are copacetic, unless they aren't
651      thisStat = 0
652  
653 <    !! Fortran's version of a cast:
444 <    FF_uses_RF = use_RF_c
653 >    electrostaticSummationMethod = thisESM
654  
655      !! init_FF is called *after* all of the atom types have been
656      !! defined in atype_module using the new_atype subroutine.
# Line 450 | Line 659 | contains
659      !! interactions are used by the force field.    
660  
661      FF_uses_DirectionalAtoms = .false.
453 –    FF_uses_LennardJones = .false.
454 –    FF_uses_Electrostatics = .false.
455 –    FF_uses_Charges = .false.    
662      FF_uses_Dipoles = .false.
457 –    FF_uses_Sticky = .false.
458 –    FF_uses_StickyPower = .false.
663      FF_uses_GayBerne = .false.
664      FF_uses_EAM = .false.
461 –    FF_uses_Shapes = .false.
462 –    FF_uses_FLARB = .false.
665  
666      call getMatchingElementList(atypes, "is_Directional", .true., &
667           nMatches, MatchList)
668      if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
669  
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 –
670      call getMatchingElementList(atypes, "is_Dipole", .true., &
671           nMatches, MatchList)
672 <    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
672 >    if (nMatches .gt. 0) FF_uses_Dipoles = .true.
673      
674      call getMatchingElementList(atypes, "is_GayBerne", .true., &
675           nMatches, MatchList)
676 <    if (nMatches .gt. 0) then
518 <       FF_uses_GayBerne = .true.
519 <       FF_uses_DirectionalAtoms = .true.
520 <    endif
676 >    if (nMatches .gt. 0) FF_uses_GayBerne = .true.
677  
678      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
679      if (nMatches .gt. 0) FF_uses_EAM = .true.
680  
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
681  
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)
682      haveSaneForceField = .true.
683  
684 <    !! check to make sure the FF_uses_RF setting makes sense
684 >    !! check to make sure the reaction field setting makes sense
685  
686 <    if (FF_uses_dipoles) then
687 <       if (FF_uses_RF) then
686 >    if (FF_uses_Dipoles) then
687 >       if (electrostaticSummationMethod == REACTION_FIELD) then
688            dielect = getDielect()
689            call initialize_rf(dielect)
690         endif
691      else
692 <       if (FF_uses_RF) then          
692 >       if (electrostaticSummationMethod == REACTION_FIELD) then
693            write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
694            thisStat = -1
695            haveSaneForceField = .false.
# Line 552 | Line 697 | contains
697         endif
698      endif
699  
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 –
700      if (FF_uses_EAM) then
701         call init_EAM_FF(my_status)
702         if (my_status /= 0) then
# Line 581 | Line 716 | contains
716         endif
717      endif
718  
584 –    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
585 –    endif
586 –
719      if (.not. haveNeighborList) then
720         !! Create neighbor lists
721         call expandNeighborList(nLocal, my_status)
# Line 617 | Line 749 | contains
749  
750      !! Stress Tensor
751      real( kind = dp), dimension(9) :: tau  
752 <    real ( kind = dp ) :: pot
752 >    real ( kind = dp ),dimension(POT_ARRAY_SIZE) :: pot
753      logical ( kind = 2) :: do_pot_c, do_stress_c
754      logical :: do_pot
755      logical :: do_stress
756      logical :: in_switching_region
757   #ifdef IS_MPI
758 <    real( kind = DP ) :: pot_local
758 >    real( kind = DP ), dimension(POT_ARRAY_SIZE) :: pot_local
759      integer :: nAtomsInRow
760      integer :: nAtomsInCol
761      integer :: nprocs
# Line 647 | Line 779 | contains
779      integer :: localError
780      integer :: propPack_i, propPack_j
781      integer :: loopStart, loopEnd, loop
782 <    integer :: iMap
783 <    real(kind=dp) :: listSkin = 1.0  
782 >    integer :: iHash
783 >  
784  
785      !! initialize local variables  
786  
# Line 739 | Line 871 | contains
871   #endif
872         outer: do i = istart, iend
873  
742 – #ifdef IS_MPI
743 –             me_i = atid_row(i)
744 – #else
745 –             me_i = atid(i)
746 – #endif
747 –
874            if (update_nlist) point(i) = nlist + 1
875  
876            n_in_i = groupStartRow(i+1) - groupStartRow(i)
# Line 779 | Line 905 | contains
905               me_j = atid(j)
906               call get_interatomic_vector(q_group(:,i), &
907                    q_group(:,j), d_grp, rgrpsq)
908 < #endif
908 > #endif      
909  
910 <             if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
910 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
911                  if (update_nlist) then
912                     nlist = nlist + 1
913  
# Line 902 | Line 1028 | contains
1028                  endif
1029               end if
1030            enddo
1031 +
1032         enddo outer
1033  
1034         if (update_nlist) then
# Line 961 | Line 1088 | contains
1088  
1089      if (do_pot) then
1090         ! scatter/gather pot_row into the members of my column
1091 <       call scatter(pot_Row, pot_Temp, plan_atom_row)
1092 <
1091 >       do i = 1,POT_ARRAY_SIZE
1092 >          call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1093 >       end do
1094         ! scatter/gather pot_local into all other procs
1095         ! add resultant to get total pot
1096         do i = 1, nlocal
1097 <          pot_local = pot_local + pot_Temp(i)
1097 >          pot_local(1:POT_ARRAY_SIZE) = pot_local(1:POT_ARRAY_SIZE) &
1098 >               + pot_Temp(1:POT_ARRAY_SIZE,i)
1099         enddo
1100  
1101         pot_Temp = 0.0_DP
1102 <
1103 <       call scatter(pot_Col, pot_Temp, plan_atom_col)
1102 >       do i = 1,POT_ARRAY_SIZE
1103 >          call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1104 >       end do
1105         do i = 1, nlocal
1106 <          pot_local = pot_local + pot_Temp(i)
1106 >          pot_local(1:POT_ARRAY_SIZE) = pot_local(1:POT_ARRAY_SIZE)&
1107 >               + pot_Temp(1:POT_ARRAY_SIZE,i)
1108         enddo
1109  
1110      endif
# Line 981 | Line 1112 | contains
1112  
1113      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1114  
1115 <       if (FF_uses_RF .and. SIM_uses_RF) then
1115 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1116  
1117   #ifdef IS_MPI
1118            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 999 | Line 1130 | contains
1130   #else
1131               me_i = atid(i)
1132   #endif
1133 <             iMap = InteractionMap(me_i, me_j)%InteractionHash
1133 >             iHash = InteractionHash(me_i,me_j)
1134              
1135 <             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1135 >             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1136  
1137                  mu_i = getDipoleMoment(me_i)
1138  
# Line 1012 | Line 1143 | contains
1143                  !! potential and torques:
1144                  call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1145   #ifdef IS_MPI
1146 <                pot_local = pot_local + rfpot
1146 >                pot_local(RF_POT) = pot_local(RF_POT) + rfpot
1147   #else
1148 <                pot = pot + rfpot
1148 >                pot(RF_POT) = pot(RF_POT) + rfpot
1149  
1150   #endif
1151               endif
# Line 1026 | Line 1157 | contains
1157   #ifdef IS_MPI
1158  
1159      if (do_pot) then
1160 <       pot = pot + pot_local
1160 >       pot(1:POT_ARRAY_SIZE) = pot(1:POT_ARRAY_SIZE) &
1161 >            + pot_local(1:POT_ARRAY_SIZE)
1162         !! we assume the c code will do the allreduce to get the total potential
1163         !! we could do it right here if we needed to...
1164      endif
# Line 1052 | Line 1184 | contains
1184    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1185         eFrame, A, f, t, pot, vpair, fpair)
1186  
1187 <    real( kind = dp ) :: pot, vpair, sw
1187 >    real( kind = dp ) :: vpair, sw
1188 >    real( kind = dp ), dimension(POT_ARRAY_SIZE) :: pot
1189      real( kind = dp ), dimension(3) :: fpair
1190      real( kind = dp ), dimension(nLocal)   :: mfact
1191      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 1065 | Line 1198 | contains
1198      real ( kind = dp ), intent(inout) :: rijsq
1199      real ( kind = dp )                :: r
1200      real ( kind = dp ), intent(inout) :: d(3)
1068 –    real ( kind = dp ) :: ebalance
1201      integer :: me_i, me_j
1202  
1203 <    integer :: iMap
1203 >    integer :: iHash
1204  
1205      r = sqrt(rijsq)
1206      vpair = 0.0d0
# Line 1082 | Line 1214 | contains
1214      me_j = atid(j)
1215   #endif
1216  
1217 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1217 >    iHash = InteractionHash(me_i, me_j)
1218  
1219 <    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1220 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1219 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1220 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot(LJ_POT), f, do_pot)
1221      endif
1222  
1223 <    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1223 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1224         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1225 <            pot, eFrame, f, t, do_pot)
1225 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1226  
1227 <       if (FF_uses_RF .and. SIM_uses_RF) then
1227 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1228  
1229            ! CHECK ME (RF needs to know about all electrostatic types)
1230            call accumulate_rf(i, j, r, eFrame, sw)
# Line 1101 | Line 1233 | contains
1233  
1234      endif
1235  
1236 <    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1236 >    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1237         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1238 <            pot, A, f, t, do_pot)
1238 >            pot(STICKY_POT), A, f, t, do_pot)
1239      endif
1240  
1241 <    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1241 >    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1242         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1243 <            pot, A, f, t, do_pot)
1243 >            pot(STICKYPOWER_POT), A, f, t, do_pot)
1244      endif
1245  
1246 <    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1246 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1247         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1248 <            pot, A, f, t, do_pot)
1248 >            pot(GAYBERNE_POT), A, f, t, do_pot)
1249      endif
1250      
1251 <    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1251 >    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1252   !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1253 < !           pot, A, f, t, do_pot)
1253 > !           pot(GAYBERNE_LJ_POT), A, f, t, do_pot)
1254      endif
1255  
1256 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1257 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1256 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1257 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot(EAM_POT), f, &
1258              do_pot)
1259      endif
1260  
1261 <    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1261 >    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1262         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1263 <            pot, A, f, t, do_pot)
1263 >            pot(SHAPE_POT), A, f, t, do_pot)
1264      endif
1265  
1266 <    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1266 >    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1267         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1268 <            pot, A, f, t, do_pot)
1268 >            pot(SHAPE_LJ_POT), A, f, t, do_pot)
1269      endif
1270      
1271    end subroutine do_pair
# Line 1141 | Line 1273 | contains
1273    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1274         do_pot, do_stress, eFrame, A, f, t, pot)
1275  
1276 <    real( kind = dp ) :: pot, sw
1276 >    real( kind = dp ) :: sw
1277 >    real( kind = dp ), dimension(POT_ARRAY_SIZE) :: pot
1278      real( kind = dp ), dimension(9,nLocal) :: eFrame
1279      real (kind=dp), dimension(9,nLocal) :: A
1280      real (kind=dp), dimension(3,nLocal) :: f
# Line 1153 | Line 1286 | contains
1286      real ( kind = dp )                :: r, rc
1287      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1288  
1289 <    integer :: me_i, me_j, iMap
1289 >    integer :: me_i, me_j, iHash
1290  
1291 +    r = sqrt(rijsq)
1292 +
1293   #ifdef IS_MPI  
1294      me_i = atid_row(i)
1295      me_j = atid_col(j)  
# Line 1163 | Line 1298 | contains
1298      me_j = atid(j)  
1299   #endif
1300  
1301 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1301 >    iHash = InteractionHash(me_i, me_j)
1302  
1303 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1303 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1304              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1305      endif
1306      
# Line 1174 | Line 1309 | contains
1309  
1310    subroutine do_preforce(nlocal,pot)
1311      integer :: nlocal
1312 <    real( kind = dp ) :: pot
1312 >    real( kind = dp ),dimension(POT_ARRAY_SIZE) :: pot
1313  
1314      if (FF_uses_EAM .and. SIM_uses_EAM) then
1315 <       call calc_EAM_preforce_Frho(nlocal,pot)
1315 >       call calc_EAM_preforce_Frho(nlocal,pot(EAM_POT))
1316      endif
1317  
1318  
# Line 1362 | Line 1497 | contains
1497  
1498    function FF_UsesDirectionalAtoms() result(doesit)
1499      logical :: doesit
1500 <    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
1500 >    doesit = FF_uses_DirectionalAtoms
1501    end function FF_UsesDirectionalAtoms
1502  
1503    function FF_RequiresPrepairCalc() result(doesit)
# Line 1374 | Line 1507 | contains
1507  
1508    function FF_RequiresPostpairCalc() result(doesit)
1509      logical :: doesit
1510 <    doesit = FF_uses_RF
1510 >    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1511    end function FF_RequiresPostpairCalc
1512  
1513   #ifdef PROFILE

Diff Legend

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