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 2405 by chrisfen, Tue Nov 1 19:24:57 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.65 2005-11-01 19:24:54 chrisfen Exp $, $Date: 2005-11-01 19:24:54 $, $Name: not supported by cvs2svn $, $Revision: 1.65 $
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.
85    logical, save :: haveSaneForceField = .false.
86 <  logical, save :: haveInteractionMap = .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 <  !!!GO AWAY---------
120 <  !!!!!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
108 <  !public :: addInteraction
109 <  !public :: setInteractionHash
110 <  !public :: getInteractionHash
111 <  public :: createInteractionMap
112 <  public :: createRcuts
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 134 | Line 118 | module doForces
118    real :: forceTimeInitial, forceTimeFinal
119    integer :: nLoops
120   #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
121    
122 <  type(Interaction), dimension(:,:),allocatable :: InteractionMap
123 <  
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 +  integer, dimension(:), allocatable, target :: groupToGtypeRow
130 +  integer, dimension(:), pointer :: groupToGtypeCol => null()
131 +
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 <
151 <  subroutine createInteractionMap(status)
147 >  subroutine createInteractionHash(status)
148      integer :: nAtypes
149 <    integer :: status
149 >    integer, intent(out) :: status
150      integer :: i
151      integer :: j
152 <    integer :: ihash
153 <    real(kind=dp) :: myRcut
158 < ! Test Types
152 >    integer :: iHash
153 >    !! Test Types
154      logical :: i_is_LJ
155      logical :: i_is_Elect
156      logical :: i_is_Sticky
# Line 170 | Line 165 | contains
165      logical :: j_is_GB
166      logical :: j_is_EAM
167      logical :: j_is_Shape
168 <    
169 <    
168 >    real(kind=dp) :: myRcut
169 >
170 >    status = 0  
171 >
172      if (.not. associated(atypes)) then
173 <       call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
173 >       call handleError("atype", "atypes was not present before call of createInteractionHash!")
174         status = -1
175         return
176      endif
# Line 185 | Line 182 | contains
182         return
183      end if
184  
185 <    if (.not. allocated(InteractionMap)) then
186 <       allocate(InteractionMap(nAtypes,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_LennardJones", i_is_LJ)
# Line 240 | Line 247 | contains
247            if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
248  
249  
250 <          InteractionMap(i,j)%InteractionHash = iHash
251 <          InteractionMap(j,i)%InteractionHash = iHash
250 >          InteractionHash(i,j) = iHash
251 >          InteractionHash(j,i) = iHash
252  
253         end do
254  
255      end do
249 –  end subroutine createInteractionMap
256  
257 < ! 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.
258 <  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
257 >    haveInteractionHash = .true.
258 >  end subroutine createInteractionHash
259  
260 <    if(.not. allocated(InteractionMap)) return
260 >  subroutine createGtypeCutoffMap(stat)
261  
262 <    nAtypes = getSize(atypes)
263 < ! If we pass a default rcut, set all atypes to that cutoff distance
264 <    if(present(defaultRList)) then
265 <       InteractionMap(:,:)%rList = defaultRList
266 <       InteractionMap(:,:)%rListSq = defaultRList*defaultRList
267 <       haveRlist = .true.
268 <       return
269 <    end if
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 <    do map_i = 1,nAtypes
273 <       do map_j = map_i,nAtypes
274 <          iMap = InteractionMap(map_i, map_j)%InteractionHash
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 <          if ( iand(iMap, LJ_PAIR).ne.0 ) then
307 < !            thisRCut = getLJCutOff(map_i,map_j)
308 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
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            
280 –          if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
281 – !            thisRCut = getElectrostaticCutOff(map_i,map_j)
282 –             if (thisRcut > actualCutoff) actualCutoff = thisRcut
283 –          endif
340            
341 <          if ( iand(iMap, STICKY_PAIR).ne.0 ) then
342 < !             thisRCut = getStickyCutOff(map_i,map_j)
343 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
288 <           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 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
342 >             biggestAtypeCutoff = atypeMaxCutoff(i)
343 >          endif
344  
345 +       endif
346 +    enddo
347 +  
348  
349 < !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
350 < !!$  subroutine setRlistDF( this_rlist )
351 < !!$
352 < !!$   real(kind=dp) :: this_rlist
353 < !!$
354 < !!$    rlist = this_rlist
355 < !!$    rlistsq = rlist * rlist
356 < !!$
357 < !!$    haveRlist = .true.
358 < !!$
359 < !!$  end subroutine setRlistDF
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()
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()
571      SIM_uses_EAM = SimUsesEAM()
350 –    SIM_uses_Shapes = SimUsesShapes()
351 –    SIM_uses_FLARB = SimUsesFLARB()
352 –    SIM_uses_RF = SimUsesRF()
572      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
573      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
574      SIM_uses_PBC = SimUsesPBC()
# Line 366 | Line 585 | contains
585  
586      error = 0
587  
588 <    if (.not. haveInteractionMap) 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 createInteractionMap(myStatus)
374 <
598 >    if (.not. haveGtypeCutoffMap) then        
599 >       myStatus = 0      
600 >       call createGtypeCutoffMap(myStatus)      
601         if (myStatus .ne. 0) then
602 <          write(default_error, *) 'createInteractionMap failed in doForces!'
602 >          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
603            error = -1
604            return
605         endif
# Line 383 | 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 412 | 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
418 <
643 >    integer, intent(in) :: thisESM
644      integer, intent(out) :: thisStat  
645      integer :: my_status, nMatches
646      integer, pointer :: MatchList(:) => null()
422 –    real(kind=dp) :: rcut, rrf, rt, dielect
647  
648      !! assume things are copacetic, unless they aren't
649      thisStat = 0
650  
651 <    !! Fortran's version of a cast:
428 <    FF_uses_RF = use_RF_c
651 >    electrostaticSummationMethod = thisESM
652  
653      !! init_FF is called *after* all of the atom types have been
654      !! defined in atype_module using the new_atype subroutine.
# Line 434 | Line 657 | contains
657      !! interactions are used by the force field.    
658  
659      FF_uses_DirectionalAtoms = .false.
437 –    FF_uses_LennardJones = .false.
438 –    FF_uses_Electrostatics = .false.
439 –    FF_uses_Charges = .false.    
660      FF_uses_Dipoles = .false.
441 –    FF_uses_Sticky = .false.
442 –    FF_uses_StickyPower = .false.
661      FF_uses_GayBerne = .false.
662      FF_uses_EAM = .false.
445 –    FF_uses_Shapes = .false.
446 –    FF_uses_FLARB = .false.
663  
664      call getMatchingElementList(atypes, "is_Directional", .true., &
665           nMatches, MatchList)
666      if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
667  
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 –
668      call getMatchingElementList(atypes, "is_Dipole", .true., &
669           nMatches, MatchList)
670 <    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
670 >    if (nMatches .gt. 0) FF_uses_Dipoles = .true.
671      
672      call getMatchingElementList(atypes, "is_GayBerne", .true., &
673           nMatches, MatchList)
674 <    if (nMatches .gt. 0) then
502 <       FF_uses_GayBerne = .true.
503 <       FF_uses_DirectionalAtoms = .true.
504 <    endif
674 >    if (nMatches .gt. 0) FF_uses_GayBerne = .true.
675  
676      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
677      if (nMatches .gt. 0) FF_uses_EAM = .true.
508 –
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
678  
516 –    call getMatchingElementList(atypes, "is_FLARB", .true., &
517 –         nMatches, MatchList)
518 –    if (nMatches .gt. 0) FF_uses_FLARB = .true.
679  
520 –    !! Assume sanity (for the sake of argument)
680      haveSaneForceField = .true.
681  
523 –    !! check to make sure the FF_uses_RF setting makes sense
524 –
525 –    if (FF_uses_dipoles) then
526 –       if (FF_uses_RF) then
527 –          dielect = getDielect()
528 –          call initialize_rf(dielect)
529 –       endif
530 –    else
531 –       if (FF_uses_RF) then          
532 –          write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
533 –          thisStat = -1
534 –          haveSaneForceField = .false.
535 –          return
536 –       endif
537 –    endif
538 –
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 –
682      if (FF_uses_EAM) then
683         call init_EAM_FF(my_status)
684         if (my_status /= 0) then
# Line 556 | Line 689 | contains
689         end if
690      endif
691  
559 –    if (FF_uses_GayBerne) then
560 –       call check_gb_pair_FF(my_status)
561 –       if (my_status .ne. 0) then
562 –          thisStat = -1
563 –          haveSaneForceField = .false.
564 –          return
565 –       endif
566 –    endif
567 –
568 –    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
569 –    endif
570 –
692      if (.not. haveNeighborList) then
693         !! Create neighbor lists
694         call expandNeighborList(nLocal, my_status)
# Line 601 | Line 722 | contains
722  
723      !! Stress Tensor
724      real( kind = dp), dimension(9) :: tau  
725 <    real ( kind = dp ) :: pot
725 >    real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
726      logical ( kind = 2) :: do_pot_c, do_stress_c
727      logical :: do_pot
728      logical :: do_stress
729      logical :: in_switching_region
730   #ifdef IS_MPI
731 <    real( kind = DP ) :: pot_local
731 >    real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
732      integer :: nAtomsInRow
733      integer :: nAtomsInCol
734      integer :: nprocs
# Line 622 | Line 743 | contains
743      integer :: nlist
744      real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
745      real( kind = DP ) :: sw, dswdr, swderiv, mf
746 +    real( kind = DP ) :: rVal
747      real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
748      real(kind=dp) :: rfpot, mu_i, virial
749      integer :: me_i, me_j, n_in_i, n_in_j
# Line 631 | Line 753 | contains
753      integer :: localError
754      integer :: propPack_i, propPack_j
755      integer :: loopStart, loopEnd, loop
756 <    integer :: iMap
757 <    real(kind=dp) :: listSkin = 1.0  
756 >    integer :: iHash
757 >    integer :: i1
758 >  
759  
760      !! initialize local variables  
761  
# Line 750 | Line 873 | contains
873               endif
874  
875   #ifdef IS_MPI
876 +             me_j = atid_col(j)
877               call get_interatomic_vector(q_group_Row(:,i), &
878                    q_group_Col(:,j), d_grp, rgrpsq)
879   #else
880 +             me_j = atid(j)
881               call get_interatomic_vector(q_group(:,i), &
882                    q_group(:,j), d_grp, rgrpsq)
883 < #endif
883 > #endif      
884  
885 <             if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
885 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
886                  if (update_nlist) then
887                     nlist = nlist + 1
888  
# Line 795 | Line 920 | contains
920                     inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
921  
922                        atom2 = groupListCol(jb)
923 <
924 <                      if (skipThisPair(atom1, atom2)) cycle inner
925 <
923 >    
924 >                      if (skipThisPair(atom1, atom2))  cycle inner
925 >                      
926                        if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
927                           d_atm(1:3) = d_grp(1:3)
928                           ratmsq = rgrpsq
# Line 824 | Line 949 | contains
949                        else
950   #ifdef IS_MPI                      
951                           call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
952 <                              do_pot, &
953 <                              eFrame, A, f, t, pot_local, vpair, fpair)
952 >                              do_pot, eFrame, A, f, t, pot_local, vpair, &
953 >                              fpair, d_grp, rgrp)
954   #else
955                           call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
956 <                              do_pot,  &
957 <                              eFrame, A, f, t, pot, vpair, fpair)
956 >                              do_pot, eFrame, A, f, t, pot, vpair, fpair, &
957 >                              d_grp, rgrp)
958   #endif
959  
960                           vij = vij + vpair
# Line 878 | Line 1003 | contains
1003                  endif
1004               end if
1005            enddo
1006 +
1007         enddo outer
1008  
1009         if (update_nlist) then
# Line 937 | Line 1063 | contains
1063  
1064      if (do_pot) then
1065         ! scatter/gather pot_row into the members of my column
1066 <       call scatter(pot_Row, pot_Temp, plan_atom_row)
1067 <
1066 >       do i = 1,LR_POT_TYPES
1067 >          call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1068 >       end do
1069         ! scatter/gather pot_local into all other procs
1070         ! add resultant to get total pot
1071         do i = 1, nlocal
1072 <          pot_local = pot_local + pot_Temp(i)
1072 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1073 >               + pot_Temp(1:LR_POT_TYPES,i)
1074         enddo
1075  
1076         pot_Temp = 0.0_DP
1077 <
1078 <       call scatter(pot_Col, pot_Temp, plan_atom_col)
1077 >       do i = 1,LR_POT_TYPES
1078 >          call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1079 >       end do
1080         do i = 1, nlocal
1081 <          pot_local = pot_local + pot_Temp(i)
1081 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1082 >               + pot_Temp(1:LR_POT_TYPES,i)
1083         enddo
1084  
1085      endif
1086   #endif
1087  
1088 <    if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1088 >    if (SIM_requires_postpair_calc) then
1089 >       do i = 1, nlocal            
1090 >          
1091 >          ! we loop only over the local atoms, so we don't need row and column
1092 >          ! lookups for the types
1093 >          
1094 >          me_i = atid(i)
1095 >          
1096 >          ! is the atom electrostatic?  See if it would have an
1097 >          ! electrostatic interaction with itself
1098 >          iHash = InteractionHash(me_i,me_i)
1099  
1100 <       if (FF_uses_RF .and. SIM_uses_RF) then
961 <
1100 >          if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1101   #ifdef IS_MPI
1102 <          call scatter(rf_Row,rf,plan_atom_row_3d)
1103 <          call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
965 <          do i = 1,nlocal
966 <             rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
967 <          end do
968 < #endif
969 <
970 <          do i = 1, nLocal
971 <
972 <             rfpot = 0.0_DP
973 < #ifdef IS_MPI
974 <             me_i = atid_row(i)
1102 >             call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1103 >                  t, do_pot)
1104   #else
1105 <             me_i = atid(i)
1105 >             call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1106 >                  t, do_pot)
1107   #endif
1108 <             iMap = InteractionMap(me_i, me_j)%InteractionHash
1108 >          endif
1109 >  
1110 >          
1111 > !!$          if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1112              
1113 <             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1114 <
1115 <                mu_i = getDipoleMoment(me_i)
1116 <
1117 <                !! The reaction field needs to include a self contribution
1118 <                !! to the field:
1119 <                call accumulate_self_rf(i, mu_i, eFrame)
1120 <                !! Get the reaction field contribution to the
1121 <                !! potential and torques:
1122 <                call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1113 >             ! loop over the excludes to accumulate RF stuff we've
1114 >             ! left out of the normal pair loop
1115 >            
1116 >             do i1 = 1, nSkipsForAtom(i)
1117 >                j = skipsForAtom(i, i1)
1118 >                
1119 >                ! prevent overcounting of the skips
1120 >                if (i.lt.j) then
1121 >                   call get_interatomic_vector(q(:,i), &
1122 >                        q(:,j), d_atm, ratmsq)
1123 >                   rVal = dsqrt(ratmsq)
1124 >                   call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1125 >                        in_switching_region)
1126   #ifdef IS_MPI
1127 <                pot_local = pot_local + rfpot
1127 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1128 >                        vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1129   #else
1130 <                pot = pot + rfpot
1131 <
1130 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1131 >                        vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1132   #endif
1133 <             endif
1134 <          enddo
1135 <       endif
1133 >                endif
1134 >             enddo
1135 > !!$          endif
1136 >       enddo
1137      endif
1138 <
1001 <
1138 >    
1139   #ifdef IS_MPI
1140 <
1140 >    
1141      if (do_pot) then
1142 <       pot = pot + pot_local
1143 <       !! we assume the c code will do the allreduce to get the total potential
1007 <       !! we could do it right here if we needed to...
1142 >       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1143 >            mpi_comm_world,mpi_err)            
1144      endif
1145 <
1145 >    
1146      if (do_stress) then
1147         call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1148              mpi_comm_world,mpi_err)
1149         call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1150              mpi_comm_world,mpi_err)
1151      endif
1152 <
1152 >    
1153   #else
1154 <
1154 >    
1155      if (do_stress) then
1156         tau = tau_Temp
1157         virial = virial_Temp
1158      endif
1159 <
1159 >    
1160   #endif
1161 <
1161 >    
1162    end subroutine do_force_loop
1163  
1164    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1165 <       eFrame, A, f, t, pot, vpair, fpair)
1165 >       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp)
1166  
1167 <    real( kind = dp ) :: pot, vpair, sw
1167 >    real( kind = dp ) :: vpair, sw
1168 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1169      real( kind = dp ), dimension(3) :: fpair
1170      real( kind = dp ), dimension(nLocal)   :: mfact
1171      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 1039 | Line 1176 | contains
1176      logical, intent(inout) :: do_pot
1177      integer, intent(in) :: i, j
1178      real ( kind = dp ), intent(inout) :: rijsq
1179 <    real ( kind = dp )                :: r
1179 >    real ( kind = dp ), intent(inout) :: r_grp
1180      real ( kind = dp ), intent(inout) :: d(3)
1181 <    real ( kind = dp ) :: ebalance
1181 >    real ( kind = dp ), intent(inout) :: d_grp(3)
1182 >    real ( kind = dp ) :: r
1183      integer :: me_i, me_j
1184  
1185 <    integer :: iMap
1185 >    integer :: iHash
1186  
1187      r = sqrt(rijsq)
1188      vpair = 0.0d0
# Line 1058 | Line 1196 | contains
1196      me_j = atid(j)
1197   #endif
1198  
1199 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1200 <
1201 <    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1202 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1199 >    iHash = InteractionHash(me_i, me_j)
1200 >    
1201 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1202 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1203 >            pot(VDW_POT), f, do_pot)
1204      endif
1205 <
1206 <    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1205 >    
1206 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1207         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1208 <            pot, eFrame, f, t, do_pot)
1070 <
1071 <       if (FF_uses_RF .and. SIM_uses_RF) then
1072 <
1073 <          ! CHECK ME (RF needs to know about all electrostatic types)
1074 <          call accumulate_rf(i, j, r, eFrame, sw)
1075 <          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1076 <       endif
1077 <
1208 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1209      endif
1210 <
1211 <    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1210 >    
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
1215 >    
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
1220 >    
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)
1230 >    
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      
1235 >    
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)
1108 <    endif
1109 <
1110 <    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1111 <       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1112 <            pot, A, f, t, do_pot)
1238 >            pot(VDW_POT), A, f, t, do_pot)
1239      endif
1240      
1241 +    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1242 +       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1243 +            pot(VDW_POT), A, f, t, do_pot)
1244 +    endif
1245 +    
1246    end subroutine do_pair
1247  
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 1129 | 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)
1270      me_j = atid_col(j)  
# Line 1139 | 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 1150 | 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 1239 | Line 1373 | contains
1373      pot_Col = 0.0_dp
1374      pot_Temp = 0.0_dp
1375  
1242 –    rf_Row = 0.0_dp
1243 –    rf_Col = 0.0_dp
1244 –    rf_Temp = 0.0_dp
1245 –
1376   #endif
1377  
1378      if (FF_uses_EAM .and. SIM_uses_EAM) then
# Line 1338 | Line 1468 | contains
1468  
1469    function FF_UsesDirectionalAtoms() result(doesit)
1470      logical :: doesit
1471 <    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
1471 >    doesit = FF_uses_DirectionalAtoms
1472    end function FF_UsesDirectionalAtoms
1473  
1474    function FF_RequiresPrepairCalc() result(doesit)
# Line 1348 | Line 1476 | contains
1476      doesit = FF_uses_EAM
1477    end function FF_RequiresPrepairCalc
1478  
1351 –  function FF_RequiresPostpairCalc() result(doesit)
1352 –    logical :: doesit
1353 –    doesit = FF_uses_RF
1354 –  end function FF_RequiresPostpairCalc
1355 –
1479   #ifdef PROFILE
1480    function getforcetime() result(totalforcetime)
1481      real(kind=dp) :: totalforcetime

Diff Legend

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