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 2269 by chuckv, Tue Aug 9 19:40:56 2005 UTC vs.
Revision 2411 by chrisfen, Wed Nov 2 21:01:21 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.27 2005-08-09 19:40:56 chuckv Exp $, $Date: 2005-08-09 19:40:56 $, $Name: not supported by cvs2svn $, $Revision: 1.27 $
48 > !! @version $Id: doForces.F90,v 1.67 2005-11-02 21:01:18 chrisfen Exp $, $Date: 2005-11-02 21:01:18 $, $Name: not supported by cvs2svn $, $Revision: 1.67 $
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 +  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 :: createGroupCutoffs
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 133 | Line 119 | module doForces
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 < !! Cuttoffs in OOPSE are handled on a Group-Group pair basis.
126 < ! Largest cutoff for atypes for all potentials
127 <  real(kind=dp), dimension(:), allocatable :: atypeMaxCuttoff
142 < ! Largest cutoff for groups
143 <  real(kind=dp), dimension(:), allocatable :: groupMaxCutoff
144 < ! Group to Gtype transformation Map
145 <  integer,dimension(:), allocatable :: groupToGtype
146 < ! Group Type Max Cutoff
147 <  real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff
148 < ! GroupType definition
149 <  type ::gtype
150 <     real(kind=dp) :: rcut ! Group Cutoff
151 <     real(kind=dp) :: rcutsq ! Group Cutoff Squared
152 <     real(kind=dp) :: rlistsq ! List cutoff Squared    
153 <  end type gtype
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(gtype), dimension(:,:), allocatable :: gtypeCutoffMap
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 <
160 <  subroutine createInteractionMap(status)
147 >  subroutine createInteractionHash(status)
148      integer :: nAtypes
149      integer, intent(out) :: status
150      integer :: i
151      integer :: j
152 <    integer :: ihash
166 <    real(kind=dp) :: myRcut
152 >    integer :: iHash
153      !! Test Types
154      logical :: i_is_LJ
155      logical :: i_is_Elect
# Line 179 | Line 165 | contains
165      logical :: j_is_GB
166      logical :: j_is_EAM
167      logical :: j_is_Shape
168 <    
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 createDefaultInteractionHash!")
173 >       call handleError("atype", "atypes was not present before call of createInteractionHash!")
174         status = -1
175         return
176      endif
# Line 197 | Line 184 | contains
184  
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 257 | Line 254 | contains
254  
255      end do
256  
257 <    haveInteractionMap = .true.
258 <  end subroutine createInteractionMap
257 >    haveInteractionHash = .true.
258 >  end subroutine createInteractionHash
259  
260 <  subroutine createGroupCutoffs(skinThickness,defaultrList,stat)
264 <    real(kind=dp), intent(in), optional :: defaultRList
265 <    real(kind-dp), intent(in), :: skinThickenss
266 <  ! Query each potential and return the cutoff for that potential. We
267 <  ! build the neighbor list based on the largest cutoff value for that
268 <  ! atype. Each potential can decide whether to calculate the force for
269 <  ! that atype based upon it's own cutoff.
270 <  
260 >  subroutine createGtypeCutoffMap(stat)
261  
262 <    real(kind=dp), intent(in), optional :: defaultRCut, defaultSkinThickness
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 :: iMap
273 <    integer :: map_i,map_j
274 <    real(kind=dp) :: thisRCut = 0.0_dp
275 <    real(kind=dp) :: actualCutoff = 0.0_dp
276 <    integer, intent(out) :: stat
277 <    integer :: nAtypes
278 <    integer :: myStatus
281 <
282 <    stat = 0
283 <    if (.not. haveInteractionMap) then
284 <
285 <       call createInteractionMap(myStatus)
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, *) 'createInteractionMap failed in doForces!'
284 >          write(default_error, *) 'createInteractionHash failed in doForces!'
285            stat = -1
286            return
287         endif
288      endif
289 <
289 > #ifdef IS_MPI
290 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
291 >    nGroupsInCol = getNgroupsInCol(plan_group_col)
292 > #endif
293      nAtypes = getSize(atypes)
294 <    !! If we pass a default rcut, set all atypes to that cutoff distance
295 <    if(present(defaultRList)) then
296 <       InteractionMap(:,:)%rCut = defaultRCut
297 <       InteractionMap(:,:)%rCutSq = defaultRCut*defaultRCut
298 <       InteractionMap(:,:)%rListSq = (defaultRCut+defaultSkinThickness)**2
299 <       haveRlist = .true.
300 <       return
301 <    end if
302 <
303 <    do map_i = 1,nAtypes
304 <       do map_j = map_i,nAtypes
306 <          iMap = InteractionMap(map_i, map_j)%InteractionHash
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            
313          if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
314             ! thisRCut = getElectrostaticCutOff(map_i,map_j)
315             if (thisRcut > actualCutoff) actualCutoff = thisRcut
316          endif
340            
341 <          if ( iand(iMap, STICKY_PAIR).ne.0 ) then
342 <             ! thisRCut = getStickyCutOff(map_i,map_j)
343 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
321 <           endif
322 <          
323 <           if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
324 <              ! thisRCut = getStickyPowerCutOff(map_i,map_j)
325 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
326 <           endif
327 <          
328 <           if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
329 <              ! thisRCut = getGayberneCutOff(map_i,map_j)
330 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
331 <           endif
332 <          
333 <           if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
334 < !              thisRCut = getGaybrneLJCutOff(map_i,map_j)
335 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
336 <           endif
337 <          
338 <           if ( iand(iMap, EAM_PAIR).ne.0 ) then      
339 < !              thisRCut = getEAMCutOff(map_i,map_j)
340 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
341 <           endif
342 <          
343 <           if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
344 < !              thisRCut = getShapeCutOff(map_i,map_j)
345 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
346 <           endif
347 <          
348 <           if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
349 < !              thisRCut = getShapeLJCutOff(map_i,map_j)
350 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
351 <           endif
352 <           InteractionMap(map_i, map_j)%rCut = actualCutoff
353 <           InteractionMap(map_i, map_j)%rCutSq = actualCutoff * actualCutoff
354 <           InteractionMap(map_i, map_j)%rListSq = (actualCutoff + skinThickness)**2
341 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
342 >             biggestAtypeCutoff = atypeMaxCutoff(i)
343 >          endif
344  
345 <           InteractionMap(map_j, map_i)%rCut = InteractionMap(map_i, map_j)%rCut
346 <           InteractionMap(map_j, map_i)%rCutSq = InteractionMap(map_i, map_j)%rCutSq
347 <           InteractionMap(map_j, map_i)%rListSq = InteractionMap(map_i, map_j)%rListSq
348 <        end do
349 <     end do
350 <     ! now the groups
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 <     haveRlist = .true.
384 <   end subroutine createGroupCutoffs
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()
370    SIM_uses_LennardJones = SimUsesLennardJones()
371    SIM_uses_Electrostatics = SimUsesElectrostatics()
372    SIM_uses_Charges = SimUsesCharges()
373    SIM_uses_Dipoles = SimUsesDipoles()
374    SIM_uses_Sticky = SimUsesSticky()
375    SIM_uses_StickyPower = SimUsesStickyPower()
376    SIM_uses_GayBerne = SimUsesGayBerne()
571      SIM_uses_EAM = SimUsesEAM()
378    SIM_uses_Shapes = SimUsesShapes()
379    SIM_uses_FLARB = SimUsesFLARB()
380    SIM_uses_RF = SimUsesRF()
572      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
573      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
574      SIM_uses_PBC = SimUsesPBC()
# Line 394 | Line 585 | contains
585  
586      error = 0
587  
588 <    if (.not. haveInteractionMap) then
398 <      
588 >    if (.not. haveInteractionHash) then      
589         myStatus = 0      
590 <       call createInteractionMap(myStatus)
401 <      
590 >       call createInteractionHash(myStatus)      
591         if (myStatus .ne. 0) then
592 <          write(default_error, *) 'createInteractionMap failed in doForces!'
592 >          write(default_error, *) 'createInteractionHash failed in doForces!'
593            error = -1
594            return
595         endif
596      endif
597  
598 +    if (.not. haveGtypeCutoffMap) then        
599 +       myStatus = 0      
600 +       call createGtypeCutoffMap(myStatus)      
601 +       if (myStatus .ne. 0) then
602 +          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
603 +          error = -1
604 +          return
605 +       endif
606 +    endif
607 +
608      if (.not. haveSIMvariables) then
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 439 | 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
445 <
643 >    integer, intent(in) :: thisESM
644      integer, intent(out) :: thisStat  
645      integer :: my_status, nMatches
646      integer, pointer :: MatchList(:) => null()
449    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:
455 <    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 461 | Line 657 | contains
657      !! interactions are used by the force field.    
658  
659      FF_uses_DirectionalAtoms = .false.
464    FF_uses_LennardJones = .false.
465    FF_uses_Electrostatics = .false.
466    FF_uses_Charges = .false.    
660      FF_uses_Dipoles = .false.
468    FF_uses_Sticky = .false.
469    FF_uses_StickyPower = .false.
661      FF_uses_GayBerne = .false.
662      FF_uses_EAM = .false.
472    FF_uses_Shapes = .false.
473    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  
479    call getMatchingElementList(atypes, "is_LennardJones", .true., &
480         nMatches, MatchList)
481    if (nMatches .gt. 0) FF_uses_LennardJones = .true.
482
483    call getMatchingElementList(atypes, "is_Electrostatic", .true., &
484         nMatches, MatchList)
485    if (nMatches .gt. 0) then
486       FF_uses_Electrostatics = .true.
487    endif
488
489    call getMatchingElementList(atypes, "is_Charge", .true., &
490         nMatches, MatchList)
491    if (nMatches .gt. 0) then
492       FF_uses_Charges = .true.  
493       FF_uses_Electrostatics = .true.
494    endif
495
668      call getMatchingElementList(atypes, "is_Dipole", .true., &
669           nMatches, MatchList)
670 <    if (nMatches .gt. 0) then
499 <       FF_uses_Dipoles = .true.
500 <       FF_uses_Electrostatics = .true.
501 <       FF_uses_DirectionalAtoms = .true.
502 <    endif
503 <
504 <    call getMatchingElementList(atypes, "is_Quadrupole", .true., &
505 <         nMatches, MatchList)
506 <    if (nMatches .gt. 0) then
507 <       FF_uses_Quadrupoles = .true.
508 <       FF_uses_Electrostatics = .true.
509 <       FF_uses_DirectionalAtoms = .true.
510 <    endif
511 <
512 <    call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
513 <         MatchList)
514 <    if (nMatches .gt. 0) then
515 <       FF_uses_Sticky = .true.
516 <       FF_uses_DirectionalAtoms = .true.
517 <    endif
518 <
519 <    call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
520 <         MatchList)
521 <    if (nMatches .gt. 0) then
522 <       FF_uses_StickyPower = .true.
523 <       FF_uses_DirectionalAtoms = .true.
524 <    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
529 <       FF_uses_GayBerne = .true.
530 <       FF_uses_DirectionalAtoms = .true.
531 <    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.
678  
536    call getMatchingElementList(atypes, "is_Shape", .true., &
537         nMatches, MatchList)
538    if (nMatches .gt. 0) then
539       FF_uses_Shapes = .true.
540       FF_uses_DirectionalAtoms = .true.
541    endif
679  
543    call getMatchingElementList(atypes, "is_FLARB", .true., &
544         nMatches, MatchList)
545    if (nMatches .gt. 0) FF_uses_FLARB = .true.
546
547    !! Assume sanity (for the sake of argument)
680      haveSaneForceField = .true.
681  
550    !! check to make sure the FF_uses_RF setting makes sense
551
552    if (FF_uses_dipoles) then
553       if (FF_uses_RF) then
554          dielect = getDielect()
555          call initialize_rf(dielect)
556       endif
557    else
558       if (FF_uses_RF) then          
559          write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
560          thisStat = -1
561          haveSaneForceField = .false.
562          return
563       endif
564    endif
565
566    !sticky module does not contain check_sticky_FF anymore
567    !if (FF_uses_sticky) then
568    !   call check_sticky_FF(my_status)
569    !   if (my_status /= 0) then
570    !      thisStat = -1
571    !      haveSaneForceField = .false.
572    !      return
573    !   end if
574    !endif
575
682      if (FF_uses_EAM) then
683         call init_EAM_FF(my_status)
684         if (my_status /= 0) then
# Line 583 | Line 689 | contains
689         end if
690      endif
691  
586    if (FF_uses_GayBerne) then
587       call check_gb_pair_FF(my_status)
588       if (my_status .ne. 0) then
589          thisStat = -1
590          haveSaneForceField = .false.
591          return
592       endif
593    endif
594
595    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
596    endif
597
692      if (.not. haveNeighborList) then
693         !! Create neighbor lists
694         call expandNeighborList(nLocal, my_status)
# Line 628 | 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 649 | 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 658 | 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 846 | contains
846   #endif
847         outer: do i = istart, iend
848  
753 #ifdef IS_MPI
754             me_i = atid_row(i)
755 #else
756             me_i = atid(i)
757 #endif
758
849            if (update_nlist) point(i) = nlist + 1
850  
851            n_in_i = groupStartRow(i+1) - groupStartRow(i)
# Line 790 | Line 880 | contains
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 812 | Line 902 | contains
902  
903                     list(nlist) = j
904                  endif
905 +                
906 +                if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
907  
908 <                if (loop .eq. PAIR_LOOP) then
909 <                   vij = 0.0d0
910 <                   fij(1:3) = 0.0d0
911 <                endif
912 <
913 <                call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
914 <                     in_switching_region)
915 <
916 <                n_in_j = groupStartCol(j+1) - groupStartCol(j)
917 <
918 <                do ia = groupStartRow(i), groupStartRow(i+1)-1
919 <
920 <                   atom1 = groupListRow(ia)
921 <
922 <                   inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
923 <
924 <                      atom2 = groupListCol(jb)
925 <
926 <                      if (skipThisPair(atom1, atom2)) cycle inner
927 <
928 <                      if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
929 <                         d_atm(1:3) = d_grp(1:3)
930 <                         ratmsq = rgrpsq
931 <                      else
908 >                   if (loop .eq. PAIR_LOOP) then
909 >                      vij = 0.0d0
910 >                      fij(1:3) = 0.0d0
911 >                   endif
912 >                  
913 >                   call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
914 >                        in_switching_region)
915 >                  
916 >                   n_in_j = groupStartCol(j+1) - groupStartCol(j)
917 >                  
918 >                   do ia = groupStartRow(i), groupStartRow(i+1)-1
919 >                      
920 >                      atom1 = groupListRow(ia)
921 >                      
922 >                      inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
923 >                        
924 >                         atom2 = groupListCol(jb)
925 >                        
926 >                         if (skipThisPair(atom1, atom2))  cycle inner
927 >                        
928 >                         if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
929 >                            d_atm(1:3) = d_grp(1:3)
930 >                            ratmsq = rgrpsq
931 >                         else
932   #ifdef IS_MPI
933 <                         call get_interatomic_vector(q_Row(:,atom1), &
934 <                              q_Col(:,atom2), d_atm, ratmsq)
933 >                            call get_interatomic_vector(q_Row(:,atom1), &
934 >                                 q_Col(:,atom2), d_atm, ratmsq)
935   #else
936 <                         call get_interatomic_vector(q(:,atom1), &
937 <                              q(:,atom2), d_atm, ratmsq)
938 < #endif
939 <                      endif
940 <
941 <                      if (loop .eq. PREPAIR_LOOP) then
936 >                            call get_interatomic_vector(q(:,atom1), &
937 >                                 q(:,atom2), d_atm, ratmsq)
938 > #endif
939 >                         endif
940 >                        
941 >                         if (loop .eq. PREPAIR_LOOP) then
942   #ifdef IS_MPI                      
943 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
944 <                              rgrpsq, d_grp, do_pot, do_stress, &
945 <                              eFrame, A, f, t, pot_local)
943 >                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
944 >                                 rgrpsq, d_grp, do_pot, do_stress, &
945 >                                 eFrame, A, f, t, pot_local)
946   #else
947 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
948 <                              rgrpsq, d_grp, do_pot, do_stress, &
949 <                              eFrame, A, f, t, pot)
947 >                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
948 >                                 rgrpsq, d_grp, do_pot, do_stress, &
949 >                                 eFrame, A, f, t, pot)
950   #endif                                              
951 <                      else
951 >                         else
952   #ifdef IS_MPI                      
953 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
954 <                              do_pot, &
955 <                              eFrame, A, f, t, pot_local, vpair, fpair)
953 >                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
954 >                                 do_pot, eFrame, A, f, t, pot_local, vpair, &
955 >                                 fpair, d_grp, rgrp)
956   #else
957 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
958 <                              do_pot,  &
959 <                              eFrame, A, f, t, pot, vpair, fpair)
957 >                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
958 >                                 do_pot, eFrame, A, f, t, pot, vpair, fpair, &
959 >                                 d_grp, rgrp)
960   #endif
961 +                            vij = vij + vpair
962 +                            fij(1:3) = fij(1:3) + fpair(1:3)
963 +                         endif
964 +                      enddo inner
965 +                   enddo
966  
967 <                         vij = vij + vpair
968 <                         fij(1:3) = fij(1:3) + fpair(1:3)
969 <                      endif
970 <                   enddo inner
971 <                enddo
972 <
973 <                if (loop .eq. PAIR_LOOP) then
974 <                   if (in_switching_region) then
975 <                      swderiv = vij*dswdr/rgrp
976 <                      fij(1) = fij(1) + swderiv*d_grp(1)
880 <                      fij(2) = fij(2) + swderiv*d_grp(2)
881 <                      fij(3) = fij(3) + swderiv*d_grp(3)
882 <
883 <                      do ia=groupStartRow(i), groupStartRow(i+1)-1
884 <                         atom1=groupListRow(ia)
885 <                         mf = mfactRow(atom1)
967 >                   if (loop .eq. PAIR_LOOP) then
968 >                      if (in_switching_region) then
969 >                         swderiv = vij*dswdr/rgrp
970 >                         fij(1) = fij(1) + swderiv*d_grp(1)
971 >                         fij(2) = fij(2) + swderiv*d_grp(2)
972 >                         fij(3) = fij(3) + swderiv*d_grp(3)
973 >                        
974 >                         do ia=groupStartRow(i), groupStartRow(i+1)-1
975 >                            atom1=groupListRow(ia)
976 >                            mf = mfactRow(atom1)
977   #ifdef IS_MPI
978 <                         f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
979 <                         f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
980 <                         f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
978 >                            f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
979 >                            f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
980 >                            f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
981   #else
982 <                         f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
983 <                         f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
984 <                         f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
982 >                            f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
983 >                            f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
984 >                            f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
985   #endif
986 <                      enddo
987 <
988 <                      do jb=groupStartCol(j), groupStartCol(j+1)-1
989 <                         atom2=groupListCol(jb)
990 <                         mf = mfactCol(atom2)
986 >                         enddo
987 >                        
988 >                         do jb=groupStartCol(j), groupStartCol(j+1)-1
989 >                            atom2=groupListCol(jb)
990 >                            mf = mfactCol(atom2)
991   #ifdef IS_MPI
992 <                         f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
993 <                         f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
994 <                         f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
992 >                            f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
993 >                            f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
994 >                            f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
995   #else
996 <                         f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
997 <                         f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
998 <                         f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
996 >                            f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
997 >                            f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
998 >                            f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
999   #endif
1000 <                      enddo
1001 <                   endif
1000 >                         enddo
1001 >                      endif
1002  
1003 <                   if (do_stress) call add_stress_tensor(d_grp, fij)
1003 >                      if (do_stress) call add_stress_tensor(d_grp, fij)
1004 >                   endif
1005                  endif
1006 <             end if
1006 >             endif
1007            enddo
1008 +          
1009         enddo outer
1010  
1011         if (update_nlist) then
# Line 972 | Line 1065 | contains
1065  
1066      if (do_pot) then
1067         ! scatter/gather pot_row into the members of my column
1068 <       call scatter(pot_Row, pot_Temp, plan_atom_row)
1069 <
1068 >       do i = 1,LR_POT_TYPES
1069 >          call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1070 >       end do
1071         ! scatter/gather pot_local into all other procs
1072         ! add resultant to get total pot
1073         do i = 1, nlocal
1074 <          pot_local = pot_local + pot_Temp(i)
1074 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1075 >               + pot_Temp(1:LR_POT_TYPES,i)
1076         enddo
1077  
1078         pot_Temp = 0.0_DP
1079 <
1080 <       call scatter(pot_Col, pot_Temp, plan_atom_col)
1079 >       do i = 1,LR_POT_TYPES
1080 >          call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1081 >       end do
1082         do i = 1, nlocal
1083 <          pot_local = pot_local + pot_Temp(i)
1083 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1084 >               + pot_Temp(1:LR_POT_TYPES,i)
1085         enddo
1086  
1087      endif
1088   #endif
1089  
1090 <    if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1090 >    if (SIM_requires_postpair_calc) then
1091 >       do i = 1, nlocal            
1092 >          
1093 >          ! we loop only over the local atoms, so we don't need row and column
1094 >          ! lookups for the types
1095 >          
1096 >          me_i = atid(i)
1097 >          
1098 >          ! is the atom electrostatic?  See if it would have an
1099 >          ! electrostatic interaction with itself
1100 >          iHash = InteractionHash(me_i,me_i)
1101  
1102 <       if (FF_uses_RF .and. SIM_uses_RF) then
996 <
1102 >          if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1103   #ifdef IS_MPI
1104 <          call scatter(rf_Row,rf,plan_atom_row_3d)
1105 <          call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
1000 <          do i = 1,nlocal
1001 <             rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
1002 <          end do
1003 < #endif
1004 <
1005 <          do i = 1, nLocal
1006 <
1007 <             rfpot = 0.0_DP
1008 < #ifdef IS_MPI
1009 <             me_i = atid_row(i)
1104 >             call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1105 >                  t, do_pot)
1106   #else
1107 <             me_i = atid(i)
1107 >             call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1108 >                  t, do_pot)
1109   #endif
1110 <             iMap = InteractionHash(me_i,me_j)
1110 >          endif
1111 >  
1112 >          
1113 >          if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1114              
1115 <             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1116 <
1117 <                mu_i = getDipoleMoment(me_i)
1118 <
1119 <                !! The reaction field needs to include a self contribution
1120 <                !! to the field:
1121 <                call accumulate_self_rf(i, mu_i, eFrame)
1122 <                !! Get the reaction field contribution to the
1123 <                !! potential and torques:
1124 <                call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1115 >             ! loop over the excludes to accumulate RF stuff we've
1116 >             ! left out of the normal pair loop
1117 >            
1118 >             do i1 = 1, nSkipsForAtom(i)
1119 >                j = skipsForAtom(i, i1)
1120 >                
1121 >                ! prevent overcounting of the skips
1122 >                if (i.lt.j) then
1123 >                   call get_interatomic_vector(q(:,i), &
1124 >                        q(:,j), d_atm, ratmsq)
1125 >                   rVal = dsqrt(ratmsq)
1126 >                   call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1127 >                        in_switching_region)
1128   #ifdef IS_MPI
1129 <                pot_local = pot_local + rfpot
1129 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1130 >                        vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1131   #else
1132 <                pot = pot + rfpot
1133 <
1132 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1133 >                        vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1134   #endif
1135 <             endif
1136 <          enddo
1137 <       endif
1135 >                endif
1136 >             enddo
1137 >          endif
1138 >       enddo
1139      endif
1140 <
1036 <
1140 >    
1141   #ifdef IS_MPI
1142 <
1142 >    
1143      if (do_pot) then
1144 <       pot = pot + pot_local
1145 <       !! we assume the c code will do the allreduce to get the total potential
1042 <       !! we could do it right here if we needed to...
1144 >       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1145 >            mpi_comm_world,mpi_err)            
1146      endif
1147 <
1147 >    
1148      if (do_stress) then
1149         call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1150              mpi_comm_world,mpi_err)
1151         call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1152              mpi_comm_world,mpi_err)
1153      endif
1154 <
1154 >    
1155   #else
1156 <
1156 >    
1157      if (do_stress) then
1158         tau = tau_Temp
1159         virial = virial_Temp
1160      endif
1161 <
1161 >    
1162   #endif
1163 <
1163 >    
1164    end subroutine do_force_loop
1165  
1166    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1167 <       eFrame, A, f, t, pot, vpair, fpair)
1167 >       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp)
1168  
1169 <    real( kind = dp ) :: pot, vpair, sw
1169 >    real( kind = dp ) :: vpair, sw
1170 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1171      real( kind = dp ), dimension(3) :: fpair
1172      real( kind = dp ), dimension(nLocal)   :: mfact
1173      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 1074 | Line 1178 | contains
1178      logical, intent(inout) :: do_pot
1179      integer, intent(in) :: i, j
1180      real ( kind = dp ), intent(inout) :: rijsq
1181 <    real ( kind = dp )                :: r
1181 >    real ( kind = dp ), intent(inout) :: r_grp
1182      real ( kind = dp ), intent(inout) :: d(3)
1183 <    real ( kind = dp ) :: ebalance
1183 >    real ( kind = dp ), intent(inout) :: d_grp(3)
1184 >    real ( kind = dp ) :: r
1185      integer :: me_i, me_j
1186  
1187 <    integer :: iMap
1187 >    integer :: iHash
1188  
1189      r = sqrt(rijsq)
1190      vpair = 0.0d0
# Line 1093 | Line 1198 | contains
1198      me_j = atid(j)
1199   #endif
1200  
1201 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1202 <
1203 <    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1204 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1201 >    iHash = InteractionHash(me_i, me_j)
1202 >    
1203 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1204 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1205 >            pot(VDW_POT), f, do_pot)
1206      endif
1207 <
1208 <    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1207 >    
1208 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1209         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1210 <            pot, eFrame, f, t, do_pot)
1105 <
1106 <       if (FF_uses_RF .and. SIM_uses_RF) then
1107 <
1108 <          ! CHECK ME (RF needs to know about all electrostatic types)
1109 <          call accumulate_rf(i, j, r, eFrame, sw)
1110 <          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1111 <       endif
1112 <
1210 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1211      endif
1212 <
1213 <    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1212 >    
1213 >    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1214         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1215 <            pot, A, f, t, do_pot)
1215 >            pot(HB_POT), A, f, t, do_pot)
1216      endif
1217 <
1218 <    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1217 >    
1218 >    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1219         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1220 <            pot, A, f, t, do_pot)
1220 >            pot(HB_POT), A, f, t, do_pot)
1221      endif
1222 <
1223 <    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1222 >    
1223 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1224         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1225 <            pot, A, f, t, do_pot)
1225 >            pot(VDW_POT), A, f, t, do_pot)
1226      endif
1227      
1228 <    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1229 < !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1230 < !           pot, A, f, t, do_pot)
1228 >    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1229 >       call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1230 >            pot(VDW_POT), A, f, t, do_pot)
1231      endif
1232 <
1233 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1234 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1235 <            do_pot)
1232 >    
1233 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1234 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1235 >            pot(METALLIC_POT), f, do_pot)
1236      endif
1237 <
1238 <    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1141 <       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1142 <            pot, A, f, t, do_pot)
1143 <    endif
1144 <
1145 <    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1237 >    
1238 >    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1239         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1240 <            pot, A, f, t, do_pot)
1240 >            pot(VDW_POT), A, f, t, do_pot)
1241      endif
1242      
1243 +    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1244 +       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1245 +            pot(VDW_POT), A, f, t, do_pot)
1246 +    endif
1247 +    
1248    end subroutine do_pair
1249  
1250    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1251         do_pot, do_stress, eFrame, A, f, t, pot)
1252  
1253 <    real( kind = dp ) :: pot, sw
1253 >    real( kind = dp ) :: sw
1254 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1255      real( kind = dp ), dimension(9,nLocal) :: eFrame
1256      real (kind=dp), dimension(9,nLocal) :: A
1257      real (kind=dp), dimension(3,nLocal) :: f
# Line 1164 | Line 1263 | contains
1263      real ( kind = dp )                :: r, rc
1264      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1265  
1266 <    integer :: me_i, me_j, iMap
1266 >    integer :: me_i, me_j, iHash
1267  
1268 +    r = sqrt(rijsq)
1269 +
1270   #ifdef IS_MPI  
1271      me_i = atid_row(i)
1272      me_j = atid_col(j)  
# Line 1174 | Line 1275 | contains
1275      me_j = atid(j)  
1276   #endif
1277  
1278 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1278 >    iHash = InteractionHash(me_i, me_j)
1279  
1280 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1280 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1281              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1282      endif
1283      
# Line 1185 | Line 1286 | contains
1286  
1287    subroutine do_preforce(nlocal,pot)
1288      integer :: nlocal
1289 <    real( kind = dp ) :: pot
1289 >    real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1290  
1291      if (FF_uses_EAM .and. SIM_uses_EAM) then
1292 <       call calc_EAM_preforce_Frho(nlocal,pot)
1292 >       call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1293      endif
1294  
1295  
# Line 1274 | Line 1375 | contains
1375      pot_Col = 0.0_dp
1376      pot_Temp = 0.0_dp
1377  
1277    rf_Row = 0.0_dp
1278    rf_Col = 0.0_dp
1279    rf_Temp = 0.0_dp
1280
1378   #endif
1379  
1380      if (FF_uses_EAM .and. SIM_uses_EAM) then
1381         call clean_EAM()
1382      endif
1383  
1287    rf = 0.0_dp
1384      tau_Temp = 0.0_dp
1385      virial_Temp = 0.0_dp
1386    end subroutine zero_work_arrays
# Line 1373 | Line 1469 | contains
1469  
1470    function FF_UsesDirectionalAtoms() result(doesit)
1471      logical :: doesit
1472 <    doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1377 <         FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1378 <         FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1472 >    doesit = FF_uses_DirectionalAtoms
1473    end function FF_UsesDirectionalAtoms
1474  
1475    function FF_RequiresPrepairCalc() result(doesit)
# Line 1383 | Line 1477 | contains
1477      doesit = FF_uses_EAM
1478    end function FF_RequiresPrepairCalc
1479  
1386  function FF_RequiresPostpairCalc() result(doesit)
1387    logical :: doesit
1388    doesit = FF_uses_RF
1389  end function FF_RequiresPostpairCalc
1390
1480   #ifdef PROFILE
1481    function getforcetime() result(totalforcetime)
1482      real(kind=dp) :: totalforcetime

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines