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

Comparing trunk/OOPSE-3.0/src/UseTheForce/doForces.F90 (file contents):
Revision 2231 by chrisfen, Wed May 18 18:31:40 2005 UTC vs.
Revision 2361 by gezelter, Wed Oct 12 21:00:50 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.18 2005-05-18 18:31:40 chrisfen Exp $, $Date: 2005-05-18 18:31:40 $, $Name: not supported by cvs2svn $, $Revision: 1.18 $
48 > !! @version $Id: doForces.F90,v 1.57 2005-10-12 21:00:45 gezelter Exp $, $Date: 2005-10-12 21:00:45 $, $Name: not supported by cvs2svn $, $Revision: 1.57 $
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  
80 –  logical, save :: haveRlist = .false.
84    logical, save :: haveNeighborList = .false.
85    logical, save :: haveSIMvariables = .false.
83 –  logical, save :: havePropertyMap = .false.
86    logical, save :: haveSaneForceField = .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
87 –  logical, save :: FF_uses_LennardJones
88 –  logical, save :: FF_uses_Electrostatics
89 –  logical, save :: FF_uses_Charges
93    logical, save :: FF_uses_Dipoles
91 –  logical, save :: FF_uses_Quadrupoles
92 –  logical, save :: FF_uses_Sticky
93 –  logical, save :: FF_uses_StickyPower
94    logical, save :: FF_uses_GayBerne
95    logical, save :: FF_uses_EAM
96 –  logical, save :: FF_uses_Shapes
97 –  logical, save :: FF_uses_FLARB
98 –  logical, save :: FF_uses_RF
96  
97    logical, save :: SIM_uses_DirectionalAtoms
101 –  logical, save :: SIM_uses_LennardJones
102 –  logical, save :: SIM_uses_Electrostatics
103 –  logical, save :: SIM_uses_Charges
104 –  logical, save :: SIM_uses_Dipoles
105 –  logical, save :: SIM_uses_Quadrupoles
106 –  logical, save :: SIM_uses_Sticky
107 –  logical, save :: SIM_uses_StickyPower
108 –  logical, save :: SIM_uses_GayBerne
98    logical, save :: SIM_uses_EAM
110 –  logical, save :: SIM_uses_Shapes
111 –  logical, save :: SIM_uses_FLARB
112 –  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
116 –  logical, save :: SIM_uses_molecular_cutoffs
102  
103 <  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
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 127 | Line 119 | module doForces
119    real :: forceTimeInitial, forceTimeFinal
120    integer :: nLoops
121   #endif
122 +  
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 <  type :: Properties
131 <     logical :: is_Directional   = .false.
133 <     logical :: is_LennardJones  = .false.
134 <     logical :: is_Electrostatic = .false.
135 <     logical :: is_Charge        = .false.
136 <     logical :: is_Dipole        = .false.
137 <     logical :: is_Quadrupole    = .false.
138 <     logical :: is_Sticky        = .false.
139 <     logical :: is_StickyPower   = .false.
140 <     logical :: is_GayBerne      = .false.
141 <     logical :: is_EAM           = .false.
142 <     logical :: is_Shape         = .false.
143 <     logical :: is_FLARB         = .false.
144 <  end type Properties
130 >  integer, dimension(:), allocatable, target :: groupToGtypeRow
131 >  integer, dimension(:), pointer :: groupToGtypeCol => null()
132  
133 <  type(Properties), dimension(:),allocatable :: PropertyMap
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 <  subroutine setRlistDF( this_rlist )
151 <
152 <    real(kind=dp) :: this_rlist
153 <
154 <    rlist = this_rlist
155 <    rlistsq = rlist * rlist
156 <
157 <    haveRlist = .true.
158 <
159 <  end subroutine setRlistDF
160 <
161 <  subroutine createPropertyMap(status)
148 >  subroutine createInteractionHash(status)
149      integer :: nAtypes
150 <    integer :: status
150 >    integer, intent(out) :: status
151      integer :: i
152 <    logical :: thisProperty
153 <    real (kind=DP) :: thisDPproperty
152 >    integer :: j
153 >    integer :: iHash
154 >    !! Test Types
155 >    logical :: i_is_LJ
156 >    logical :: i_is_Elect
157 >    logical :: i_is_Sticky
158 >    logical :: i_is_StickyP
159 >    logical :: i_is_GB
160 >    logical :: i_is_EAM
161 >    logical :: i_is_Shape
162 >    logical :: j_is_LJ
163 >    logical :: j_is_Elect
164 >    logical :: j_is_Sticky
165 >    logical :: j_is_StickyP
166 >    logical :: j_is_GB
167 >    logical :: j_is_EAM
168 >    logical :: j_is_Shape
169 >    real(kind=dp) :: myRcut
170  
171 <    status = 0
171 >    status = 0  
172  
173 +    if (.not. associated(atypes)) then
174 +       call handleError("atype", "atypes was not present before call of createInteractionHash!")
175 +       status = -1
176 +       return
177 +    endif
178 +    
179      nAtypes = getSize(atypes)
180 <
180 >    
181      if (nAtypes == 0) then
182         status = -1
183         return
184      end if
185  
186 <    if (.not. allocated(PropertyMap)) then
187 <       allocate(PropertyMap(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_Directional", thisProperty)
202 <       PropertyMap(i)%is_Directional = thisProperty
201 >       call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
202 >       call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
203 >       call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
204 >       call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
205 >       call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
206 >       call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
207 >       call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
208  
209 <       call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
186 <       PropertyMap(i)%is_LennardJones = thisProperty
209 >       do j = i, nAtypes
210  
211 <       call getElementProperty(atypes, i, "is_Electrostatic", thisProperty)
212 <       PropertyMap(i)%is_Electrostatic = thisProperty
211 >          iHash = 0
212 >          myRcut = 0.0_dp
213  
214 <       call getElementProperty(atypes, i, "is_Charge", thisProperty)
215 <       PropertyMap(i)%is_Charge = thisProperty
214 >          call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ)
215 >          call getElementProperty(atypes, j, "is_Electrostatic", j_is_Elect)
216 >          call getElementProperty(atypes, j, "is_Sticky", j_is_Sticky)
217 >          call getElementProperty(atypes, j, "is_StickyPower", j_is_StickyP)
218 >          call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
219 >          call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
220 >          call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
221  
222 <       call getElementProperty(atypes, i, "is_Dipole", thisProperty)
223 <       PropertyMap(i)%is_Dipole = thisProperty
222 >          if (i_is_LJ .and. j_is_LJ) then
223 >             iHash = ior(iHash, LJ_PAIR)            
224 >          endif
225 >          
226 >          if (i_is_Elect .and. j_is_Elect) then
227 >             iHash = ior(iHash, ELECTROSTATIC_PAIR)
228 >          endif
229 >          
230 >          if (i_is_Sticky .and. j_is_Sticky) then
231 >             iHash = ior(iHash, STICKY_PAIR)
232 >          endif
233  
234 <       call getElementProperty(atypes, i, "is_Quadrupole", thisProperty)
235 <       PropertyMap(i)%is_Quadrupole = thisProperty
234 >          if (i_is_StickyP .and. j_is_StickyP) then
235 >             iHash = ior(iHash, STICKYPOWER_PAIR)
236 >          endif
237  
238 <       call getElementProperty(atypes, i, "is_Sticky", thisProperty)
239 <       PropertyMap(i)%is_Sticky = thisProperty
240 <      
203 <       call getElementProperty(atypes, i, "is_StickyPower", thisProperty)
204 <       PropertyMap(i)%is_StickyPower = thisProperty
238 >          if (i_is_EAM .and. j_is_EAM) then
239 >             iHash = ior(iHash, EAM_PAIR)
240 >          endif
241  
242 <       call getElementProperty(atypes, i, "is_GayBerne", thisProperty)
243 <       PropertyMap(i)%is_GayBerne = thisProperty
242 >          if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
243 >          if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
244 >          if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
245  
246 <       call getElementProperty(atypes, i, "is_EAM", thisProperty)
247 <       PropertyMap(i)%is_EAM = thisProperty
246 >          if (i_is_Shape .and. j_is_Shape) iHash = ior(iHash, SHAPE_PAIR)
247 >          if (i_is_Shape .and. j_is_LJ) iHash = ior(iHash, SHAPE_LJ)
248 >          if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
249  
212 –       call getElementProperty(atypes, i, "is_Shape", thisProperty)
213 –       PropertyMap(i)%is_Shape = thisProperty
250  
251 <       call getElementProperty(atypes, i, "is_FLARB", thisProperty)
252 <       PropertyMap(i)%is_FLARB = thisProperty
251 >          InteractionHash(i,j) = iHash
252 >          InteractionHash(j,i) = iHash
253 >
254 >       end do
255 >
256      end do
257  
258 <    havePropertyMap = .true.
258 >    haveInteractionHash = .true.
259 >  end subroutine createInteractionHash
260  
261 <  end subroutine createPropertyMap
261 >  subroutine createGtypeCutoffMap(stat)
262  
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 +    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 +
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 +          
342 +          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
343 +             biggestAtypeCutoff = atypeMaxCutoff(i)
344 +          endif
345 +
346 +       endif
347 +    enddo
348 +  
349 +
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()
225 –    SIM_uses_LennardJones = SimUsesLennardJones()
226 –    SIM_uses_Electrostatics = SimUsesElectrostatics()
227 –    SIM_uses_Charges = SimUsesCharges()
228 –    SIM_uses_Dipoles = SimUsesDipoles()
229 –    SIM_uses_Sticky = SimUsesSticky()
230 –    SIM_uses_StickyPower = SimUsesStickyPower()
231 –    SIM_uses_GayBerne = SimUsesGayBerne()
572      SIM_uses_EAM = SimUsesEAM()
233 –    SIM_uses_Shapes = SimUsesShapes()
234 –    SIM_uses_FLARB = SimUsesFLARB()
235 –    SIM_uses_RF = SimUsesRF()
573      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
574      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
575      SIM_uses_PBC = SimUsesPBC()
# Line 249 | Line 586 | contains
586  
587      error = 0
588  
589 <    if (.not. havePropertyMap) 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 createPropertyMap(myStatus)
257 <
599 >    if (.not. haveGtypeCutoffMap) then        
600 >       myStatus = 0      
601 >       call createGtypeCutoffMap(myStatus)      
602         if (myStatus .ne. 0) then
603 <          write(default_error, *) 'createPropertyMap failed in doForces!'
603 >          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
604            error = -1
605            return
606         endif
# Line 266 | 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 295 | 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
301 <
644 >    integer, intent(in) :: thisESM
645      integer, intent(out) :: thisStat  
646      integer :: my_status, nMatches
647      integer, pointer :: MatchList(:) => null()
# Line 307 | Line 650 | contains
650      !! assume things are copacetic, unless they aren't
651      thisStat = 0
652  
653 <    !! Fortran's version of a cast:
311 <    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 317 | Line 659 | contains
659      !! interactions are used by the force field.    
660  
661      FF_uses_DirectionalAtoms = .false.
320 –    FF_uses_LennardJones = .false.
321 –    FF_uses_Electrostatics = .false.
322 –    FF_uses_Charges = .false.    
662      FF_uses_Dipoles = .false.
324 –    FF_uses_Sticky = .false.
325 –    FF_uses_StickyPower = .false.
663      FF_uses_GayBerne = .false.
664      FF_uses_EAM = .false.
328 –    FF_uses_Shapes = .false.
329 –    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  
335 –    call getMatchingElementList(atypes, "is_LennardJones", .true., &
336 –         nMatches, MatchList)
337 –    if (nMatches .gt. 0) FF_uses_LennardJones = .true.
338 –
339 –    call getMatchingElementList(atypes, "is_Electrostatic", .true., &
340 –         nMatches, MatchList)
341 –    if (nMatches .gt. 0) then
342 –       FF_uses_Electrostatics = .true.
343 –    endif
344 –
345 –    call getMatchingElementList(atypes, "is_Charge", .true., &
346 –         nMatches, MatchList)
347 –    if (nMatches .gt. 0) then
348 –       FF_uses_Charges = .true.  
349 –       FF_uses_Electrostatics = .true.
350 –    endif
351 –
670      call getMatchingElementList(atypes, "is_Dipole", .true., &
671           nMatches, MatchList)
672 <    if (nMatches .gt. 0) then
355 <       FF_uses_Dipoles = .true.
356 <       FF_uses_Electrostatics = .true.
357 <       FF_uses_DirectionalAtoms = .true.
358 <    endif
359 <
360 <    call getMatchingElementList(atypes, "is_Quadrupole", .true., &
361 <         nMatches, MatchList)
362 <    if (nMatches .gt. 0) then
363 <       FF_uses_Quadrupoles = .true.
364 <       FF_uses_Electrostatics = .true.
365 <       FF_uses_DirectionalAtoms = .true.
366 <    endif
367 <
368 <    call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
369 <         MatchList)
370 <    if (nMatches .gt. 0) then
371 <       FF_uses_Sticky = .true.
372 <       FF_uses_DirectionalAtoms = .true.
373 <    endif
374 <
375 <    call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
376 <         MatchList)
377 <    if (nMatches .gt. 0) then
378 <       FF_uses_StickyPower = .true.
379 <       FF_uses_DirectionalAtoms = .true.
380 <    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
385 <       FF_uses_GayBerne = .true.
386 <       FF_uses_DirectionalAtoms = .true.
387 <    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  
392 –    call getMatchingElementList(atypes, "is_Shape", .true., &
393 –         nMatches, MatchList)
394 –    if (nMatches .gt. 0) then
395 –       FF_uses_Shapes = .true.
396 –       FF_uses_DirectionalAtoms = .true.
397 –    endif
681  
399 –    call getMatchingElementList(atypes, "is_FLARB", .true., &
400 –         nMatches, MatchList)
401 –    if (nMatches .gt. 0) FF_uses_FLARB = .true.
402 –
403 –    !! 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 419 | Line 697 | contains
697         endif
698      endif
699  
422 –    !sticky module does not contain check_sticky_FF anymore
423 –    !if (FF_uses_sticky) then
424 –    !   call check_sticky_FF(my_status)
425 –    !   if (my_status /= 0) then
426 –    !      thisStat = -1
427 –    !      haveSaneForceField = .false.
428 –    !      return
429 –    !   end if
430 –    !endif
431 –
700      if (FF_uses_EAM) then
701         call init_EAM_FF(my_status)
702         if (my_status /= 0) then
# Line 448 | Line 716 | contains
716         endif
717      endif
718  
451 –    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
452 –    endif
453 –
719      if (.not. haveNeighborList) then
720         !! Create neighbor lists
721         call expandNeighborList(nLocal, my_status)
# Line 484 | 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 514 | Line 779 | contains
779      integer :: localError
780      integer :: propPack_i, propPack_j
781      integer :: loopStart, loopEnd, loop
782 +    integer :: iHash
783 +  
784  
518 –    real(kind=dp) :: listSkin = 1.0  
519 –
785      !! initialize local variables  
786  
787   #ifdef IS_MPI
# Line 633 | 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 < rlistsq) then
910 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
911                  if (update_nlist) then
912                     nlist = nlist + 1
913  
# Line 761 | Line 1028 | contains
1028                  endif
1029               end if
1030            enddo
1031 +
1032         enddo outer
1033  
1034         if (update_nlist) then
# Line 820 | 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 840 | 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 858 | Line 1130 | contains
1130   #else
1131               me_i = atid(i)
1132   #endif
1133 +             iHash = InteractionHash(me_i,me_j)
1134 +            
1135 +             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1136  
862 –             if (PropertyMap(me_i)%is_Dipole) then
863 –
1137                  mu_i = getDipoleMoment(me_i)
1138  
1139                  !! The reaction field needs to include a self contribution
# Line 870 | 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 884 | 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 910 | 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 923 | Line 1198 | contains
1198      real ( kind = dp ), intent(inout) :: rijsq
1199      real ( kind = dp )                :: r
1200      real ( kind = dp ), intent(inout) :: d(3)
926 –    real ( kind = dp ) :: ebalance
1201      integer :: me_i, me_j
1202  
1203 +    integer :: iHash
1204 +
1205      r = sqrt(rijsq)
1206      vpair = 0.0d0
1207      fpair(1:3) = 0.0d0
# Line 938 | Line 1214 | contains
1214      me_j = atid(j)
1215   #endif
1216  
1217 <    !    write(*,*) i, j, me_i, me_j
942 <
943 <    if (FF_uses_LennardJones .and. SIM_uses_LennardJones) then
1217 >    iHash = InteractionHash(me_i, me_j)
1218  
1219 <       if ( PropertyMap(me_i)%is_LennardJones .and. &
1220 <            PropertyMap(me_j)%is_LennardJones ) then
1221 <          call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
948 <       endif
949 <
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 (FF_uses_Electrostatics .and. SIM_uses_Electrostatics) then
1224 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1225 >       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1226 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1227  
1228 <       if (PropertyMap(me_i)%is_Electrostatic .and. &
955 <            PropertyMap(me_j)%is_Electrostatic) then
956 <          call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
957 <               pot, eFrame, f, t, do_pot, ebalance)
958 <       endif
1228 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1229  
1230 <       if (FF_uses_dipoles .and. SIM_uses_dipoles) then      
1231 <          if ( PropertyMap(me_i)%is_Dipole .and. &
1232 <               PropertyMap(me_j)%is_Dipole) then
963 <             if (FF_uses_RF .and. SIM_uses_RF) then
964 <                call accumulate_rf(i, j, r, eFrame, sw)
965 <                call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
966 <             endif
967 <          endif
1230 >          ! CHECK ME (RF needs to know about all electrostatic types)
1231 >          call accumulate_rf(i, j, r, eFrame, sw)
1232 >          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1233         endif
1234 +
1235      endif
1236  
1237 +    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1238 +       call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1239 +            pot(HB_POT), A, f, t, do_pot)
1240 +    endif
1241  
1242 <    if (FF_uses_Sticky .and. SIM_uses_sticky) 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(HB_POT), A, f, t, do_pot)
1245 >    endif
1246  
1247 <       if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
1248 <          call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1249 <               pot, A, f, t, do_pot)
977 <       endif
978 <
1247 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1248 >       call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1249 >            pot(VDW_POT), A, f, t, do_pot)
1250      endif
980 –
981 –    if (FF_uses_StickyPower .and. SIM_uses_stickypower) then
982 –       if ( PropertyMap(me_i)%is_StickyPower .and. &
983 –            PropertyMap(me_j)%is_StickyPower) then
984 –          call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
985 –               pot, A, f, t, do_pot, ebalance)
986 –       endif
987 –    endif
1251      
1252 <    if (FF_uses_GayBerne .and. SIM_uses_GayBerne) then
1253 <
1254 <       if ( PropertyMap(me_i)%is_GayBerne .and. &
992 <            PropertyMap(me_j)%is_GayBerne) then
993 <          call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
994 <               pot, A, f, t, do_pot)
995 <       endif
996 <
1252 >    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1253 > !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1254 > !           pot(VDW_POT), A, f, t, do_pot)
1255      endif
1256  
1257 <    if (FF_uses_EAM .and. SIM_uses_EAM) then
1258 <
1259 <       if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
1002 <          call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1003 <               do_pot)
1004 <       endif
1005 <
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 <
1263 <    !    write(*,*) PropertyMap(me_i)%is_Shape,PropertyMap(me_j)%is_Shape
1264 <
1011 <    if (FF_uses_Shapes .and. SIM_uses_Shapes) then
1012 <       if ( PropertyMap(me_i)%is_Shape .and. &
1013 <            PropertyMap(me_j)%is_Shape ) then
1014 <          call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1015 <               pot, A, f, t, do_pot)
1016 <       endif
1017 <       if ( (PropertyMap(me_i)%is_Shape .and. &
1018 <            PropertyMap(me_j)%is_LennardJones) .or. &
1019 <            (PropertyMap(me_i)%is_LennardJones .and. &
1020 <            PropertyMap(me_j)%is_Shape) ) then
1021 <          call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1022 <               pot, A, f, t, do_pot)
1023 <       endif
1262 >    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1263 >       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1264 >            pot(VDW_POT), A, f, t, do_pot)
1265      endif
1266  
1267 +    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1268 +       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1269 +            pot(VDW_POT), A, f, t, do_pot)
1270 +    endif
1271 +    
1272    end subroutine do_pair
1273  
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 1040 | Line 1287 | contains
1287      real ( kind = dp )                :: r, rc
1288      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1289  
1290 <    logical :: is_EAM_i, is_EAM_j
1290 >    integer :: me_i, me_j, iHash
1291  
1045 –    integer :: me_i, me_j
1046 –
1047 –
1292      r = sqrt(rijsq)
1049 –    if (SIM_uses_molecular_cutoffs) then
1050 –       rc = sqrt(rcijsq)
1051 –    else
1052 –       rc = r
1053 –    endif
1293  
1055 –
1294   #ifdef IS_MPI  
1295      me_i = atid_row(i)
1296      me_j = atid_col(j)  
# Line 1061 | Line 1299 | contains
1299      me_j = atid(j)  
1300   #endif
1301  
1302 <    if (FF_uses_EAM .and. SIM_uses_EAM) then
1302 >    iHash = InteractionHash(me_i, me_j)
1303  
1304 <       if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1304 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1305              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1068 –
1306      endif
1307 <
1307 >    
1308    end subroutine do_prepair
1309  
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 1261 | Line 1498 | contains
1498  
1499    function FF_UsesDirectionalAtoms() result(doesit)
1500      logical :: doesit
1501 <    doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1265 <         FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1266 <         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 1273 | 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