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

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2262 by chuckv, Sun Jul 3 20:53:43 2005 UTC vs.
Revision 2367 by kdaily, Fri Oct 14 15:44:37 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.23 2005-07-03 20:53:43 chuckv Exp $, $Date: 2005-07-03 20:53:43 $, $Name: not supported by cvs2svn $, $Revision: 1.23 $
48 > !! @version $Id: doForces.F90,v 1.58 2005-10-14 15:44:37 kdaily Exp $, $Date: 2005-10-14 15:44:37 $, $Name: not supported by cvs2svn $, $Revision: 1.58 $
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 :: status
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 <    
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 185 | 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 240 | 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
249 –  end subroutine createInteractionMap
257  
258 < ! 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.
259 <  subroutine createRcuts(defaultRList)
253 <    real(kind=dp), intent(in), optional :: defaultRList
254 <    integer :: iMap
255 <    integer :: map_i,map_j
256 <    real(kind=dp) :: thisRCut = 0.0_dp
257 <    real(kind=dp) :: actualCutoff = 0.0_dp
258 <    integer :: nAtypes
258 >    haveInteractionHash = .true.
259 >  end subroutine createInteractionHash
260  
261 <    if(.not. allocated(InteractionMap)) return
261 >  subroutine createGtypeCutoffMap(stat)
262  
263 <    nAtypes = getSize(atypes)
264 < ! If we pass a default rcut, set all atypes to that cutoff distance
265 <    if(present(defaultRList)) then
266 <       InteractionMap(:,:)%rList = defaultRList
267 <       InteractionMap(:,:)%rListSq = defaultRList*defaultRList
268 <       haveRlist = .true.
269 <       return
270 <    end if
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 <    do map_i = 1,nAtypes
274 <       do map_j = map_i,nAtypes
275 <          iMap = InteractionMap(map_i, map_j)%InteractionHash
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, *) 'createInteractionHash failed in doForces!'
286 >          stat = -1
287 >          return
288 >       endif
289 >    endif
290 > #ifdef IS_MPI
291 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
292 >    nGroupsInCol = getNgroupsInCol(plan_group_col)
293 > #endif
294 >    nAtypes = getSize(atypes)
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
289 <          
290 <           if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
291 < !              thisRCut = getStickyPowerCutOff(map_i,map_j)
292 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
293 <           endif
294 <          
295 <           if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
296 < !              thisRCut = getGayberneCutOff(map_i,map_j)
297 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
298 <           endif
299 <          
300 <           if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
301 < !              thisRCut = getGaybrneLJCutOff(map_i,map_j)
302 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
303 <           endif
304 <          
305 <           if ( iand(iMap, EAM_PAIR).ne.0 ) then      
306 < !              thisRCut = getEAMCutOff(map_i,map_j)
307 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
308 <           endif
309 <          
310 <           if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
311 < !              thisRCut = getShapeCutOff(map_i,map_j)
312 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
313 <           endif
314 <          
315 <           if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
316 < !              thisRCut = getShapeLJCutOff(map_i,map_j)
317 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
318 <           endif
319 <           InteractionMap(map_i, map_j)%rList = actualCutoff
320 <           InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
321 <        end do
322 <     end do
323 <          haveRlist = .true.
324 <  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()
342 –    SIM_uses_LennardJones = SimUsesLennardJones()
343 –    SIM_uses_Electrostatics = SimUsesElectrostatics()
344 –    SIM_uses_Charges = SimUsesCharges()
345 –    SIM_uses_Dipoles = SimUsesDipoles()
346 –    SIM_uses_Sticky = SimUsesSticky()
347 –    SIM_uses_StickyPower = SimUsesStickyPower()
348 –    SIM_uses_GayBerne = SimUsesGayBerne()
572      SIM_uses_EAM = SimUsesEAM()
350 –    SIM_uses_Shapes = SimUsesShapes()
351 –    SIM_uses_FLARB = SimUsesFLARB()
352 –    SIM_uses_RF = SimUsesRF()
573      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
574      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
575      SIM_uses_PBC = SimUsesPBC()
# Line 366 | 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)
374 <
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 383 | 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 412 | 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
418 <
644 >    integer, intent(in) :: thisESM
645      integer, intent(out) :: thisStat  
646      integer :: my_status, nMatches
647      integer, pointer :: MatchList(:) => null()
# Line 424 | Line 650 | contains
650      !! assume things are copacetic, unless they aren't
651      thisStat = 0
652  
653 <    !! Fortran's version of a cast:
428 <    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 434 | Line 659 | contains
659      !! interactions are used by the force field.    
660  
661      FF_uses_DirectionalAtoms = .false.
437 –    FF_uses_LennardJones = .false.
438 –    FF_uses_Electrostatics = .false.
439 –    FF_uses_Charges = .false.    
662      FF_uses_Dipoles = .false.
441 –    FF_uses_Sticky = .false.
442 –    FF_uses_StickyPower = .false.
663      FF_uses_GayBerne = .false.
664      FF_uses_EAM = .false.
445 –    FF_uses_Shapes = .false.
446 –    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  
452 –    call getMatchingElementList(atypes, "is_LennardJones", .true., &
453 –         nMatches, MatchList)
454 –    if (nMatches .gt. 0) FF_uses_LennardJones = .true.
455 –
456 –    call getMatchingElementList(atypes, "is_Electrostatic", .true., &
457 –         nMatches, MatchList)
458 –    if (nMatches .gt. 0) then
459 –       FF_uses_Electrostatics = .true.
460 –    endif
461 –
462 –    call getMatchingElementList(atypes, "is_Charge", .true., &
463 –         nMatches, MatchList)
464 –    if (nMatches .gt. 0) then
465 –       FF_uses_Charges = .true.  
466 –       FF_uses_Electrostatics = .true.
467 –    endif
468 –
670      call getMatchingElementList(atypes, "is_Dipole", .true., &
671           nMatches, MatchList)
672 <    if (nMatches .gt. 0) then
472 <       FF_uses_Dipoles = .true.
473 <       FF_uses_Electrostatics = .true.
474 <       FF_uses_DirectionalAtoms = .true.
475 <    endif
476 <
477 <    call getMatchingElementList(atypes, "is_Quadrupole", .true., &
478 <         nMatches, MatchList)
479 <    if (nMatches .gt. 0) then
480 <       FF_uses_Quadrupoles = .true.
481 <       FF_uses_Electrostatics = .true.
482 <       FF_uses_DirectionalAtoms = .true.
483 <    endif
484 <
485 <    call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
486 <         MatchList)
487 <    if (nMatches .gt. 0) then
488 <       FF_uses_Sticky = .true.
489 <       FF_uses_DirectionalAtoms = .true.
490 <    endif
491 <
492 <    call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
493 <         MatchList)
494 <    if (nMatches .gt. 0) then
495 <       FF_uses_StickyPower = .true.
496 <       FF_uses_DirectionalAtoms = .true.
497 <    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
502 <       FF_uses_GayBerne = .true.
503 <       FF_uses_DirectionalAtoms = .true.
504 <    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  
509 –    call getMatchingElementList(atypes, "is_Shape", .true., &
510 –         nMatches, MatchList)
511 –    if (nMatches .gt. 0) then
512 –       FF_uses_Shapes = .true.
513 –       FF_uses_DirectionalAtoms = .true.
514 –    endif
681  
516 –    call getMatchingElementList(atypes, "is_FLARB", .true., &
517 –         nMatches, MatchList)
518 –    if (nMatches .gt. 0) FF_uses_FLARB = .true.
519 –
520 –    !! 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 536 | Line 697 | contains
697         endif
698      endif
699  
539 –    !sticky module does not contain check_sticky_FF anymore
540 –    !if (FF_uses_sticky) then
541 –    !   call check_sticky_FF(my_status)
542 –    !   if (my_status /= 0) then
543 –    !      thisStat = -1
544 –    !      haveSaneForceField = .false.
545 –    !      return
546 –    !   end if
547 –    !endif
548 –
700      if (FF_uses_EAM) then
701         call init_EAM_FF(my_status)
702         if (my_status /= 0) then
# Line 565 | Line 716 | contains
716         endif
717      endif
718  
568 –    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
569 –    endif
570 –
719      if (.not. haveNeighborList) then
720         !! Create neighbor lists
721         call expandNeighborList(nLocal, my_status)
# Line 601 | Line 749 | contains
749  
750      !! Stress Tensor
751      real( kind = dp), dimension(9) :: tau  
752 <    real ( kind = dp ) :: pot
752 >    real ( kind = dp ),dimension(LR_POT_TYPES) :: 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(LR_POT_TYPES) :: pot_local
759      integer :: nAtomsInRow
760      integer :: nAtomsInCol
761      integer :: nprocs
# Line 631 | 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 750 | Line 898 | contains
898               endif
899  
900   #ifdef IS_MPI
901 +             me_j = atid_col(j)
902               call get_interatomic_vector(q_group_Row(:,i), &
903                    q_group_Col(:,j), d_grp, rgrpsq)
904   #else
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 878 | Line 1028 | contains
1028                  endif
1029               end if
1030            enddo
1031 +
1032         enddo outer
1033  
1034         if (update_nlist) then
# Line 937 | 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,LR_POT_TYPES
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:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1098 >               + pot_Temp(1:LR_POT_TYPES,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,LR_POT_TYPES
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:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1107 >               + pot_Temp(1:LR_POT_TYPES,i)
1108         enddo
1109  
1110      endif
# Line 957 | 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 975 | 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 988 | 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(ELECTROSTATIC_POT) = pot_local(ELECTROSTATIC_POT) + rfpot
1147   #else
1148 <                pot = pot + rfpot
1148 >                pot(ELECTROSTATIC_POT) = pot(ELECTROSTATIC_POT) + rfpot
1149  
1150   #endif
1151               endif
# Line 1002 | Line 1157 | contains
1157   #ifdef IS_MPI
1158  
1159      if (do_pot) then
1160 <       pot = pot + pot_local
1160 >       pot(1:LR_POT_TYPES) = pot(1:LR_POT_TYPES) &
1161 >            + pot_local(1:LR_POT_TYPES)
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 1028 | 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(LR_POT_TYPES) :: pot
1189      real( kind = dp ), dimension(3) :: fpair
1190      real( kind = dp ), dimension(nLocal)   :: mfact
1191      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 1041 | Line 1198 | contains
1198      real ( kind = dp ), intent(inout) :: rijsq
1199      real ( kind = dp )                :: r
1200      real ( kind = dp ), intent(inout) :: d(3)
1044 –    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 1058 | 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, &
1221 >            pot(VDW_POT), f, do_pot)
1222      endif
1223  
1224 <    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1224 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1225         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1226 <            pot, eFrame, f, t, do_pot)
1226 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1227  
1228 <       if (FF_uses_RF .and. SIM_uses_RF) then
1228 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1229  
1230            ! CHECK ME (RF needs to know about all electrostatic types)
1231            call accumulate_rf(i, j, r, eFrame, sw)
# Line 1077 | Line 1234 | contains
1234  
1235      endif
1236  
1237 <    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1237 >    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1238         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1239 <            pot, A, f, t, do_pot)
1239 >            pot(HB_POT), A, f, t, do_pot)
1240      endif
1241  
1242 <    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1242 >    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1243         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1244 <            pot, A, f, t, do_pot)
1244 >            pot(HB_POT), A, f, t, do_pot)
1245      endif
1246  
1247 <    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1247 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1248         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1249 <            pot, A, f, t, do_pot)
1249 >            pot(VDW_POT), A, f, t, do_pot)
1250      endif
1251      
1252 <    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1253 < !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1254 < !           pot, A, f, t, do_pot)
1252 >    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1253 >       call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1254 >            pot(VDW_POT), A, f, t, do_pot)
1255      endif
1256  
1257 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1258 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1259 <            do_pot)
1257 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1258 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1259 >            pot(METALLIC_POT), f, do_pot)
1260      endif
1261  
1262 <    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1262 >    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1263         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1264 <            pot, A, f, t, do_pot)
1264 >            pot(VDW_POT), A, f, t, do_pot)
1265      endif
1266  
1267 <    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1267 >    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1268         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1269 <            pot, A, f, t, do_pot)
1269 >            pot(VDW_POT), A, f, t, do_pot)
1270      endif
1271      
1272    end subroutine do_pair
# Line 1117 | Line 1274 | contains
1274    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1275         do_pot, do_stress, eFrame, A, f, t, pot)
1276  
1277 <    real( kind = dp ) :: pot, sw
1277 >    real( kind = dp ) :: sw
1278 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1279      real( kind = dp ), dimension(9,nLocal) :: eFrame
1280      real (kind=dp), dimension(9,nLocal) :: A
1281      real (kind=dp), dimension(3,nLocal) :: f
# Line 1129 | Line 1287 | contains
1287      real ( kind = dp )                :: r, rc
1288      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1289  
1290 <    integer :: me_i, me_j, iMap
1290 >    integer :: me_i, me_j, iHash
1291  
1292 +    r = sqrt(rijsq)
1293 +
1294   #ifdef IS_MPI  
1295      me_i = atid_row(i)
1296      me_j = atid_col(j)  
# Line 1139 | Line 1299 | contains
1299      me_j = atid(j)  
1300   #endif
1301  
1302 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1302 >    iHash = InteractionHash(me_i, me_j)
1303  
1304 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1304 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1305              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1306      endif
1307      
# Line 1150 | Line 1310 | contains
1310  
1311    subroutine do_preforce(nlocal,pot)
1312      integer :: nlocal
1313 <    real( kind = dp ) :: pot
1313 >    real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1314  
1315      if (FF_uses_EAM .and. SIM_uses_EAM) then
1316 <       call calc_EAM_preforce_Frho(nlocal,pot)
1316 >       call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1317      endif
1318  
1319  
# Line 1338 | Line 1498 | contains
1498  
1499    function FF_UsesDirectionalAtoms() result(doesit)
1500      logical :: doesit
1501 <    doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1342 <         FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1343 <         FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1501 >    doesit = FF_uses_DirectionalAtoms
1502    end function FF_UsesDirectionalAtoms
1503  
1504    function FF_RequiresPrepairCalc() result(doesit)
# Line 1350 | Line 1508 | contains
1508  
1509    function FF_RequiresPostpairCalc() result(doesit)
1510      logical :: doesit
1511 <    doesit = FF_uses_RF
1511 >    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1512    end function FF_RequiresPostpairCalc
1513  
1514   #ifdef PROFILE

Diff Legend

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