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

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2259 by gezelter, Mon Jun 27 21:01:36 2005 UTC vs.
Revision 3397 by chuckv, Tue May 27 16:39:06 2008 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.20 2005-06-27 21:01:30 gezelter Exp $, $Date: 2005-06-27 21:01:30 $, $Name: not supported by cvs2svn $, $Revision: 1.20 $
48 > !! @version $Id: doForces.F90,v 1.93 2008-05-27 16:39:05 chuckv Exp $, $Date: 2008-05-27 16:39:05 $, $Name: not supported by cvs2svn $, $Revision: 1.93 $
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
65 +  use MetalNonMetal
66 +  use suttonchen
67    use status
68   #ifdef IS_MPI
69    use mpiSimulation
# Line 72 | Line 73 | module doForces
73    PRIVATE
74  
75   #define __FORTRAN90
76 < #include "UseTheForce/fSwitchingFunction.h"
76 > #include "UseTheForce/fCutoffPolicy.h"
77   #include "UseTheForce/DarkSide/fInteractionMap.h"
78 + #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
79  
80    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
81    INTEGER, PARAMETER:: PAIR_LOOP    = 2
82  
81  logical, save :: haveRlist = .false.
83    logical, save :: haveNeighborList = .false.
84    logical, save :: haveSIMvariables = .false.
84  logical, save :: havePropertyMap = .false.
85    logical, save :: haveSaneForceField = .false.
86 +  logical, save :: haveInteractionHash = .false.
87 +  logical, save :: haveGtypeCutoffMap = .false.
88 +  logical, save :: haveDefaultCutoffs = .false.
89 +  logical, save :: haveSkinThickness = .false.
90 +  logical, save :: haveElectrostaticSummationMethod = .false.
91 +  logical, save :: haveCutoffPolicy = .false.
92 +  logical, save :: VisitCutoffsAfterComputing = .false.
93 +  logical, save :: do_box_dipole = .false.
94  
95    logical, save :: FF_uses_DirectionalAtoms
88  logical, save :: FF_uses_LennardJones
89  logical, save :: FF_uses_Electrostatics
90  logical, save :: FF_uses_Charges
96    logical, save :: FF_uses_Dipoles
92  logical, save :: FF_uses_Quadrupoles
93  logical, save :: FF_uses_Sticky
94  logical, save :: FF_uses_StickyPower
97    logical, save :: FF_uses_GayBerne
98    logical, save :: FF_uses_EAM
99 <  logical, save :: FF_uses_Shapes
100 <  logical, save :: FF_uses_FLARB
101 <  logical, save :: FF_uses_RF
99 >  logical, save :: FF_uses_SC
100 >  logical, save :: FF_uses_MNM
101 >
102  
103    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
104    logical, save :: SIM_uses_EAM
105 <  logical, save :: SIM_uses_Shapes
106 <  logical, save :: SIM_uses_FLARB
113 <  logical, save :: SIM_uses_RF
105 >  logical, save :: SIM_uses_SC
106 >  logical, save :: SIM_uses_MNM
107    logical, save :: SIM_requires_postpair_calc
108    logical, save :: SIM_requires_prepair_calc
109    logical, save :: SIM_uses_PBC
110 <  logical, save :: SIM_uses_molecular_cutoffs
110 >  logical, save :: SIM_uses_AtomicVirial
111  
112 <  real(kind=dp), save :: rlist, rlistsq
112 >  integer, save :: electrostaticSummationMethod
113 >  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
114  
115 +  real(kind=dp), save :: defaultRcut, defaultRsw, largestRcut
116 +  real(kind=dp), save :: skinThickness
117 +  logical, save :: defaultDoShiftPot
118 +  logical, save :: defaultDoShiftFrc
119 +
120    public :: init_FF
121 +  public :: setCutoffs
122 +  public :: cWasLame
123 +  public :: setElectrostaticMethod
124 +  public :: setBoxDipole
125 +  public :: getBoxDipole
126 +  public :: setCutoffPolicy
127 +  public :: setSkinThickness
128    public :: do_force_loop
123  public :: setRlistDF
129  
130   #ifdef PROFILE
131    public :: getforcetime
# Line 128 | Line 133 | module doForces
133    real :: forceTimeInitial, forceTimeFinal
134    integer :: nLoops
135   #endif
136 +  
137 +  !! Variables for cutoff mapping and interaction mapping
138 +  ! Bit hash to determine pair-pair interactions.
139 +  integer, dimension(:,:), allocatable :: InteractionHash
140 +  real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
141 +  real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
142 +  real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
143  
144 <  type :: Properties
145 <     logical :: is_Directional   = .false.
134 <     logical :: is_LennardJones  = .false.
135 <     logical :: is_Electrostatic = .false.
136 <     logical :: is_Charge        = .false.
137 <     logical :: is_Dipole        = .false.
138 <     logical :: is_Quadrupole    = .false.
139 <     logical :: is_Sticky        = .false.
140 <     logical :: is_StickyPower   = .false.
141 <     logical :: is_GayBerne      = .false.
142 <     logical :: is_EAM           = .false.
143 <     logical :: is_Shape         = .false.
144 <     logical :: is_FLARB         = .false.
145 <  end type Properties
144 >  integer, dimension(:), allocatable, target :: groupToGtypeRow
145 >  integer, dimension(:), pointer :: groupToGtypeCol => null()
146  
147 <  type(Properties), dimension(:),allocatable :: PropertyMap
147 >  real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
148 >  real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
149 >  type ::gtypeCutoffs
150 >     real(kind=dp) :: rcut
151 >     real(kind=dp) :: rcutsq
152 >     real(kind=dp) :: rlistsq
153 >  end type gtypeCutoffs
154 >  type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
155  
156 +  real(kind=dp), dimension(3) :: boxDipole
157 +
158   contains
159  
160 <  subroutine setRlistDF( this_rlist )
152 <
153 <    real(kind=dp) :: this_rlist
154 <
155 <    rlist = this_rlist
156 <    rlistsq = rlist * rlist
157 <
158 <    haveRlist = .true.
159 <
160 <  end subroutine setRlistDF
161 <
162 <  subroutine createPropertyMap(status)
160 >  subroutine createInteractionHash()
161      integer :: nAtypes
164    integer :: status
162      integer :: i
163 <    logical :: thisProperty
164 <    real (kind=DP) :: thisDPproperty
163 >    integer :: j
164 >    integer :: iHash
165 >    !! Test Types
166 >    logical :: i_is_LJ
167 >    logical :: i_is_Elect
168 >    logical :: i_is_Sticky
169 >    logical :: i_is_StickyP
170 >    logical :: i_is_GB
171 >    logical :: i_is_EAM
172 >    logical :: i_is_Shape
173 >    logical :: i_is_SC
174 >    logical :: j_is_LJ
175 >    logical :: j_is_Elect
176 >    logical :: j_is_Sticky
177 >    logical :: j_is_StickyP
178 >    logical :: j_is_GB
179 >    logical :: j_is_EAM
180 >    logical :: j_is_Shape
181 >    logical :: j_is_SC
182 >    real(kind=dp) :: myRcut
183  
184 <    status = 0
185 <
184 >    if (.not. associated(atypes)) then
185 >       call handleError("doForces", "atypes was not present before call of createInteractionHash!")
186 >       return
187 >    endif
188 >    
189      nAtypes = getSize(atypes)
190 <
190 >    
191      if (nAtypes == 0) then
192 <       status = -1
192 >       call handleError("doForces", "nAtypes was zero during call of createInteractionHash!")
193         return
194      end if
195  
196 <    if (.not. allocated(PropertyMap)) then
197 <       allocate(PropertyMap(nAtypes))
196 >    if (.not. allocated(InteractionHash)) then
197 >       allocate(InteractionHash(nAtypes,nAtypes))
198 >    else
199 >       deallocate(InteractionHash)
200 >       allocate(InteractionHash(nAtypes,nAtypes))
201      endif
202  
203 +    if (.not. allocated(atypeMaxCutoff)) then
204 +       allocate(atypeMaxCutoff(nAtypes))
205 +    else
206 +       deallocate(atypeMaxCutoff)
207 +       allocate(atypeMaxCutoff(nAtypes))
208 +    endif
209 +        
210      do i = 1, nAtypes
211 <       call getElementProperty(atypes, i, "is_Directional", thisProperty)
212 <       PropertyMap(i)%is_Directional = thisProperty
211 >       call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
212 >       call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
213 >       call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
214 >       call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
215 >       call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
216 >       call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
217 >       call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
218 >       call getElementProperty(atypes, i, "is_SC", i_is_SC)
219  
220 <       call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
187 <       PropertyMap(i)%is_LennardJones = thisProperty
220 >       do j = i, nAtypes
221  
222 <       call getElementProperty(atypes, i, "is_Electrostatic", thisProperty)
223 <       PropertyMap(i)%is_Electrostatic = thisProperty
222 >          iHash = 0
223 >          myRcut = 0.0_dp
224  
225 <       call getElementProperty(atypes, i, "is_Charge", thisProperty)
226 <       PropertyMap(i)%is_Charge = thisProperty
225 >          call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ)
226 >          call getElementProperty(atypes, j, "is_Electrostatic", j_is_Elect)
227 >          call getElementProperty(atypes, j, "is_Sticky", j_is_Sticky)
228 >          call getElementProperty(atypes, j, "is_StickyPower", j_is_StickyP)
229 >          call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
230 >          call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
231 >          call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
232 >          call getElementProperty(atypes, j, "is_SC", j_is_SC)
233  
234 <       call getElementProperty(atypes, i, "is_Dipole", thisProperty)
235 <       PropertyMap(i)%is_Dipole = thisProperty
234 >          if (i_is_LJ .and. j_is_LJ) then
235 >             iHash = ior(iHash, LJ_PAIR)            
236 >          endif
237 >          
238 >          if (i_is_Elect .and. j_is_Elect) then
239 >             iHash = ior(iHash, ELECTROSTATIC_PAIR)
240 >          endif
241 >          
242 >          if (i_is_Sticky .and. j_is_Sticky) then
243 >             iHash = ior(iHash, STICKY_PAIR)
244 >          endif
245  
246 <       call getElementProperty(atypes, i, "is_Quadrupole", thisProperty)
247 <       PropertyMap(i)%is_Quadrupole = thisProperty
246 >          if (i_is_StickyP .and. j_is_StickyP) then
247 >             iHash = ior(iHash, STICKYPOWER_PAIR)
248 >          endif
249  
250 <       call getElementProperty(atypes, i, "is_Sticky", thisProperty)
251 <       PropertyMap(i)%is_Sticky = thisProperty
252 <      
204 <       call getElementProperty(atypes, i, "is_StickyPower", thisProperty)
205 <       PropertyMap(i)%is_StickyPower = thisProperty
250 >          if (i_is_EAM .and. j_is_EAM) then
251 >             iHash = ior(iHash, EAM_PAIR)
252 >          endif
253  
254 <       call getElementProperty(atypes, i, "is_GayBerne", thisProperty)
255 <       PropertyMap(i)%is_GayBerne = thisProperty
254 >          if (i_is_SC .and. j_is_SC) then
255 >             iHash = ior(iHash, SC_PAIR)
256 >          endif
257  
258 <       call getElementProperty(atypes, i, "is_EAM", thisProperty)
259 <       PropertyMap(i)%is_EAM = thisProperty
258 >          if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
259 >          if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
260 >          if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
261 >        
262 >          if ((i_is_EAM.or.i_is_SC).and.(.not.(j_is_EAM.or.j_is_SC))) iHash = ior(iHash, MNM_PAIR)
263 >          if ((j_is_EAM.or.j_is_SC).and.(.not.(i_is_EAM.or.i_is_SC))) iHash = ior(iHash, MNM_PAIR)
264  
265 <       call getElementProperty(atypes, i, "is_Shape", thisProperty)
266 <       PropertyMap(i)%is_Shape = thisProperty
265 >          if (i_is_Shape .and. j_is_Shape) iHash = ior(iHash, SHAPE_PAIR)
266 >          if (i_is_Shape .and. j_is_LJ) iHash = ior(iHash, SHAPE_LJ)
267 >          if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
268  
269 <       call getElementProperty(atypes, i, "is_FLARB", thisProperty)
270 <       PropertyMap(i)%is_FLARB = thisProperty
269 >
270 >          InteractionHash(i,j) = iHash
271 >          InteractionHash(j,i) = iHash
272 >
273 >       end do
274 >
275      end do
276  
277 <    havePropertyMap = .true.
277 >    haveInteractionHash = .true.
278 >  end subroutine createInteractionHash
279  
280 <  end subroutine createPropertyMap
280 >  subroutine createGtypeCutoffMap()
281  
282 <  subroutine setSimVariables()
283 <    SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
284 <    SIM_uses_LennardJones = SimUsesLennardJones()
285 <    SIM_uses_Electrostatics = SimUsesElectrostatics()
286 <    SIM_uses_Charges = SimUsesCharges()
287 <    SIM_uses_Dipoles = SimUsesDipoles()
288 <    SIM_uses_Sticky = SimUsesSticky()
289 <    SIM_uses_StickyPower = SimUsesStickyPower()
290 <    SIM_uses_GayBerne = SimUsesGayBerne()
233 <    SIM_uses_EAM = SimUsesEAM()
234 <    SIM_uses_Shapes = SimUsesShapes()
235 <    SIM_uses_FLARB = SimUsesFLARB()
236 <    SIM_uses_RF = SimUsesRF()
237 <    SIM_requires_postpair_calc = SimRequiresPostpairCalc()
238 <    SIM_requires_prepair_calc = SimRequiresPrepairCalc()
239 <    SIM_uses_PBC = SimUsesPBC()
282 >    logical :: i_is_LJ
283 >    logical :: i_is_Elect
284 >    logical :: i_is_Sticky
285 >    logical :: i_is_StickyP
286 >    logical :: i_is_GB
287 >    logical :: i_is_EAM
288 >    logical :: i_is_Shape
289 >    logical :: i_is_SC
290 >    logical :: GtypeFound
291  
292 <    haveSIMvariables = .true.
292 >    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
293 >    integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
294 >    integer :: nGroupsInRow
295 >    integer :: nGroupsInCol
296 >    integer :: nGroupTypesRow,nGroupTypesCol
297 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol
298 >    real(kind=dp) :: biggestAtypeCutoff
299  
300 <    return
301 <  end subroutine setSimVariables
300 >    if (.not. haveInteractionHash) then
301 >       call createInteractionHash()      
302 >    endif
303 > #ifdef IS_MPI
304 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
305 >    nGroupsInCol = getNgroupsInCol(plan_group_col)
306 > #endif
307 >    nAtypes = getSize(atypes)
308 > ! Set all of the initial cutoffs to zero.
309 >    atypeMaxCutoff = 0.0_dp
310 >    do i = 1, nAtypes
311 >       if (SimHasAtype(i)) then    
312 >          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
313 >          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
314 >          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
315 >          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
316 >          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
317 >          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
318 >          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
319 >          call getElementProperty(atypes, i, "is_SC", i_is_SC)
320 >
321 >          if (haveDefaultCutoffs) then
322 >             atypeMaxCutoff(i) = defaultRcut
323 >          else
324 >             if (i_is_LJ) then          
325 >                thisRcut = getSigma(i) * 2.5_dp
326 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
327 >             endif
328 >             if (i_is_Elect) then
329 >                thisRcut = defaultRcut
330 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
331 >             endif
332 >             if (i_is_Sticky) then
333 >                thisRcut = getStickyCut(i)
334 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
335 >             endif
336 >             if (i_is_StickyP) then
337 >                thisRcut = getStickyPowerCut(i)
338 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
339 >             endif
340 >             if (i_is_GB) then
341 >                thisRcut = getGayBerneCut(i)
342 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
343 >             endif
344 >             if (i_is_EAM) then
345 >                thisRcut = getEAMCut(i)
346 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
347 >             endif
348 >             if (i_is_Shape) then
349 >                thisRcut = getShapeCut(i)
350 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
351 >             endif
352 >             if (i_is_SC) then
353 >                thisRcut = getSCCut(i)
354 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
355 >             endif
356 >          endif
357 >                    
358 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
359 >             biggestAtypeCutoff = atypeMaxCutoff(i)
360 >          endif
361  
362 <  subroutine doReadyCheck(error)
363 <    integer, intent(out) :: error
362 >       endif
363 >    enddo
364 >    
365 >    istart = 1
366 >    jstart = 1
367 > #ifdef IS_MPI
368 >    iend = nGroupsInRow
369 >    jend = nGroupsInCol
370 > #else
371 >    iend = nGroups
372 >    jend = nGroups
373 > #endif
374 >    
375 >    !! allocate the groupToGtype and gtypeMaxCutoff here.
376 >    if(.not.allocated(groupToGtypeRow)) then
377 >     !  allocate(groupToGtype(iend))
378 >       allocate(groupToGtypeRow(iend))
379 >    else
380 >       deallocate(groupToGtypeRow)
381 >       allocate(groupToGtypeRow(iend))
382 >    endif
383 >    if(.not.allocated(groupMaxCutoffRow)) then
384 >       allocate(groupMaxCutoffRow(iend))
385 >    else
386 >       deallocate(groupMaxCutoffRow)
387 >       allocate(groupMaxCutoffRow(iend))
388 >    end if
389  
390 <    integer :: myStatus
390 >    if(.not.allocated(gtypeMaxCutoffRow)) then
391 >       allocate(gtypeMaxCutoffRow(iend))
392 >    else
393 >       deallocate(gtypeMaxCutoffRow)
394 >       allocate(gtypeMaxCutoffRow(iend))
395 >    endif
396  
251    error = 0
397  
398 <    if (.not. havePropertyMap) then
398 > #ifdef IS_MPI
399 >       ! We only allocate new storage if we are in MPI because Ncol /= Nrow
400 >    if(.not.associated(groupToGtypeCol)) then
401 >       allocate(groupToGtypeCol(jend))
402 >    else
403 >       deallocate(groupToGtypeCol)
404 >       allocate(groupToGtypeCol(jend))
405 >    end if
406  
407 <       myStatus = 0
407 >    if(.not.associated(groupMaxCutoffCol)) then
408 >       allocate(groupMaxCutoffCol(jend))
409 >    else
410 >       deallocate(groupMaxCutoffCol)
411 >       allocate(groupMaxCutoffCol(jend))
412 >    end if
413 >    if(.not.associated(gtypeMaxCutoffCol)) then
414 >       allocate(gtypeMaxCutoffCol(jend))
415 >    else
416 >       deallocate(gtypeMaxCutoffCol)      
417 >       allocate(gtypeMaxCutoffCol(jend))
418 >    end if
419  
420 <       call createPropertyMap(myStatus)
420 >       groupMaxCutoffCol = 0.0_dp
421 >       gtypeMaxCutoffCol = 0.0_dp
422  
423 <       if (myStatus .ne. 0) then
424 <          write(default_error, *) 'createPropertyMap failed in doForces!'
425 <          error = -1
426 <          return
423 > #endif
424 >       groupMaxCutoffRow = 0.0_dp
425 >       gtypeMaxCutoffRow = 0.0_dp
426 >
427 >
428 >    !! first we do a single loop over the cutoff groups to find the
429 >    !! largest cutoff for any atypes present in this group.  We also
430 >    !! create gtypes at this point.
431 >    
432 >    tol = 1.0e-6_dp
433 >    nGroupTypesRow = 0
434 >    nGroupTypesCol = 0
435 >    do i = istart, iend      
436 >       n_in_i = groupStartRow(i+1) - groupStartRow(i)
437 >       groupMaxCutoffRow(i) = 0.0_dp
438 >       do ia = groupStartRow(i), groupStartRow(i+1)-1
439 >          atom1 = groupListRow(ia)
440 > #ifdef IS_MPI
441 >          me_i = atid_row(atom1)
442 > #else
443 >          me_i = atid(atom1)
444 > #endif          
445 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
446 >             groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
447 >          endif          
448 >       enddo
449 >       if (nGroupTypesRow.eq.0) then
450 >          nGroupTypesRow = nGroupTypesRow + 1
451 >          gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
452 >          groupToGtypeRow(i) = nGroupTypesRow
453 >       else
454 >          GtypeFound = .false.
455 >          do g = 1, nGroupTypesRow
456 >             if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
457 >                groupToGtypeRow(i) = g
458 >                GtypeFound = .true.
459 >             endif
460 >          enddo
461 >          if (.not.GtypeFound) then            
462 >             nGroupTypesRow = nGroupTypesRow + 1
463 >             gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
464 >             groupToGtypeRow(i) = nGroupTypesRow
465 >          endif
466         endif
467 +    enddo    
468 +
469 + #ifdef IS_MPI
470 +    do j = jstart, jend      
471 +       n_in_j = groupStartCol(j+1) - groupStartCol(j)
472 +       groupMaxCutoffCol(j) = 0.0_dp
473 +       do ja = groupStartCol(j), groupStartCol(j+1)-1
474 +          atom1 = groupListCol(ja)
475 +
476 +          me_j = atid_col(atom1)
477 +
478 +          if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
479 +             groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
480 +          endif          
481 +       enddo
482 +
483 +       if (nGroupTypesCol.eq.0) then
484 +          nGroupTypesCol = nGroupTypesCol + 1
485 +          gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
486 +          groupToGtypeCol(j) = nGroupTypesCol
487 +       else
488 +          GtypeFound = .false.
489 +          do g = 1, nGroupTypesCol
490 +             if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
491 +                groupToGtypeCol(j) = g
492 +                GtypeFound = .true.
493 +             endif
494 +          enddo
495 +          if (.not.GtypeFound) then            
496 +             nGroupTypesCol = nGroupTypesCol + 1
497 +             gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
498 +             groupToGtypeCol(j) = nGroupTypesCol
499 +          endif
500 +       endif
501 +    enddo    
502 +
503 + #else
504 + ! Set pointers to information we just found
505 +    nGroupTypesCol = nGroupTypesRow
506 +    groupToGtypeCol => groupToGtypeRow
507 +    gtypeMaxCutoffCol => gtypeMaxCutoffRow
508 +    groupMaxCutoffCol => groupMaxCutoffRow
509 + #endif
510 +
511 +    !! allocate the gtypeCutoffMap here.
512 +    allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
513 +    !! then we do a double loop over all the group TYPES to find the cutoff
514 +    !! map between groups of two types
515 +    tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
516 +
517 +    do i = 1, nGroupTypesRow      
518 +       do j = 1, nGroupTypesCol
519 +      
520 +          select case(cutoffPolicy)
521 +          case(TRADITIONAL_CUTOFF_POLICY)
522 +             thisRcut = tradRcut
523 +          case(MIX_CUTOFF_POLICY)
524 +             thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
525 +          case(MAX_CUTOFF_POLICY)
526 +             thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
527 +          case default
528 +             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
529 +             return
530 +          end select
531 +          gtypeCutoffMap(i,j)%rcut = thisRcut
532 +          
533 +          if (thisRcut.gt.largestRcut) largestRcut = thisRcut
534 +
535 +          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
536 +
537 +          if (.not.haveSkinThickness) then
538 +             skinThickness = 1.0_dp
539 +          endif
540 +
541 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
542 +
543 +          ! sanity check
544 +
545 +          if (haveDefaultCutoffs) then
546 +             if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
547 +                call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
548 +             endif
549 +          endif
550 +       enddo
551 +    enddo
552 +
553 +    if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
554 +    if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
555 +    if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
556 + #ifdef IS_MPI
557 +    if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
558 +    if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
559 + #endif
560 +    groupMaxCutoffCol => null()
561 +    gtypeMaxCutoffCol => null()
562 +    
563 +    haveGtypeCutoffMap = .true.
564 +   end subroutine createGtypeCutoffMap
565 +
566 +   subroutine setCutoffs(defRcut, defRsw, defSP, defSF)
567 +
568 +     real(kind=dp),intent(in) :: defRcut, defRsw
569 +     logical, intent(in) :: defSP, defSF
570 +     character(len = statusMsgSize) :: errMsg
571 +     integer :: localError
572 +
573 +     defaultRcut = defRcut
574 +     defaultRsw = defRsw
575 +    
576 +     defaultDoShiftPot = defSP
577 +     defaultDoShiftFrc = defSF
578 +
579 +     if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
580 +        if (defaultDoShiftFrc) then
581 +           write(errMsg, *) &
582 +                'cutoffRadius and switchingRadius are set to the', newline &
583 +                // tab, 'same value.  OOPSE will use shifted force', newline &
584 +                // tab, 'potentials instead of switching functions.'
585 +          
586 +           call handleInfo("setCutoffs", errMsg)
587 +        else
588 +           write(errMsg, *) &
589 +                'cutoffRadius and switchingRadius are set to the', newline &
590 +                // tab, 'same value.  OOPSE will use shifted', newline &
591 +                // tab, 'potentials instead of switching functions.'
592 +          
593 +           call handleInfo("setCutoffs", errMsg)
594 +          
595 +           defaultDoShiftPot = .true.
596 +        endif
597 +                
598 +     endif
599 +    
600 +     localError = 0
601 +     call setLJDefaultCutoff( defaultRcut, defaultDoShiftPot, &
602 +          defaultDoShiftFrc )
603 +     call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
604 +     call setCutoffEAM( defaultRcut )
605 +     call setCutoffSC( defaultRcut )
606 +     call setMnMDefaultCutoff( defaultRcut, defaultDoShiftPot, &
607 +          defaultDoShiftFrc )
608 +     call set_switch(defaultRsw, defaultRcut)
609 +     call setHmatDangerousRcutValue(defaultRcut)
610 +        
611 +     haveDefaultCutoffs = .true.
612 +     haveGtypeCutoffMap = .false.
613 +
614 +   end subroutine setCutoffs
615 +
616 +   subroutine cWasLame()
617 +    
618 +     VisitCutoffsAfterComputing = .true.
619 +     return
620 +    
621 +   end subroutine cWasLame
622 +  
623 +   subroutine setCutoffPolicy(cutPolicy)
624 +    
625 +     integer, intent(in) :: cutPolicy
626 +    
627 +     cutoffPolicy = cutPolicy
628 +     haveCutoffPolicy = .true.
629 +     haveGtypeCutoffMap = .false.
630 +    
631 +   end subroutine setCutoffPolicy
632 +    
633 +   subroutine setBoxDipole()
634 +
635 +     do_box_dipole = .true.
636 +    
637 +   end subroutine setBoxDipole
638 +
639 +   subroutine getBoxDipole( box_dipole )
640 +
641 +     real(kind=dp), intent(inout), dimension(3) :: box_dipole
642 +
643 +     box_dipole = boxDipole
644 +
645 +   end subroutine getBoxDipole
646 +
647 +   subroutine setElectrostaticMethod( thisESM )
648 +
649 +     integer, intent(in) :: thisESM
650 +
651 +     electrostaticSummationMethod = thisESM
652 +     haveElectrostaticSummationMethod = .true.
653 +    
654 +   end subroutine setElectrostaticMethod
655 +
656 +   subroutine setSkinThickness( thisSkin )
657 +    
658 +     real(kind=dp), intent(in) :: thisSkin
659 +    
660 +     skinThickness = thisSkin
661 +     haveSkinThickness = .true.    
662 +     haveGtypeCutoffMap = .false.
663 +    
664 +   end subroutine setSkinThickness
665 +      
666 +   subroutine setSimVariables()
667 +     SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
668 +     SIM_uses_EAM = SimUsesEAM()
669 +     SIM_requires_postpair_calc = SimRequiresPostpairCalc()
670 +     SIM_requires_prepair_calc = SimRequiresPrepairCalc()
671 +     SIM_uses_PBC = SimUsesPBC()
672 +     SIM_uses_SC = SimUsesSC()
673 +     SIM_uses_AtomicVirial = SimUsesAtomicVirial()
674 +
675 +     haveSIMvariables = .true.
676 +    
677 +     return
678 +   end subroutine setSimVariables
679 +
680 +  subroutine doReadyCheck(error)
681 +    integer, intent(out) :: error
682 +    integer :: myStatus
683 +
684 +    error = 0
685 +
686 +    if (.not. haveInteractionHash) then      
687 +       call createInteractionHash()      
688      endif
689  
690 <    if (.not. haveSIMvariables) then
691 <       call setSimVariables()
690 >    if (.not. haveGtypeCutoffMap) then        
691 >       call createGtypeCutoffMap()      
692      endif
693  
694 <    if (.not. haveRlist) then
695 <       write(default_error, *) 'rList has not been set in doForces!'
696 <       error = -1
697 <       return
694 >    if (VisitCutoffsAfterComputing) then
695 >       call set_switch(largestRcut, largestRcut)      
696 >       call setHmatDangerousRcutValue(largestRcut)
697 >       call setCutoffEAM(largestRcut)
698 >       call setCutoffSC(largestRcut)
699 >       VisitCutoffsAfterComputing = .false.
700      endif
701  
702 +    if (.not. haveSIMvariables) then
703 +       call setSimVariables()
704 +    endif
705 +
706      if (.not. haveNeighborList) then
707         write(default_error, *) 'neighbor list has not been initialized in doForces!'
708         error = -1
709         return
710      end if
711 <
711 >    
712      if (.not. haveSaneForceField) then
713         write(default_error, *) 'Force Field is not sane in doForces!'
714         error = -1
715         return
716      end if
717 <
717 >    
718   #ifdef IS_MPI
719      if (.not. isMPISimSet()) then
720         write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
# Line 296 | Line 726 | contains
726    end subroutine doReadyCheck
727  
728  
729 <  subroutine init_FF(use_RF_c, thisStat)
729 >  subroutine init_FF(thisStat)
730  
301    logical, intent(in) :: use_RF_c
302
731      integer, intent(out) :: thisStat  
732      integer :: my_status, nMatches
733      integer, pointer :: MatchList(:) => null()
306    real(kind=dp) :: rcut, rrf, rt, dielect
734  
735      !! assume things are copacetic, unless they aren't
736      thisStat = 0
737  
311    !! Fortran's version of a cast:
312    FF_uses_RF = use_RF_c
313
738      !! init_FF is called *after* all of the atom types have been
739      !! defined in atype_module using the new_atype subroutine.
740      !!
# Line 318 | Line 742 | contains
742      !! interactions are used by the force field.    
743  
744      FF_uses_DirectionalAtoms = .false.
321    FF_uses_LennardJones = .false.
322    FF_uses_Electrostatics = .false.
323    FF_uses_Charges = .false.    
745      FF_uses_Dipoles = .false.
325    FF_uses_Sticky = .false.
326    FF_uses_StickyPower = .false.
746      FF_uses_GayBerne = .false.
747      FF_uses_EAM = .false.
748 <    FF_uses_Shapes = .false.
330 <    FF_uses_FLARB = .false.
748 >    FF_uses_SC = .false.
749  
750      call getMatchingElementList(atypes, "is_Directional", .true., &
751           nMatches, MatchList)
752      if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
753  
336    call getMatchingElementList(atypes, "is_LennardJones", .true., &
337         nMatches, MatchList)
338    if (nMatches .gt. 0) FF_uses_LennardJones = .true.
339
340    call getMatchingElementList(atypes, "is_Electrostatic", .true., &
341         nMatches, MatchList)
342    if (nMatches .gt. 0) then
343       FF_uses_Electrostatics = .true.
344    endif
345
346    call getMatchingElementList(atypes, "is_Charge", .true., &
347         nMatches, MatchList)
348    if (nMatches .gt. 0) then
349       FF_uses_Charges = .true.  
350       FF_uses_Electrostatics = .true.
351    endif
352
754      call getMatchingElementList(atypes, "is_Dipole", .true., &
755           nMatches, MatchList)
756 <    if (nMatches .gt. 0) then
356 <       FF_uses_Dipoles = .true.
357 <       FF_uses_Electrostatics = .true.
358 <       FF_uses_DirectionalAtoms = .true.
359 <    endif
360 <
361 <    call getMatchingElementList(atypes, "is_Quadrupole", .true., &
362 <         nMatches, MatchList)
363 <    if (nMatches .gt. 0) then
364 <       FF_uses_Quadrupoles = .true.
365 <       FF_uses_Electrostatics = .true.
366 <       FF_uses_DirectionalAtoms = .true.
367 <    endif
368 <
369 <    call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
370 <         MatchList)
371 <    if (nMatches .gt. 0) then
372 <       FF_uses_Sticky = .true.
373 <       FF_uses_DirectionalAtoms = .true.
374 <    endif
375 <
376 <    call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
377 <         MatchList)
378 <    if (nMatches .gt. 0) then
379 <       FF_uses_StickyPower = .true.
380 <       FF_uses_DirectionalAtoms = .true.
381 <    endif
756 >    if (nMatches .gt. 0) FF_uses_Dipoles = .true.
757      
758      call getMatchingElementList(atypes, "is_GayBerne", .true., &
759           nMatches, MatchList)
760 <    if (nMatches .gt. 0) then
386 <       FF_uses_GayBerne = .true.
387 <       FF_uses_DirectionalAtoms = .true.
388 <    endif
760 >    if (nMatches .gt. 0) FF_uses_GayBerne = .true.
761  
762      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
763      if (nMatches .gt. 0) FF_uses_EAM = .true.
764  
765 <    call getMatchingElementList(atypes, "is_Shape", .true., &
766 <         nMatches, MatchList)
395 <    if (nMatches .gt. 0) then
396 <       FF_uses_Shapes = .true.
397 <       FF_uses_DirectionalAtoms = .true.
398 <    endif
765 >    call getMatchingElementList(atypes, "is_SC", .true., nMatches, MatchList)
766 >    if (nMatches .gt. 0) FF_uses_SC = .true.
767  
400    call getMatchingElementList(atypes, "is_FLARB", .true., &
401         nMatches, MatchList)
402    if (nMatches .gt. 0) FF_uses_FLARB = .true.
768  
404    !! Assume sanity (for the sake of argument)
769      haveSaneForceField = .true.
770  
407    !! check to make sure the FF_uses_RF setting makes sense
408
409    if (FF_uses_dipoles) then
410       if (FF_uses_RF) then
411          dielect = getDielect()
412          call initialize_rf(dielect)
413       endif
414    else
415       if (FF_uses_RF) then          
416          write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
417          thisStat = -1
418          haveSaneForceField = .false.
419          return
420       endif
421    endif
422
423    !sticky module does not contain check_sticky_FF anymore
424    !if (FF_uses_sticky) then
425    !   call check_sticky_FF(my_status)
426    !   if (my_status /= 0) then
427    !      thisStat = -1
428    !      haveSaneForceField = .false.
429    !      return
430    !   end if
431    !endif
432
771      if (FF_uses_EAM) then
772         call init_EAM_FF(my_status)
773         if (my_status /= 0) then
# Line 440 | Line 778 | contains
778         end if
779      endif
780  
443    if (FF_uses_GayBerne) then
444       call check_gb_pair_FF(my_status)
445       if (my_status .ne. 0) then
446          thisStat = -1
447          haveSaneForceField = .false.
448          return
449       endif
450    endif
451
452    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
453    endif
454
781      if (.not. haveNeighborList) then
782         !! Create neighbor lists
783         call expandNeighborList(nLocal, my_status)
# Line 468 | Line 794 | contains
794  
795    !! Does force loop over i,j pairs. Calls do_pair to calculates forces.
796    !------------------------------------------------------------->
797 <  subroutine do_force_loop(q, q_group, A, eFrame, f, t, tau, pot, &
797 >  subroutine do_force_loop(q, q_group, A, eFrame, f, t, tau, pot, particle_pot,&
798         do_pot_c, do_stress_c, error)
799      !! Position array provided by C, dimensioned by getNlocal
800      real ( kind = dp ), dimension(3, nLocal) :: q
# Line 485 | Line 811 | contains
811  
812      !! Stress Tensor
813      real( kind = dp), dimension(9) :: tau  
814 <    real ( kind = dp ) :: pot
814 >    real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
815 >    real( kind = dp ), dimension(nLocal) :: particle_pot
816      logical ( kind = 2) :: do_pot_c, do_stress_c
817      logical :: do_pot
818      logical :: do_stress
819      logical :: in_switching_region
820   #ifdef IS_MPI
821 <    real( kind = DP ) :: pot_local
821 >    real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
822      integer :: nAtomsInRow
823      integer :: nAtomsInCol
824      integer :: nprocs
# Line 504 | Line 831 | contains
831      integer :: istart, iend
832      integer :: ia, jb, atom1, atom2
833      integer :: nlist
834 <    real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
834 >    real( kind = DP ) :: ratmsq, rgrpsq, rgrp, rag, vpair, vij
835      real( kind = DP ) :: sw, dswdr, swderiv, mf
836 <    real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
837 <    real(kind=dp) :: rfpot, mu_i, virial
836 >    real( kind = DP ) :: rVal
837 >    real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij, fg, dag
838 >    real(kind=dp) :: rfpot, mu_i
839 >    real(kind=dp):: rCut
840      integer :: me_i, me_j, n_in_i, n_in_j
841      logical :: is_dp_i
842      integer :: neighborListSize
# Line 515 | Line 844 | contains
844      integer :: localError
845      integer :: propPack_i, propPack_j
846      integer :: loopStart, loopEnd, loop
847 +    integer :: iHash
848 +    integer :: i1
849  
850 <    real(kind=dp) :: listSkin = 1.0  
850 >    !! the variables for the box dipole moment
851 > #ifdef IS_MPI
852 >    integer :: pChgCount_local
853 >    integer :: nChgCount_local
854 >    real(kind=dp) :: pChg_local
855 >    real(kind=dp) :: nChg_local
856 >    real(kind=dp), dimension(3) :: pChgPos_local
857 >    real(kind=dp), dimension(3) :: nChgPos_local
858 >    real(kind=dp), dimension(3) :: dipVec_local
859 > #endif
860 >    integer :: pChgCount
861 >    integer :: nChgCount
862 >    real(kind=dp) :: pChg
863 >    real(kind=dp) :: nChg
864 >    real(kind=dp) :: chg_value
865 >    real(kind=dp), dimension(3) :: pChgPos
866 >    real(kind=dp), dimension(3) :: nChgPos
867 >    real(kind=dp), dimension(3) :: dipVec
868 >    real(kind=dp), dimension(3) :: chgVec
869  
870 +    !! initialize box dipole variables
871 +    if (do_box_dipole) then
872 + #ifdef IS_MPI
873 +       pChg_local = 0.0_dp
874 +       nChg_local = 0.0_dp
875 +       pChgCount_local = 0
876 +       nChgCount_local = 0
877 +       do i=1, 3
878 +          pChgPos_local = 0.0_dp
879 +          nChgPos_local = 0.0_dp
880 +          dipVec_local = 0.0_dp
881 +       enddo
882 + #endif
883 +       pChg = 0.0_dp
884 +       nChg = 0.0_dp
885 +       pChgCount = 0
886 +       nChgCount = 0
887 +       chg_value = 0.0_dp
888 +      
889 +       do i=1, 3
890 +          pChgPos(i) = 0.0_dp
891 +          nChgPos(i) = 0.0_dp
892 +          dipVec(i) = 0.0_dp
893 +          chgVec(i) = 0.0_dp
894 +          boxDipole(i) = 0.0_dp
895 +       enddo
896 +    endif
897 +
898      !! initialize local variables  
899  
900   #ifdef IS_MPI
# Line 580 | Line 957 | contains
957         ! (but only on the first time through):
958         if (loop .eq. loopStart) then
959   #ifdef IS_MPI
960 <          call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
960 >          call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
961                 update_nlist)
962   #else
963 <          call checkNeighborList(nGroups, q_group, listSkin, &
963 >          call checkNeighborList(nGroups, q_group, skinThickness, &
964                 update_nlist)
965   #endif
966         endif
# Line 634 | Line 1011 | contains
1011               endif
1012  
1013   #ifdef IS_MPI
1014 +             me_j = atid_col(j)
1015               call get_interatomic_vector(q_group_Row(:,i), &
1016                    q_group_Col(:,j), d_grp, rgrpsq)
1017   #else
1018 +             me_j = atid(j)
1019               call get_interatomic_vector(q_group(:,i), &
1020                    q_group(:,j), d_grp, rgrpsq)
1021 < #endif
1021 > #endif      
1022  
1023 <             if (rgrpsq < rlistsq) then
1023 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
1024                  if (update_nlist) then
1025                     nlist = nlist + 1
1026  
# Line 661 | Line 1040 | contains
1040  
1041                     list(nlist) = j
1042                  endif
1043 +                
1044 +                if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
1045  
1046 <                if (loop .eq. PAIR_LOOP) then
1047 <                   vij = 0.0d0
1048 <                   fij(1:3) = 0.0d0
1049 <                endif
1050 <
1051 <                call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
1052 <                     in_switching_region)
1053 <
1054 <                n_in_j = groupStartCol(j+1) - groupStartCol(j)
1055 <
1056 <                do ia = groupStartRow(i), groupStartRow(i+1)-1
1057 <
1058 <                   atom1 = groupListRow(ia)
1059 <
1060 <                   inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
1061 <
1062 <                      atom2 = groupListCol(jb)
1063 <
1064 <                      if (skipThisPair(atom1, atom2)) cycle inner
1065 <
1066 <                      if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
1067 <                         d_atm(1:3) = d_grp(1:3)
1068 <                         ratmsq = rgrpsq
1069 <                      else
1070 < #ifdef IS_MPI
1071 <                         call get_interatomic_vector(q_Row(:,atom1), &
1072 <                              q_Col(:,atom2), d_atm, ratmsq)
1073 < #else
1074 <                         call get_interatomic_vector(q(:,atom1), &
1075 <                              q(:,atom2), d_atm, ratmsq)
1076 < #endif
696 <                      endif
697 <
698 <                      if (loop .eq. PREPAIR_LOOP) then
699 < #ifdef IS_MPI                      
700 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
701 <                              rgrpsq, d_grp, do_pot, do_stress, &
702 <                              eFrame, A, f, t, pot_local)
1046 >                   rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
1047 >                   if (loop .eq. PAIR_LOOP) then
1048 >                      vij = 0.0_dp
1049 >                      fij(1) = 0.0_dp
1050 >                      fij(2) = 0.0_dp
1051 >                      fij(3) = 0.0_dp
1052 >                   endif
1053 >                  
1054 >                   call get_switch(rgrpsq, sw, dswdr,rgrp, in_switching_region)
1055 >                  
1056 >                   n_in_j = groupStartCol(j+1) - groupStartCol(j)
1057 >                  
1058 >                   do ia = groupStartRow(i), groupStartRow(i+1)-1
1059 >                      
1060 >                      atom1 = groupListRow(ia)
1061 >                      
1062 >                      inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
1063 >                        
1064 >                         atom2 = groupListCol(jb)
1065 >                        
1066 >                         if (skipThisPair(atom1, atom2))  cycle inner
1067 >                        
1068 >                         if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
1069 >                            d_atm(1) = d_grp(1)
1070 >                            d_atm(2) = d_grp(2)
1071 >                            d_atm(3) = d_grp(3)
1072 >                            ratmsq = rgrpsq
1073 >                         else
1074 > #ifdef IS_MPI
1075 >                            call get_interatomic_vector(q_Row(:,atom1), &
1076 >                                 q_Col(:,atom2), d_atm, ratmsq)
1077   #else
1078 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1079 <                              rgrpsq, d_grp, do_pot, do_stress, &
1080 <                              eFrame, A, f, t, pot)
1078 >                            call get_interatomic_vector(q(:,atom1), &
1079 >                                 q(:,atom2), d_atm, ratmsq)
1080 > #endif
1081 >                         endif
1082 >                        
1083 >                         if (loop .eq. PREPAIR_LOOP) then
1084 > #ifdef IS_MPI                      
1085 >                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1086 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1087 >                                 eFrame, A, f, t, pot_local)
1088 > #else
1089 >                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1090 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1091 >                                 eFrame, A, f, t, pot)
1092   #endif                                              
1093 <                      else
1093 >                         else
1094   #ifdef IS_MPI                      
1095 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1096 <                              do_pot, &
1097 <                              eFrame, A, f, t, pot_local, vpair, fpair)
1095 >                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1096 >                                 do_pot, eFrame, A, f, t, pot_local, vpair, &
1097 >                                 fpair, d_grp, rgrp, rCut)
1098 >                            ! particle_pot will be accumulated from row & column
1099 >                            ! arrays later
1100   #else
1101 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1102 <                              do_pot,  &
1103 <                              eFrame, A, f, t, pot, vpair, fpair)
1101 >                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1102 >                                 do_pot, eFrame, A, f, t, pot, vpair, &
1103 >                                 fpair, d_grp, rgrp, rCut)
1104 >                            particle_pot(atom1) = particle_pot(atom1) + vpair*sw
1105 >                            particle_pot(atom2) = particle_pot(atom2) + vpair*sw
1106   #endif
1107 +                            vij = vij + vpair
1108 +                            fij(1) = fij(1) + fpair(1)
1109 +                            fij(2) = fij(2) + fpair(2)
1110 +                            fij(3) = fij(3) + fpair(3)
1111 +                            if (do_stress) then
1112 +                               call add_stress_tensor(d_atm, fpair, tau)
1113 +                            endif
1114 +                         endif
1115 +                      enddo inner
1116 +                   enddo
1117  
1118 <                         vij = vij + vpair
1119 <                         fij(1:3) = fij(1:3) + fpair(1:3)
1120 <                      endif
1121 <                   enddo inner
1122 <                enddo
1123 <
1124 <                if (loop .eq. PAIR_LOOP) then
1125 <                   if (in_switching_region) then
1126 <                      swderiv = vij*dswdr/rgrp
1127 <                      fij(1) = fij(1) + swderiv*d_grp(1)
1128 <                      fij(2) = fij(2) + swderiv*d_grp(2)
1129 <                      fij(3) = fij(3) + swderiv*d_grp(3)
1130 <
1131 <                      do ia=groupStartRow(i), groupStartRow(i+1)-1
1132 <                         atom1=groupListRow(ia)
1133 <                         mf = mfactRow(atom1)
1134 < #ifdef IS_MPI
1135 <                         f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1136 <                         f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1137 <                         f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1138 < #else
1139 <                         f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1140 <                         f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1141 <                         f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1118 >                   if (loop .eq. PAIR_LOOP) then
1119 >                      if (in_switching_region) then
1120 >                         swderiv = vij*dswdr/rgrp
1121 >                         fg = swderiv*d_grp
1122 >
1123 >                         fij(1) = fij(1) + fg(1)
1124 >                         fij(2) = fij(2) + fg(2)
1125 >                         fij(3) = fij(3) + fg(3)
1126 >                        
1127 >                         if (do_stress .and. (n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
1128 >                            call add_stress_tensor(d_atm, fg, tau)
1129 >                         endif  
1130 >                        
1131 >                         do ia=groupStartRow(i), groupStartRow(i+1)-1
1132 >                            atom1=groupListRow(ia)
1133 >                            mf = mfactRow(atom1)
1134 >                            ! fg is the force on atom ia due to cutoff group's
1135 >                            ! presence in switching region
1136 >                            fg = swderiv*d_grp*mf
1137 > #ifdef IS_MPI
1138 >                            f_Row(1,atom1) = f_Row(1,atom1) + fg(1)
1139 >                            f_Row(2,atom1) = f_Row(2,atom1) + fg(2)
1140 >                            f_Row(3,atom1) = f_Row(3,atom1) + fg(3)
1141 > #else
1142 >                            f(1,atom1) = f(1,atom1) + fg(1)
1143 >                            f(2,atom1) = f(2,atom1) + fg(2)
1144 >                            f(3,atom1) = f(3,atom1) + fg(3)
1145   #endif
1146 <                      enddo
1147 <
1148 <                      do jb=groupStartCol(j), groupStartCol(j+1)-1
1149 <                         atom2=groupListCol(jb)
748 <                         mf = mfactCol(atom2)
1146 >                            if (n_in_i .gt. 1) then
1147 >                               if (do_stress.and.SIM_uses_AtomicVirial) then
1148 >                                  ! find the distance between the atom and the center of
1149 >                                  ! the cutoff group:
1150   #ifdef IS_MPI
1151 <                         f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1152 <                         f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
752 <                         f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1151 >                                  call get_interatomic_vector(q_Row(:,atom1), &
1152 >                                       q_group_Row(:,i), dag, rag)
1153   #else
1154 <                         f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1155 <                         f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
756 <                         f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1154 >                                  call get_interatomic_vector(q(:,atom1), &
1155 >                                       q_group(:,i), dag, rag)
1156   #endif
1157 <                      enddo
1157 >                                  call add_stress_tensor(dag,fg,tau)
1158 >                               endif
1159 >                            endif
1160 >                         enddo
1161 >                        
1162 >                         do jb=groupStartCol(j), groupStartCol(j+1)-1
1163 >                            atom2=groupListCol(jb)
1164 >                            mf = mfactCol(atom2)
1165 >                            ! fg is the force on atom jb due to cutoff group's
1166 >                            ! presence in switching region
1167 >                            fg = -swderiv*d_grp*mf
1168 > #ifdef IS_MPI
1169 >                            f_Col(1,atom2) = f_Col(1,atom2) + fg(1)
1170 >                            f_Col(2,atom2) = f_Col(2,atom2) + fg(2)
1171 >                            f_Col(3,atom2) = f_Col(3,atom2) + fg(3)
1172 > #else
1173 >                            f(1,atom2) = f(1,atom2) + fg(1)
1174 >                            f(2,atom2) = f(2,atom2) + fg(2)
1175 >                            f(3,atom2) = f(3,atom2) + fg(3)
1176 > #endif
1177 >                            if (n_in_j .gt. 1) then
1178 >                               if (do_stress.and.SIM_uses_AtomicVirial) then
1179 >                                  ! find the distance between the atom and the center of
1180 >                                  ! the cutoff group:
1181 > #ifdef IS_MPI
1182 >                                  call get_interatomic_vector(q_Col(:,atom2), &
1183 >                                       q_group_Col(:,j), dag, rag)
1184 > #else
1185 >                                  call get_interatomic_vector(q(:,atom2), &
1186 >                                       q_group(:,j), dag, rag)
1187 > #endif
1188 >                                  call add_stress_tensor(dag,fg,tau)                              
1189 >                               endif
1190 >                            endif                            
1191 >                         enddo
1192 >                      endif
1193 >                      !if (do_stress.and.(.not.SIM_uses_AtomicVirial)) then
1194 >                      !   call add_stress_tensor(d_grp, fij, tau)
1195 >                      !endif
1196                     endif
760
761                   if (do_stress) call add_stress_tensor(d_grp, fij)
1197                  endif
1198 <             end if
1198 >             endif
1199            enddo
1200 +          
1201         enddo outer
1202  
1203         if (update_nlist) then
# Line 779 | Line 1215 | contains
1215         endif
1216  
1217         if (loop .eq. PREPAIR_LOOP) then
1218 + #ifdef IS_MPI
1219 +          call do_preforce(nlocal, pot_local)
1220 + #else
1221            call do_preforce(nlocal, pot)
1222 + #endif
1223         endif
1224  
1225      enddo
# Line 821 | Line 1261 | contains
1261  
1262      if (do_pot) then
1263         ! scatter/gather pot_row into the members of my column
1264 <       call scatter(pot_Row, pot_Temp, plan_atom_row)
1265 <
1264 >       do i = 1,LR_POT_TYPES
1265 >          call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1266 >       end do
1267         ! scatter/gather pot_local into all other procs
1268         ! add resultant to get total pot
1269         do i = 1, nlocal
1270 <          pot_local = pot_local + pot_Temp(i)
1270 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1271 >               + pot_Temp(1:LR_POT_TYPES,i)
1272         enddo
1273  
1274 +       do i = 1,LR_POT_TYPES
1275 +          particle_pot(1:nlocal) = particle_pot(1:nlocal) + pot_Temp(i,1:nlocal)
1276 +       enddo
1277 +
1278         pot_Temp = 0.0_DP
1279  
1280 <       call scatter(pot_Col, pot_Temp, plan_atom_col)
1280 >       do i = 1,LR_POT_TYPES
1281 >          call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1282 >       end do
1283 >
1284         do i = 1, nlocal
1285 <          pot_local = pot_local + pot_Temp(i)
1285 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1286 >               + pot_Temp(1:LR_POT_TYPES,i)
1287         enddo
1288  
1289 <    endif
1290 < #endif
1289 >       do i = 1,LR_POT_TYPES
1290 >          particle_pot(1:nlocal) = particle_pot(1:nlocal) + pot_Temp(i,1:nlocal)
1291 >       enddo
1292  
842    if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1293  
1294 <       if (FF_uses_RF .and. SIM_uses_RF) then
845 <
846 < #ifdef IS_MPI
847 <          call scatter(rf_Row,rf,plan_atom_row_3d)
848 <          call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
849 <          do i = 1,nlocal
850 <             rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
851 <          end do
1294 >    endif
1295   #endif
1296  
1297 <          do i = 1, nLocal
1297 >    if (SIM_requires_postpair_calc) then
1298 >       do i = 1, nlocal            
1299 >          
1300 >          ! we loop only over the local atoms, so we don't need row and column
1301 >          ! lookups for the types
1302 >          
1303 >          me_i = atid(i)
1304 >          
1305 >          ! is the atom electrostatic?  See if it would have an
1306 >          ! electrostatic interaction with itself
1307 >          iHash = InteractionHash(me_i,me_i)
1308  
1309 <             rfpot = 0.0_DP
1309 >          if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1310   #ifdef IS_MPI
1311 <             me_i = atid_row(i)
1311 >             call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1312 >                  t, do_pot)
1313   #else
1314 <             me_i = atid(i)
1314 >             call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1315 >                  t, do_pot)
1316   #endif
1317 <
1318 <             if (PropertyMap(me_i)%is_Dipole) then
1319 <
1320 <                mu_i = getDipoleMoment(me_i)
1321 <
1322 <                !! The reaction field needs to include a self contribution
1323 <                !! to the field:
1324 <                call accumulate_self_rf(i, mu_i, eFrame)
1325 <                !! Get the reaction field contribution to the
1326 <                !! potential and torques:
1327 <                call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1328 < #ifdef IS_MPI
1329 <                pot_local = pot_local + rfpot
1317 >          endif
1318 >  
1319 >          
1320 >          if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1321 >            
1322 >             ! loop over the excludes to accumulate RF stuff we've
1323 >             ! left out of the normal pair loop
1324 >            
1325 >             do i1 = 1, nSkipsForAtom(i)
1326 >                j = skipsForAtom(i, i1)
1327 >                
1328 >                ! prevent overcounting of the skips
1329 >                if (i.lt.j) then
1330 >                   call get_interatomic_vector(q(:,i), q(:,j), d_atm, ratmsq)
1331 >                   rVal = sqrt(ratmsq)
1332 >                   call get_switch(ratmsq, sw, dswdr, rVal,in_switching_region)
1333 > #ifdef IS_MPI
1334 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1335 >                        vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1336   #else
1337 <                pot = pot + rfpot
1337 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1338 >                        vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1339 > #endif
1340 >                endif
1341 >             enddo
1342 >          endif
1343  
1344 +          if (do_box_dipole) then
1345 + #ifdef IS_MPI
1346 +             call accumulate_box_dipole(i, eFrame, q(:,i), pChg_local, &
1347 +                  nChg_local, pChgPos_local, nChgPos_local, dipVec_local, &
1348 +                  pChgCount_local, nChgCount_local)
1349 + #else
1350 +             call accumulate_box_dipole(i, eFrame, q(:,i), pChg, nChg, &
1351 +                  pChgPos, nChgPos, dipVec, pChgCount, nChgCount)
1352   #endif
1353 <             endif
1354 <          enddo
881 <       endif
1353 >          endif
1354 >       enddo
1355      endif
1356  
884
1357   #ifdef IS_MPI
886
1358      if (do_pot) then
1359 <       pot = pot + pot_local
1360 <       !! we assume the c code will do the allreduce to get the total potential
1361 <       !! we could do it right here if we needed to...
1359 > #ifdef SINGLE_PRECISION
1360 >       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_real,mpi_sum, &
1361 >            mpi_comm_world,mpi_err)            
1362 > #else
1363 >       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision, &
1364 >            mpi_sum, mpi_comm_world,mpi_err)            
1365 > #endif
1366      endif
1367 +        
1368 +    if (do_box_dipole) then
1369  
1370 <    if (do_stress) then
1371 <       call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1372 <            mpi_comm_world,mpi_err)
1373 <       call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1374 <            mpi_comm_world,mpi_err)
1375 <    endif
1376 <
1370 > #ifdef SINGLE_PRECISION
1371 >       call mpi_allreduce(pChg_local, pChg, 1, mpi_real, mpi_sum, &
1372 >            mpi_comm_world, mpi_err)
1373 >       call mpi_allreduce(nChg_local, nChg, 1, mpi_real, mpi_sum, &
1374 >            mpi_comm_world, mpi_err)
1375 >       call mpi_allreduce(pChgCount_local, pChgCount, 1, mpi_integer, mpi_sum,&
1376 >            mpi_comm_world, mpi_err)
1377 >       call mpi_allreduce(nChgCount_local, nChgCount, 1, mpi_integer, mpi_sum,&
1378 >            mpi_comm_world, mpi_err)
1379 >       call mpi_allreduce(pChgPos_local, pChgPos, 3, mpi_real, mpi_sum, &
1380 >            mpi_comm_world, mpi_err)
1381 >       call mpi_allreduce(nChgPos_local, nChgPos, 3, mpi_real, mpi_sum, &
1382 >            mpi_comm_world, mpi_err)
1383 >       call mpi_allreduce(dipVec_local, dipVec, 3, mpi_real, mpi_sum, &
1384 >            mpi_comm_world, mpi_err)
1385   #else
1386 +       call mpi_allreduce(pChg_local, pChg, 1, mpi_double_precision, mpi_sum, &
1387 +            mpi_comm_world, mpi_err)
1388 +       call mpi_allreduce(nChg_local, nChg, 1, mpi_double_precision, mpi_sum, &
1389 +            mpi_comm_world, mpi_err)
1390 +       call mpi_allreduce(pChgCount_local, pChgCount, 1, mpi_integer,&
1391 +            mpi_sum, mpi_comm_world, mpi_err)
1392 +       call mpi_allreduce(nChgCount_local, nChgCount, 1, mpi_integer,&
1393 +            mpi_sum, mpi_comm_world, mpi_err)
1394 +       call mpi_allreduce(pChgPos_local, pChgPos, 3, mpi_double_precision, &
1395 +            mpi_sum, mpi_comm_world, mpi_err)
1396 +       call mpi_allreduce(nChgPos_local, nChgPos, 3, mpi_double_precision, &
1397 +            mpi_sum, mpi_comm_world, mpi_err)
1398 +       call mpi_allreduce(dipVec_local, dipVec, 3, mpi_double_precision, &
1399 +            mpi_sum, mpi_comm_world, mpi_err)
1400 + #endif
1401  
902    if (do_stress) then
903       tau = tau_Temp
904       virial = virial_Temp
1402      endif
1403 <
1403 >    
1404   #endif
1405  
1406 +    if (do_box_dipole) then
1407 +       ! first load the accumulated dipole moment (if dipoles were present)
1408 +       boxDipole(1) = dipVec(1)
1409 +       boxDipole(2) = dipVec(2)
1410 +       boxDipole(3) = dipVec(3)
1411 +
1412 +       ! now include the dipole moment due to charges
1413 +       ! use the lesser of the positive and negative charge totals
1414 +       if (nChg .le. pChg) then
1415 +          chg_value = nChg
1416 +       else
1417 +          chg_value = pChg
1418 +       endif
1419 +      
1420 +       ! find the average positions
1421 +       if (pChgCount .gt. 0 .and. nChgCount .gt. 0) then
1422 +          pChgPos = pChgPos / pChgCount
1423 +          nChgPos = nChgPos / nChgCount
1424 +       endif
1425 +
1426 +       ! dipole is from the negative to the positive (physics notation)
1427 +       chgVec(1) = pChgPos(1) - nChgPos(1)
1428 +       chgVec(2) = pChgPos(2) - nChgPos(2)
1429 +       chgVec(3) = pChgPos(3) - nChgPos(3)
1430 +
1431 +       boxDipole(1) = boxDipole(1) + chgVec(1) * chg_value
1432 +       boxDipole(2) = boxDipole(2) + chgVec(2) * chg_value
1433 +       boxDipole(3) = boxDipole(3) + chgVec(3) * chg_value
1434 +
1435 +    endif
1436 +
1437    end subroutine do_force_loop
1438  
1439    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1440 <       eFrame, A, f, t, pot, vpair, fpair)
1440 >       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1441  
1442 <    real( kind = dp ) :: pot, vpair, sw
1442 >    real( kind = dp ) :: vpair, sw
1443 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1444 >    real( kind = dp ), dimension(nLocal) :: particle_pot
1445      real( kind = dp ), dimension(3) :: fpair
1446      real( kind = dp ), dimension(nLocal)   :: mfact
1447      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 922 | Line 1452 | contains
1452      logical, intent(inout) :: do_pot
1453      integer, intent(in) :: i, j
1454      real ( kind = dp ), intent(inout) :: rijsq
1455 <    real ( kind = dp )                :: r
1455 >    real ( kind = dp ), intent(inout) :: r_grp
1456      real ( kind = dp ), intent(inout) :: d(3)
1457 <    real ( kind = dp ) :: ebalance
1457 >    real ( kind = dp ), intent(inout) :: d_grp(3)
1458 >    real ( kind = dp ), intent(inout) :: rCut
1459 >    real ( kind = dp ) :: r, pair_pot
1460 >    real ( kind = dp ) :: a_k, b_k, c_k, d_k, dx
1461      integer :: me_i, me_j
1462 +    integer :: k
1463  
1464 <    integer :: iMap
1464 >    integer :: iHash
1465  
1466      r = sqrt(rijsq)
1467 <    vpair = 0.0d0
1468 <    fpair(1:3) = 0.0d0
1467 >    
1468 >    vpair = 0.0_dp
1469 >    fpair(1:3) = 0.0_dp
1470  
1471   #ifdef IS_MPI
1472      me_i = atid_row(i)
# Line 941 | Line 1476 | contains
1476      me_j = atid(j)
1477   #endif
1478  
1479 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1480 <
1481 <    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1482 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1479 >    iHash = InteractionHash(me_i, me_j)
1480 >    
1481 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1482 >       call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1483 >            pot(VDW_POT), f, do_pot)
1484      endif
1485 <
1486 <    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1487 <       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1488 <            pot, eFrame, f, t, do_pot)
953 <
954 <       if (FF_uses_dipoles .and. SIM_uses_dipoles) then                
955 <          if ( PropertyMap(me_i)%is_Dipole .and. &
956 <               PropertyMap(me_j)%is_Dipole) then
957 <             if (FF_uses_RF .and. SIM_uses_RF) then
958 <                call accumulate_rf(i, j, r, eFrame, sw)
959 <                call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
960 <             endif
961 <          endif
962 <       endif
1485 >    
1486 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1487 >       call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1488 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1489      endif
1490 <
1491 <    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1490 >    
1491 >    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1492         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1493 <            pot, A, f, t, do_pot)
1493 >            pot(HB_POT), A, f, t, do_pot)
1494      endif
1495 <
1496 <    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1495 >    
1496 >    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1497         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1498 <            pot, A, f, t, do_pot)
1498 >            pot(HB_POT), A, f, t, do_pot)
1499      endif
1500 <
1501 <    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1500 >    
1501 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1502         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1503 <            pot, A, f, t, do_pot)
1503 >            pot(VDW_POT), A, f, t, do_pot)
1504      endif
1505      
1506 <    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1507 <       call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1508 <            pot, A, f, t, do_pot)
1506 >    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1507 >       call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1508 >            pot(VDW_POT), A, f, t, do_pot)
1509      endif
1510 <
1511 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1512 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1513 <            do_pot)
1510 >    
1511 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1512 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1513 >            pot(METALLIC_POT), f, do_pot)
1514      endif
1515 <
1516 <    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1515 >    
1516 >    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1517         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1518 <            pot, A, f, t, do_pot)
1518 >            pot(VDW_POT), A, f, t, do_pot)
1519      endif
1520 <
1521 <    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1520 >    
1521 >    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1522         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1523 <            pot, A, f, t, do_pot)
1523 >            pot(VDW_POT), A, f, t, do_pot)
1524      endif
1525 <    
1525 >
1526 >    if ( iand(iHash, SC_PAIR).ne.0 ) then      
1527 >       call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1528 >            pot(METALLIC_POT), f, do_pot)
1529 >    endif
1530 >    
1531 >    if ( iand(iHash, MNM_PAIR).ne.0 ) then      
1532 >       call do_mnm_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1533 >            pot(VDW_POT), A, f, t, do_pot)
1534 >    endif
1535 >
1536    end subroutine do_pair
1537  
1538 <  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1538 >  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1539         do_pot, do_stress, eFrame, A, f, t, pot)
1540  
1541 <    real( kind = dp ) :: pot, sw
1541 >    real( kind = dp ) :: sw
1542 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1543      real( kind = dp ), dimension(9,nLocal) :: eFrame
1544      real (kind=dp), dimension(9,nLocal) :: A
1545      real (kind=dp), dimension(3,nLocal) :: f
# Line 1010 | Line 1547 | contains
1547  
1548      logical, intent(inout) :: do_pot, do_stress
1549      integer, intent(in) :: i, j
1550 <    real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1550 >    real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1551      real ( kind = dp )                :: r, rc
1552      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1553  
1554 <    integer :: me_i, me_j, iMap
1554 >    integer :: me_i, me_j, iHash
1555  
1556 +    r = sqrt(rijsq)
1557 +    
1558   #ifdef IS_MPI  
1559      me_i = atid_row(i)
1560      me_j = atid_col(j)  
# Line 1024 | Line 1563 | contains
1563      me_j = atid(j)  
1564   #endif
1565  
1566 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1566 >    iHash = InteractionHash(me_i, me_j)
1567  
1568 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1569 <            call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1568 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1569 >            call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1570      endif
1571 +
1572 +    if ( iand(iHash, SC_PAIR).ne.0 ) then      
1573 +            call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1574 +    endif
1575      
1576    end subroutine do_prepair
1577  
1578  
1579    subroutine do_preforce(nlocal,pot)
1580      integer :: nlocal
1581 <    real( kind = dp ) :: pot
1581 >    real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1582  
1583      if (FF_uses_EAM .and. SIM_uses_EAM) then
1584 <       call calc_EAM_preforce_Frho(nlocal,pot)
1584 >       call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1585      endif
1586 <
1587 <
1586 >    if (FF_uses_SC .and. SIM_uses_SC) then
1587 >       call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1588 >    endif
1589    end subroutine do_preforce
1590  
1591  
# Line 1053 | Line 1597 | contains
1597      real( kind = dp ) :: d(3), scaled(3)
1598      integer i
1599  
1600 <    d(1:3) = q_j(1:3) - q_i(1:3)
1600 >    d(1) = q_j(1) - q_i(1)
1601 >    d(2) = q_j(2) - q_i(2)
1602 >    d(3) = q_j(3) - q_i(3)
1603  
1604      ! Wrap back into periodic box if necessary
1605      if ( SIM_uses_PBC ) then
1606  
1607         if( .not.boxIsOrthorhombic ) then
1608            ! calc the scaled coordinates.
1609 +          ! scaled = matmul(HmatInv, d)
1610  
1611 <          scaled = matmul(HmatInv, d)
1612 <
1611 >          scaled(1) = HmatInv(1,1)*d(1) + HmatInv(1,2)*d(2) + HmatInv(1,3)*d(3)
1612 >          scaled(2) = HmatInv(2,1)*d(1) + HmatInv(2,2)*d(2) + HmatInv(2,3)*d(3)
1613 >          scaled(3) = HmatInv(3,1)*d(1) + HmatInv(3,2)*d(2) + HmatInv(3,3)*d(3)
1614 >          
1615            ! wrap the scaled coordinates
1616  
1617 <          scaled = scaled  - anint(scaled)
1617 >          scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1618 >          scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1619 >          scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1620  
1070
1621            ! calc the wrapped real coordinates from the wrapped scaled
1622            ! coordinates
1623 +          ! d = matmul(Hmat,scaled)
1624 +          d(1)= Hmat(1,1)*scaled(1) + Hmat(1,2)*scaled(2) + Hmat(1,3)*scaled(3)
1625 +          d(2)= Hmat(2,1)*scaled(1) + Hmat(2,2)*scaled(2) + Hmat(2,3)*scaled(3)
1626 +          d(3)= Hmat(3,1)*scaled(1) + Hmat(3,2)*scaled(2) + Hmat(3,3)*scaled(3)
1627  
1074          d = matmul(Hmat,scaled)
1075
1628         else
1629            ! calc the scaled coordinates.
1630  
1631 <          do i = 1, 3
1632 <             scaled(i) = d(i) * HmatInv(i,i)
1631 >          scaled(1) = d(1) * HmatInv(1,1)
1632 >          scaled(2) = d(2) * HmatInv(2,2)
1633 >          scaled(3) = d(3) * HmatInv(3,3)
1634 >          
1635 >          ! wrap the scaled coordinates
1636 >          
1637 >          scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1638 >          scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1639 >          scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1640  
1641 <             ! wrap the scaled coordinates
1641 >          ! calc the wrapped real coordinates from the wrapped scaled
1642 >          ! coordinates
1643  
1644 <             scaled(i) = scaled(i) - anint(scaled(i))
1644 >          d(1) = scaled(1)*Hmat(1,1)
1645 >          d(2) = scaled(2)*Hmat(2,2)
1646 >          d(3) = scaled(3)*Hmat(3,3)
1647  
1086             ! calc the wrapped real coordinates from the wrapped scaled
1087             ! coordinates
1088
1089             d(i) = scaled(i)*Hmat(i,i)
1090          enddo
1648         endif
1649  
1650      endif
1651  
1652 <    r_sq = dot_product(d,d)
1652 >    r_sq = d(1)*d(1) + d(2)*d(2) + d(3)*d(3)
1653  
1654    end subroutine get_interatomic_vector
1655  
# Line 1124 | Line 1681 | contains
1681      pot_Col = 0.0_dp
1682      pot_Temp = 0.0_dp
1683  
1127    rf_Row = 0.0_dp
1128    rf_Col = 0.0_dp
1129    rf_Temp = 0.0_dp
1130
1684   #endif
1685  
1686      if (FF_uses_EAM .and. SIM_uses_EAM) then
1687         call clean_EAM()
1688      endif
1689  
1137    rf = 0.0_dp
1138    tau_Temp = 0.0_dp
1139    virial_Temp = 0.0_dp
1690    end subroutine zero_work_arrays
1691  
1692    function skipThisPair(atom1, atom2) result(skip_it)
# Line 1223 | Line 1773 | contains
1773  
1774    function FF_UsesDirectionalAtoms() result(doesit)
1775      logical :: doesit
1776 <    doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1227 <         FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1228 <         FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1776 >    doesit = FF_uses_DirectionalAtoms
1777    end function FF_UsesDirectionalAtoms
1778  
1779    function FF_RequiresPrepairCalc() result(doesit)
1780      logical :: doesit
1781 <    doesit = FF_uses_EAM
1781 >    doesit = FF_uses_EAM .or. FF_uses_SC
1782    end function FF_RequiresPrepairCalc
1783  
1236  function FF_RequiresPostpairCalc() result(doesit)
1237    logical :: doesit
1238    doesit = FF_uses_RF
1239  end function FF_RequiresPostpairCalc
1240
1784   #ifdef PROFILE
1785    function getforcetime() result(totalforcetime)
1786      real(kind=dp) :: totalforcetime
# Line 1247 | Line 1790 | contains
1790  
1791    !! This cleans componets of force arrays belonging only to fortran
1792  
1793 <  subroutine add_stress_tensor(dpair, fpair)
1793 >  subroutine add_stress_tensor(dpair, fpair, tau)
1794  
1795      real( kind = dp ), dimension(3), intent(in) :: dpair, fpair
1796 +    real( kind = dp ), dimension(9), intent(inout) :: tau
1797  
1798      ! because the d vector is the rj - ri vector, and
1799      ! because fx, fy, fz are the force on atom i, we need a
1800      ! negative sign here:  
1801  
1802 <    tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1)
1803 <    tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2)
1804 <    tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3)
1805 <    tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1)
1806 <    tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2)
1807 <    tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3)
1808 <    tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1)
1809 <    tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2)
1810 <    tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3)
1802 >    tau(1) = tau(1) - dpair(1) * fpair(1)
1803 >    tau(2) = tau(2) - dpair(1) * fpair(2)
1804 >    tau(3) = tau(3) - dpair(1) * fpair(3)
1805 >    tau(4) = tau(4) - dpair(2) * fpair(1)
1806 >    tau(5) = tau(5) - dpair(2) * fpair(2)
1807 >    tau(6) = tau(6) - dpair(2) * fpair(3)
1808 >    tau(7) = tau(7) - dpair(3) * fpair(1)
1809 >    tau(8) = tau(8) - dpair(3) * fpair(2)
1810 >    tau(9) = tau(9) - dpair(3) * fpair(3)
1811  
1268    virial_Temp = virial_Temp + &
1269         (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1270
1812    end subroutine add_stress_tensor
1813  
1814   end module doForces

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines