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 2280 by gezelter, Thu Sep 1 20:17:55 2005 UTC vs.
Revision 2461 by gezelter, Mon Nov 21 22:59:02 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.34 2005-09-01 20:17:55 gezelter Exp $, $Date: 2005-09-01 20:17:55 $, $Name: not supported by cvs2svn $, $Revision: 1.34 $
48 > !! @version $Id: doForces.F90,v 1.69 2005-11-21 22:58:35 gezelter Exp $, $Date: 2005-11-21 22:58:35 $, $Name: not supported by cvs2svn $, $Revision: 1.69 $
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          
305 >       if (SimHasAtype(i)) then    
306            call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
307            call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
308            call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
# 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 <          if (i_is_LJ) then
315 <             thisRcut = getSigma(i) * 2.5_dp
316 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
317 <          endif
318 <          if (i_is_Elect) then
319 <             thisRcut = defaultRcut
320 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
321 <          endif
322 <          if (i_is_Sticky) then
323 <             thisRcut = getStickyCut(i)
324 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
325 <          endif
326 <          if (i_is_StickyP) then
327 <             thisRcut = getStickyPowerCut(i)
328 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
329 <          endif
330 <          if (i_is_GB) then
331 <             thisRcut = getGayBerneCut(i)
332 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
333 <          endif
334 <          if (i_is_EAM) then
335 <             thisRcut = getEAMCut(i)
336 <             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
312 <             thisRcut = getShapeCut(i)
313 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
314 <          endif
315 <          
347 >                    
348            if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
349               biggestAtypeCutoff = atypeMaxCutoff(i)
350            endif
351 +
352         endif
353      enddo
321 –  
322 –    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 <
364 >    
365      !! allocate the groupToGtype and gtypeMaxCutoff here.
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      
333 –    !! first we do a single loop over the cutoff groups to find the largest cutoff for any atypes
334 –    !! present in this group.   We also create gtypes at this point.
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 344 | 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 >
440 >       if (nGroupTypesRow.eq.0) then
441 >          nGroupTypesRow = nGroupTypesRow + 1
442 >          gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
443 >          groupToGtypeRow(i) = nGroupTypesRow
444         else
445 <          do g = 1, nGroupTypes
446 <             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).gt.tol) then
447 <                nGroupTypes = nGroupTypes + 1
448 <                gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
449 <                groupToGtype(i) = nGroupTypes
361 <             else
362 <                groupToGtype(i) = g
445 >          GtypeFound = .false.
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 +             nGroupTypesRow = nGroupTypesRow + 1
454 +             gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
455 +             groupToGtypeRow(i) = nGroupTypesRow
456 +          endif
457         endif
458 <    enddo
367 <    
368 <    !! allocate the gtypeCutoffMap here.
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(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)
514 <             case(MIX_CUTOFF_POLICY)
515 <                thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
516 <             case(MAX_CUTOFF_POLICY)
517 <                thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
518 <             case default
519 <                call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
520 <                return
521 <          end select      
522 <         gtypeCutoffMap(i,j)%rcut = thisRcut
523 <         gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
524 <         skin = defaultRlist - defaultRcut
525 <         gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
512 >          case(TRADITIONAL_CUTOFF_POLICY)
513 >             thisRcut = tradRcut
514 >          case(MIX_CUTOFF_POLICY)
515 >             thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
516 >          case(MAX_CUTOFF_POLICY)
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
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     end subroutine createGtypeCutoffMap
556  
557 <   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
398 <     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
399 <     integer, intent(in) :: cutPolicy
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 <     defaultRlist = defRlist
566 <     cutoffPolicy = cutPolicy
567 <   end subroutine setDefaultCutoffs
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 <   subroutine setCutoffPolicy(cutPolicy)
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  
606 <   end subroutine setCutoffPolicy
414 <    
606 >     call createGtypeCutoffMap()
607      
608 <  subroutine setSimVariables()
609 <    SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
610 <    SIM_uses_EAM = SimUsesEAM()
419 <    SIM_uses_RF = SimUsesRF()
420 <    SIM_requires_postpair_calc = SimRequiresPostpairCalc()
421 <    SIM_requires_prepair_calc = SimRequiresPrepairCalc()
422 <    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 +    
626 +     call createGtypeCutoffMap()
627 +    
628 +   end subroutine setSkinThickness
629 +      
630 +   subroutine setSimVariables()
631 +     SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
632 +     SIM_uses_EAM = SimUsesEAM()
633 +     SIM_requires_postpair_calc = SimRequiresPostpairCalc()
634 +     SIM_requires_prepair_calc = SimRequiresPrepairCalc()
635 +     SIM_uses_PBC = SimUsesPBC()
636 +    
637 +     haveSIMvariables = .true.
638 +    
639 +     return
640 +   end subroutine setSimVariables
641 +
642    subroutine doReadyCheck(error)
643      integer, intent(out) :: error
644  
# Line 434 | Line 647 | contains
647      error = 0
648  
649      if (.not. haveInteractionHash) then      
650 <       myStatus = 0      
438 <       call createInteractionHash(myStatus)      
439 <       if (myStatus .ne. 0) then
440 <          write(default_error, *) 'createInteractionHash failed in doForces!'
441 <          error = -1
442 <          return
443 <       endif
650 >       call createInteractionHash()      
651      endif
652  
653      if (.not. haveGtypeCutoffMap) then        
654 <       myStatus = 0      
448 <       call createGtypeCutoffMap(myStatus)      
449 <       if (myStatus .ne. 0) then
450 <          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
451 <          error = -1
452 <          return
453 <       endif
654 >       call createGtypeCutoffMap()      
655      endif
656  
657 +
658 +    if (VisitCutoffsAfterComputing) then
659 +       call set_switch(GROUP_SWITCH, largestRcut, largestRcut)      
660 +    endif
661 +
662 +
663      if (.not. haveSIMvariables) then
664         call setSimVariables()
665      endif
666  
667 <    if (.not. haveRlist) then
668 <       write(default_error, *) 'rList has not been set in doForces!'
669 <       error = -1
670 <       return
671 <    endif
667 >  !  if (.not. haveRlist) then
668 >  !     write(default_error, *) 'rList has not been set in doForces!'
669 >  !     error = -1
670 >  !     return
671 >  !  endif
672  
673      if (.not. haveNeighborList) then
674         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 486 | Line 693 | contains
693    end subroutine doReadyCheck
694  
695  
696 <  subroutine init_FF(use_RF, use_UW, use_DW, thisStat)
696 >  subroutine init_FF(thisStat)
697  
491 –    logical, intent(in) :: use_RF
492 –    logical, intent(in) :: use_UW
493 –    logical, intent(in) :: use_DW
698      integer, intent(out) :: thisStat  
699      integer :: my_status, nMatches
496 –    integer :: corrMethod
700      integer, pointer :: MatchList(:) => null()
498 –    real(kind=dp) :: rcut, rrf, rt, dielect
701  
702      !! assume things are copacetic, unless they aren't
703      thisStat = 0
704  
503 –    !! Fortran's version of a cast:
504 –    FF_uses_RF = use_RF
505 –
506 –    !! set the electrostatic correction method
507 –    if (use_UW) then
508 –       corrMethod = 1
509 –    elseif (use_DW) then
510 –       corrMethod = 2
511 –    else
512 –       corrMethod = 0
513 –    endif
514 –    
705      !! init_FF is called *after* all of the atom types have been
706      !! defined in atype_module using the new_atype subroutine.
707      !!
# Line 541 | Line 731 | contains
731  
732      haveSaneForceField = .true.
733  
544 –    !! check to make sure the FF_uses_RF setting makes sense
545 –
546 –    if (FF_uses_Dipoles) then
547 –       if (FF_uses_RF) then
548 –          dielect = getDielect()
549 –          call initialize_rf(dielect)
550 –       endif
551 –    else
552 –       if (FF_uses_RF) then          
553 –          write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
554 –          thisStat = -1
555 –          haveSaneForceField = .false.
556 –          return
557 –       endif
558 –    endif
559 –
734      if (FF_uses_EAM) then
735         call init_EAM_FF(my_status)
736         if (my_status /= 0) then
# Line 567 | Line 741 | contains
741         end if
742      endif
743  
570 –    if (FF_uses_GayBerne) then
571 –       call check_gb_pair_FF(my_status)
572 –       if (my_status .ne. 0) then
573 –          thisStat = -1
574 –          haveSaneForceField = .false.
575 –          return
576 –       endif
577 –    endif
578 –
744      if (.not. haveNeighborList) then
745         !! Create neighbor lists
746         call expandNeighborList(nLocal, my_status)
# Line 609 | Line 774 | contains
774  
775      !! Stress Tensor
776      real( kind = dp), dimension(9) :: tau  
777 <    real ( kind = dp ) :: pot
777 >    real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
778      logical ( kind = 2) :: do_pot_c, do_stress_c
779      logical :: do_pot
780      logical :: do_stress
781      logical :: in_switching_region
782   #ifdef IS_MPI
783 <    real( kind = DP ) :: pot_local
783 >    real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
784      integer :: nAtomsInRow
785      integer :: nAtomsInCol
786      integer :: nprocs
# Line 630 | Line 795 | contains
795      integer :: nlist
796      real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
797      real( kind = DP ) :: sw, dswdr, swderiv, mf
798 +    real( kind = DP ) :: rVal
799      real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
800      real(kind=dp) :: rfpot, mu_i, virial
801 +    real(kind=dp):: rCut
802      integer :: me_i, me_j, n_in_i, n_in_j
803      logical :: is_dp_i
804      integer :: neighborListSize
# Line 640 | Line 807 | contains
807      integer :: propPack_i, propPack_j
808      integer :: loopStart, loopEnd, loop
809      integer :: iHash
810 <    real(kind=dp) :: listSkin = 1.0  
810 >    integer :: i1
811 >  
812  
813      !! initialize local variables  
814  
# Line 704 | Line 872 | contains
872         ! (but only on the first time through):
873         if (loop .eq. loopStart) then
874   #ifdef IS_MPI
875 <          call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
875 >          call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
876                 update_nlist)
877   #else
878 <          call checkNeighborList(nGroups, q_group, listSkin, &
878 >          call checkNeighborList(nGroups, q_group, skinThickness, &
879                 update_nlist)
880   #endif
881         endif
# Line 765 | Line 933 | contains
933               me_j = atid(j)
934               call get_interatomic_vector(q_group(:,i), &
935                    q_group(:,j), d_grp, rgrpsq)
936 < #endif
936 > #endif      
937  
938 <             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
938 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
939                  if (update_nlist) then
940                     nlist = nlist + 1
941  
# Line 787 | Line 955 | contains
955  
956                     list(nlist) = j
957                  endif
958 +
959  
960 <                if (loop .eq. PAIR_LOOP) then
961 <                   vij = 0.0d0
793 <                   fij(1:3) = 0.0d0
794 <                endif
960 >                
961 >                if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
962  
963 <                call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
964 <                     in_switching_region)
965 <
966 <                n_in_j = groupStartCol(j+1) - groupStartCol(j)
967 <
968 <                do ia = groupStartRow(i), groupStartRow(i+1)-1
969 <
970 <                   atom1 = groupListRow(ia)
971 <
972 <                   inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
973 <
974 <                      atom2 = groupListCol(jb)
975 <
976 <                      if (skipThisPair(atom1, atom2)) cycle inner
977 <
978 <                      if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
979 <                         d_atm(1:3) = d_grp(1:3)
980 <                         ratmsq = rgrpsq
981 <                      else
982 < #ifdef IS_MPI
983 <                         call get_interatomic_vector(q_Row(:,atom1), &
984 <                              q_Col(:,atom2), d_atm, ratmsq)
963 >                   rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
964 >                   if (loop .eq. PAIR_LOOP) then
965 >                      vij = 0.0d0
966 >                      fij(1:3) = 0.0d0
967 >                   endif
968 >                  
969 >                   call get_switch(rgrpsq, sw, dswdr, rgrp, &
970 >                        group_switch, in_switching_region)
971 >                  
972 >                   n_in_j = groupStartCol(j+1) - groupStartCol(j)
973 >                  
974 >                   do ia = groupStartRow(i), groupStartRow(i+1)-1
975 >                      
976 >                      atom1 = groupListRow(ia)
977 >                      
978 >                      inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
979 >                        
980 >                         atom2 = groupListCol(jb)
981 >                        
982 >                         if (skipThisPair(atom1, atom2))  cycle inner
983 >                        
984 >                         if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
985 >                            d_atm(1:3) = d_grp(1:3)
986 >                            ratmsq = rgrpsq
987 >                         else
988 > #ifdef IS_MPI
989 >                            call get_interatomic_vector(q_Row(:,atom1), &
990 >                                 q_Col(:,atom2), d_atm, ratmsq)
991   #else
992 <                         call get_interatomic_vector(q(:,atom1), &
993 <                              q(:,atom2), d_atm, ratmsq)
992 >                            call get_interatomic_vector(q(:,atom1), &
993 >                                 q(:,atom2), d_atm, ratmsq)
994   #endif
995 <                      endif
996 <
997 <                      if (loop .eq. PREPAIR_LOOP) then
995 >                         endif
996 >                        
997 >                         if (loop .eq. PREPAIR_LOOP) then
998   #ifdef IS_MPI                      
999 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1000 <                              rgrpsq, d_grp, do_pot, do_stress, &
1001 <                              eFrame, A, f, t, pot_local)
999 >                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1000 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1001 >                                 eFrame, A, f, t, pot_local)
1002   #else
1003 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1004 <                              rgrpsq, d_grp, do_pot, do_stress, &
1005 <                              eFrame, A, f, t, pot)
1003 >                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1004 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1005 >                                 eFrame, A, f, t, pot)
1006   #endif                                              
1007 <                      else
1007 >                         else
1008   #ifdef IS_MPI                      
1009 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1010 <                              do_pot, &
1011 <                              eFrame, A, f, t, pot_local, vpair, fpair)
1009 >                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1010 >                                 do_pot, eFrame, A, f, t, pot_local, vpair, &
1011 >                                 fpair, d_grp, rgrp, rCut)
1012   #else
1013 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1014 <                              do_pot,  &
1015 <                              eFrame, A, f, t, pot, vpair, fpair)
1013 >                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1014 >                                 do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1015 >                                 d_grp, rgrp, rCut)
1016   #endif
1017 +                            vij = vij + vpair
1018 +                            fij(1:3) = fij(1:3) + fpair(1:3)
1019 +                         endif
1020 +                      enddo inner
1021 +                   enddo
1022  
1023 <                         vij = vij + vpair
1024 <                         fij(1:3) = fij(1:3) + fpair(1:3)
1025 <                      endif
1026 <                   enddo inner
1027 <                enddo
1028 <
1029 <                if (loop .eq. PAIR_LOOP) then
1030 <                   if (in_switching_region) then
1031 <                      swderiv = vij*dswdr/rgrp
1032 <                      fij(1) = fij(1) + swderiv*d_grp(1)
855 <                      fij(2) = fij(2) + swderiv*d_grp(2)
856 <                      fij(3) = fij(3) + swderiv*d_grp(3)
857 <
858 <                      do ia=groupStartRow(i), groupStartRow(i+1)-1
859 <                         atom1=groupListRow(ia)
860 <                         mf = mfactRow(atom1)
1023 >                   if (loop .eq. PAIR_LOOP) then
1024 >                      if (in_switching_region) then
1025 >                         swderiv = vij*dswdr/rgrp
1026 >                         fij(1) = fij(1) + swderiv*d_grp(1)
1027 >                         fij(2) = fij(2) + swderiv*d_grp(2)
1028 >                         fij(3) = fij(3) + swderiv*d_grp(3)
1029 >                        
1030 >                         do ia=groupStartRow(i), groupStartRow(i+1)-1
1031 >                            atom1=groupListRow(ia)
1032 >                            mf = mfactRow(atom1)
1033   #ifdef IS_MPI
1034 <                         f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1035 <                         f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1036 <                         f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1034 >                            f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1035 >                            f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1036 >                            f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1037   #else
1038 <                         f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1039 <                         f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1040 <                         f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1038 >                            f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1039 >                            f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1040 >                            f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1041   #endif
1042 <                      enddo
1043 <
1044 <                      do jb=groupStartCol(j), groupStartCol(j+1)-1
1045 <                         atom2=groupListCol(jb)
1046 <                         mf = mfactCol(atom2)
1042 >                         enddo
1043 >                        
1044 >                         do jb=groupStartCol(j), groupStartCol(j+1)-1
1045 >                            atom2=groupListCol(jb)
1046 >                            mf = mfactCol(atom2)
1047   #ifdef IS_MPI
1048 <                         f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1049 <                         f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1050 <                         f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1048 >                            f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1049 >                            f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1050 >                            f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1051   #else
1052 <                         f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1053 <                         f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1054 <                         f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1052 >                            f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1053 >                            f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1054 >                            f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1055   #endif
1056 <                      enddo
1057 <                   endif
1056 >                         enddo
1057 >                      endif
1058  
1059 <                   if (do_stress) call add_stress_tensor(d_grp, fij)
1059 >                      if (do_stress) call add_stress_tensor(d_grp, fij)
1060 >                   endif
1061                  endif
1062 <             end if
1062 >             endif
1063            enddo
1064 +          
1065         enddo outer
1066  
1067         if (update_nlist) then
# Line 947 | Line 1121 | contains
1121  
1122      if (do_pot) then
1123         ! scatter/gather pot_row into the members of my column
1124 <       call scatter(pot_Row, pot_Temp, plan_atom_row)
1125 <
1124 >       do i = 1,LR_POT_TYPES
1125 >          call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1126 >       end do
1127         ! scatter/gather pot_local into all other procs
1128         ! add resultant to get total pot
1129         do i = 1, nlocal
1130 <          pot_local = pot_local + pot_Temp(i)
1130 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1131 >               + pot_Temp(1:LR_POT_TYPES,i)
1132         enddo
1133  
1134         pot_Temp = 0.0_DP
1135 <
1136 <       call scatter(pot_Col, pot_Temp, plan_atom_col)
1135 >       do i = 1,LR_POT_TYPES
1136 >          call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1137 >       end do
1138         do i = 1, nlocal
1139 <          pot_local = pot_local + pot_Temp(i)
1139 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1140 >               + pot_Temp(1:LR_POT_TYPES,i)
1141         enddo
1142  
1143      endif
1144   #endif
1145  
1146 <    if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1146 >    if (SIM_requires_postpair_calc) then
1147 >       do i = 1, nlocal            
1148 >          
1149 >          ! we loop only over the local atoms, so we don't need row and column
1150 >          ! lookups for the types
1151 >          
1152 >          me_i = atid(i)
1153 >          
1154 >          ! is the atom electrostatic?  See if it would have an
1155 >          ! electrostatic interaction with itself
1156 >          iHash = InteractionHash(me_i,me_i)
1157  
1158 <       if (FF_uses_RF .and. SIM_uses_RF) then
971 <
1158 >          if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1159   #ifdef IS_MPI
1160 <          call scatter(rf_Row,rf,plan_atom_row_3d)
1161 <          call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
975 <          do i = 1,nlocal
976 <             rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
977 <          end do
978 < #endif
979 <
980 <          do i = 1, nLocal
981 <
982 <             rfpot = 0.0_DP
983 < #ifdef IS_MPI
984 <             me_i = atid_row(i)
1160 >             call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1161 >                  t, do_pot)
1162   #else
1163 <             me_i = atid(i)
1163 >             call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1164 >                  t, do_pot)
1165   #endif
1166 <             iHash = InteractionHash(me_i,me_j)
1166 >          endif
1167 >  
1168 >          
1169 >          if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1170              
1171 <             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1172 <
1173 <                mu_i = getDipoleMoment(me_i)
1174 <
1175 <                !! The reaction field needs to include a self contribution
1176 <                !! to the field:
1177 <                call accumulate_self_rf(i, mu_i, eFrame)
1178 <                !! Get the reaction field contribution to the
1179 <                !! potential and torques:
1180 <                call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1171 >             ! loop over the excludes to accumulate RF stuff we've
1172 >             ! left out of the normal pair loop
1173 >            
1174 >             do i1 = 1, nSkipsForAtom(i)
1175 >                j = skipsForAtom(i, i1)
1176 >                
1177 >                ! prevent overcounting of the skips
1178 >                if (i.lt.j) then
1179 >                   call get_interatomic_vector(q(:,i), &
1180 >                        q(:,j), d_atm, ratmsq)
1181 >                   rVal = dsqrt(ratmsq)
1182 >                   call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1183 >                        in_switching_region)
1184   #ifdef IS_MPI
1185 <                pot_local = pot_local + rfpot
1185 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1186 >                        vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1187   #else
1188 <                pot = pot + rfpot
1189 <
1188 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1189 >                        vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1190   #endif
1191 <             endif
1192 <          enddo
1193 <       endif
1191 >                endif
1192 >             enddo
1193 >          endif
1194 >       enddo
1195      endif
1196 <
1011 <
1196 >    
1197   #ifdef IS_MPI
1198 <
1198 >    
1199      if (do_pot) then
1200 <       pot = pot + pot_local
1201 <       !! we assume the c code will do the allreduce to get the total potential
1017 <       !! we could do it right here if we needed to...
1200 >       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1201 >            mpi_comm_world,mpi_err)            
1202      endif
1203 <
1203 >    
1204      if (do_stress) then
1205         call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1206              mpi_comm_world,mpi_err)
1207         call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1208              mpi_comm_world,mpi_err)
1209      endif
1210 <
1210 >    
1211   #else
1212 <
1212 >    
1213      if (do_stress) then
1214         tau = tau_Temp
1215         virial = virial_Temp
1216      endif
1217 <
1217 >    
1218   #endif
1219 <
1219 >    
1220    end subroutine do_force_loop
1221  
1222    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1223 <       eFrame, A, f, t, pot, vpair, fpair)
1223 >       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1224  
1225 <    real( kind = dp ) :: pot, vpair, sw
1225 >    real( kind = dp ) :: vpair, sw
1226 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1227      real( kind = dp ), dimension(3) :: fpair
1228      real( kind = dp ), dimension(nLocal)   :: mfact
1229      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 1049 | Line 1234 | contains
1234      logical, intent(inout) :: do_pot
1235      integer, intent(in) :: i, j
1236      real ( kind = dp ), intent(inout) :: rijsq
1237 <    real ( kind = dp )                :: r
1237 >    real ( kind = dp ), intent(inout) :: r_grp
1238      real ( kind = dp ), intent(inout) :: d(3)
1239 <    real ( kind = dp ) :: ebalance
1239 >    real ( kind = dp ), intent(inout) :: d_grp(3)
1240 >    real ( kind = dp ), intent(inout) :: rCut
1241 >    real ( kind = dp ) :: r
1242      integer :: me_i, me_j
1243  
1244      integer :: iHash
# Line 1069 | Line 1256 | contains
1256   #endif
1257  
1258      iHash = InteractionHash(me_i, me_j)
1259 <
1259 >    
1260      if ( iand(iHash, LJ_PAIR).ne.0 ) then
1261 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1262 <    endif
1076 <
1077 <    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1078 <       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1079 <            pot, eFrame, f, t, do_pot, corrMethod)
1080 <
1081 <       if (FF_uses_RF .and. SIM_uses_RF) then
1082 <
1083 <          ! CHECK ME (RF needs to know about all electrostatic types)
1084 <          call accumulate_rf(i, j, r, eFrame, sw)
1085 <          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1086 <       endif
1087 <
1261 >       call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1262 >            pot(VDW_POT), f, do_pot)
1263      endif
1264 <
1264 >    
1265 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1266 >       call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1267 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1268 >    endif
1269 >    
1270      if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1271         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1272 <            pot, A, f, t, do_pot)
1272 >            pot(HB_POT), A, f, t, do_pot)
1273      endif
1274 <
1274 >    
1275      if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1276         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1277 <            pot, A, f, t, do_pot)
1277 >            pot(HB_POT), A, f, t, do_pot)
1278      endif
1279 <
1279 >    
1280      if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1281         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1282 <            pot, A, f, t, do_pot)
1282 >            pot(VDW_POT), A, f, t, do_pot)
1283      endif
1284      
1285      if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1286 < !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1287 < !           pot, A, f, t, do_pot)
1286 >       call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1287 >            pot(VDW_POT), A, f, t, do_pot)
1288      endif
1289 <
1289 >    
1290      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1291 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1292 <            do_pot)
1291 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1292 >            pot(METALLIC_POT), f, do_pot)
1293      endif
1294 <
1294 >    
1295      if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1296         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1297 <            pot, A, f, t, do_pot)
1297 >            pot(VDW_POT), A, f, t, do_pot)
1298      endif
1299 <
1299 >    
1300      if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1301         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1302 <            pot, A, f, t, do_pot)
1302 >            pot(VDW_POT), A, f, t, do_pot)
1303      endif
1304 +
1305 +    if ( iand(iHash, SC_PAIR).ne.0 ) then      
1306 +       call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1307 +            pot(METALLIC_POT), f, do_pot)
1308 +    endif
1309 +
1310      
1311 +    
1312    end subroutine do_pair
1313  
1314 <  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1314 >  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1315         do_pot, do_stress, eFrame, A, f, t, pot)
1316  
1317 <    real( kind = dp ) :: pot, sw
1317 >    real( kind = dp ) :: sw
1318 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1319      real( kind = dp ), dimension(9,nLocal) :: eFrame
1320      real (kind=dp), dimension(9,nLocal) :: A
1321      real (kind=dp), dimension(3,nLocal) :: f
# Line 1135 | Line 1323 | contains
1323  
1324      logical, intent(inout) :: do_pot, do_stress
1325      integer, intent(in) :: i, j
1326 <    real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1326 >    real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1327      real ( kind = dp )                :: r, rc
1328      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1329  
1330      integer :: me_i, me_j, iHash
1331  
1332 +    r = sqrt(rijsq)
1333 +
1334   #ifdef IS_MPI  
1335      me_i = atid_row(i)
1336      me_j = atid_col(j)  
# Line 1152 | Line 1342 | contains
1342      iHash = InteractionHash(me_i, me_j)
1343  
1344      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1345 <            call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1345 >            call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1346      endif
1347 +
1348 +    if ( iand(iHash, SC_PAIR).ne.0 ) then      
1349 +            call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1350 +    endif
1351      
1352    end subroutine do_prepair
1353  
1354  
1355    subroutine do_preforce(nlocal,pot)
1356      integer :: nlocal
1357 <    real( kind = dp ) :: pot
1357 >    real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1358  
1359      if (FF_uses_EAM .and. SIM_uses_EAM) then
1360 <       call calc_EAM_preforce_Frho(nlocal,pot)
1360 >       call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1361      endif
1362 +    if (FF_uses_SC .and. SIM_uses_SC) then
1363 +       call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1364 +    endif
1365  
1366  
1367    end subroutine do_preforce
# Line 1249 | Line 1446 | contains
1446      pot_Col = 0.0_dp
1447      pot_Temp = 0.0_dp
1448  
1252 –    rf_Row = 0.0_dp
1253 –    rf_Col = 0.0_dp
1254 –    rf_Temp = 0.0_dp
1255 –
1449   #endif
1450  
1451      if (FF_uses_EAM .and. SIM_uses_EAM) then
1452         call clean_EAM()
1453      endif
1454  
1262 –    rf = 0.0_dp
1455      tau_Temp = 0.0_dp
1456      virial_Temp = 0.0_dp
1457    end subroutine zero_work_arrays
# Line 1353 | Line 1545 | contains
1545  
1546    function FF_RequiresPrepairCalc() result(doesit)
1547      logical :: doesit
1548 <    doesit = FF_uses_EAM
1548 >    doesit = FF_uses_EAM .or. FF_uses_SC &
1549 >         .or. FF_uses_MEAM
1550    end function FF_RequiresPrepairCalc
1551  
1359 –  function FF_RequiresPostpairCalc() result(doesit)
1360 –    logical :: doesit
1361 –    doesit = FF_uses_RF
1362 –  end function FF_RequiresPostpairCalc
1363 –
1552   #ifdef PROFILE
1553    function getforcetime() result(totalforcetime)
1554      real(kind=dp) :: totalforcetime

Diff Legend

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