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 2284 by gezelter, Wed Sep 7 19:18:54 2005 UTC vs.
Revision 2512 by gezelter, Thu Dec 15 21:43:16 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.37 2005-09-07 19:18:54 gezelter Exp $, $Date: 2005-09-07 19:18:54 $, $Name: not supported by cvs2svn $, $Revision: 1.37 $
48 > !! @version $Id: doForces.F90,v 1.71 2005-12-15 21:43:16 gezelter Exp $, $Date: 2005-12-15 21:43:16 $, $Name: not supported by cvs2svn $, $Revision: 1.71 $
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          
# 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 256 | Line 281 | contains
281      logical :: i_is_GB
282      logical :: i_is_EAM
283      logical :: i_is_Shape
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  
265 –    stat = 0
294      if (.not. haveInteractionHash) then
295 <       call createInteractionHash(myStatus)      
268 <       if (myStatus .ne. 0) then
269 <          write(default_error, *) 'createInteractionHash failed in doForces!'
270 <          stat = -1
271 <          return
272 <       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 284 | 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
338 <          endif
339 <          if (i_is_Shape) then
340 <             thisRcut = getShapeCut(i)
341 <             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 <          
347 >                    
348            if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
349               biggestAtypeCutoff = atypeMaxCutoff(i)
350            endif
351 +
352         endif
353      enddo
322 –  
323 –    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 351 | 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)
437 <          endif
435 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
436 >             groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
437 >          endif          
438         enddo
439 <       if (nGroupTypes.eq.0) then
440 <          nGroupTypes = nGroupTypes + 1
441 <          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
442 <          groupToGtype(i) = nGroupTypes
439 >       if (nGroupTypesRow.eq.0) then
440 >          nGroupTypesRow = nGroupTypesRow + 1
441 >          gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
442 >          groupToGtypeRow(i) = nGroupTypesRow
443         else
444 <          do g = 1, nGroupTypes
445 <             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).gt.tol) then
446 <                nGroupTypes = nGroupTypes + 1
447 <                gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
448 <                groupToGtype(i) = nGroupTypes
368 <             else
369 <                groupToGtype(i) = g
444 >          GtypeFound = .false.
445 >          do g = 1, nGroupTypesRow
446 >             if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
447 >                groupToGtypeRow(i) = g
448 >                GtypeFound = .true.
449               endif
450            enddo
451 +          if (.not.GtypeFound) then            
452 +             nGroupTypesRow = nGroupTypesRow + 1
453 +             gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
454 +             groupToGtypeRow(i) = nGroupTypesRow
455 +          endif
456         endif
457 <    enddo
458 <    
457 >    enddo    
458 >
459 > #ifdef IS_MPI
460 >    do j = jstart, jend      
461 >       n_in_j = groupStartCol(j+1) - groupStartCol(j)
462 >       groupMaxCutoffCol(j) = 0.0_dp
463 >       do ja = groupStartCol(j), groupStartCol(j+1)-1
464 >          atom1 = groupListCol(ja)
465 >
466 >          me_j = atid_col(atom1)
467 >
468 >          if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
469 >             groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
470 >          endif          
471 >       enddo
472 >
473 >       if (nGroupTypesCol.eq.0) then
474 >          nGroupTypesCol = nGroupTypesCol + 1
475 >          gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
476 >          groupToGtypeCol(j) = nGroupTypesCol
477 >       else
478 >          GtypeFound = .false.
479 >          do g = 1, nGroupTypesCol
480 >             if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
481 >                groupToGtypeCol(j) = g
482 >                GtypeFound = .true.
483 >             endif
484 >          enddo
485 >          if (.not.GtypeFound) then            
486 >             nGroupTypesCol = nGroupTypesCol + 1
487 >             gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
488 >             groupToGtypeCol(j) = nGroupTypesCol
489 >          endif
490 >       endif
491 >    enddo    
492 >
493 > #else
494 > ! Set pointers to information we just found
495 >    nGroupTypesCol = nGroupTypesRow
496 >    groupToGtypeCol => groupToGtypeRow
497 >    gtypeMaxCutoffCol => gtypeMaxCutoffRow
498 >    groupMaxCutoffCol => groupMaxCutoffRow
499 > #endif
500 >
501      !! allocate the gtypeCutoffMap here.
502 <    allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
502 >    allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
503      !! then we do a double loop over all the group TYPES to find the cutoff
504      !! map between groups of two types
505 <    
506 <    do i = 1, nGroupTypes
507 <       do j = 1, nGroupTypes
505 >    tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
506 >
507 >    do i = 1, nGroupTypesRow      
508 >       do j = 1, nGroupTypesCol
509        
510            select case(cutoffPolicy)
511            case(TRADITIONAL_CUTOFF_POLICY)
512 <             thisRcut = maxval(gtypeMaxCutoff)
512 >             thisRcut = tradRcut
513            case(MIX_CUTOFF_POLICY)
514 <             thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
514 >             thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
515            case(MAX_CUTOFF_POLICY)
516 <             thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
516 >             thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
517            case default
518               call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
519               return
520            end select
521            gtypeCutoffMap(i,j)%rcut = thisRcut
522 +          
523 +          if (thisRcut.gt.largestRcut) largestRcut = thisRcut
524 +
525            gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
396 –          skin = defaultRlist - defaultRcut
397 –          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
526  
527 +          if (.not.haveSkinThickness) then
528 +             skinThickness = 1.0_dp
529 +          endif
530 +
531 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
532 +
533 +          ! sanity check
534 +
535 +          if (haveDefaultCutoffs) then
536 +             if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
537 +                call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
538 +             endif
539 +          endif
540         enddo
541      enddo
542 +
543 +    if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
544 +    if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
545 +    if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
546 + #ifdef IS_MPI
547 +    if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
548 +    if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
549 + #endif
550 +    groupMaxCutoffCol => null()
551 +    gtypeMaxCutoffCol => null()
552      
553      haveGtypeCutoffMap = .true.
554 <    
404 <  end subroutine createGtypeCutoffMap
405 <  
406 <  subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
407 <    real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
408 <    integer, intent(in) :: cutPolicy
409 <    
410 <    defaultRcut = defRcut
411 <    defaultRsw = defRsw
412 <    defaultRlist = defRlist
413 <    cutoffPolicy = cutPolicy
414 <  end subroutine setDefaultCutoffs
415 <  
416 <  subroutine setCutoffPolicy(cutPolicy)
554 >   end subroutine createGtypeCutoffMap
555  
556 +   subroutine setCutoffs(defRcut, defRsw)
557 +
558 +     real(kind=dp),intent(in) :: defRcut, defRsw
559 +     character(len = statusMsgSize) :: errMsg
560 +     integer :: localError
561 +
562 +     defaultRcut = defRcut
563 +     defaultRsw = defRsw
564 +    
565 +     defaultDoShift = .false.
566 +     if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
567 +        
568 +        write(errMsg, *) &
569 +             'cutoffRadius and switchingRadius are set to the same', newline &
570 +             // tab, 'value.  OOPSE will use shifted ', newline &
571 +             // tab, 'potentials instead of switching functions.'
572 +        
573 +        call handleInfo("setCutoffs", errMsg)
574 +        
575 +        defaultDoShift = .true.
576 +        
577 +     endif
578 +
579 +     localError = 0
580 +     call setLJDefaultCutoff( defaultRcut, defaultDoShift )
581 +     call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
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 +     haveGtypeCutoffMap = .false.
591 +   end subroutine setCutoffs
592 +
593 +   subroutine cWasLame()
594 +    
595 +     VisitCutoffsAfterComputing = .true.
596 +     return
597 +    
598 +   end subroutine cWasLame
599 +  
600 +   subroutine setCutoffPolicy(cutPolicy)
601 +    
602       integer, intent(in) :: cutPolicy
603 +    
604       cutoffPolicy = cutPolicy
605 <     call createGtypeCutoffMap()
606 <
422 <   end subroutine setCutoffPolicy
423 <    
605 >     haveCutoffPolicy = .true.
606 >     haveGtypeCutoffMap = .false.
607      
608 <  subroutine setSimVariables()
609 <    SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
610 <    SIM_uses_EAM = SimUsesEAM()
428 <    SIM_uses_RF = SimUsesRF()
429 <    SIM_requires_postpair_calc = SimRequiresPostpairCalc()
430 <    SIM_requires_prepair_calc = SimRequiresPrepairCalc()
431 <    SIM_uses_PBC = SimUsesPBC()
608 >   end subroutine setCutoffPolicy
609 >  
610 >   subroutine setElectrostaticMethod( thisESM )
611  
612 <    haveSIMvariables = .true.
612 >     integer, intent(in) :: thisESM
613  
614 <    return
615 <  end subroutine setSimVariables
614 >     electrostaticSummationMethod = thisESM
615 >     haveElectrostaticSummationMethod = .true.
616 >    
617 >   end subroutine setElectrostaticMethod
618  
619 +   subroutine setSkinThickness( thisSkin )
620 +    
621 +     real(kind=dp), intent(in) :: thisSkin
622 +    
623 +     skinThickness = thisSkin
624 +     haveSkinThickness = .true.    
625 +     haveGtypeCutoffMap = .false.
626 +    
627 +   end subroutine setSkinThickness
628 +      
629 +   subroutine setSimVariables()
630 +     SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
631 +     SIM_uses_EAM = SimUsesEAM()
632 +     SIM_requires_postpair_calc = SimRequiresPostpairCalc()
633 +     SIM_requires_prepair_calc = SimRequiresPrepairCalc()
634 +     SIM_uses_PBC = SimUsesPBC()
635 +    
636 +     haveSIMvariables = .true.
637 +    
638 +     return
639 +   end subroutine setSimVariables
640 +
641    subroutine doReadyCheck(error)
642      integer, intent(out) :: error
643  
# Line 443 | Line 646 | contains
646      error = 0
647  
648      if (.not. haveInteractionHash) then      
649 <       myStatus = 0      
447 <       call createInteractionHash(myStatus)      
448 <       if (myStatus .ne. 0) then
449 <          write(default_error, *) 'createInteractionHash failed in doForces!'
450 <          error = -1
451 <          return
452 <       endif
649 >       call createInteractionHash()      
650      endif
651  
652      if (.not. haveGtypeCutoffMap) then        
653 <       myStatus = 0      
457 <       call createGtypeCutoffMap(myStatus)      
458 <       if (myStatus .ne. 0) then
459 <          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
460 <          error = -1
461 <          return
462 <       endif
653 >       call createGtypeCutoffMap()      
654      endif
655  
656 +
657 +    if (VisitCutoffsAfterComputing) then
658 +       call set_switch(GROUP_SWITCH, largestRcut, largestRcut)      
659 +    endif
660 +
661 +
662      if (.not. haveSIMvariables) then
663         call setSimVariables()
664      endif
# Line 495 | Line 692 | contains
692    end subroutine doReadyCheck
693  
694  
695 <  subroutine init_FF(use_RF, use_UW, use_DW, thisStat)
695 >  subroutine init_FF(thisStat)
696  
500 –    logical, intent(in) :: use_RF
501 –    logical, intent(in) :: use_UW
502 –    logical, intent(in) :: use_DW
697      integer, intent(out) :: thisStat  
698      integer :: my_status, nMatches
505 –    integer :: corrMethod
699      integer, pointer :: MatchList(:) => null()
507 –    real(kind=dp) :: rcut, rrf, rt, dielect
700  
701      !! assume things are copacetic, unless they aren't
702      thisStat = 0
703  
512 –    !! Fortran's version of a cast:
513 –    FF_uses_RF = use_RF
514 –
515 –    !! set the electrostatic correction method
516 –    if (use_UW) then
517 –       corrMethod = 1
518 –    elseif (use_DW) then
519 –       corrMethod = 2
520 –    else
521 –       corrMethod = 0
522 –    endif
523 –    
704      !! init_FF is called *after* all of the atom types have been
705      !! defined in atype_module using the new_atype subroutine.
706      !!
# Line 549 | Line 729 | contains
729  
730  
731      haveSaneForceField = .true.
552 –
553 –    !! check to make sure the FF_uses_RF setting makes sense
554 –
555 –    if (FF_uses_Dipoles) then
556 –       if (FF_uses_RF) then
557 –          dielect = getDielect()
558 –          call initialize_rf(dielect)
559 –       endif
560 –    else
561 –       if (FF_uses_RF) then          
562 –          write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
563 –          thisStat = -1
564 –          haveSaneForceField = .false.
565 –          return
566 –       endif
567 –    endif
732  
733      if (FF_uses_EAM) then
734         call init_EAM_FF(my_status)
# Line 576 | Line 740 | contains
740         end if
741      endif
742  
579 –    if (FF_uses_GayBerne) then
580 –       call check_gb_pair_FF(my_status)
581 –       if (my_status .ne. 0) then
582 –          thisStat = -1
583 –          haveSaneForceField = .false.
584 –          return
585 –       endif
586 –    endif
587 –
743      if (.not. haveNeighborList) then
744         !! Create neighbor lists
745         call expandNeighborList(nLocal, my_status)
# Line 618 | Line 773 | contains
773  
774      !! Stress Tensor
775      real( kind = dp), dimension(9) :: tau  
776 <    real ( kind = dp ) :: pot
776 >    real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
777      logical ( kind = 2) :: do_pot_c, do_stress_c
778      logical :: do_pot
779      logical :: do_stress
780      logical :: in_switching_region
781   #ifdef IS_MPI
782 <    real( kind = DP ) :: pot_local
782 >    real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
783      integer :: nAtomsInRow
784      integer :: nAtomsInCol
785      integer :: nprocs
# Line 639 | Line 794 | contains
794      integer :: nlist
795      real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
796      real( kind = DP ) :: sw, dswdr, swderiv, mf
797 +    real( kind = DP ) :: rVal
798      real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
799      real(kind=dp) :: rfpot, mu_i, virial
800 +    real(kind=dp):: rCut
801      integer :: me_i, me_j, n_in_i, n_in_j
802      logical :: is_dp_i
803      integer :: neighborListSize
# Line 649 | Line 806 | contains
806      integer :: propPack_i, propPack_j
807      integer :: loopStart, loopEnd, loop
808      integer :: iHash
809 <    real(kind=dp) :: listSkin = 1.0  
809 >    integer :: i1
810 >  
811  
812      !! initialize local variables  
813  
# Line 713 | Line 871 | contains
871         ! (but only on the first time through):
872         if (loop .eq. loopStart) then
873   #ifdef IS_MPI
874 <          call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
874 >          call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
875                 update_nlist)
876   #else
877 <          call checkNeighborList(nGroups, q_group, listSkin, &
877 >          call checkNeighborList(nGroups, q_group, skinThickness, &
878                 update_nlist)
879   #endif
880         endif
# Line 774 | Line 932 | contains
932               me_j = atid(j)
933               call get_interatomic_vector(q_group(:,i), &
934                    q_group(:,j), d_grp, rgrpsq)
935 < #endif
935 > #endif      
936  
937 <             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
937 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
938                  if (update_nlist) then
939                     nlist = nlist + 1
940  
# Line 796 | Line 954 | contains
954  
955                     list(nlist) = j
956                  endif
957 +
958  
959 <                if (loop .eq. PAIR_LOOP) then
960 <                   vij = 0.0d0
802 <                   fij(1:3) = 0.0d0
803 <                endif
959 >                
960 >                if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
961  
962 <                call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
963 <                     in_switching_region)
964 <
965 <                n_in_j = groupStartCol(j+1) - groupStartCol(j)
966 <
967 <                do ia = groupStartRow(i), groupStartRow(i+1)-1
968 <
969 <                   atom1 = groupListRow(ia)
970 <
971 <                   inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
972 <
973 <                      atom2 = groupListCol(jb)
974 <
975 <                      if (skipThisPair(atom1, atom2)) cycle inner
976 <
977 <                      if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
978 <                         d_atm(1:3) = d_grp(1:3)
979 <                         ratmsq = rgrpsq
980 <                      else
981 < #ifdef IS_MPI
982 <                         call get_interatomic_vector(q_Row(:,atom1), &
983 <                              q_Col(:,atom2), d_atm, ratmsq)
984 < #else
985 <                         call get_interatomic_vector(q(:,atom1), &
986 <                              q(:,atom2), d_atm, ratmsq)
987 < #endif
988 <                      endif
989 <
833 <                      if (loop .eq. PREPAIR_LOOP) then
834 < #ifdef IS_MPI                      
835 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
836 <                              rgrpsq, d_grp, do_pot, do_stress, &
837 <                              eFrame, A, f, t, pot_local)
962 >                   rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
963 >                   if (loop .eq. PAIR_LOOP) then
964 >                      vij = 0.0d0
965 >                      fij(1:3) = 0.0d0
966 >                   endif
967 >                  
968 >                   call get_switch(rgrpsq, sw, dswdr, rgrp, &
969 >                        group_switch, in_switching_region)
970 >                  
971 >                   n_in_j = groupStartCol(j+1) - groupStartCol(j)
972 >                  
973 >                   do ia = groupStartRow(i), groupStartRow(i+1)-1
974 >                      
975 >                      atom1 = groupListRow(ia)
976 >                      
977 >                      inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
978 >                        
979 >                         atom2 = groupListCol(jb)
980 >                        
981 >                         if (skipThisPair(atom1, atom2))  cycle inner
982 >                        
983 >                         if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
984 >                            d_atm(1:3) = d_grp(1:3)
985 >                            ratmsq = rgrpsq
986 >                         else
987 > #ifdef IS_MPI
988 >                            call get_interatomic_vector(q_Row(:,atom1), &
989 >                                 q_Col(:,atom2), d_atm, ratmsq)
990   #else
991 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
992 <                              rgrpsq, d_grp, do_pot, do_stress, &
993 <                              eFrame, A, f, t, pot)
991 >                            call get_interatomic_vector(q(:,atom1), &
992 >                                 q(:,atom2), d_atm, ratmsq)
993 > #endif
994 >                         endif
995 >                        
996 >                         if (loop .eq. PREPAIR_LOOP) then
997 > #ifdef IS_MPI                      
998 >                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
999 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1000 >                                 eFrame, A, f, t, pot_local)
1001 > #else
1002 >                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1003 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1004 >                                 eFrame, A, f, t, pot)
1005   #endif                                              
1006 <                      else
1006 >                         else
1007   #ifdef IS_MPI                      
1008 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1009 <                              do_pot, &
1010 <                              eFrame, A, f, t, pot_local, vpair, fpair)
1008 >                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1009 >                                 do_pot, eFrame, A, f, t, pot_local, vpair, &
1010 >                                 fpair, d_grp, rgrp, rCut)
1011   #else
1012 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1013 <                              do_pot,  &
1014 <                              eFrame, A, f, t, pot, vpair, fpair)
1012 >                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1013 >                                 do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1014 >                                 d_grp, rgrp, rCut)
1015   #endif
1016 <
1017 <                         vij = vij + vpair
1018 <                         fij(1:3) = fij(1:3) + fpair(1:3)
1019 <                      endif
1020 <                   enddo inner
858 <                enddo
1016 >                            vij = vij + vpair
1017 >                            fij(1:3) = fij(1:3) + fpair(1:3)
1018 >                         endif
1019 >                      enddo inner
1020 >                   enddo
1021  
1022 <                if (loop .eq. PAIR_LOOP) then
1023 <                   if (in_switching_region) then
1024 <                      swderiv = vij*dswdr/rgrp
1025 <                      fij(1) = fij(1) + swderiv*d_grp(1)
1026 <                      fij(2) = fij(2) + swderiv*d_grp(2)
1027 <                      fij(3) = fij(3) + swderiv*d_grp(3)
1028 <
1029 <                      do ia=groupStartRow(i), groupStartRow(i+1)-1
1030 <                         atom1=groupListRow(ia)
1031 <                         mf = mfactRow(atom1)
1022 >                   if (loop .eq. PAIR_LOOP) then
1023 >                      if (in_switching_region) then
1024 >                         swderiv = vij*dswdr/rgrp
1025 >                         fij(1) = fij(1) + swderiv*d_grp(1)
1026 >                         fij(2) = fij(2) + swderiv*d_grp(2)
1027 >                         fij(3) = fij(3) + swderiv*d_grp(3)
1028 >                        
1029 >                         do ia=groupStartRow(i), groupStartRow(i+1)-1
1030 >                            atom1=groupListRow(ia)
1031 >                            mf = mfactRow(atom1)
1032   #ifdef IS_MPI
1033 <                         f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1034 <                         f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1035 <                         f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1033 >                            f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1034 >                            f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1035 >                            f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1036   #else
1037 <                         f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1038 <                         f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1039 <                         f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1037 >                            f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1038 >                            f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1039 >                            f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1040   #endif
1041 <                      enddo
1042 <
1043 <                      do jb=groupStartCol(j), groupStartCol(j+1)-1
1044 <                         atom2=groupListCol(jb)
1045 <                         mf = mfactCol(atom2)
1041 >                         enddo
1042 >                        
1043 >                         do jb=groupStartCol(j), groupStartCol(j+1)-1
1044 >                            atom2=groupListCol(jb)
1045 >                            mf = mfactCol(atom2)
1046   #ifdef IS_MPI
1047 <                         f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1048 <                         f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1049 <                         f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1047 >                            f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1048 >                            f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1049 >                            f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1050   #else
1051 <                         f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1052 <                         f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1053 <                         f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1051 >                            f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1052 >                            f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1053 >                            f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1054   #endif
1055 <                      enddo
1056 <                   endif
1055 >                         enddo
1056 >                      endif
1057  
1058 <                   if (do_stress) call add_stress_tensor(d_grp, fij)
1058 >                      if (do_stress) call add_stress_tensor(d_grp, fij)
1059 >                   endif
1060                  endif
1061 <             end if
1061 >             endif
1062            enddo
1063 +          
1064         enddo outer
1065  
1066         if (update_nlist) then
# Line 956 | Line 1120 | contains
1120  
1121      if (do_pot) then
1122         ! scatter/gather pot_row into the members of my column
1123 <       call scatter(pot_Row, pot_Temp, plan_atom_row)
1124 <
1123 >       do i = 1,LR_POT_TYPES
1124 >          call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1125 >       end do
1126         ! scatter/gather pot_local into all other procs
1127         ! add resultant to get total pot
1128         do i = 1, nlocal
1129 <          pot_local = pot_local + pot_Temp(i)
1129 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1130 >               + pot_Temp(1:LR_POT_TYPES,i)
1131         enddo
1132  
1133         pot_Temp = 0.0_DP
1134 <
1135 <       call scatter(pot_Col, pot_Temp, plan_atom_col)
1134 >       do i = 1,LR_POT_TYPES
1135 >          call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1136 >       end do
1137         do i = 1, nlocal
1138 <          pot_local = pot_local + pot_Temp(i)
1138 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1139 >               + pot_Temp(1:LR_POT_TYPES,i)
1140         enddo
1141  
1142      endif
1143   #endif
1144  
1145 <    if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1145 >    if (SIM_requires_postpair_calc) then
1146 >       do i = 1, nlocal            
1147 >          
1148 >          ! we loop only over the local atoms, so we don't need row and column
1149 >          ! lookups for the types
1150 >          
1151 >          me_i = atid(i)
1152 >          
1153 >          ! is the atom electrostatic?  See if it would have an
1154 >          ! electrostatic interaction with itself
1155 >          iHash = InteractionHash(me_i,me_i)
1156  
1157 <       if (FF_uses_RF .and. SIM_uses_RF) then
980 <
1157 >          if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1158   #ifdef IS_MPI
1159 <          call scatter(rf_Row,rf,plan_atom_row_3d)
1160 <          call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
984 <          do i = 1,nlocal
985 <             rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
986 <          end do
987 < #endif
988 <
989 <          do i = 1, nLocal
990 <
991 <             rfpot = 0.0_DP
992 < #ifdef IS_MPI
993 <             me_i = atid_row(i)
1159 >             call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1160 >                  t, do_pot)
1161   #else
1162 <             me_i = atid(i)
1162 >             call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1163 >                  t, do_pot)
1164   #endif
1165 <             iHash = InteractionHash(me_i,me_j)
1165 >          endif
1166 >  
1167 >          
1168 >          if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1169              
1170 <             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1171 <
1172 <                mu_i = getDipoleMoment(me_i)
1173 <
1174 <                !! The reaction field needs to include a self contribution
1175 <                !! to the field:
1176 <                call accumulate_self_rf(i, mu_i, eFrame)
1177 <                !! Get the reaction field contribution to the
1178 <                !! potential and torques:
1179 <                call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1170 >             ! loop over the excludes to accumulate RF stuff we've
1171 >             ! left out of the normal pair loop
1172 >            
1173 >             do i1 = 1, nSkipsForAtom(i)
1174 >                j = skipsForAtom(i, i1)
1175 >                
1176 >                ! prevent overcounting of the skips
1177 >                if (i.lt.j) then
1178 >                   call get_interatomic_vector(q(:,i), &
1179 >                        q(:,j), d_atm, ratmsq)
1180 >                   rVal = dsqrt(ratmsq)
1181 >                   call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1182 >                        in_switching_region)
1183   #ifdef IS_MPI
1184 <                pot_local = pot_local + rfpot
1184 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1185 >                        vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1186   #else
1187 <                pot = pot + rfpot
1188 <
1187 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1188 >                        vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1189   #endif
1190 <             endif
1191 <          enddo
1192 <       endif
1190 >                endif
1191 >             enddo
1192 >          endif
1193 >       enddo
1194      endif
1195 <
1020 <
1195 >    
1196   #ifdef IS_MPI
1197 <
1197 >    
1198      if (do_pot) then
1199 <       pot = pot + pot_local
1200 <       !! we assume the c code will do the allreduce to get the total potential
1026 <       !! we could do it right here if we needed to...
1199 >       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1200 >            mpi_comm_world,mpi_err)            
1201      endif
1202 <
1202 >    
1203      if (do_stress) then
1204         call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1205              mpi_comm_world,mpi_err)
1206         call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1207              mpi_comm_world,mpi_err)
1208      endif
1209 <
1209 >    
1210   #else
1211 <
1211 >    
1212      if (do_stress) then
1213         tau = tau_Temp
1214         virial = virial_Temp
1215      endif
1216 <
1216 >    
1217   #endif
1218 <
1218 >    
1219    end subroutine do_force_loop
1220  
1221    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1222 <       eFrame, A, f, t, pot, vpair, fpair)
1222 >       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1223  
1224 <    real( kind = dp ) :: pot, vpair, sw
1224 >    real( kind = dp ) :: vpair, sw
1225 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1226      real( kind = dp ), dimension(3) :: fpair
1227      real( kind = dp ), dimension(nLocal)   :: mfact
1228      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 1058 | Line 1233 | contains
1233      logical, intent(inout) :: do_pot
1234      integer, intent(in) :: i, j
1235      real ( kind = dp ), intent(inout) :: rijsq
1236 <    real ( kind = dp )                :: r
1236 >    real ( kind = dp ), intent(inout) :: r_grp
1237      real ( kind = dp ), intent(inout) :: d(3)
1238 <    real ( kind = dp ) :: ebalance
1238 >    real ( kind = dp ), intent(inout) :: d_grp(3)
1239 >    real ( kind = dp ), intent(inout) :: rCut
1240 >    real ( kind = dp ) :: r
1241      integer :: me_i, me_j
1242  
1243      integer :: iHash
# Line 1078 | Line 1255 | contains
1255   #endif
1256  
1257      iHash = InteractionHash(me_i, me_j)
1258 <
1258 >    
1259      if ( iand(iHash, LJ_PAIR).ne.0 ) then
1260 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1260 >       call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1261 >            pot(VDW_POT), f, do_pot)
1262      endif
1263 <
1263 >    
1264      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1265 <       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1266 <            pot, eFrame, f, t, do_pot, corrMethod)
1089 <
1090 <       if (FF_uses_RF .and. SIM_uses_RF) then
1091 <
1092 <          ! CHECK ME (RF needs to know about all electrostatic types)
1093 <          call accumulate_rf(i, j, r, eFrame, sw)
1094 <          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1095 <       endif
1096 <
1265 >       call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1266 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1267      endif
1268 <
1268 >    
1269      if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1270         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1271 <            pot, A, f, t, do_pot)
1271 >            pot(HB_POT), A, f, t, do_pot)
1272      endif
1273 <
1273 >    
1274      if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1275         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1276 <            pot, A, f, t, do_pot)
1276 >            pot(HB_POT), A, f, t, do_pot)
1277      endif
1278 <
1278 >    
1279      if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1280         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1281 <            pot, A, f, t, do_pot)
1281 >            pot(VDW_POT), A, f, t, do_pot)
1282      endif
1283      
1284      if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1285 < !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1286 < !           pot, A, f, t, do_pot)
1285 >       call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1286 >            pot(VDW_POT), A, f, t, do_pot)
1287      endif
1288 <
1288 >    
1289      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1290 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1291 <            do_pot)
1290 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1291 >            pot(METALLIC_POT), f, do_pot)
1292      endif
1293 <
1293 >    
1294      if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1295         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1296 <            pot, A, f, t, do_pot)
1296 >            pot(VDW_POT), A, f, t, do_pot)
1297      endif
1298 <
1298 >    
1299      if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1300         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1301 <            pot, A, f, t, do_pot)
1301 >            pot(VDW_POT), A, f, t, do_pot)
1302      endif
1303 +
1304 +    if ( iand(iHash, SC_PAIR).ne.0 ) then      
1305 +       call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1306 +            pot(METALLIC_POT), f, do_pot)
1307 +    endif
1308 +
1309      
1310 +    
1311    end subroutine do_pair
1312  
1313 <  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1313 >  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1314         do_pot, do_stress, eFrame, A, f, t, pot)
1315  
1316 <    real( kind = dp ) :: pot, sw
1316 >    real( kind = dp ) :: sw
1317 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1318      real( kind = dp ), dimension(9,nLocal) :: eFrame
1319      real (kind=dp), dimension(9,nLocal) :: A
1320      real (kind=dp), dimension(3,nLocal) :: f
# Line 1144 | Line 1322 | contains
1322  
1323      logical, intent(inout) :: do_pot, do_stress
1324      integer, intent(in) :: i, j
1325 <    real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1325 >    real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1326      real ( kind = dp )                :: r, rc
1327      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1328  
1329      integer :: me_i, me_j, iHash
1330  
1331 +    r = sqrt(rijsq)
1332 +
1333   #ifdef IS_MPI  
1334      me_i = atid_row(i)
1335      me_j = atid_col(j)  
# Line 1161 | Line 1341 | contains
1341      iHash = InteractionHash(me_i, me_j)
1342  
1343      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1344 <            call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1344 >            call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1345      endif
1346 +
1347 +    if ( iand(iHash, SC_PAIR).ne.0 ) then      
1348 +            call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1349 +    endif
1350      
1351    end subroutine do_prepair
1352  
1353  
1354    subroutine do_preforce(nlocal,pot)
1355      integer :: nlocal
1356 <    real( kind = dp ) :: pot
1356 >    real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1357  
1358      if (FF_uses_EAM .and. SIM_uses_EAM) then
1359 <       call calc_EAM_preforce_Frho(nlocal,pot)
1359 >       call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1360      endif
1361 +    if (FF_uses_SC .and. SIM_uses_SC) then
1362 +       call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1363 +    endif
1364  
1365  
1366    end subroutine do_preforce
# Line 1258 | Line 1445 | contains
1445      pot_Col = 0.0_dp
1446      pot_Temp = 0.0_dp
1447  
1261 –    rf_Row = 0.0_dp
1262 –    rf_Col = 0.0_dp
1263 –    rf_Temp = 0.0_dp
1264 –
1448   #endif
1449  
1450      if (FF_uses_EAM .and. SIM_uses_EAM) then
1451         call clean_EAM()
1452      endif
1453  
1271 –    rf = 0.0_dp
1454      tau_Temp = 0.0_dp
1455      virial_Temp = 0.0_dp
1456    end subroutine zero_work_arrays
# Line 1362 | Line 1544 | contains
1544  
1545    function FF_RequiresPrepairCalc() result(doesit)
1546      logical :: doesit
1547 <    doesit = FF_uses_EAM
1547 >    doesit = FF_uses_EAM .or. FF_uses_SC &
1548 >         .or. FF_uses_MEAM
1549    end function FF_RequiresPrepairCalc
1550  
1368 –  function FF_RequiresPostpairCalc() result(doesit)
1369 –    logical :: doesit
1370 –    doesit = FF_uses_RF
1371 –  end function FF_RequiresPostpairCalc
1372 –
1551   #ifdef PROFILE
1552    function getforcetime() result(totalforcetime)
1553      real(kind=dp) :: totalforcetime

Diff Legend

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