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 2270 by gezelter, Tue Aug 9 22:33:37 2005 UTC vs.
Revision 2540 by chuckv, Mon Jan 9 22:22:35 2006 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.28 2005-08-09 22:33:37 gezelter Exp $, $Date: 2005-08-09 22:33:37 $, $Name: not supported by cvs2svn $, $Revision: 1.28 $
48 > !! @version $Id: doForces.F90,v 1.75 2006-01-09 22:22:35 chuckv Exp $, $Date: 2006-01-09 22:22:35 $, $Name: not supported by cvs2svn $, $Revision: 1.75 $
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 73 | Line 73 | module doForces
73  
74   #define __FORTRAN90
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
82    INTEGER, PARAMETER:: PAIR_LOOP    = 2
83  
# Line 83 | Line 86 | module doForces
86    logical, save :: haveSaneForceField = .false.
87    logical, save :: haveInteractionHash = .false.
88    logical, save :: haveGtypeCutoffMap = .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 :: 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 :: setCutoffs
120 +  public :: cWasLame
121 +  public :: setElectrostaticMethod
122 +  public :: setCutoffPolicy
123 +  public :: setSkinThickness
124    public :: do_force_loop
102 –  public :: createInteractionHash
103 –  public :: createGtypeCutoffMap
125  
126   #ifdef PROFILE
127    public :: getforcetime
# Line 113 | 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
148       real(kind=dp) :: rlistsq
149    end type gtypeCutoffs
150    type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
151 <  
151 >
152   contains
153  
154 <  subroutine createInteractionHash(status)
154 >  subroutine createInteractionHash()
155      integer :: nAtypes
130 –    integer, intent(out) :: status
156      integer :: i
157      integer :: j
158      integer :: iHash
# Line 139 | 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 146 | Line 173 | contains
173      logical :: j_is_GB
174      logical :: j_is_EAM
175      logical :: j_is_Shape
176 <    
177 <    status = 0  
176 >    logical :: j_is_SC
177 >    logical :: j_is_MEAM
178 >    real(kind=dp) :: myRcut
179  
180      if (.not. associated(atypes)) then
181 <       call handleError("atype", "atypes was not present before call of createInteractionHash!")
154 <       status = -1
181 >       call handleError("doForces", "atypes was not present before call of createInteractionHash!")
182         return
183      endif
184      
185      nAtypes = getSize(atypes)
186      
187      if (nAtypes == 0) then
188 <       status = -1
188 >       call handleError("doForces", "nAtypes was zero during call of createInteractionHash!")
189         return
190      end if
191  
192      if (.not. allocated(InteractionHash)) then
193         allocate(InteractionHash(nAtypes,nAtypes))
194 +    else
195 +       deallocate(InteractionHash)
196 +       allocate(InteractionHash(nAtypes,nAtypes))
197      endif
198  
199      if (.not. allocated(atypeMaxCutoff)) then
200         allocate(atypeMaxCutoff(nAtypes))
201 +    else
202 +       deallocate(atypeMaxCutoff)
203 +       allocate(atypeMaxCutoff(nAtypes))
204      endif
205          
206      do i = 1, nAtypes
# Line 178 | 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 <
215 <       if (i_is_LJ) then
183 <          thisCut = getDefaultLJCutoff(i)
184 <          if (thisCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisCut
185 <       endif
186 <
214 >       call getElementProperty(atypes, i, "is_SC", i_is_SC)
215 >       call getElementProperty(atypes, i, "is_MEAM", i_is_MEAM)
216  
188 –
217         do j = i, nAtypes
218  
219            iHash = 0
# Line 198 | 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 219 | 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 238 | Line 272 | contains
272      haveInteractionHash = .true.
273    end subroutine createInteractionHash
274  
275 <  subroutine createGtypeCutoffMap(defaultRcut, defaultSkinThickness, stat)
275 >  subroutine createGtypeCutoffMap()
276  
277 <    real(kind=dp), intent(in), optional :: defaultRCut, defaultSkinThickness
278 <    integer, intent(out) :: stat
277 >    logical :: i_is_LJ
278 >    logical :: i_is_Elect
279 >    logical :: i_is_Sticky
280 >    logical :: i_is_StickyP
281 >    logical :: i_is_GB
282 >    logical :: i_is_EAM
283 >    logical :: i_is_Shape
284 >    logical :: i_is_SC
285 >    logical :: GtypeFound
286  
287 <    integer :: myStatus, nAtypes
287 >    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
288 >    integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
289 >    integer :: nGroupsInRow
290 >    integer :: nGroupsInCol
291 >    integer :: nGroupTypesRow,nGroupTypesCol
292 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol
293 >    real(kind=dp) :: biggestAtypeCutoff
294  
248 –    stat = 0
295      if (.not. haveInteractionHash) then
296 <       call createInteractionHash(myStatus)      
251 <       if (myStatus .ne. 0) then
252 <          write(default_error, *) 'createInteractionHash failed in doForces!'
253 <          stat = -1
254 <          return
255 <       endif
296 >       call createInteractionHash()      
297      endif
298 <
298 > #ifdef IS_MPI
299 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
300 >    nGroupsInCol = getNgroupsInCol(plan_group_col)
301 > #endif
302      nAtypes = getSize(atypes)
303 <
303 > ! Set all of the initial cutoffs to zero.
304 >    atypeMaxCutoff = 0.0_dp
305      do i = 1, nAtypes
306 <      
307 <       atypeMaxCutoff(i) =
306 >       if (SimHasAtype(i)) then    
307 >          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
308 >          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
309 >          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
310 >          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
311 >          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
312 >          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
313 >          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
314 >          call getElementProperty(atypes, i, "is_SC", i_is_SC)
315 >
316 >          if (haveDefaultCutoffs) then
317 >             atypeMaxCutoff(i) = defaultRcut
318 >          else
319 >             if (i_is_LJ) then          
320 >                thisRcut = getSigma(i) * 2.5_dp
321 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
322 >             endif
323 >             if (i_is_Elect) then
324 >                thisRcut = defaultRcut
325 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
326 >             endif
327 >             if (i_is_Sticky) then
328 >                thisRcut = getStickyCut(i)
329 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
330 >             endif
331 >             if (i_is_StickyP) then
332 >                thisRcut = getStickyPowerCut(i)
333 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
334 >             endif
335 >             if (i_is_GB) then
336 >                thisRcut = getGayBerneCut(i)
337 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
338 >             endif
339 >             if (i_is_EAM) then
340 >                thisRcut = getEAMCut(i)
341 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
342 >             endif
343 >             if (i_is_Shape) then
344 >                thisRcut = getShapeCut(i)
345 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
346 >             endif
347 >             if (i_is_SC) then
348 >                thisRcut = getSCCut(i)
349 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
350 >             endif
351 >          endif
352 >                    
353 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
354 >             biggestAtypeCutoff = atypeMaxCutoff(i)
355 >          endif
356  
357 +       endif
358 +    enddo
359      
360 +    istart = 1
361 +    jstart = 1
362 + #ifdef IS_MPI
363 +    iend = nGroupsInRow
364 +    jend = nGroupsInCol
365 + #else
366 +    iend = nGroups
367 +    jend = nGroups
368 + #endif
369 +    
370 +    !! allocate the groupToGtype and gtypeMaxCutoff here.
371 +    if(.not.allocated(groupToGtypeRow)) then
372 +     !  allocate(groupToGtype(iend))
373 +       allocate(groupToGtypeRow(iend))
374 +    else
375 +       deallocate(groupToGtypeRow)
376 +       allocate(groupToGtypeRow(iend))
377 +    endif
378 +    if(.not.allocated(groupMaxCutoffRow)) then
379 +       allocate(groupMaxCutoffRow(iend))
380 +    else
381 +       deallocate(groupMaxCutoffRow)
382 +       allocate(groupMaxCutoffRow(iend))
383 +    end if
384  
385 +    if(.not.allocated(gtypeMaxCutoffRow)) then
386 +       allocate(gtypeMaxCutoffRow(iend))
387 +    else
388 +       deallocate(gtypeMaxCutoffRow)
389 +       allocate(gtypeMaxCutoffRow(iend))
390 +    endif
391  
392  
393 <     haveGtypeCutoffMap = .true.
394 <   end subroutine createGtypeCutoffMap
393 > #ifdef IS_MPI
394 >       ! We only allocate new storage if we are in MPI because Ncol /= Nrow
395 >    if(.not.associated(groupToGtypeCol)) then
396 >       allocate(groupToGtypeCol(jend))
397 >    else
398 >       deallocate(groupToGtypeCol)
399 >       allocate(groupToGtypeCol(jend))
400 >    end if
401  
402 <  subroutine setSimVariables()
403 <    SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
404 <    SIM_uses_EAM = SimUsesEAM()
405 <    SIM_uses_RF = SimUsesRF()
406 <    SIM_requires_postpair_calc = SimRequiresPostpairCalc()
407 <    SIM_requires_prepair_calc = SimRequiresPrepairCalc()
408 <    SIM_uses_PBC = SimUsesPBC()
402 >    if(.not.associated(groupMaxCutoffCol)) then
403 >       allocate(groupMaxCutoffCol(jend))
404 >    else
405 >       deallocate(groupMaxCutoffCol)
406 >       allocate(groupMaxCutoffCol(jend))
407 >    end if
408 >    if(.not.associated(gtypeMaxCutoffCol)) then
409 >       allocate(gtypeMaxCutoffCol(jend))
410 >    else
411 >       deallocate(gtypeMaxCutoffCol)      
412 >       allocate(gtypeMaxCutoffCol(jend))
413 >    end if
414  
415 <    haveSIMvariables = .true.
415 >       groupMaxCutoffCol = 0.0_dp
416 >       gtypeMaxCutoffCol = 0.0_dp
417  
418 <    return
419 <  end subroutine setSimVariables
418 > #endif
419 >       groupMaxCutoffRow = 0.0_dp
420 >       gtypeMaxCutoffRow = 0.0_dp
421  
422 +
423 +    !! first we do a single loop over the cutoff groups to find the
424 +    !! largest cutoff for any atypes present in this group.  We also
425 +    !! create gtypes at this point.
426 +    
427 +    tol = 1.0d-6
428 +    nGroupTypesRow = 0
429 +    nGroupTypesCol = 0
430 +    do i = istart, iend      
431 +       n_in_i = groupStartRow(i+1) - groupStartRow(i)
432 +       groupMaxCutoffRow(i) = 0.0_dp
433 +       do ia = groupStartRow(i), groupStartRow(i+1)-1
434 +          atom1 = groupListRow(ia)
435 + #ifdef IS_MPI
436 +          me_i = atid_row(atom1)
437 + #else
438 +          me_i = atid(atom1)
439 + #endif          
440 +          if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
441 +             groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
442 +          endif          
443 +       enddo
444 +       if (nGroupTypesRow.eq.0) then
445 +          nGroupTypesRow = nGroupTypesRow + 1
446 +          gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
447 +          groupToGtypeRow(i) = nGroupTypesRow
448 +       else
449 +          GtypeFound = .false.
450 +          do g = 1, nGroupTypesRow
451 +             if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
452 +                groupToGtypeRow(i) = g
453 +                GtypeFound = .true.
454 +             endif
455 +          enddo
456 +          if (.not.GtypeFound) then            
457 +             nGroupTypesRow = nGroupTypesRow + 1
458 +             gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
459 +             groupToGtypeRow(i) = nGroupTypesRow
460 +          endif
461 +       endif
462 +    enddo    
463 +
464 + #ifdef IS_MPI
465 +    do j = jstart, jend      
466 +       n_in_j = groupStartCol(j+1) - groupStartCol(j)
467 +       groupMaxCutoffCol(j) = 0.0_dp
468 +       do ja = groupStartCol(j), groupStartCol(j+1)-1
469 +          atom1 = groupListCol(ja)
470 +
471 +          me_j = atid_col(atom1)
472 +
473 +          if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
474 +             groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
475 +          endif          
476 +       enddo
477 +
478 +       if (nGroupTypesCol.eq.0) then
479 +          nGroupTypesCol = nGroupTypesCol + 1
480 +          gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
481 +          groupToGtypeCol(j) = nGroupTypesCol
482 +       else
483 +          GtypeFound = .false.
484 +          do g = 1, nGroupTypesCol
485 +             if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
486 +                groupToGtypeCol(j) = g
487 +                GtypeFound = .true.
488 +             endif
489 +          enddo
490 +          if (.not.GtypeFound) then            
491 +             nGroupTypesCol = nGroupTypesCol + 1
492 +             gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
493 +             groupToGtypeCol(j) = nGroupTypesCol
494 +          endif
495 +       endif
496 +    enddo    
497 +
498 + #else
499 + ! Set pointers to information we just found
500 +    nGroupTypesCol = nGroupTypesRow
501 +    groupToGtypeCol => groupToGtypeRow
502 +    gtypeMaxCutoffCol => gtypeMaxCutoffRow
503 +    groupMaxCutoffCol => groupMaxCutoffRow
504 + #endif
505 +
506 +    !! allocate the gtypeCutoffMap here.
507 +    allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
508 +    !! then we do a double loop over all the group TYPES to find the cutoff
509 +    !! map between groups of two types
510 +    tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
511 +
512 +    do i = 1, nGroupTypesRow      
513 +       do j = 1, nGroupTypesCol
514 +      
515 +          select case(cutoffPolicy)
516 +          case(TRADITIONAL_CUTOFF_POLICY)
517 +             thisRcut = tradRcut
518 +          case(MIX_CUTOFF_POLICY)
519 +             thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
520 +          case(MAX_CUTOFF_POLICY)
521 +             thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
522 +          case default
523 +             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
524 +             return
525 +          end select
526 +          gtypeCutoffMap(i,j)%rcut = thisRcut
527 +          
528 +          if (thisRcut.gt.largestRcut) largestRcut = thisRcut
529 +
530 +          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
531 +
532 +          if (.not.haveSkinThickness) then
533 +             skinThickness = 1.0_dp
534 +          endif
535 +
536 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
537 +
538 +          ! sanity check
539 +
540 +          if (haveDefaultCutoffs) then
541 +             if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
542 +                call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
543 +             endif
544 +          endif
545 +       enddo
546 +    enddo
547 +
548 +    if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
549 +    if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
550 +    if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
551 + #ifdef IS_MPI
552 +    if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
553 +    if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
554 + #endif
555 +    groupMaxCutoffCol => null()
556 +    gtypeMaxCutoffCol => null()
557 +    
558 +    haveGtypeCutoffMap = .true.
559 +   end subroutine createGtypeCutoffMap
560 +
561 +   subroutine setCutoffs(defRcut, defRsw)
562 +
563 +     real(kind=dp),intent(in) :: defRcut, defRsw
564 +     character(len = statusMsgSize) :: errMsg
565 +     integer :: localError
566 +
567 +     defaultRcut = defRcut
568 +     defaultRsw = defRsw
569 +    
570 +     defaultDoShift = .false.
571 +     if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
572 +        
573 +        write(errMsg, *) &
574 +             'cutoffRadius and switchingRadius are set to the same', newline &
575 +             // tab, 'value.  OOPSE will use shifted ', newline &
576 +             // tab, 'potentials instead of switching functions.'
577 +        
578 +        call handleInfo("setCutoffs", errMsg)
579 +        
580 +        defaultDoShift = .true.
581 +        
582 +     endif
583 +
584 +     localError = 0
585 +     call setLJDefaultCutoff( defaultRcut, defaultDoShift )
586 +     call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
587 +     call setCutoffEAM( defaultRcut, localError)
588 +     if (localError /= 0) then
589 +       write(errMsg, *) 'An error has occured in setting the EAM cutoff'
590 +       call handleError("setCutoffs", errMsg)
591 +     end if
592 +     call set_switch(GROUP_SWITCH, defaultRsw, defaultRcut)
593 +
594 +     haveDefaultCutoffs = .true.
595 +     haveGtypeCutoffMap = .false.
596 +   end subroutine setCutoffs
597 +
598 +   subroutine cWasLame()
599 +    
600 +     VisitCutoffsAfterComputing = .true.
601 +     return
602 +    
603 +   end subroutine cWasLame
604 +  
605 +   subroutine setCutoffPolicy(cutPolicy)
606 +    
607 +     integer, intent(in) :: cutPolicy
608 +    
609 +     cutoffPolicy = cutPolicy
610 +     haveCutoffPolicy = .true.
611 +     haveGtypeCutoffMap = .false.
612 +    
613 +   end subroutine setCutoffPolicy
614 +  
615 +   subroutine setElectrostaticMethod( thisESM )
616 +
617 +     integer, intent(in) :: thisESM
618 +
619 +     electrostaticSummationMethod = thisESM
620 +     haveElectrostaticSummationMethod = .true.
621 +    
622 +   end subroutine setElectrostaticMethod
623 +
624 +   subroutine setSkinThickness( thisSkin )
625 +    
626 +     real(kind=dp), intent(in) :: thisSkin
627 +    
628 +     skinThickness = thisSkin
629 +     haveSkinThickness = .true.    
630 +     haveGtypeCutoffMap = .false.
631 +    
632 +   end subroutine setSkinThickness
633 +      
634 +   subroutine setSimVariables()
635 +     SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
636 +     SIM_uses_EAM = SimUsesEAM()
637 +     SIM_requires_postpair_calc = SimRequiresPostpairCalc()
638 +     SIM_requires_prepair_calc = SimRequiresPrepairCalc()
639 +     SIM_uses_PBC = SimUsesPBC()
640 +     SIM_uses_SC = SimUsesSC()
641 +    
642 +     haveSIMvariables = .true.
643 +    
644 +     return
645 +   end subroutine setSimVariables
646 +
647    subroutine doReadyCheck(error)
648      integer, intent(out) :: error
649  
# Line 289 | Line 652 | contains
652      error = 0
653  
654      if (.not. haveInteractionHash) then      
655 <       myStatus = 0      
293 <       call createInteractionHash(myStatus)      
294 <       if (myStatus .ne. 0) then
295 <          write(default_error, *) 'createInteractionHash failed in doForces!'
296 <          error = -1
297 <          return
298 <       endif
655 >       call createInteractionHash()      
656      endif
657  
658      if (.not. haveGtypeCutoffMap) then        
659 <       myStatus = 0      
303 <       call createGtypeCutoffMap(myStatus)      
304 <       if (myStatus .ne. 0) then
305 <          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
306 <          error = -1
307 <          return
308 <       endif
659 >       call createGtypeCutoffMap()      
660      endif
661  
662 +
663 +    if (VisitCutoffsAfterComputing) then
664 +       call set_switch(GROUP_SWITCH, largestRcut, largestRcut)      
665 +    endif
666 +
667 +
668      if (.not. haveSIMvariables) then
669         call setSimVariables()
670      endif
671  
672 <    if (.not. haveRlist) then
673 <       write(default_error, *) 'rList has not been set in doForces!'
674 <       error = -1
675 <       return
676 <    endif
672 >  !  if (.not. haveRlist) then
673 >  !     write(default_error, *) 'rList has not been set in doForces!'
674 >  !     error = -1
675 >  !     return
676 >  !  endif
677  
678      if (.not. haveNeighborList) then
679         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 341 | Line 698 | contains
698    end subroutine doReadyCheck
699  
700  
701 <  subroutine init_FF(use_RF_c, thisStat)
701 >  subroutine init_FF(thisStat)
702  
346 –    logical, intent(in) :: use_RF_c
347 –
703      integer, intent(out) :: thisStat  
704      integer :: my_status, nMatches
705      integer, pointer :: MatchList(:) => null()
351 –    real(kind=dp) :: rcut, rrf, rt, dielect
706  
707      !! assume things are copacetic, unless they aren't
708      thisStat = 0
709  
356 –    !! Fortran's version of a cast:
357 –    FF_uses_RF = use_RF_c
358 –
710      !! init_FF is called *after* all of the atom types have been
711      !! defined in atype_module using the new_atype subroutine.
712      !!
# Line 366 | Line 717 | contains
717      FF_uses_Dipoles = .false.
718      FF_uses_GayBerne = .false.
719      FF_uses_EAM = .false.
720 +    FF_uses_SC = .false.
721  
722      call getMatchingElementList(atypes, "is_Directional", .true., &
723           nMatches, MatchList)
# Line 382 | Line 734 | contains
734      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
735      if (nMatches .gt. 0) FF_uses_EAM = .true.
736  
737 +    call getMatchingElementList(atypes, "is_SC", .true., nMatches, MatchList)
738 +    if (nMatches .gt. 0) FF_uses_SC = .true.
739  
386 –    haveSaneForceField = .true.
740  
741 <    !! check to make sure the FF_uses_RF setting makes sense
389 <
390 <    if (FF_uses_Dipoles) then
391 <       if (FF_uses_RF) then
392 <          dielect = getDielect()
393 <          call initialize_rf(dielect)
394 <       endif
395 <    else
396 <       if (FF_uses_RF) then          
397 <          write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
398 <          thisStat = -1
399 <          haveSaneForceField = .false.
400 <          return
401 <       endif
402 <    endif
741 >    haveSaneForceField = .true.
742  
743      if (FF_uses_EAM) then
744         call init_EAM_FF(my_status)
# Line 411 | Line 750 | contains
750         end if
751      endif
752  
414 –    if (FF_uses_GayBerne) then
415 –       call check_gb_pair_FF(my_status)
416 –       if (my_status .ne. 0) then
417 –          thisStat = -1
418 –          haveSaneForceField = .false.
419 –          return
420 –       endif
421 –    endif
422 –
753      if (.not. haveNeighborList) then
754         !! Create neighbor lists
755         call expandNeighborList(nLocal, my_status)
# Line 453 | Line 783 | contains
783  
784      !! Stress Tensor
785      real( kind = dp), dimension(9) :: tau  
786 <    real ( kind = dp ) :: pot
786 >    real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
787      logical ( kind = 2) :: do_pot_c, do_stress_c
788      logical :: do_pot
789      logical :: do_stress
790      logical :: in_switching_region
791   #ifdef IS_MPI
792 <    real( kind = DP ) :: pot_local
792 >    real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
793      integer :: nAtomsInRow
794      integer :: nAtomsInCol
795      integer :: nprocs
# Line 474 | Line 804 | contains
804      integer :: nlist
805      real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
806      real( kind = DP ) :: sw, dswdr, swderiv, mf
807 +    real( kind = DP ) :: rVal
808      real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
809      real(kind=dp) :: rfpot, mu_i, virial
810 +    real(kind=dp):: rCut
811      integer :: me_i, me_j, n_in_i, n_in_j
812      logical :: is_dp_i
813      integer :: neighborListSize
# Line 484 | Line 816 | contains
816      integer :: propPack_i, propPack_j
817      integer :: loopStart, loopEnd, loop
818      integer :: iHash
819 <    real(kind=dp) :: listSkin = 1.0  
819 >    integer :: i1
820 >  
821  
822      !! initialize local variables  
823  
# Line 548 | Line 881 | contains
881         ! (but only on the first time through):
882         if (loop .eq. loopStart) then
883   #ifdef IS_MPI
884 <          call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
884 >          call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
885                 update_nlist)
886   #else
887 <          call checkNeighborList(nGroups, q_group, listSkin, &
887 >          call checkNeighborList(nGroups, q_group, skinThickness, &
888                 update_nlist)
889   #endif
890         endif
# Line 575 | Line 908 | contains
908   #endif
909         outer: do i = istart, iend
910  
578 – #ifdef IS_MPI
579 –             me_i = atid_row(i)
580 – #else
581 –             me_i = atid(i)
582 – #endif
583 –
911            if (update_nlist) point(i) = nlist + 1
912  
913            n_in_i = groupStartRow(i+1) - groupStartRow(i)
# Line 615 | Line 942 | contains
942               me_j = atid(j)
943               call get_interatomic_vector(q_group(:,i), &
944                    q_group(:,j), d_grp, rgrpsq)
945 < #endif
945 > #endif      
946  
947 <             if (rgrpsq < InteractionHash(me_i,me_j)%rListsq) then
947 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
948                  if (update_nlist) then
949                     nlist = nlist + 1
950  
# Line 637 | Line 964 | contains
964  
965                     list(nlist) = j
966                  endif
967 +
968  
969 <                if (loop .eq. PAIR_LOOP) then
970 <                   vij = 0.0d0
643 <                   fij(1:3) = 0.0d0
644 <                endif
969 >                
970 >                if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
971  
972 <                call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
973 <                     in_switching_region)
974 <
975 <                n_in_j = groupStartCol(j+1) - groupStartCol(j)
976 <
977 <                do ia = groupStartRow(i), groupStartRow(i+1)-1
978 <
979 <                   atom1 = groupListRow(ia)
980 <
981 <                   inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
982 <
983 <                      atom2 = groupListCol(jb)
984 <
985 <                      if (skipThisPair(atom1, atom2)) cycle inner
986 <
987 <                      if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
988 <                         d_atm(1:3) = d_grp(1:3)
989 <                         ratmsq = rgrpsq
990 <                      else
991 < #ifdef IS_MPI
992 <                         call get_interatomic_vector(q_Row(:,atom1), &
993 <                              q_Col(:,atom2), d_atm, ratmsq)
972 >                   rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
973 >                   if (loop .eq. PAIR_LOOP) then
974 >                      vij = 0.0d0
975 >                      fij(1:3) = 0.0d0
976 >                   endif
977 >                  
978 >                   call get_switch(rgrpsq, sw, dswdr, rgrp, &
979 >                        group_switch, in_switching_region)
980 >                  
981 >                   n_in_j = groupStartCol(j+1) - groupStartCol(j)
982 >                  
983 >                   do ia = groupStartRow(i), groupStartRow(i+1)-1
984 >                      
985 >                      atom1 = groupListRow(ia)
986 >                      
987 >                      inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
988 >                        
989 >                         atom2 = groupListCol(jb)
990 >                        
991 >                         if (skipThisPair(atom1, atom2))  cycle inner
992 >                        
993 >                         if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
994 >                            d_atm(1:3) = d_grp(1:3)
995 >                            ratmsq = rgrpsq
996 >                         else
997 > #ifdef IS_MPI
998 >                            call get_interatomic_vector(q_Row(:,atom1), &
999 >                                 q_Col(:,atom2), d_atm, ratmsq)
1000   #else
1001 <                         call get_interatomic_vector(q(:,atom1), &
1002 <                              q(:,atom2), d_atm, ratmsq)
1001 >                            call get_interatomic_vector(q(:,atom1), &
1002 >                                 q(:,atom2), d_atm, ratmsq)
1003   #endif
1004 <                      endif
1005 <
1006 <                      if (loop .eq. PREPAIR_LOOP) then
1004 >                         endif
1005 >                        
1006 >                         if (loop .eq. PREPAIR_LOOP) then
1007   #ifdef IS_MPI                      
1008 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1009 <                              rgrpsq, d_grp, do_pot, do_stress, &
1010 <                              eFrame, A, f, t, pot_local)
1008 >                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1009 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1010 >                                 eFrame, A, f, t, pot_local)
1011   #else
1012 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1013 <                              rgrpsq, d_grp, do_pot, do_stress, &
1014 <                              eFrame, A, f, t, pot)
1012 >                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1013 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1014 >                                 eFrame, A, f, t, pot)
1015   #endif                                              
1016 <                      else
1016 >                         else
1017   #ifdef IS_MPI                      
1018 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1019 <                              do_pot, &
1020 <                              eFrame, A, f, t, pot_local, vpair, fpair)
1018 >                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1019 >                                 do_pot, eFrame, A, f, t, pot_local, vpair, &
1020 >                                 fpair, d_grp, rgrp, rCut)
1021   #else
1022 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1023 <                              do_pot,  &
1024 <                              eFrame, A, f, t, pot, vpair, fpair)
1025 < #endif
1022 >                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1023 >                                 do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1024 >                                 d_grp, rgrp, rCut)
1025 > #endif
1026 >                            vij = vij + vpair
1027 >                            fij(1:3) = fij(1:3) + fpair(1:3)
1028 >                         endif
1029 >                      enddo inner
1030 >                   enddo
1031  
1032 <                         vij = vij + vpair
1033 <                         fij(1:3) = fij(1:3) + fpair(1:3)
1034 <                      endif
1035 <                   enddo inner
1036 <                enddo
1037 <
1038 <                if (loop .eq. PAIR_LOOP) then
1039 <                   if (in_switching_region) then
1040 <                      swderiv = vij*dswdr/rgrp
1041 <                      fij(1) = fij(1) + swderiv*d_grp(1)
705 <                      fij(2) = fij(2) + swderiv*d_grp(2)
706 <                      fij(3) = fij(3) + swderiv*d_grp(3)
707 <
708 <                      do ia=groupStartRow(i), groupStartRow(i+1)-1
709 <                         atom1=groupListRow(ia)
710 <                         mf = mfactRow(atom1)
1032 >                   if (loop .eq. PAIR_LOOP) then
1033 >                      if (in_switching_region) then
1034 >                         swderiv = vij*dswdr/rgrp
1035 >                         fij(1) = fij(1) + swderiv*d_grp(1)
1036 >                         fij(2) = fij(2) + swderiv*d_grp(2)
1037 >                         fij(3) = fij(3) + swderiv*d_grp(3)
1038 >                        
1039 >                         do ia=groupStartRow(i), groupStartRow(i+1)-1
1040 >                            atom1=groupListRow(ia)
1041 >                            mf = mfactRow(atom1)
1042   #ifdef IS_MPI
1043 <                         f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1044 <                         f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1045 <                         f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1043 >                            f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1044 >                            f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1045 >                            f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1046   #else
1047 <                         f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1048 <                         f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1049 <                         f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1047 >                            f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1048 >                            f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1049 >                            f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1050   #endif
1051 <                      enddo
1052 <
1053 <                      do jb=groupStartCol(j), groupStartCol(j+1)-1
1054 <                         atom2=groupListCol(jb)
1055 <                         mf = mfactCol(atom2)
1051 >                         enddo
1052 >                        
1053 >                         do jb=groupStartCol(j), groupStartCol(j+1)-1
1054 >                            atom2=groupListCol(jb)
1055 >                            mf = mfactCol(atom2)
1056   #ifdef IS_MPI
1057 <                         f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1058 <                         f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1059 <                         f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1057 >                            f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1058 >                            f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1059 >                            f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1060   #else
1061 <                         f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1062 <                         f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1063 <                         f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1061 >                            f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1062 >                            f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1063 >                            f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1064   #endif
1065 <                      enddo
1066 <                   endif
1065 >                         enddo
1066 >                      endif
1067  
1068 <                   if (do_stress) call add_stress_tensor(d_grp, fij)
1068 >                      if (do_stress) call add_stress_tensor(d_grp, fij)
1069 >                   endif
1070                  endif
1071 <             end if
1071 >             endif
1072            enddo
1073 +          
1074         enddo outer
1075  
1076         if (update_nlist) then
# Line 797 | Line 1130 | contains
1130  
1131      if (do_pot) then
1132         ! scatter/gather pot_row into the members of my column
1133 <       call scatter(pot_Row, pot_Temp, plan_atom_row)
1134 <
1133 >       do i = 1,LR_POT_TYPES
1134 >          call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1135 >       end do
1136         ! scatter/gather pot_local into all other procs
1137         ! add resultant to get total pot
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         pot_Temp = 0.0_DP
1144 <
1145 <       call scatter(pot_Col, pot_Temp, plan_atom_col)
1144 >       do i = 1,LR_POT_TYPES
1145 >          call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1146 >       end do
1147         do i = 1, nlocal
1148 <          pot_local = pot_local + pot_Temp(i)
1148 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1149 >               + pot_Temp(1:LR_POT_TYPES,i)
1150         enddo
1151  
1152      endif
1153   #endif
1154  
1155 <    if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1155 >    if (SIM_requires_postpair_calc) then
1156 >       do i = 1, nlocal            
1157 >          
1158 >          ! we loop only over the local atoms, so we don't need row and column
1159 >          ! lookups for the types
1160 >          
1161 >          me_i = atid(i)
1162 >          
1163 >          ! is the atom electrostatic?  See if it would have an
1164 >          ! electrostatic interaction with itself
1165 >          iHash = InteractionHash(me_i,me_i)
1166  
1167 <       if (FF_uses_RF .and. SIM_uses_RF) then
821 <
822 < #ifdef IS_MPI
823 <          call scatter(rf_Row,rf,plan_atom_row_3d)
824 <          call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
825 <          do i = 1,nlocal
826 <             rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
827 <          end do
828 < #endif
829 <
830 <          do i = 1, nLocal
831 <
832 <             rfpot = 0.0_DP
1167 >          if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1168   #ifdef IS_MPI
1169 <             me_i = atid_row(i)
1169 >             call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1170 >                  t, do_pot)
1171   #else
1172 <             me_i = atid(i)
1172 >             call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1173 >                  t, do_pot)
1174   #endif
1175 <             iHash = InteractionHash(me_i,me_j)
1175 >          endif
1176 >  
1177 >          
1178 >          if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1179              
1180 <             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1181 <
1182 <                mu_i = getDipoleMoment(me_i)
1183 <
1184 <                !! The reaction field needs to include a self contribution
1185 <                !! to the field:
1186 <                call accumulate_self_rf(i, mu_i, eFrame)
1187 <                !! Get the reaction field contribution to the
1188 <                !! potential and torques:
1189 <                call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1180 >             ! loop over the excludes to accumulate RF stuff we've
1181 >             ! left out of the normal pair loop
1182 >            
1183 >             do i1 = 1, nSkipsForAtom(i)
1184 >                j = skipsForAtom(i, i1)
1185 >                
1186 >                ! prevent overcounting of the skips
1187 >                if (i.lt.j) then
1188 >                   call get_interatomic_vector(q(:,i), &
1189 >                        q(:,j), d_atm, ratmsq)
1190 >                   rVal = dsqrt(ratmsq)
1191 >                   call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1192 >                        in_switching_region)
1193   #ifdef IS_MPI
1194 <                pot_local = pot_local + rfpot
1194 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1195 >                        vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1196   #else
1197 <                pot = pot + rfpot
1198 <
1197 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1198 >                        vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1199   #endif
1200 <             endif
1201 <          enddo
1202 <       endif
1200 >                endif
1201 >             enddo
1202 >          endif
1203 >       enddo
1204      endif
1205 <
861 <
1205 >    
1206   #ifdef IS_MPI
1207 <
1207 >    
1208      if (do_pot) then
1209 <       pot = pot + pot_local
1210 <       !! we assume the c code will do the allreduce to get the total potential
867 <       !! we could do it right here if we needed to...
1209 >       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1210 >            mpi_comm_world,mpi_err)            
1211      endif
1212 <
1212 >    
1213      if (do_stress) then
1214         call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1215              mpi_comm_world,mpi_err)
1216         call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1217              mpi_comm_world,mpi_err)
1218      endif
1219 <
1219 >    
1220   #else
1221 <
1221 >    
1222      if (do_stress) then
1223         tau = tau_Temp
1224         virial = virial_Temp
1225      endif
1226 <
1226 >    
1227   #endif
1228 <
1228 >    
1229    end subroutine do_force_loop
1230  
1231    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1232 <       eFrame, A, f, t, pot, vpair, fpair)
1232 >       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1233  
1234 <    real( kind = dp ) :: pot, vpair, sw
1234 >    real( kind = dp ) :: vpair, sw
1235 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1236      real( kind = dp ), dimension(3) :: fpair
1237      real( kind = dp ), dimension(nLocal)   :: mfact
1238      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 899 | Line 1243 | contains
1243      logical, intent(inout) :: do_pot
1244      integer, intent(in) :: i, j
1245      real ( kind = dp ), intent(inout) :: rijsq
1246 <    real ( kind = dp )                :: r
1246 >    real ( kind = dp ), intent(inout) :: r_grp
1247      real ( kind = dp ), intent(inout) :: d(3)
1248 <    real ( kind = dp ) :: ebalance
1248 >    real ( kind = dp ), intent(inout) :: d_grp(3)
1249 >    real ( kind = dp ), intent(inout) :: rCut
1250 >    real ( kind = dp ) :: r
1251      integer :: me_i, me_j
1252  
1253      integer :: iHash
# Line 919 | Line 1265 | contains
1265   #endif
1266  
1267      iHash = InteractionHash(me_i, me_j)
1268 <
1268 >    
1269      if ( iand(iHash, LJ_PAIR).ne.0 ) then
1270 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1270 >       call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1271 >            pot(VDW_POT), f, do_pot)
1272      endif
1273 <
1274 <    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1275 <       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1276 <            pot, eFrame, f, t, do_pot)
930 <
931 <       if (FF_uses_RF .and. SIM_uses_RF) then
932 <
933 <          ! CHECK ME (RF needs to know about all electrostatic types)
934 <          call accumulate_rf(i, j, r, eFrame, sw)
935 <          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
936 <       endif
937 <
1273 >    
1274 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1275 >       call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1276 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1277      endif
1278 <
1278 >    
1279      if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1280         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1281 <            pot, A, f, t, do_pot)
1281 >            pot(HB_POT), A, f, t, do_pot)
1282      endif
1283 <
1283 >    
1284      if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1285         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1286 <            pot, A, f, t, do_pot)
1286 >            pot(HB_POT), A, f, t, do_pot)
1287      endif
1288 <
1288 >    
1289      if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1290         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1291 <            pot, A, f, t, do_pot)
1291 >            pot(VDW_POT), A, f, t, do_pot)
1292      endif
1293      
1294      if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1295 < !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1296 < !           pot, A, f, t, do_pot)
1295 >       call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1296 >            pot(VDW_POT), A, f, t, do_pot)
1297      endif
1298 <
1298 >    
1299      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1300 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1301 <            do_pot)
1300 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1301 >            pot(METALLIC_POT), f, do_pot)
1302      endif
1303 <
1303 >    
1304      if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1305         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1306 <            pot, A, f, t, do_pot)
1306 >            pot(VDW_POT), A, f, t, do_pot)
1307      endif
1308 <
1308 >    
1309      if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1310         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1311 <            pot, A, f, t, do_pot)
1311 >            pot(VDW_POT), A, f, t, do_pot)
1312      endif
1313 +
1314 +    if ( iand(iHash, SC_PAIR).ne.0 ) then      
1315 +       call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1316 +            pot(METALLIC_POT), f, do_pot)
1317 +    endif
1318 +
1319      
1320 +    
1321    end subroutine do_pair
1322  
1323 <  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1323 >  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1324         do_pot, do_stress, eFrame, A, f, t, pot)
1325  
1326 <    real( kind = dp ) :: pot, sw
1326 >    real( kind = dp ) :: sw
1327 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1328      real( kind = dp ), dimension(9,nLocal) :: eFrame
1329      real (kind=dp), dimension(9,nLocal) :: A
1330      real (kind=dp), dimension(3,nLocal) :: f
# Line 985 | Line 1332 | contains
1332  
1333      logical, intent(inout) :: do_pot, do_stress
1334      integer, intent(in) :: i, j
1335 <    real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1335 >    real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1336      real ( kind = dp )                :: r, rc
1337      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1338  
1339      integer :: me_i, me_j, iHash
1340  
1341 +    r = sqrt(rijsq)
1342 +
1343   #ifdef IS_MPI  
1344      me_i = atid_row(i)
1345      me_j = atid_col(j)  
# Line 1002 | Line 1351 | contains
1351      iHash = InteractionHash(me_i, me_j)
1352  
1353      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1354 <            call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1354 >            call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1355      endif
1356 +
1357 +    if ( iand(iHash, SC_PAIR).ne.0 ) then      
1358 +            call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1359 +    endif
1360      
1361    end subroutine do_prepair
1362  
1363  
1364    subroutine do_preforce(nlocal,pot)
1365      integer :: nlocal
1366 <    real( kind = dp ) :: pot
1366 >    real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1367  
1368      if (FF_uses_EAM .and. SIM_uses_EAM) then
1369 <       call calc_EAM_preforce_Frho(nlocal,pot)
1369 >       call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1370      endif
1371 +    if (FF_uses_SC .and. SIM_uses_SC) then
1372 +       call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1373 +    endif
1374  
1375  
1376    end subroutine do_preforce
# Line 1099 | Line 1455 | contains
1455      pot_Col = 0.0_dp
1456      pot_Temp = 0.0_dp
1457  
1102 –    rf_Row = 0.0_dp
1103 –    rf_Col = 0.0_dp
1104 –    rf_Temp = 0.0_dp
1105 –
1458   #endif
1459  
1460      if (FF_uses_EAM .and. SIM_uses_EAM) then
1461         call clean_EAM()
1462      endif
1463  
1112 –    rf = 0.0_dp
1464      tau_Temp = 0.0_dp
1465      virial_Temp = 0.0_dp
1466    end subroutine zero_work_arrays
# Line 1203 | Line 1554 | contains
1554  
1555    function FF_RequiresPrepairCalc() result(doesit)
1556      logical :: doesit
1557 <    doesit = FF_uses_EAM
1557 >    doesit = FF_uses_EAM .or. FF_uses_SC &
1558 >         .or. FF_uses_MEAM
1559    end function FF_RequiresPrepairCalc
1560  
1209 –  function FF_RequiresPostpairCalc() result(doesit)
1210 –    logical :: doesit
1211 –    doesit = FF_uses_RF
1212 –  end function FF_RequiresPostpairCalc
1213 –
1561   #ifdef PROFILE
1562    function getforcetime() result(totalforcetime)
1563      real(kind=dp) :: totalforcetime

Diff Legend

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