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 2259 by gezelter, Mon Jun 27 21:01:36 2005 UTC vs.
Revision 2381 by chrisfen, Tue Oct 18 15:01:42 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.20 2005-06-27 21:01:30 gezelter Exp $, $Date: 2005-06-27 21:01:30 $, $Name: not supported by cvs2svn $, $Revision: 1.20 $
48 > !! @version $Id: doForces.F90,v 1.60 2005-10-18 15:01:42 chrisfen Exp $, $Date: 2005-10-18 15:01:42 $, $Name: not supported by cvs2svn $, $Revision: 1.60 $
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
62 <  use gb_pair
61 >  use gayberne
62    use shapes
63    use vector_class
64    use eam
# Line 73 | Line 72 | module doForces
72  
73   #define __FORTRAN90
74   #include "UseTheForce/fSwitchingFunction.h"
75 + #include "UseTheForce/fCutoffPolicy.h"
76   #include "UseTheForce/DarkSide/fInteractionMap.h"
77 + #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
78  
79 +
80    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
81    INTEGER, PARAMETER:: PAIR_LOOP    = 2
82  
81 –  logical, save :: haveRlist = .false.
83    logical, save :: haveNeighborList = .false.
84    logical, save :: haveSIMvariables = .false.
84 –  logical, save :: havePropertyMap = .false.
85    logical, save :: haveSaneForceField = .false.
86 +  logical, save :: haveInteractionHash = .false.
87 +  logical, save :: haveGtypeCutoffMap = .false.
88 +  logical, save :: haveDefaultCutoffs = .false.
89 +  logical, save :: haveRlist = .false.
90  
91    logical, save :: FF_uses_DirectionalAtoms
88 –  logical, save :: FF_uses_LennardJones
89 –  logical, save :: FF_uses_Electrostatics
90 –  logical, save :: FF_uses_Charges
92    logical, save :: FF_uses_Dipoles
92 –  logical, save :: FF_uses_Quadrupoles
93 –  logical, save :: FF_uses_Sticky
94 –  logical, save :: FF_uses_StickyPower
93    logical, save :: FF_uses_GayBerne
94    logical, save :: FF_uses_EAM
97 –  logical, save :: FF_uses_Shapes
98 –  logical, save :: FF_uses_FLARB
99 –  logical, save :: FF_uses_RF
95  
96    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
97    logical, save :: SIM_uses_EAM
111 –  logical, save :: SIM_uses_Shapes
112 –  logical, save :: SIM_uses_FLARB
113 –  logical, save :: SIM_uses_RF
98    logical, save :: SIM_requires_postpair_calc
99    logical, save :: SIM_requires_prepair_calc
100    logical, save :: SIM_uses_PBC
117 –  logical, save :: SIM_uses_molecular_cutoffs
101  
102 <  real(kind=dp), save :: rlist, rlistsq
102 >  integer, save :: electrostaticSummationMethod
103  
104    public :: init_FF
105 +  public :: setDefaultCutoffs
106    public :: do_force_loop
107 <  public :: setRlistDF
107 >  public :: createInteractionHash
108 >  public :: createGtypeCutoffMap
109 >  public :: getStickyCut
110 >  public :: getStickyPowerCut
111 >  public :: getGayBerneCut
112 >  public :: getEAMCut
113 >  public :: getShapeCut
114  
115   #ifdef PROFILE
116    public :: getforcetime
# Line 128 | Line 118 | module doForces
118    real :: forceTimeInitial, forceTimeFinal
119    integer :: nLoops
120   #endif
121 +  
122 +  !! Variables for cutoff mapping and interaction mapping
123 +  ! Bit hash to determine pair-pair interactions.
124 +  integer, dimension(:,:), allocatable :: InteractionHash
125 +  real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
126 +  real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
127 +  real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
128  
129 <  type :: Properties
130 <     logical :: is_Directional   = .false.
134 <     logical :: is_LennardJones  = .false.
135 <     logical :: is_Electrostatic = .false.
136 <     logical :: is_Charge        = .false.
137 <     logical :: is_Dipole        = .false.
138 <     logical :: is_Quadrupole    = .false.
139 <     logical :: is_Sticky        = .false.
140 <     logical :: is_StickyPower   = .false.
141 <     logical :: is_GayBerne      = .false.
142 <     logical :: is_EAM           = .false.
143 <     logical :: is_Shape         = .false.
144 <     logical :: is_FLARB         = .false.
145 <  end type Properties
129 >  integer, dimension(:), allocatable, target :: groupToGtypeRow
130 >  integer, dimension(:), pointer :: groupToGtypeCol => null()
131  
132 <  type(Properties), dimension(:),allocatable :: PropertyMap
132 >  real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
133 >  real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
134 >  type ::gtypeCutoffs
135 >     real(kind=dp) :: rcut
136 >     real(kind=dp) :: rcutsq
137 >     real(kind=dp) :: rlistsq
138 >  end type gtypeCutoffs
139 >  type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
140  
141 +  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
142 +  real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
143 +  real(kind=dp),save :: listSkin
144 +  
145   contains
146  
147 <  subroutine setRlistDF( this_rlist )
152 <
153 <    real(kind=dp) :: this_rlist
154 <
155 <    rlist = this_rlist
156 <    rlistsq = rlist * rlist
157 <
158 <    haveRlist = .true.
159 <
160 <  end subroutine setRlistDF
161 <
162 <  subroutine createPropertyMap(status)
147 >  subroutine createInteractionHash(status)
148      integer :: nAtypes
149 <    integer :: status
149 >    integer, intent(out) :: status
150      integer :: i
151 <    logical :: thisProperty
152 <    real (kind=DP) :: thisDPproperty
151 >    integer :: j
152 >    integer :: iHash
153 >    !! Test Types
154 >    logical :: i_is_LJ
155 >    logical :: i_is_Elect
156 >    logical :: i_is_Sticky
157 >    logical :: i_is_StickyP
158 >    logical :: i_is_GB
159 >    logical :: i_is_EAM
160 >    logical :: i_is_Shape
161 >    logical :: j_is_LJ
162 >    logical :: j_is_Elect
163 >    logical :: j_is_Sticky
164 >    logical :: j_is_StickyP
165 >    logical :: j_is_GB
166 >    logical :: j_is_EAM
167 >    logical :: j_is_Shape
168 >    real(kind=dp) :: myRcut
169  
170 <    status = 0
170 >    status = 0  
171  
172 +    if (.not. associated(atypes)) then
173 +       call handleError("atype", "atypes was not present before call of createInteractionHash!")
174 +       status = -1
175 +       return
176 +    endif
177 +    
178      nAtypes = getSize(atypes)
179 <
179 >    
180      if (nAtypes == 0) then
181         status = -1
182         return
183      end if
184  
185 <    if (.not. allocated(PropertyMap)) then
186 <       allocate(PropertyMap(nAtypes))
185 >    if (.not. allocated(InteractionHash)) then
186 >       allocate(InteractionHash(nAtypes,nAtypes))
187 >    else
188 >       deallocate(InteractionHash)
189 >       allocate(InteractionHash(nAtypes,nAtypes))
190      endif
191  
192 +    if (.not. allocated(atypeMaxCutoff)) then
193 +       allocate(atypeMaxCutoff(nAtypes))
194 +    else
195 +       deallocate(atypeMaxCutoff)
196 +       allocate(atypeMaxCutoff(nAtypes))
197 +    endif
198 +        
199      do i = 1, nAtypes
200 <       call getElementProperty(atypes, i, "is_Directional", thisProperty)
201 <       PropertyMap(i)%is_Directional = thisProperty
200 >       call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
201 >       call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
202 >       call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
203 >       call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
204 >       call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
205 >       call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
206 >       call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
207  
208 <       call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
187 <       PropertyMap(i)%is_LennardJones = thisProperty
208 >       do j = i, nAtypes
209  
210 <       call getElementProperty(atypes, i, "is_Electrostatic", thisProperty)
211 <       PropertyMap(i)%is_Electrostatic = thisProperty
210 >          iHash = 0
211 >          myRcut = 0.0_dp
212  
213 <       call getElementProperty(atypes, i, "is_Charge", thisProperty)
214 <       PropertyMap(i)%is_Charge = thisProperty
213 >          call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ)
214 >          call getElementProperty(atypes, j, "is_Electrostatic", j_is_Elect)
215 >          call getElementProperty(atypes, j, "is_Sticky", j_is_Sticky)
216 >          call getElementProperty(atypes, j, "is_StickyPower", j_is_StickyP)
217 >          call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
218 >          call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
219 >          call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
220  
221 <       call getElementProperty(atypes, i, "is_Dipole", thisProperty)
222 <       PropertyMap(i)%is_Dipole = thisProperty
221 >          if (i_is_LJ .and. j_is_LJ) then
222 >             iHash = ior(iHash, LJ_PAIR)            
223 >          endif
224 >          
225 >          if (i_is_Elect .and. j_is_Elect) then
226 >             iHash = ior(iHash, ELECTROSTATIC_PAIR)
227 >          endif
228 >          
229 >          if (i_is_Sticky .and. j_is_Sticky) then
230 >             iHash = ior(iHash, STICKY_PAIR)
231 >          endif
232  
233 <       call getElementProperty(atypes, i, "is_Quadrupole", thisProperty)
234 <       PropertyMap(i)%is_Quadrupole = thisProperty
233 >          if (i_is_StickyP .and. j_is_StickyP) then
234 >             iHash = ior(iHash, STICKYPOWER_PAIR)
235 >          endif
236  
237 <       call getElementProperty(atypes, i, "is_Sticky", thisProperty)
238 <       PropertyMap(i)%is_Sticky = thisProperty
239 <      
204 <       call getElementProperty(atypes, i, "is_StickyPower", thisProperty)
205 <       PropertyMap(i)%is_StickyPower = thisProperty
237 >          if (i_is_EAM .and. j_is_EAM) then
238 >             iHash = ior(iHash, EAM_PAIR)
239 >          endif
240  
241 <       call getElementProperty(atypes, i, "is_GayBerne", thisProperty)
242 <       PropertyMap(i)%is_GayBerne = thisProperty
241 >          if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
242 >          if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
243 >          if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
244  
245 <       call getElementProperty(atypes, i, "is_EAM", thisProperty)
246 <       PropertyMap(i)%is_EAM = thisProperty
245 >          if (i_is_Shape .and. j_is_Shape) iHash = ior(iHash, SHAPE_PAIR)
246 >          if (i_is_Shape .and. j_is_LJ) iHash = ior(iHash, SHAPE_LJ)
247 >          if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
248  
213 –       call getElementProperty(atypes, i, "is_Shape", thisProperty)
214 –       PropertyMap(i)%is_Shape = thisProperty
249  
250 <       call getElementProperty(atypes, i, "is_FLARB", thisProperty)
251 <       PropertyMap(i)%is_FLARB = thisProperty
218 <    end do
250 >          InteractionHash(i,j) = iHash
251 >          InteractionHash(j,i) = iHash
252  
253 <    havePropertyMap = .true.
253 >       end do
254  
255 <  end subroutine createPropertyMap
255 >    end do
256  
257 +    haveInteractionHash = .true.
258 +  end subroutine createInteractionHash
259 +
260 +  subroutine createGtypeCutoffMap(stat)
261 +
262 +    integer, intent(out), optional :: stat
263 +    logical :: i_is_LJ
264 +    logical :: i_is_Elect
265 +    logical :: i_is_Sticky
266 +    logical :: i_is_StickyP
267 +    logical :: i_is_GB
268 +    logical :: i_is_EAM
269 +    logical :: i_is_Shape
270 +    logical :: GtypeFound
271 +
272 +    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
273 +    integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
274 +    integer :: nGroupsInRow
275 +    integer :: nGroupsInCol
276 +    integer :: nGroupTypesRow,nGroupTypesCol
277 +    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin
278 +    real(kind=dp) :: biggestAtypeCutoff
279 +
280 +    stat = 0
281 +    if (.not. haveInteractionHash) then
282 +       call createInteractionHash(myStatus)      
283 +       if (myStatus .ne. 0) then
284 +          write(default_error, *) 'createInteractionHash failed in doForces!'
285 +          stat = -1
286 +          return
287 +       endif
288 +    endif
289 + #ifdef IS_MPI
290 +    nGroupsInRow = getNgroupsInRow(plan_group_row)
291 +    nGroupsInCol = getNgroupsInCol(plan_group_col)
292 + #endif
293 +    nAtypes = getSize(atypes)
294 + ! Set all of the initial cutoffs to zero.
295 +    atypeMaxCutoff = 0.0_dp
296 +    do i = 1, nAtypes
297 +       if (SimHasAtype(i)) then    
298 +          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
299 +          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
300 +          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
301 +          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
302 +          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
303 +          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
304 +          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
305 +          
306 +
307 +          if (haveDefaultCutoffs) then
308 +             atypeMaxCutoff(i) = defaultRcut
309 +          else
310 +             if (i_is_LJ) then          
311 +                thisRcut = getSigma(i) * 2.5_dp
312 +                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
313 +             endif
314 +             if (i_is_Elect) then
315 +                thisRcut = defaultRcut
316 +                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
317 +             endif
318 +             if (i_is_Sticky) then
319 +                thisRcut = getStickyCut(i)
320 +                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
321 +             endif
322 +             if (i_is_StickyP) then
323 +                thisRcut = getStickyPowerCut(i)
324 +                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
325 +             endif
326 +             if (i_is_GB) then
327 +                thisRcut = getGayBerneCut(i)
328 +                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
329 +             endif
330 +             if (i_is_EAM) then
331 +                thisRcut = getEAMCut(i)
332 +                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
333 +             endif
334 +             if (i_is_Shape) then
335 +                thisRcut = getShapeCut(i)
336 +                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
337 +             endif
338 +          endif
339 +          
340 +          
341 +          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
342 +             biggestAtypeCutoff = atypeMaxCutoff(i)
343 +          endif
344 +
345 +       endif
346 +    enddo
347 +  
348 +
349 +    
350 +    istart = 1
351 +    jstart = 1
352 + #ifdef IS_MPI
353 +    iend = nGroupsInRow
354 +    jend = nGroupsInCol
355 + #else
356 +    iend = nGroups
357 +    jend = nGroups
358 + #endif
359 +    
360 +    !! allocate the groupToGtype and gtypeMaxCutoff here.
361 +    if(.not.allocated(groupToGtypeRow)) then
362 +     !  allocate(groupToGtype(iend))
363 +       allocate(groupToGtypeRow(iend))
364 +    else
365 +       deallocate(groupToGtypeRow)
366 +       allocate(groupToGtypeRow(iend))
367 +    endif
368 +    if(.not.allocated(groupMaxCutoffRow)) then
369 +       allocate(groupMaxCutoffRow(iend))
370 +    else
371 +       deallocate(groupMaxCutoffRow)
372 +       allocate(groupMaxCutoffRow(iend))
373 +    end if
374 +
375 +    if(.not.allocated(gtypeMaxCutoffRow)) then
376 +       allocate(gtypeMaxCutoffRow(iend))
377 +    else
378 +       deallocate(gtypeMaxCutoffRow)
379 +       allocate(gtypeMaxCutoffRow(iend))
380 +    endif
381 +
382 +
383 + #ifdef IS_MPI
384 +       ! We only allocate new storage if we are in MPI because Ncol /= Nrow
385 +    if(.not.associated(groupToGtypeCol)) then
386 +       allocate(groupToGtypeCol(jend))
387 +    else
388 +       deallocate(groupToGtypeCol)
389 +       allocate(groupToGtypeCol(jend))
390 +    end if
391 +
392 +    if(.not.associated(groupToGtypeCol)) then
393 +       allocate(groupToGtypeCol(jend))
394 +    else
395 +       deallocate(groupToGtypeCol)
396 +       allocate(groupToGtypeCol(jend))
397 +    end if
398 +    if(.not.associated(gtypeMaxCutoffCol)) then
399 +       allocate(gtypeMaxCutoffCol(jend))
400 +    else
401 +       deallocate(gtypeMaxCutoffCol)      
402 +       allocate(gtypeMaxCutoffCol(jend))
403 +    end if
404 +
405 +       groupMaxCutoffCol = 0.0_dp
406 +       gtypeMaxCutoffCol = 0.0_dp
407 +
408 + #endif
409 +       groupMaxCutoffRow = 0.0_dp
410 +       gtypeMaxCutoffRow = 0.0_dp
411 +
412 +
413 +    !! first we do a single loop over the cutoff groups to find the
414 +    !! largest cutoff for any atypes present in this group.  We also
415 +    !! create gtypes at this point.
416 +    
417 +    tol = 1.0d-6
418 +    nGroupTypesRow = 0
419 +
420 +    do i = istart, iend      
421 +       n_in_i = groupStartRow(i+1) - groupStartRow(i)
422 +       groupMaxCutoffRow(i) = 0.0_dp
423 +       do ia = groupStartRow(i), groupStartRow(i+1)-1
424 +          atom1 = groupListRow(ia)
425 + #ifdef IS_MPI
426 +          me_i = atid_row(atom1)
427 + #else
428 +          me_i = atid(atom1)
429 + #endif          
430 +          if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
431 +             groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
432 +          endif          
433 +       enddo
434 +
435 +       if (nGroupTypesRow.eq.0) then
436 +          nGroupTypesRow = nGroupTypesRow + 1
437 +          gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
438 +          groupToGtypeRow(i) = nGroupTypesRow
439 +       else
440 +          GtypeFound = .false.
441 +          do g = 1, nGroupTypesRow
442 +             if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
443 +                groupToGtypeRow(i) = g
444 +                GtypeFound = .true.
445 +             endif
446 +          enddo
447 +          if (.not.GtypeFound) then            
448 +             nGroupTypesRow = nGroupTypesRow + 1
449 +             gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
450 +             groupToGtypeRow(i) = nGroupTypesRow
451 +          endif
452 +       endif
453 +    enddo    
454 +
455 + #ifdef IS_MPI
456 +    do j = jstart, jend      
457 +       n_in_j = groupStartCol(j+1) - groupStartCol(j)
458 +       groupMaxCutoffCol(j) = 0.0_dp
459 +       do ja = groupStartCol(j), groupStartCol(j+1)-1
460 +          atom1 = groupListCol(ja)
461 +
462 +          me_j = atid_col(atom1)
463 +
464 +          if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
465 +             groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
466 +          endif          
467 +       enddo
468 +
469 +       if (nGroupTypesCol.eq.0) then
470 +          nGroupTypesCol = nGroupTypesCol + 1
471 +          gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
472 +          groupToGtypeCol(j) = nGroupTypesCol
473 +       else
474 +          GtypeFound = .false.
475 +          do g = 1, nGroupTypesCol
476 +             if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
477 +                groupToGtypeCol(j) = g
478 +                GtypeFound = .true.
479 +             endif
480 +          enddo
481 +          if (.not.GtypeFound) then            
482 +             nGroupTypesCol = nGroupTypesCol + 1
483 +             gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
484 +             groupToGtypeCol(j) = nGroupTypesCol
485 +          endif
486 +       endif
487 +    enddo    
488 +
489 + #else
490 + ! Set pointers to information we just found
491 +    nGroupTypesCol = nGroupTypesRow
492 +    groupToGtypeCol => groupToGtypeRow
493 +    gtypeMaxCutoffCol => gtypeMaxCutoffRow
494 +    groupMaxCutoffCol => groupMaxCutoffRow
495 + #endif
496 +
497 +
498 +
499 +
500 +
501 +    !! allocate the gtypeCutoffMap here.
502 +    allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
503 +    !! then we do a double loop over all the group TYPES to find the cutoff
504 +    !! map between groups of two types
505 +    tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
506 +
507 +    do i = 1, nGroupTypesRow
508 +       do j = 1, nGroupTypesCol
509 +      
510 +          select case(cutoffPolicy)
511 +          case(TRADITIONAL_CUTOFF_POLICY)
512 +             thisRcut = tradRcut
513 +          case(MIX_CUTOFF_POLICY)
514 +             thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
515 +          case(MAX_CUTOFF_POLICY)
516 +             thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
517 +          case default
518 +             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
519 +             return
520 +          end select
521 +          gtypeCutoffMap(i,j)%rcut = thisRcut
522 +          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
523 +          skin = defaultRlist - defaultRcut
524 +          listSkin = skin ! set neighbor list skin thickness
525 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
526 +
527 +          ! sanity check
528 +
529 +          if (haveDefaultCutoffs) then
530 +             if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
531 +                call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
532 +             endif
533 +          endif
534 +       enddo
535 +    enddo
536 +    if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
537 +    if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
538 +    if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
539 + #ifdef IS_MPI
540 +    if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
541 +    if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
542 + #endif
543 +    groupMaxCutoffCol => null()
544 +    gtypeMaxCutoffCol => null()
545 +    
546 +    haveGtypeCutoffMap = .true.
547 +   end subroutine createGtypeCutoffMap
548 +
549 +   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
550 +     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
551 +     integer, intent(in) :: cutPolicy
552 +
553 +     defaultRcut = defRcut
554 +     defaultRsw = defRsw
555 +     defaultRlist = defRlist
556 +     cutoffPolicy = cutPolicy
557 +
558 +     haveDefaultCutoffs = .true.
559 +   end subroutine setDefaultCutoffs
560 +
561 +   subroutine setCutoffPolicy(cutPolicy)
562 +
563 +     integer, intent(in) :: cutPolicy
564 +     cutoffPolicy = cutPolicy
565 +     call createGtypeCutoffMap()
566 +   end subroutine setCutoffPolicy
567 +    
568 +    
569    subroutine setSimVariables()
570      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
226 –    SIM_uses_LennardJones = SimUsesLennardJones()
227 –    SIM_uses_Electrostatics = SimUsesElectrostatics()
228 –    SIM_uses_Charges = SimUsesCharges()
229 –    SIM_uses_Dipoles = SimUsesDipoles()
230 –    SIM_uses_Sticky = SimUsesSticky()
231 –    SIM_uses_StickyPower = SimUsesStickyPower()
232 –    SIM_uses_GayBerne = SimUsesGayBerne()
571      SIM_uses_EAM = SimUsesEAM()
234 –    SIM_uses_Shapes = SimUsesShapes()
235 –    SIM_uses_FLARB = SimUsesFLARB()
236 –    SIM_uses_RF = SimUsesRF()
572      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
573      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
574      SIM_uses_PBC = SimUsesPBC()
# Line 250 | Line 585 | contains
585  
586      error = 0
587  
588 <    if (.not. havePropertyMap) then
588 >    if (.not. haveInteractionHash) then      
589 >       myStatus = 0      
590 >       call createInteractionHash(myStatus)      
591 >       if (myStatus .ne. 0) then
592 >          write(default_error, *) 'createInteractionHash failed in doForces!'
593 >          error = -1
594 >          return
595 >       endif
596 >    endif
597  
598 <       myStatus = 0
599 <
600 <       call createPropertyMap(myStatus)
258 <
598 >    if (.not. haveGtypeCutoffMap) then        
599 >       myStatus = 0      
600 >       call createGtypeCutoffMap(myStatus)      
601         if (myStatus .ne. 0) then
602 <          write(default_error, *) 'createPropertyMap failed in doForces!'
602 >          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
603            error = -1
604            return
605         endif
# Line 267 | Line 609 | contains
609         call setSimVariables()
610      endif
611  
612 <    if (.not. haveRlist) then
613 <       write(default_error, *) 'rList has not been set in doForces!'
614 <       error = -1
615 <       return
616 <    endif
612 >  !  if (.not. haveRlist) then
613 >  !     write(default_error, *) 'rList has not been set in doForces!'
614 >  !     error = -1
615 >  !     return
616 >  !  endif
617  
618      if (.not. haveNeighborList) then
619         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 296 | Line 638 | contains
638    end subroutine doReadyCheck
639  
640  
641 <  subroutine init_FF(use_RF_c, thisStat)
641 >  subroutine init_FF(thisESM, thisStat)
642  
643 <    logical, intent(in) :: use_RF_c
302 <
643 >    integer, intent(in) :: thisESM
644      integer, intent(out) :: thisStat  
645      integer :: my_status, nMatches
646      integer, pointer :: MatchList(:) => null()
# Line 308 | Line 649 | contains
649      !! assume things are copacetic, unless they aren't
650      thisStat = 0
651  
652 <    !! Fortran's version of a cast:
312 <    FF_uses_RF = use_RF_c
652 >    electrostaticSummationMethod = thisESM
653  
654      !! init_FF is called *after* all of the atom types have been
655      !! defined in atype_module using the new_atype subroutine.
# Line 318 | Line 658 | contains
658      !! interactions are used by the force field.    
659  
660      FF_uses_DirectionalAtoms = .false.
321 –    FF_uses_LennardJones = .false.
322 –    FF_uses_Electrostatics = .false.
323 –    FF_uses_Charges = .false.    
661      FF_uses_Dipoles = .false.
325 –    FF_uses_Sticky = .false.
326 –    FF_uses_StickyPower = .false.
662      FF_uses_GayBerne = .false.
663      FF_uses_EAM = .false.
329 –    FF_uses_Shapes = .false.
330 –    FF_uses_FLARB = .false.
664  
665      call getMatchingElementList(atypes, "is_Directional", .true., &
666           nMatches, MatchList)
667      if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
668  
336 –    call getMatchingElementList(atypes, "is_LennardJones", .true., &
337 –         nMatches, MatchList)
338 –    if (nMatches .gt. 0) FF_uses_LennardJones = .true.
339 –
340 –    call getMatchingElementList(atypes, "is_Electrostatic", .true., &
341 –         nMatches, MatchList)
342 –    if (nMatches .gt. 0) then
343 –       FF_uses_Electrostatics = .true.
344 –    endif
345 –
346 –    call getMatchingElementList(atypes, "is_Charge", .true., &
347 –         nMatches, MatchList)
348 –    if (nMatches .gt. 0) then
349 –       FF_uses_Charges = .true.  
350 –       FF_uses_Electrostatics = .true.
351 –    endif
352 –
669      call getMatchingElementList(atypes, "is_Dipole", .true., &
670           nMatches, MatchList)
671 <    if (nMatches .gt. 0) then
356 <       FF_uses_Dipoles = .true.
357 <       FF_uses_Electrostatics = .true.
358 <       FF_uses_DirectionalAtoms = .true.
359 <    endif
360 <
361 <    call getMatchingElementList(atypes, "is_Quadrupole", .true., &
362 <         nMatches, MatchList)
363 <    if (nMatches .gt. 0) then
364 <       FF_uses_Quadrupoles = .true.
365 <       FF_uses_Electrostatics = .true.
366 <       FF_uses_DirectionalAtoms = .true.
367 <    endif
368 <
369 <    call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
370 <         MatchList)
371 <    if (nMatches .gt. 0) then
372 <       FF_uses_Sticky = .true.
373 <       FF_uses_DirectionalAtoms = .true.
374 <    endif
375 <
376 <    call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
377 <         MatchList)
378 <    if (nMatches .gt. 0) then
379 <       FF_uses_StickyPower = .true.
380 <       FF_uses_DirectionalAtoms = .true.
381 <    endif
671 >    if (nMatches .gt. 0) FF_uses_Dipoles = .true.
672      
673      call getMatchingElementList(atypes, "is_GayBerne", .true., &
674           nMatches, MatchList)
675 <    if (nMatches .gt. 0) then
386 <       FF_uses_GayBerne = .true.
387 <       FF_uses_DirectionalAtoms = .true.
388 <    endif
675 >    if (nMatches .gt. 0) FF_uses_GayBerne = .true.
676  
677      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
678      if (nMatches .gt. 0) FF_uses_EAM = .true.
679  
393 –    call getMatchingElementList(atypes, "is_Shape", .true., &
394 –         nMatches, MatchList)
395 –    if (nMatches .gt. 0) then
396 –       FF_uses_Shapes = .true.
397 –       FF_uses_DirectionalAtoms = .true.
398 –    endif
680  
400 –    call getMatchingElementList(atypes, "is_FLARB", .true., &
401 –         nMatches, MatchList)
402 –    if (nMatches .gt. 0) FF_uses_FLARB = .true.
403 –
404 –    !! Assume sanity (for the sake of argument)
681      haveSaneForceField = .true.
682  
407 –    !! check to make sure the FF_uses_RF setting makes sense
408 –
409 –    if (FF_uses_dipoles) then
410 –       if (FF_uses_RF) then
411 –          dielect = getDielect()
412 –          call initialize_rf(dielect)
413 –       endif
414 –    else
415 –       if (FF_uses_RF) then          
416 –          write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
417 –          thisStat = -1
418 –          haveSaneForceField = .false.
419 –          return
420 –       endif
421 –    endif
422 –
423 –    !sticky module does not contain check_sticky_FF anymore
424 –    !if (FF_uses_sticky) then
425 –    !   call check_sticky_FF(my_status)
426 –    !   if (my_status /= 0) then
427 –    !      thisStat = -1
428 –    !      haveSaneForceField = .false.
429 –    !      return
430 –    !   end if
431 –    !endif
432 –
683      if (FF_uses_EAM) then
684         call init_EAM_FF(my_status)
685         if (my_status /= 0) then
# Line 440 | Line 690 | contains
690         end if
691      endif
692  
443 –    if (FF_uses_GayBerne) then
444 –       call check_gb_pair_FF(my_status)
445 –       if (my_status .ne. 0) then
446 –          thisStat = -1
447 –          haveSaneForceField = .false.
448 –          return
449 –       endif
450 –    endif
451 –
452 –    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
453 –    endif
454 –
693      if (.not. haveNeighborList) then
694         !! Create neighbor lists
695         call expandNeighborList(nLocal, my_status)
# Line 485 | Line 723 | contains
723  
724      !! Stress Tensor
725      real( kind = dp), dimension(9) :: tau  
726 <    real ( kind = dp ) :: pot
726 >    real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
727      logical ( kind = 2) :: do_pot_c, do_stress_c
728      logical :: do_pot
729      logical :: do_stress
730      logical :: in_switching_region
731   #ifdef IS_MPI
732 <    real( kind = DP ) :: pot_local
732 >    real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
733      integer :: nAtomsInRow
734      integer :: nAtomsInCol
735      integer :: nprocs
# Line 515 | Line 753 | contains
753      integer :: localError
754      integer :: propPack_i, propPack_j
755      integer :: loopStart, loopEnd, loop
756 +    integer :: iHash
757 +  
758  
519 –    real(kind=dp) :: listSkin = 1.0  
520 –
759      !! initialize local variables  
760  
761   #ifdef IS_MPI
# Line 634 | Line 872 | contains
872               endif
873  
874   #ifdef IS_MPI
875 +             me_j = atid_col(j)
876               call get_interatomic_vector(q_group_Row(:,i), &
877                    q_group_Col(:,j), d_grp, rgrpsq)
878   #else
879 +             me_j = atid(j)
880               call get_interatomic_vector(q_group(:,i), &
881                    q_group(:,j), d_grp, rgrpsq)
882 < #endif
882 > #endif      
883  
884 <             if (rgrpsq < rlistsq) then
884 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
885                  if (update_nlist) then
886                     nlist = nlist + 1
887  
# Line 762 | Line 1002 | contains
1002                  endif
1003               end if
1004            enddo
1005 +
1006         enddo outer
1007  
1008         if (update_nlist) then
# Line 821 | Line 1062 | contains
1062  
1063      if (do_pot) then
1064         ! scatter/gather pot_row into the members of my column
1065 <       call scatter(pot_Row, pot_Temp, plan_atom_row)
1066 <
1065 >       do i = 1,LR_POT_TYPES
1066 >          call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1067 >       end do
1068         ! scatter/gather pot_local into all other procs
1069         ! add resultant to get total pot
1070         do i = 1, nlocal
1071 <          pot_local = pot_local + pot_Temp(i)
1071 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1072 >               + pot_Temp(1:LR_POT_TYPES,i)
1073         enddo
1074  
1075         pot_Temp = 0.0_DP
1076 <
1077 <       call scatter(pot_Col, pot_Temp, plan_atom_col)
1076 >       do i = 1,LR_POT_TYPES
1077 >          call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1078 >       end do
1079         do i = 1, nlocal
1080 <          pot_local = pot_local + pot_Temp(i)
1080 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1081 >               + pot_Temp(1:LR_POT_TYPES,i)
1082         enddo
1083  
1084      endif
# Line 841 | Line 1086 | contains
1086  
1087      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1088  
1089 <       if (FF_uses_RF .and. SIM_uses_RF) then
1089 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1090  
1091   #ifdef IS_MPI
1092            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 859 | Line 1104 | contains
1104   #else
1105               me_i = atid(i)
1106   #endif
1107 <
1108 <             if (PropertyMap(me_i)%is_Dipole) then
1107 >             iHash = InteractionHash(me_i,me_j)
1108 >            
1109 >             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1110  
1111                  mu_i = getDipoleMoment(me_i)
1112  
# Line 871 | Line 1117 | contains
1117                  !! potential and torques:
1118                  call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1119   #ifdef IS_MPI
1120 <                pot_local = pot_local + rfpot
1120 >                pot_local(ELECTROSTATIC_POT) = pot_local(ELECTROSTATIC_POT) + rfpot
1121   #else
1122 <                pot = pot + rfpot
1122 >                pot(ELECTROSTATIC_POT) = pot(ELECTROSTATIC_POT) + rfpot
1123  
1124   #endif
1125               endif
# Line 885 | Line 1131 | contains
1131   #ifdef IS_MPI
1132  
1133      if (do_pot) then
1134 <       pot = pot + pot_local
1134 >       pot(1:LR_POT_TYPES) = pot(1:LR_POT_TYPES) &
1135 >            + pot_local(1:LR_POT_TYPES)
1136         !! we assume the c code will do the allreduce to get the total potential
1137         !! we could do it right here if we needed to...
1138      endif
# Line 911 | Line 1158 | contains
1158    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1159         eFrame, A, f, t, pot, vpair, fpair)
1160  
1161 <    real( kind = dp ) :: pot, vpair, sw
1161 >    real( kind = dp ) :: vpair, sw
1162 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1163      real( kind = dp ), dimension(3) :: fpair
1164      real( kind = dp ), dimension(nLocal)   :: mfact
1165      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 924 | Line 1172 | contains
1172      real ( kind = dp ), intent(inout) :: rijsq
1173      real ( kind = dp )                :: r
1174      real ( kind = dp ), intent(inout) :: d(3)
927 –    real ( kind = dp ) :: ebalance
1175      integer :: me_i, me_j
1176  
1177 <    integer :: iMap
1177 >    integer :: iHash
1178  
1179      r = sqrt(rijsq)
1180      vpair = 0.0d0
# Line 941 | Line 1188 | contains
1188      me_j = atid(j)
1189   #endif
1190  
1191 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1191 >    iHash = InteractionHash(me_i, me_j)
1192  
1193 <    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1194 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1193 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1194 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1195 >            pot(VDW_POT), f, do_pot)
1196      endif
1197  
1198 <    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1198 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1199         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1200 <            pot, eFrame, f, t, do_pot)
1200 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1201  
1202 <       if (FF_uses_dipoles .and. SIM_uses_dipoles) then                
1203 <          if ( PropertyMap(me_i)%is_Dipole .and. &
1204 <               PropertyMap(me_j)%is_Dipole) then
1205 <             if (FF_uses_RF .and. SIM_uses_RF) then
1206 <                call accumulate_rf(i, j, r, eFrame, sw)
959 <                call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
960 <             endif
961 <          endif
1202 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1203 >
1204 >          ! CHECK ME (RF needs to know about all electrostatic types)
1205 >          call accumulate_rf(i, j, r, eFrame, sw)
1206 >          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1207         endif
1208 +
1209      endif
1210  
1211 <    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1211 >    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1212         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1213 <            pot, A, f, t, do_pot)
1213 >            pot(HB_POT), A, f, t, do_pot)
1214      endif
1215  
1216 <    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1216 >    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1217         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1218 <            pot, A, f, t, do_pot)
1218 >            pot(HB_POT), A, f, t, do_pot)
1219      endif
1220  
1221 <    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1221 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1222         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1223 <            pot, A, f, t, do_pot)
1223 >            pot(VDW_POT), A, f, t, do_pot)
1224      endif
1225      
1226 <    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1227 <       call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1228 <            pot, A, f, t, do_pot)
1226 >    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1227 >       call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1228 >            pot(VDW_POT), A, f, t, do_pot)
1229      endif
1230  
1231 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1232 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1233 <            do_pot)
1231 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1232 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1233 >            pot(METALLIC_POT), f, do_pot)
1234      endif
1235  
1236 <    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1236 >    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1237         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1238 <            pot, A, f, t, do_pot)
1238 >            pot(VDW_POT), A, f, t, do_pot)
1239      endif
1240  
1241 <    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1241 >    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1242         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1243 <            pot, A, f, t, do_pot)
1243 >            pot(VDW_POT), A, f, t, do_pot)
1244      endif
1245      
1246    end subroutine do_pair
# Line 1002 | Line 1248 | contains
1248    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1249         do_pot, do_stress, eFrame, A, f, t, pot)
1250  
1251 <    real( kind = dp ) :: pot, sw
1251 >    real( kind = dp ) :: sw
1252 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1253      real( kind = dp ), dimension(9,nLocal) :: eFrame
1254      real (kind=dp), dimension(9,nLocal) :: A
1255      real (kind=dp), dimension(3,nLocal) :: f
# Line 1014 | Line 1261 | contains
1261      real ( kind = dp )                :: r, rc
1262      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1263  
1264 <    integer :: me_i, me_j, iMap
1264 >    integer :: me_i, me_j, iHash
1265 >
1266 >    r = sqrt(rijsq)
1267  
1268   #ifdef IS_MPI  
1269      me_i = atid_row(i)
# Line 1024 | Line 1273 | contains
1273      me_j = atid(j)  
1274   #endif
1275  
1276 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1276 >    iHash = InteractionHash(me_i, me_j)
1277  
1278 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1278 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1279              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1280      endif
1281      
# Line 1035 | Line 1284 | contains
1284  
1285    subroutine do_preforce(nlocal,pot)
1286      integer :: nlocal
1287 <    real( kind = dp ) :: pot
1287 >    real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1288  
1289      if (FF_uses_EAM .and. SIM_uses_EAM) then
1290 <       call calc_EAM_preforce_Frho(nlocal,pot)
1290 >       call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1291      endif
1292  
1293  
# Line 1223 | Line 1472 | contains
1472  
1473    function FF_UsesDirectionalAtoms() result(doesit)
1474      logical :: doesit
1475 <    doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1227 <         FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1228 <         FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1475 >    doesit = FF_uses_DirectionalAtoms
1476    end function FF_UsesDirectionalAtoms
1477  
1478    function FF_RequiresPrepairCalc() result(doesit)
# Line 1235 | Line 1482 | contains
1482  
1483    function FF_RequiresPostpairCalc() result(doesit)
1484      logical :: doesit
1485 <    doesit = FF_uses_RF
1485 >    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1486    end function FF_RequiresPostpairCalc
1487  
1488   #ifdef PROFILE

Diff Legend

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