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 2286 by gezelter, Wed Sep 7 22:08:39 2005 UTC vs.
Revision 2503 by gezelter, Thu Dec 8 22:04:40 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.39 2005-09-07 22:08:39 gezelter Exp $, $Date: 2005-09-07 22:08:39 $, $Name: not supported by cvs2svn $, $Revision: 1.39 $
48 > !! @version $Id: doForces.F90,v 1.70 2005-12-08 22:04:40 gezelter Exp $, $Date: 2005-12-08 22:04:40 $, $Name: not supported by cvs2svn $, $Revision: 1.70 $
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 suttonchen
66    use status
67   #ifdef IS_MPI
68    use mpiSimulation
# Line 75 | Line 75 | module doForces
75   #include "UseTheForce/fSwitchingFunction.h"
76   #include "UseTheForce/fCutoffPolicy.h"
77   #include "UseTheForce/DarkSide/fInteractionMap.h"
78 + #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
79  
80  
81    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
# Line 85 | Line 86 | module doForces
86    logical, save :: haveSaneForceField = .false.
87    logical, save :: haveInteractionHash = .false.
88    logical, save :: haveGtypeCutoffMap = .false.
89 <  logical, save :: haveRlist = .false.
89 >  logical, save :: haveDefaultCutoffs = .false.
90 >  logical, save :: haveSkinThickness = .false.
91 >  logical, save :: haveElectrostaticSummationMethod = .false.
92 >  logical, save :: haveCutoffPolicy = .false.
93 >  logical, save :: VisitCutoffsAfterComputing = .false.
94  
95    logical, save :: FF_uses_DirectionalAtoms
96    logical, save :: FF_uses_Dipoles
97    logical, save :: FF_uses_GayBerne
98    logical, save :: FF_uses_EAM
99 <  logical, save :: FF_uses_RF
99 >  logical, save :: FF_uses_SC
100 >  logical, save :: FF_uses_MEAM
101 >
102  
103    logical, save :: SIM_uses_DirectionalAtoms
104    logical, save :: SIM_uses_EAM
105 <  logical, save :: SIM_uses_RF
105 >  logical, save :: SIM_uses_SC
106 >  logical, save :: SIM_uses_MEAM
107    logical, save :: SIM_requires_postpair_calc
108    logical, save :: SIM_requires_prepair_calc
109    logical, save :: SIM_uses_PBC
110  
111 <  integer, save :: corrMethod
111 >  integer, save :: electrostaticSummationMethod
112 >  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
113  
114 +  real(kind=dp), save :: defaultRcut, defaultRsw, largestRcut
115 +  real(kind=dp), save :: skinThickness
116 +  logical, save :: defaultDoShift
117 +
118    public :: init_FF
119 <  public :: setDefaultCutoffs
119 >  public :: setCutoffs
120 >  public :: cWasLame
121 >  public :: setElectrostaticMethod
122 >  public :: setCutoffPolicy
123 >  public :: setSkinThickness
124    public :: do_force_loop
108 –  public :: createInteractionHash
109 –  public :: createGtypeCutoffMap
110 –  public :: getStickyCut
111 –  public :: getStickyPowerCut
112 –  public :: getGayBerneCut
113 –  public :: getEAMCut
114 –  public :: getShapeCut
125  
126   #ifdef PROFILE
127    public :: getforcetime
# Line 124 | Line 134 | module doForces
134    ! Bit hash to determine pair-pair interactions.
135    integer, dimension(:,:), allocatable :: InteractionHash
136    real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
137 <  real(kind=dp), dimension(:), allocatable :: groupMaxCutoff
138 <  integer, dimension(:), allocatable :: groupToGtype
139 <  real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff
137 >  real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
138 >  real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
139 >
140 >  integer, dimension(:), allocatable, target :: groupToGtypeRow
141 >  integer, dimension(:), pointer :: groupToGtypeCol => null()
142 >
143 >  real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
144 >  real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
145    type ::gtypeCutoffs
146       real(kind=dp) :: rcut
147       real(kind=dp) :: rcutsq
# Line 134 | Line 149 | module doForces
149    end type gtypeCutoffs
150    type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
151  
137 –  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
138 –  real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
139 –  
152   contains
153  
154 <  subroutine createInteractionHash(status)
154 >  subroutine createInteractionHash()
155      integer :: nAtypes
144 –    integer, intent(out) :: status
156      integer :: i
157      integer :: j
158      integer :: iHash
# Line 153 | Line 164 | contains
164      logical :: i_is_GB
165      logical :: i_is_EAM
166      logical :: i_is_Shape
167 +    logical :: i_is_SC
168 +    logical :: i_is_MEAM
169      logical :: j_is_LJ
170      logical :: j_is_Elect
171      logical :: j_is_Sticky
# Line 160 | Line 173 | contains
173      logical :: j_is_GB
174      logical :: j_is_EAM
175      logical :: j_is_Shape
176 +    logical :: j_is_SC
177 +    logical :: j_is_MEAM
178      real(kind=dp) :: myRcut
179  
165 –    status = 0  
166 –
180      if (.not. associated(atypes)) then
181 <       call handleError("atype", "atypes was not present before call of createInteractionHash!")
169 <       status = -1
181 >       call handleError("doForces", "atypes was not present before call of createInteractionHash!")
182         return
183      endif
184      
185      nAtypes = getSize(atypes)
186      
187      if (nAtypes == 0) then
188 <       status = -1
188 >       call handleError("doForces", "nAtypes was zero during call of createInteractionHash!")
189         return
190      end if
191  
192      if (.not. allocated(InteractionHash)) then
193         allocate(InteractionHash(nAtypes,nAtypes))
194 +    else
195 +       deallocate(InteractionHash)
196 +       allocate(InteractionHash(nAtypes,nAtypes))
197      endif
198  
199      if (.not. allocated(atypeMaxCutoff)) then
200         allocate(atypeMaxCutoff(nAtypes))
201 +    else
202 +       deallocate(atypeMaxCutoff)
203 +       allocate(atypeMaxCutoff(nAtypes))
204      endif
205          
206      do i = 1, nAtypes
# Line 193 | Line 211 | contains
211         call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
212         call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
213         call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
214 +       call getElementProperty(atypes, i, "is_SC", i_is_SC)
215 +       call getElementProperty(atypes, i, "is_MEAM", i_is_MEAM)
216  
217         do j = i, nAtypes
218  
# Line 206 | Line 226 | contains
226            call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
227            call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
228            call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
229 +          call getElementProperty(atypes, j, "is_SC", j_is_SC)
230 +          call getElementProperty(atypes, j, "is_MEAM", j_is_MEAM)
231  
232            if (i_is_LJ .and. j_is_LJ) then
233               iHash = ior(iHash, LJ_PAIR)            
# Line 227 | Line 249 | contains
249               iHash = ior(iHash, EAM_PAIR)
250            endif
251  
252 +          if (i_is_SC .and. j_is_SC) then
253 +             iHash = ior(iHash, SC_PAIR)
254 +          endif
255 +
256            if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
257            if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
258            if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
# Line 246 | Line 272 | contains
272      haveInteractionHash = .true.
273    end subroutine createInteractionHash
274  
275 <  subroutine createGtypeCutoffMap(stat)
275 >  subroutine createGtypeCutoffMap()
276  
251 –    integer, intent(out), optional :: stat
277      logical :: i_is_LJ
278      logical :: i_is_Elect
279      logical :: i_is_Sticky
# Line 259 | Line 284 | contains
284      logical :: GtypeFound
285  
286      integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
287 <    integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
288 <    real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
287 >    integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
288 >    integer :: nGroupsInRow
289 >    integer :: nGroupsInCol
290 >    integer :: nGroupTypesRow,nGroupTypesCol
291 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol
292      real(kind=dp) :: biggestAtypeCutoff
293  
266 –    stat = 0
294      if (.not. haveInteractionHash) then
295 <       call createInteractionHash(myStatus)      
269 <       if (myStatus .ne. 0) then
270 <          write(default_error, *) 'createInteractionHash failed in doForces!'
271 <          stat = -1
272 <          return
273 <       endif
295 >       call createInteractionHash()      
296      endif
297 <
297 > #ifdef IS_MPI
298 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
299 >    nGroupsInCol = getNgroupsInCol(plan_group_col)
300 > #endif
301      nAtypes = getSize(atypes)
302 <    
302 > ! Set all of the initial cutoffs to zero.
303 >    atypeMaxCutoff = 0.0_dp
304      do i = 1, nAtypes
305         if (SimHasAtype(i)) then    
306            call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
# Line 285 | Line 311 | contains
311            call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
312            call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
313            
314 <          atypeMaxCutoff(i) = 0.0_dp
315 <          if (i_is_LJ) then
316 <             thisRcut = getSigma(i) * 2.5_dp
317 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
318 <          endif
319 <          if (i_is_Elect) then
320 <             thisRcut = defaultRcut
321 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
322 <          endif
323 <          if (i_is_Sticky) then
324 <             thisRcut = getStickyCut(i)
325 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
326 <          endif
327 <          if (i_is_StickyP) then
328 <             thisRcut = getStickyPowerCut(i)
329 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
330 <          endif
331 <          if (i_is_GB) then
332 <             thisRcut = getGayBerneCut(i)
333 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
334 <          endif
335 <          if (i_is_EAM) then
336 <             thisRcut = getEAMCut(i)
337 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
314 >
315 >          if (haveDefaultCutoffs) then
316 >             atypeMaxCutoff(i) = defaultRcut
317 >          else
318 >             if (i_is_LJ) then          
319 >                thisRcut = getSigma(i) * 2.5_dp
320 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
321 >             endif
322 >             if (i_is_Elect) then
323 >                thisRcut = defaultRcut
324 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
325 >             endif
326 >             if (i_is_Sticky) then
327 >                thisRcut = getStickyCut(i)
328 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
329 >             endif
330 >             if (i_is_StickyP) then
331 >                thisRcut = getStickyPowerCut(i)
332 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
333 >             endif
334 >             if (i_is_GB) then
335 >                thisRcut = getGayBerneCut(i)
336 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
337 >             endif
338 >             if (i_is_EAM) then
339 >                thisRcut = getEAMCut(i)
340 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
341 >             endif
342 >             if (i_is_Shape) then
343 >                thisRcut = getShapeCut(i)
344 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
345 >             endif
346            endif
347 <          if (i_is_Shape) then
314 <             thisRcut = getShapeCut(i)
315 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
316 <          endif
317 <          
347 >                    
348            if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
349               biggestAtypeCutoff = atypeMaxCutoff(i)
350            endif
351 +
352         endif
353      enddo
323 –  
324 –    nGroupTypes = 0
354      
355      istart = 1
356 +    jstart = 1
357   #ifdef IS_MPI
358      iend = nGroupsInRow
359 +    jend = nGroupsInCol
360   #else
361      iend = nGroups
362 +    jend = nGroups
363   #endif
364      
365      !! allocate the groupToGtype and gtypeMaxCutoff here.
366 <    if(.not.allocated(groupToGtype)) then
367 <       allocate(groupToGtype(iend))
368 <       allocate(groupMaxCutoff(iend))
369 <       allocate(gtypeMaxCutoff(iend))
366 >    if(.not.allocated(groupToGtypeRow)) then
367 >     !  allocate(groupToGtype(iend))
368 >       allocate(groupToGtypeRow(iend))
369 >    else
370 >       deallocate(groupToGtypeRow)
371 >       allocate(groupToGtypeRow(iend))
372      endif
373 +    if(.not.allocated(groupMaxCutoffRow)) then
374 +       allocate(groupMaxCutoffRow(iend))
375 +    else
376 +       deallocate(groupMaxCutoffRow)
377 +       allocate(groupMaxCutoffRow(iend))
378 +    end if
379 +
380 +    if(.not.allocated(gtypeMaxCutoffRow)) then
381 +       allocate(gtypeMaxCutoffRow(iend))
382 +    else
383 +       deallocate(gtypeMaxCutoffRow)
384 +       allocate(gtypeMaxCutoffRow(iend))
385 +    endif
386 +
387 +
388 + #ifdef IS_MPI
389 +       ! We only allocate new storage if we are in MPI because Ncol /= Nrow
390 +    if(.not.associated(groupToGtypeCol)) then
391 +       allocate(groupToGtypeCol(jend))
392 +    else
393 +       deallocate(groupToGtypeCol)
394 +       allocate(groupToGtypeCol(jend))
395 +    end if
396 +
397 +    if(.not.associated(groupToGtypeCol)) then
398 +       allocate(groupToGtypeCol(jend))
399 +    else
400 +       deallocate(groupToGtypeCol)
401 +       allocate(groupToGtypeCol(jend))
402 +    end if
403 +    if(.not.associated(gtypeMaxCutoffCol)) then
404 +       allocate(gtypeMaxCutoffCol(jend))
405 +    else
406 +       deallocate(gtypeMaxCutoffCol)      
407 +       allocate(gtypeMaxCutoffCol(jend))
408 +    end if
409 +
410 +       groupMaxCutoffCol = 0.0_dp
411 +       gtypeMaxCutoffCol = 0.0_dp
412 +
413 + #endif
414 +       groupMaxCutoffRow = 0.0_dp
415 +       gtypeMaxCutoffRow = 0.0_dp
416 +
417 +
418      !! first we do a single loop over the cutoff groups to find the
419      !! largest cutoff for any atypes present in this group.  We also
420      !! create gtypes at this point.
421      
422      tol = 1.0d-6
423 <    
423 >    nGroupTypesRow = 0
424 >
425      do i = istart, iend      
426         n_in_i = groupStartRow(i+1) - groupStartRow(i)
427 <       groupMaxCutoff(i) = 0.0_dp
427 >       groupMaxCutoffRow(i) = 0.0_dp
428         do ia = groupStartRow(i), groupStartRow(i+1)-1
429            atom1 = groupListRow(ia)
430   #ifdef IS_MPI
# Line 352 | Line 432 | contains
432   #else
433            me_i = atid(atom1)
434   #endif          
435 <          if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
436 <             groupMaxCutoff(i)=atypeMaxCutoff(me_i)
435 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
436 >             groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
437            endif          
438         enddo
439  
440 <       if (nGroupTypes.eq.0) then
441 <          nGroupTypes = nGroupTypes + 1
442 <          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
443 <          groupToGtype(i) = nGroupTypes
440 >       if (nGroupTypesRow.eq.0) then
441 >          nGroupTypesRow = nGroupTypesRow + 1
442 >          gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
443 >          groupToGtypeRow(i) = nGroupTypesRow
444         else
445            GtypeFound = .false.
446 <          do g = 1, nGroupTypes
447 <             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
448 <                groupToGtype(i) = g
446 >          do g = 1, nGroupTypesRow
447 >             if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
448 >                groupToGtypeRow(i) = g
449                  GtypeFound = .true.
450               endif
451            enddo
452            if (.not.GtypeFound) then            
453 <             nGroupTypes = nGroupTypes + 1
454 <             gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
455 <             groupToGtype(i) = nGroupTypes
453 >             nGroupTypesRow = nGroupTypesRow + 1
454 >             gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
455 >             groupToGtypeRow(i) = nGroupTypesRow
456            endif
457         endif
458      enddo    
459  
460 + #ifdef IS_MPI
461 +    do j = jstart, jend      
462 +       n_in_j = groupStartCol(j+1) - groupStartCol(j)
463 +       groupMaxCutoffCol(j) = 0.0_dp
464 +       do ja = groupStartCol(j), groupStartCol(j+1)-1
465 +          atom1 = groupListCol(ja)
466 +
467 +          me_j = atid_col(atom1)
468 +
469 +          if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
470 +             groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
471 +          endif          
472 +       enddo
473 +
474 +       if (nGroupTypesCol.eq.0) then
475 +          nGroupTypesCol = nGroupTypesCol + 1
476 +          gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
477 +          groupToGtypeCol(j) = nGroupTypesCol
478 +       else
479 +          GtypeFound = .false.
480 +          do g = 1, nGroupTypesCol
481 +             if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
482 +                groupToGtypeCol(j) = g
483 +                GtypeFound = .true.
484 +             endif
485 +          enddo
486 +          if (.not.GtypeFound) then            
487 +             nGroupTypesCol = nGroupTypesCol + 1
488 +             gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
489 +             groupToGtypeCol(j) = nGroupTypesCol
490 +          endif
491 +       endif
492 +    enddo    
493 +
494 + #else
495 + ! Set pointers to information we just found
496 +    nGroupTypesCol = nGroupTypesRow
497 +    groupToGtypeCol => groupToGtypeRow
498 +    gtypeMaxCutoffCol => gtypeMaxCutoffRow
499 +    groupMaxCutoffCol => groupMaxCutoffRow
500 + #endif
501 +
502      !! allocate the gtypeCutoffMap here.
503 <    allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
503 >    allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
504      !! then we do a double loop over all the group TYPES to find the cutoff
505      !! map between groups of two types
506 <    
507 <    do i = 1, nGroupTypes
508 <       do j = 1, nGroupTypes
506 >    tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
507 >
508 >    do i = 1, nGroupTypesRow      
509 >       do j = 1, nGroupTypesCol
510        
511            select case(cutoffPolicy)
512            case(TRADITIONAL_CUTOFF_POLICY)
513 <             thisRcut = maxval(gtypeMaxCutoff)
513 >             thisRcut = tradRcut
514            case(MIX_CUTOFF_POLICY)
515 <             thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
515 >             thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
516            case(MAX_CUTOFF_POLICY)
517 <             thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
517 >             thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
518            case default
519               call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
520               return
521            end select
522            gtypeCutoffMap(i,j)%rcut = thisRcut
523 +          
524 +          if (thisRcut.gt.largestRcut) largestRcut = thisRcut
525 +
526            gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
401 –          skin = defaultRlist - defaultRcut
402 –          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
527  
528 +          if (.not.haveSkinThickness) then
529 +             skinThickness = 1.0_dp
530 +          endif
531 +
532 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
533 +
534 +          ! sanity check
535 +
536 +          if (haveDefaultCutoffs) then
537 +             if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
538 +                call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
539 +             endif
540 +          endif
541         enddo
542      enddo
543 +
544 +    if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
545 +    if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
546 +    if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
547 + #ifdef IS_MPI
548 +    if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
549 +    if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
550 + #endif
551 +    groupMaxCutoffCol => null()
552 +    gtypeMaxCutoffCol => null()
553      
554      haveGtypeCutoffMap = .true.
555 <    
409 <  end subroutine createGtypeCutoffMap
410 <  
411 <  subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
412 <    real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
413 <    integer, intent(in) :: cutPolicy
414 <    
415 <    defaultRcut = defRcut
416 <    defaultRsw = defRsw
417 <    defaultRlist = defRlist
418 <    cutoffPolicy = cutPolicy
419 <  end subroutine setDefaultCutoffs
420 <  
421 <  subroutine setCutoffPolicy(cutPolicy)
555 >   end subroutine createGtypeCutoffMap
556  
557 +   subroutine setCutoffs(defRcut, defRsw)
558 +
559 +     real(kind=dp),intent(in) :: defRcut, defRsw
560 +     character(len = statusMsgSize) :: errMsg
561 +     integer :: localError
562 +
563 +     defaultRcut = defRcut
564 +     defaultRsw = defRsw
565 +    
566 +     defaultDoShift = .false.
567 +     if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
568 +        
569 +        write(errMsg, *) &
570 +             'cutoffRadius and switchingRadius are set to the same', newline &
571 +             // tab, 'value.  OOPSE will use shifted ', newline &
572 +             // tab, 'potentials instead of switching functions.'
573 +        
574 +        call handleInfo("setCutoffs", errMsg)
575 +        
576 +        defaultDoShift = .true.
577 +        
578 +     endif
579 +
580 +     localError = 0
581 +     call setLJDefaultCutoff( defaultRcut, defaultDoShift )
582 +     call setCutoffEAM( defaultRcut, localError)
583 +     if (localError /= 0) then
584 +       write(errMsg, *) 'An error has occured in setting the EAM cutoff'
585 +       call handleError("setCutoffs", errMsg)
586 +     end if
587 +     call set_switch(GROUP_SWITCH, defaultRsw, defaultRcut)
588 +    
589 +     haveDefaultCutoffs = .true.
590 +   end subroutine setCutoffs
591 +
592 +   subroutine cWasLame()
593 +    
594 +     VisitCutoffsAfterComputing = .true.
595 +     return
596 +    
597 +   end subroutine cWasLame
598 +  
599 +   subroutine setCutoffPolicy(cutPolicy)
600 +    
601       integer, intent(in) :: cutPolicy
602 +    
603       cutoffPolicy = cutPolicy
604 <     call createGtypeCutoffMap()
604 >     haveCutoffPolicy = .true.
605 >     write(*,*) 'have cutoffPolicy in F = ', cutPolicy
606  
607 <   end subroutine setCutoffPolicy
428 <    
607 >     call createGtypeCutoffMap()
608      
609 <  subroutine setSimVariables()
610 <    SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
611 <    SIM_uses_EAM = SimUsesEAM()
433 <    SIM_uses_RF = SimUsesRF()
434 <    SIM_requires_postpair_calc = SimRequiresPostpairCalc()
435 <    SIM_requires_prepair_calc = SimRequiresPrepairCalc()
436 <    SIM_uses_PBC = SimUsesPBC()
609 >   end subroutine setCutoffPolicy
610 >  
611 >   subroutine setElectrostaticMethod( thisESM )
612  
613 <    haveSIMvariables = .true.
613 >     integer, intent(in) :: thisESM
614  
615 <    return
616 <  end subroutine setSimVariables
615 >     electrostaticSummationMethod = thisESM
616 >     haveElectrostaticSummationMethod = .true.
617 >    
618 >   end subroutine setElectrostaticMethod
619  
620 +   subroutine setSkinThickness( thisSkin )
621 +    
622 +     real(kind=dp), intent(in) :: thisSkin
623 +    
624 +     skinThickness = thisSkin
625 +     haveSkinThickness = .true.
626 +    
627 +     call createGtypeCutoffMap()
628 +    
629 +   end subroutine setSkinThickness
630 +      
631 +   subroutine setSimVariables()
632 +     SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
633 +     SIM_uses_EAM = SimUsesEAM()
634 +     SIM_requires_postpair_calc = SimRequiresPostpairCalc()
635 +     SIM_requires_prepair_calc = SimRequiresPrepairCalc()
636 +     SIM_uses_PBC = SimUsesPBC()
637 +    
638 +     haveSIMvariables = .true.
639 +    
640 +     return
641 +   end subroutine setSimVariables
642 +
643    subroutine doReadyCheck(error)
644      integer, intent(out) :: error
645  
# Line 448 | Line 648 | contains
648      error = 0
649  
650      if (.not. haveInteractionHash) then      
651 <       myStatus = 0      
452 <       call createInteractionHash(myStatus)      
453 <       if (myStatus .ne. 0) then
454 <          write(default_error, *) 'createInteractionHash failed in doForces!'
455 <          error = -1
456 <          return
457 <       endif
651 >       call createInteractionHash()      
652      endif
653  
654      if (.not. haveGtypeCutoffMap) then        
655 <       myStatus = 0      
462 <       call createGtypeCutoffMap(myStatus)      
463 <       if (myStatus .ne. 0) then
464 <          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
465 <          error = -1
466 <          return
467 <       endif
655 >       call createGtypeCutoffMap()      
656      endif
657  
658 +
659 +    if (VisitCutoffsAfterComputing) then
660 +       call set_switch(GROUP_SWITCH, largestRcut, largestRcut)      
661 +    endif
662 +
663 +
664      if (.not. haveSIMvariables) then
665         call setSimVariables()
666      endif
# Line 500 | Line 694 | contains
694    end subroutine doReadyCheck
695  
696  
697 <  subroutine init_FF(use_RF, use_UW, use_DW, thisStat)
697 >  subroutine init_FF(thisStat)
698  
505 –    logical, intent(in) :: use_RF
506 –    logical, intent(in) :: use_UW
507 –    logical, intent(in) :: use_DW
699      integer, intent(out) :: thisStat  
700      integer :: my_status, nMatches
510 –    integer :: corrMethod
701      integer, pointer :: MatchList(:) => null()
512 –    real(kind=dp) :: rcut, rrf, rt, dielect
702  
703      !! assume things are copacetic, unless they aren't
704      thisStat = 0
705  
517 –    !! Fortran's version of a cast:
518 –    FF_uses_RF = use_RF
519 –
520 –    !! set the electrostatic correction method
521 –    if (use_UW) then
522 –       corrMethod = 1
523 –    elseif (use_DW) then
524 –       corrMethod = 2
525 –    else
526 –       corrMethod = 0
527 –    endif
528 –    
706      !! init_FF is called *after* all of the atom types have been
707      !! defined in atype_module using the new_atype subroutine.
708      !!
# Line 555 | Line 732 | contains
732  
733      haveSaneForceField = .true.
734  
558 –    !! check to make sure the FF_uses_RF setting makes sense
559 –
560 –    if (FF_uses_Dipoles) then
561 –       if (FF_uses_RF) then
562 –          dielect = getDielect()
563 –          call initialize_rf(dielect)
564 –       endif
565 –    else
566 –       if (FF_uses_RF) then          
567 –          write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
568 –          thisStat = -1
569 –          haveSaneForceField = .false.
570 –          return
571 –       endif
572 –    endif
573 –
735      if (FF_uses_EAM) then
736         call init_EAM_FF(my_status)
737         if (my_status /= 0) then
# Line 581 | Line 742 | contains
742         end if
743      endif
744  
584 –    if (FF_uses_GayBerne) then
585 –       call check_gb_pair_FF(my_status)
586 –       if (my_status .ne. 0) then
587 –          thisStat = -1
588 –          haveSaneForceField = .false.
589 –          return
590 –       endif
591 –    endif
592 –
745      if (.not. haveNeighborList) then
746         !! Create neighbor lists
747         call expandNeighborList(nLocal, my_status)
# Line 623 | Line 775 | contains
775  
776      !! Stress Tensor
777      real( kind = dp), dimension(9) :: tau  
778 <    real ( kind = dp ) :: pot
778 >    real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
779      logical ( kind = 2) :: do_pot_c, do_stress_c
780      logical :: do_pot
781      logical :: do_stress
782      logical :: in_switching_region
783   #ifdef IS_MPI
784 <    real( kind = DP ) :: pot_local
784 >    real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
785      integer :: nAtomsInRow
786      integer :: nAtomsInCol
787      integer :: nprocs
# Line 644 | Line 796 | contains
796      integer :: nlist
797      real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
798      real( kind = DP ) :: sw, dswdr, swderiv, mf
799 +    real( kind = DP ) :: rVal
800      real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
801      real(kind=dp) :: rfpot, mu_i, virial
802 +    real(kind=dp):: rCut
803      integer :: me_i, me_j, n_in_i, n_in_j
804      logical :: is_dp_i
805      integer :: neighborListSize
# Line 654 | Line 808 | contains
808      integer :: propPack_i, propPack_j
809      integer :: loopStart, loopEnd, loop
810      integer :: iHash
811 <    real(kind=dp) :: listSkin = 1.0  
811 >    integer :: i1
812 >  
813  
814      !! initialize local variables  
815  
# Line 718 | Line 873 | contains
873         ! (but only on the first time through):
874         if (loop .eq. loopStart) then
875   #ifdef IS_MPI
876 <          call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
876 >          call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
877                 update_nlist)
878   #else
879 <          call checkNeighborList(nGroups, q_group, listSkin, &
879 >          call checkNeighborList(nGroups, q_group, skinThickness, &
880                 update_nlist)
881   #endif
882         endif
# Line 779 | Line 934 | contains
934               me_j = atid(j)
935               call get_interatomic_vector(q_group(:,i), &
936                    q_group(:,j), d_grp, rgrpsq)
937 < #endif
937 > #endif      
938  
939 <             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
939 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
940                  if (update_nlist) then
941                     nlist = nlist + 1
942  
# Line 801 | Line 956 | contains
956  
957                     list(nlist) = j
958                  endif
959 +
960  
961 <                if (loop .eq. PAIR_LOOP) then
962 <                   vij = 0.0d0
807 <                   fij(1:3) = 0.0d0
808 <                endif
961 >                
962 >                if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
963  
964 <                call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
965 <                     in_switching_region)
966 <
967 <                n_in_j = groupStartCol(j+1) - groupStartCol(j)
968 <
969 <                do ia = groupStartRow(i), groupStartRow(i+1)-1
970 <
971 <                   atom1 = groupListRow(ia)
972 <
973 <                   inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
974 <
975 <                      atom2 = groupListCol(jb)
976 <
977 <                      if (skipThisPair(atom1, atom2)) cycle inner
978 <
979 <                      if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
980 <                         d_atm(1:3) = d_grp(1:3)
981 <                         ratmsq = rgrpsq
982 <                      else
983 < #ifdef IS_MPI
984 <                         call get_interatomic_vector(q_Row(:,atom1), &
985 <                              q_Col(:,atom2), d_atm, ratmsq)
986 < #else
987 <                         call get_interatomic_vector(q(:,atom1), &
988 <                              q(:,atom2), d_atm, ratmsq)
989 < #endif
990 <                      endif
991 <
992 <                      if (loop .eq. PREPAIR_LOOP) then
993 < #ifdef IS_MPI                      
994 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
995 <                              rgrpsq, d_grp, do_pot, do_stress, &
996 <                              eFrame, A, f, t, pot_local)
997 < #else
998 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
999 <                              rgrpsq, d_grp, do_pot, do_stress, &
1000 <                              eFrame, A, f, t, pot)
964 >                   rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
965 >                   if (loop .eq. PAIR_LOOP) then
966 >                      vij = 0.0d0
967 >                      fij(1:3) = 0.0d0
968 >                   endif
969 >                  
970 >                   call get_switch(rgrpsq, sw, dswdr, rgrp, &
971 >                        group_switch, in_switching_region)
972 >                  
973 >                   n_in_j = groupStartCol(j+1) - groupStartCol(j)
974 >                  
975 >                   do ia = groupStartRow(i), groupStartRow(i+1)-1
976 >                      
977 >                      atom1 = groupListRow(ia)
978 >                      
979 >                      inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
980 >                        
981 >                         atom2 = groupListCol(jb)
982 >                        
983 >                         if (skipThisPair(atom1, atom2))  cycle inner
984 >                        
985 >                         if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
986 >                            d_atm(1:3) = d_grp(1:3)
987 >                            ratmsq = rgrpsq
988 >                         else
989 > #ifdef IS_MPI
990 >                            call get_interatomic_vector(q_Row(:,atom1), &
991 >                                 q_Col(:,atom2), d_atm, ratmsq)
992 > #else
993 >                            call get_interatomic_vector(q(:,atom1), &
994 >                                 q(:,atom2), d_atm, ratmsq)
995 > #endif
996 >                         endif
997 >                        
998 >                         if (loop .eq. PREPAIR_LOOP) then
999 > #ifdef IS_MPI                      
1000 >                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1001 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1002 >                                 eFrame, A, f, t, pot_local)
1003 > #else
1004 >                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1005 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1006 >                                 eFrame, A, f, t, pot)
1007   #endif                                              
1008 <                      else
1008 >                         else
1009   #ifdef IS_MPI                      
1010 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1011 <                              do_pot, &
1012 <                              eFrame, A, f, t, pot_local, vpair, fpair)
1010 >                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1011 >                                 do_pot, eFrame, A, f, t, pot_local, vpair, &
1012 >                                 fpair, d_grp, rgrp, rCut)
1013   #else
1014 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1015 <                              do_pot,  &
1016 <                              eFrame, A, f, t, pot, vpair, fpair)
1014 >                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1015 >                                 do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1016 >                                 d_grp, rgrp, rCut)
1017   #endif
1018 +                            vij = vij + vpair
1019 +                            fij(1:3) = fij(1:3) + fpair(1:3)
1020 +                         endif
1021 +                      enddo inner
1022 +                   enddo
1023  
1024 <                         vij = vij + vpair
1025 <                         fij(1:3) = fij(1:3) + fpair(1:3)
1026 <                      endif
1027 <                   enddo inner
1028 <                enddo
1029 <
1030 <                if (loop .eq. PAIR_LOOP) then
1031 <                   if (in_switching_region) then
1032 <                      swderiv = vij*dswdr/rgrp
1033 <                      fij(1) = fij(1) + swderiv*d_grp(1)
869 <                      fij(2) = fij(2) + swderiv*d_grp(2)
870 <                      fij(3) = fij(3) + swderiv*d_grp(3)
871 <
872 <                      do ia=groupStartRow(i), groupStartRow(i+1)-1
873 <                         atom1=groupListRow(ia)
874 <                         mf = mfactRow(atom1)
1024 >                   if (loop .eq. PAIR_LOOP) then
1025 >                      if (in_switching_region) then
1026 >                         swderiv = vij*dswdr/rgrp
1027 >                         fij(1) = fij(1) + swderiv*d_grp(1)
1028 >                         fij(2) = fij(2) + swderiv*d_grp(2)
1029 >                         fij(3) = fij(3) + swderiv*d_grp(3)
1030 >                        
1031 >                         do ia=groupStartRow(i), groupStartRow(i+1)-1
1032 >                            atom1=groupListRow(ia)
1033 >                            mf = mfactRow(atom1)
1034   #ifdef IS_MPI
1035 <                         f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1036 <                         f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1037 <                         f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1035 >                            f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1036 >                            f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1037 >                            f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1038   #else
1039 <                         f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1040 <                         f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1041 <                         f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1039 >                            f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1040 >                            f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1041 >                            f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1042   #endif
1043 <                      enddo
1044 <
1045 <                      do jb=groupStartCol(j), groupStartCol(j+1)-1
1046 <                         atom2=groupListCol(jb)
1047 <                         mf = mfactCol(atom2)
1043 >                         enddo
1044 >                        
1045 >                         do jb=groupStartCol(j), groupStartCol(j+1)-1
1046 >                            atom2=groupListCol(jb)
1047 >                            mf = mfactCol(atom2)
1048   #ifdef IS_MPI
1049 <                         f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1050 <                         f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1051 <                         f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1049 >                            f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1050 >                            f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1051 >                            f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1052   #else
1053 <                         f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1054 <                         f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1055 <                         f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1053 >                            f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1054 >                            f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1055 >                            f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1056   #endif
1057 <                      enddo
1058 <                   endif
1057 >                         enddo
1058 >                      endif
1059  
1060 <                   if (do_stress) call add_stress_tensor(d_grp, fij)
1060 >                      if (do_stress) call add_stress_tensor(d_grp, fij)
1061 >                   endif
1062                  endif
1063 <             end if
1063 >             endif
1064            enddo
1065 +          
1066         enddo outer
1067  
1068         if (update_nlist) then
# Line 961 | Line 1122 | contains
1122  
1123      if (do_pot) then
1124         ! scatter/gather pot_row into the members of my column
1125 <       call scatter(pot_Row, pot_Temp, plan_atom_row)
1126 <
1125 >       do i = 1,LR_POT_TYPES
1126 >          call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1127 >       end do
1128         ! scatter/gather pot_local into all other procs
1129         ! add resultant to get total pot
1130         do i = 1, nlocal
1131 <          pot_local = pot_local + pot_Temp(i)
1131 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1132 >               + pot_Temp(1:LR_POT_TYPES,i)
1133         enddo
1134  
1135         pot_Temp = 0.0_DP
1136 <
1137 <       call scatter(pot_Col, pot_Temp, plan_atom_col)
1136 >       do i = 1,LR_POT_TYPES
1137 >          call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1138 >       end do
1139         do i = 1, nlocal
1140 <          pot_local = pot_local + pot_Temp(i)
1140 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1141 >               + pot_Temp(1:LR_POT_TYPES,i)
1142         enddo
1143  
1144      endif
1145   #endif
1146  
1147 <    if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1147 >    if (SIM_requires_postpair_calc) then
1148 >       do i = 1, nlocal            
1149 >          
1150 >          ! we loop only over the local atoms, so we don't need row and column
1151 >          ! lookups for the types
1152 >          
1153 >          me_i = atid(i)
1154 >          
1155 >          ! is the atom electrostatic?  See if it would have an
1156 >          ! electrostatic interaction with itself
1157 >          iHash = InteractionHash(me_i,me_i)
1158  
1159 <       if (FF_uses_RF .and. SIM_uses_RF) then
985 <
1159 >          if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1160   #ifdef IS_MPI
1161 <          call scatter(rf_Row,rf,plan_atom_row_3d)
1162 <          call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
1163 <          do i = 1,nlocal
1164 <             rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
1165 <          end do
1161 >             call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1162 >                  t, do_pot)
1163 > #else
1164 >             call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1165 >                  t, do_pot)
1166   #endif
1167 <
1168 <          do i = 1, nLocal
1169 <
1170 <             rfpot = 0.0_DP
997 < #ifdef IS_MPI
998 <             me_i = atid_row(i)
999 < #else
1000 <             me_i = atid(i)
1001 < #endif
1002 <             iHash = InteractionHash(me_i,me_j)
1167 >          endif
1168 >  
1169 >          
1170 >          if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1171              
1172 <             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1173 <
1174 <                mu_i = getDipoleMoment(me_i)
1175 <
1176 <                !! The reaction field needs to include a self contribution
1177 <                !! to the field:
1178 <                call accumulate_self_rf(i, mu_i, eFrame)
1179 <                !! Get the reaction field contribution to the
1180 <                !! potential and torques:
1181 <                call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1172 >             ! loop over the excludes to accumulate RF stuff we've
1173 >             ! left out of the normal pair loop
1174 >            
1175 >             do i1 = 1, nSkipsForAtom(i)
1176 >                j = skipsForAtom(i, i1)
1177 >                
1178 >                ! prevent overcounting of the skips
1179 >                if (i.lt.j) then
1180 >                   call get_interatomic_vector(q(:,i), &
1181 >                        q(:,j), d_atm, ratmsq)
1182 >                   rVal = dsqrt(ratmsq)
1183 >                   call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1184 >                        in_switching_region)
1185   #ifdef IS_MPI
1186 <                pot_local = pot_local + rfpot
1186 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1187 >                        vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1188   #else
1189 <                pot = pot + rfpot
1190 <
1189 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1190 >                        vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1191   #endif
1192 <             endif
1193 <          enddo
1194 <       endif
1192 >                endif
1193 >             enddo
1194 >          endif
1195 >       enddo
1196      endif
1197 <
1025 <
1197 >    
1198   #ifdef IS_MPI
1199 <
1199 >    
1200      if (do_pot) then
1201 <       pot = pot + pot_local
1202 <       !! we assume the c code will do the allreduce to get the total potential
1031 <       !! we could do it right here if we needed to...
1201 >       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1202 >            mpi_comm_world,mpi_err)            
1203      endif
1204 <
1204 >    
1205      if (do_stress) then
1206         call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1207              mpi_comm_world,mpi_err)
1208         call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1209              mpi_comm_world,mpi_err)
1210      endif
1211 <
1211 >    
1212   #else
1213 <
1213 >    
1214      if (do_stress) then
1215         tau = tau_Temp
1216         virial = virial_Temp
1217      endif
1218 <
1218 >    
1219   #endif
1220 <
1220 >    
1221    end subroutine do_force_loop
1222  
1223    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1224 <       eFrame, A, f, t, pot, vpair, fpair)
1224 >       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1225  
1226 <    real( kind = dp ) :: pot, vpair, sw
1226 >    real( kind = dp ) :: vpair, sw
1227 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1228      real( kind = dp ), dimension(3) :: fpair
1229      real( kind = dp ), dimension(nLocal)   :: mfact
1230      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 1063 | Line 1235 | contains
1235      logical, intent(inout) :: do_pot
1236      integer, intent(in) :: i, j
1237      real ( kind = dp ), intent(inout) :: rijsq
1238 <    real ( kind = dp )                :: r
1238 >    real ( kind = dp ), intent(inout) :: r_grp
1239      real ( kind = dp ), intent(inout) :: d(3)
1240 <    real ( kind = dp ) :: ebalance
1240 >    real ( kind = dp ), intent(inout) :: d_grp(3)
1241 >    real ( kind = dp ), intent(inout) :: rCut
1242 >    real ( kind = dp ) :: r
1243      integer :: me_i, me_j
1244  
1245      integer :: iHash
# Line 1083 | Line 1257 | contains
1257   #endif
1258  
1259      iHash = InteractionHash(me_i, me_j)
1260 <
1260 >    
1261      if ( iand(iHash, LJ_PAIR).ne.0 ) then
1262 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1262 >       call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1263 >            pot(VDW_POT), f, do_pot)
1264      endif
1265 <
1265 >    
1266      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1267 <       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1268 <            pot, eFrame, f, t, do_pot, corrMethod)
1094 <
1095 <       if (FF_uses_RF .and. SIM_uses_RF) then
1096 <
1097 <          ! CHECK ME (RF needs to know about all electrostatic types)
1098 <          call accumulate_rf(i, j, r, eFrame, sw)
1099 <          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1100 <       endif
1101 <
1267 >       call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1268 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1269      endif
1270 <
1270 >    
1271      if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1272         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1273 <            pot, A, f, t, do_pot)
1273 >            pot(HB_POT), A, f, t, do_pot)
1274      endif
1275 <
1275 >    
1276      if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1277         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1278 <            pot, A, f, t, do_pot)
1278 >            pot(HB_POT), A, f, t, do_pot)
1279      endif
1280 <
1280 >    
1281      if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1282         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1283 <            pot, A, f, t, do_pot)
1283 >            pot(VDW_POT), A, f, t, do_pot)
1284      endif
1285      
1286      if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1287 < !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1288 < !           pot, A, f, t, do_pot)
1287 >       call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1288 >            pot(VDW_POT), A, f, t, do_pot)
1289      endif
1290 <
1290 >    
1291      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1292 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1293 <            do_pot)
1292 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1293 >            pot(METALLIC_POT), f, do_pot)
1294      endif
1295 <
1295 >    
1296      if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1297         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1298 <            pot, A, f, t, do_pot)
1298 >            pot(VDW_POT), A, f, t, do_pot)
1299      endif
1300 <
1300 >    
1301      if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1302         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1303 <            pot, A, f, t, do_pot)
1303 >            pot(VDW_POT), A, f, t, do_pot)
1304      endif
1305 +
1306 +    if ( iand(iHash, SC_PAIR).ne.0 ) then      
1307 +       call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1308 +            pot(METALLIC_POT), f, do_pot)
1309 +    endif
1310 +
1311      
1312 +    
1313    end subroutine do_pair
1314  
1315 <  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1315 >  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1316         do_pot, do_stress, eFrame, A, f, t, pot)
1317  
1318 <    real( kind = dp ) :: pot, sw
1318 >    real( kind = dp ) :: sw
1319 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1320      real( kind = dp ), dimension(9,nLocal) :: eFrame
1321      real (kind=dp), dimension(9,nLocal) :: A
1322      real (kind=dp), dimension(3,nLocal) :: f
# Line 1149 | Line 1324 | contains
1324  
1325      logical, intent(inout) :: do_pot, do_stress
1326      integer, intent(in) :: i, j
1327 <    real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1327 >    real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1328      real ( kind = dp )                :: r, rc
1329      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1330  
1331      integer :: me_i, me_j, iHash
1332  
1333 +    r = sqrt(rijsq)
1334 +
1335   #ifdef IS_MPI  
1336      me_i = atid_row(i)
1337      me_j = atid_col(j)  
# Line 1166 | Line 1343 | contains
1343      iHash = InteractionHash(me_i, me_j)
1344  
1345      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1346 <            call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1346 >            call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1347      endif
1348 +
1349 +    if ( iand(iHash, SC_PAIR).ne.0 ) then      
1350 +            call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1351 +    endif
1352      
1353    end subroutine do_prepair
1354  
1355  
1356    subroutine do_preforce(nlocal,pot)
1357      integer :: nlocal
1358 <    real( kind = dp ) :: pot
1358 >    real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1359  
1360      if (FF_uses_EAM .and. SIM_uses_EAM) then
1361 <       call calc_EAM_preforce_Frho(nlocal,pot)
1361 >       call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1362      endif
1363 +    if (FF_uses_SC .and. SIM_uses_SC) then
1364 +       call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1365 +    endif
1366  
1367  
1368    end subroutine do_preforce
# Line 1263 | Line 1447 | contains
1447      pot_Col = 0.0_dp
1448      pot_Temp = 0.0_dp
1449  
1266 –    rf_Row = 0.0_dp
1267 –    rf_Col = 0.0_dp
1268 –    rf_Temp = 0.0_dp
1269 –
1450   #endif
1451  
1452      if (FF_uses_EAM .and. SIM_uses_EAM) then
1453         call clean_EAM()
1454      endif
1455  
1276 –    rf = 0.0_dp
1456      tau_Temp = 0.0_dp
1457      virial_Temp = 0.0_dp
1458    end subroutine zero_work_arrays
# Line 1367 | Line 1546 | contains
1546  
1547    function FF_RequiresPrepairCalc() result(doesit)
1548      logical :: doesit
1549 <    doesit = FF_uses_EAM
1549 >    doesit = FF_uses_EAM .or. FF_uses_SC &
1550 >         .or. FF_uses_MEAM
1551    end function FF_RequiresPrepairCalc
1552  
1373 –  function FF_RequiresPostpairCalc() result(doesit)
1374 –    logical :: doesit
1375 –    doesit = FF_uses_RF
1376 –  end function FF_RequiresPostpairCalc
1377 –
1553   #ifdef PROFILE
1554    function getforcetime() result(totalforcetime)
1555      real(kind=dp) :: totalforcetime

Diff Legend

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