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 2342 by chrisfen, Tue Oct 4 19:32:58 2005 UTC vs.
Revision 2512 by gezelter, Thu Dec 15 21:43:16 2005 UTC

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

Diff Legend

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