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 2267 by tim, Fri Jul 29 17:34:06 2005 UTC vs.
Revision 2532 by tim, Fri Dec 30 21:25:56 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.25 2005-07-29 17:34:06 tim Exp $, $Date: 2005-07-29 17:34:06 $, $Name: not supported by cvs2svn $, $Revision: 1.25 $
48 > !! @version $Id: doForces.F90,v 1.73 2005-12-30 21:25:56 tim Exp $, $Date: 2005-12-30 21:25:56 $, $Name: not supported by cvs2svn $, $Revision: 1.73 $
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  
81  logical, save :: haveRlist = .false.
84    logical, save :: haveNeighborList = .false.
85    logical, save :: haveSIMvariables = .false.
86    logical, save :: haveSaneForceField = .false.
87 <  logical, save :: haveInteractionMap = .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
88  logical, save :: FF_uses_LennardJones
89  logical, save :: FF_uses_Electrostatics
90  logical, save :: FF_uses_Charges
96    logical, save :: FF_uses_Dipoles
92  logical, save :: FF_uses_Quadrupoles
93  logical, save :: FF_uses_Sticky
94  logical, save :: FF_uses_StickyPower
97    logical, save :: FF_uses_GayBerne
98    logical, save :: FF_uses_EAM
99 <  logical, save :: FF_uses_Shapes
100 <  logical, save :: FF_uses_FLARB
101 <  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
102  logical, save :: SIM_uses_LennardJones
103  logical, save :: SIM_uses_Electrostatics
104  logical, save :: SIM_uses_Charges
105  logical, save :: SIM_uses_Dipoles
106  logical, save :: SIM_uses_Quadrupoles
107  logical, save :: SIM_uses_Sticky
108  logical, save :: SIM_uses_StickyPower
109  logical, save :: SIM_uses_GayBerne
104    logical, save :: SIM_uses_EAM
105 <  logical, save :: SIM_uses_Shapes
106 <  logical, save :: SIM_uses_FLARB
113 <  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
117  logical, save :: SIM_uses_molecular_cutoffs
110  
111 <  !!!GO AWAY---------
112 <  !!!!!real(kind=dp), save :: rlist, rlistsq
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
124 !  public :: setRlistDF
125  !public :: addInteraction
126  !public :: setInteractionHash
127  !public :: getInteractionHash
128  public :: createInteractionMap
129  public :: createRcuts
125  
126   #ifdef PROFILE
127    public :: getforcetime
# Line 134 | Line 129 | module doForces
129    real :: forceTimeInitial, forceTimeFinal
130    integer :: nLoops
131   #endif
137
138  type, public :: Interaction
139     integer :: InteractionHash
140     real(kind=dp) :: rList = 0.0_dp
141     real(kind=dp) :: rListSq = 0.0_dp
142  end type Interaction
132    
133 <  type(Interaction), dimension(:,:),allocatable :: InteractionMap
134 <  
133 >  !! Variables for cutoff mapping and interaction mapping
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, target :: groupMaxCutoffRow
138 >  real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
139  
140 <  
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 >
152   contains
153  
154 <
151 <  subroutine createInteractionMap(status)
154 >  subroutine createInteractionHash()
155      integer :: nAtypes
153    integer, intent(out) :: status
156      integer :: i
157      integer :: j
158 <    integer :: ihash
157 <    real(kind=dp) :: myRcut
158 >    integer :: iHash
159      !! Test Types
160      logical :: i_is_LJ
161      logical :: i_is_Elect
# Line 163 | 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 170 | Line 173 | contains
173      logical :: j_is_GB
174      logical :: j_is_EAM
175      logical :: j_is_Shape
176 <    
177 <    status = 0
178 <    
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 createDefaultInteractionMap!")
178 <       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(InteractionMap)) then
193 <       allocate(InteractionMap(nAtypes,nAtypes))
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
207         call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
# Line 198 | 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 211 | 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 232 | 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 241 | Line 262 | contains
262            if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
263  
264  
265 <          InteractionMap(i,j)%InteractionHash = iHash
266 <          InteractionMap(j,i)%InteractionHash = iHash
265 >          InteractionHash(i,j) = iHash
266 >          InteractionHash(j,i) = iHash
267  
268         end do
269  
270      end do
271  
272 <    haveInteractionMap = .true.
273 <  end subroutine createInteractionMap
272 >    haveInteractionHash = .true.
273 >  end subroutine createInteractionHash
274  
275 < ! Query each potential and return the cutoff for that potential. We build the neighbor list based on the largest cutoff value for that atype. Each potential can decide whether to calculate the force for that atype based upon it's own cutoff.
255 <  subroutine createRcuts(defaultRList,stat)
256 <    real(kind=dp), intent(in), optional :: defaultRList
257 <    integer :: iMap
258 <    integer :: map_i,map_j
259 <    real(kind=dp) :: thisRCut = 0.0_dp
260 <    real(kind=dp) :: actualCutoff = 0.0_dp
261 <    integer, intent(out) :: stat
262 <    integer :: nAtypes
263 <    integer :: myStatus
275 >  subroutine createGtypeCutoffMap()
276  
277 <    stat = 0
278 <    if (.not. haveInteractionMap) then
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 <       call createInteractionMap(myStatus)
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  
295 <       if (myStatus .ne. 0) then
296 <          write(default_error, *) 'createInteractionMap failed in doForces!'
297 <          stat = -1
298 <          return
295 >    if (.not. haveInteractionHash) then
296 >       call createInteractionHash()      
297 >    endif
298 > #ifdef IS_MPI
299 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
300 >    nGroupsInCol = getNgroupsInCol(plan_group_col)
301 > #endif
302 >    nAtypes = getSize(atypes)
303 > ! Set all of the initial cutoffs to zero.
304 >    atypeMaxCutoff = 0.0_dp
305 >    do i = 1, nAtypes
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 <    nAtypes = getSize(atypes)
393 <    !! If we pass a default rcut, set all atypes to that cutoff distance
394 <    if(present(defaultRList)) then
395 <       InteractionMap(:,:)%rList = defaultRList
396 <       InteractionMap(:,:)%rListSq = defaultRList*defaultRList
397 <       haveRlist = .true.
398 <       return
392 >
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 <    do map_i = 1,nAtypes
403 <       do map_j = map_i,nAtypes
404 <          iMap = InteractionMap(map_i, map_j)%InteractionHash
405 <          
406 <          if ( iand(iMap, LJ_PAIR).ne.0 ) then
407 <             ! thisRCut = getLJCutOff(map_i,map_j)
408 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
409 <          endif
410 <          
411 <          if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
412 <             ! thisRCut = getElectrostaticCutOff(map_i,map_j)
413 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
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 >       groupMaxCutoffCol = 0.0_dp
416 >       gtypeMaxCutoffCol = 0.0_dp
417 >
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 ( iand(iMap, STICKY_PAIR).ne.0 ) then
302 <             ! thisRCut = getStickyCutOff(map_i,map_j)
303 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
304 <           endif
305 <          
306 <           if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
307 <              ! thisRCut = getStickyPowerCutOff(map_i,map_j)
308 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
309 <           endif
310 <          
311 <           if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
312 <              ! thisRCut = getGayberneCutOff(map_i,map_j)
313 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
314 <           endif
315 <          
316 <           if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
317 < !              thisRCut = getGaybrneLJCutOff(map_i,map_j)
318 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
319 <           endif
320 <          
321 <           if ( iand(iMap, EAM_PAIR).ne.0 ) then      
322 < !              thisRCut = getEAMCutOff(map_i,map_j)
323 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
324 <           endif
325 <          
326 <           if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
327 < !              thisRCut = getShapeCutOff(map_i,map_j)
328 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
329 <           endif
330 <          
331 <           if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
332 < !              thisRCut = getShapeLJCutOff(map_i,map_j)
333 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
334 <           endif
335 <           InteractionMap(map_i, map_j)%rList = actualCutoff
336 <           InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
337 <        end do
338 <     end do
339 <     haveRlist = .true.
340 <  end subroutine createRcuts
528 >          if (thisRcut.gt.largestRcut) largestRcut = thisRcut
529  
530 +          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
531  
532 < !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
533 < !!$  subroutine setRlistDF( this_rlist )
534 < !!$
346 < !!$   real(kind=dp) :: this_rlist
347 < !!$
348 < !!$    rlist = this_rlist
349 < !!$    rlistsq = rlist * rlist
350 < !!$
351 < !!$    haveRlist = .true.
352 < !!$
353 < !!$  end subroutine setRlistDF
532 >          if (.not.haveSkinThickness) then
533 >             skinThickness = 1.0_dp
534 >          endif
535  
536 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
537  
538 <  subroutine setSimVariables()
357 <    SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
358 <    SIM_uses_LennardJones = SimUsesLennardJones()
359 <    SIM_uses_Electrostatics = SimUsesElectrostatics()
360 <    SIM_uses_Charges = SimUsesCharges()
361 <    SIM_uses_Dipoles = SimUsesDipoles()
362 <    SIM_uses_Sticky = SimUsesSticky()
363 <    SIM_uses_StickyPower = SimUsesStickyPower()
364 <    SIM_uses_GayBerne = SimUsesGayBerne()
365 <    SIM_uses_EAM = SimUsesEAM()
366 <    SIM_uses_Shapes = SimUsesShapes()
367 <    SIM_uses_FLARB = SimUsesFLARB()
368 <    SIM_uses_RF = SimUsesRF()
369 <    SIM_requires_postpair_calc = SimRequiresPostpairCalc()
370 <    SIM_requires_prepair_calc = SimRequiresPrepairCalc()
371 <    SIM_uses_PBC = SimUsesPBC()
538 >          ! sanity check
539  
540 <    haveSIMvariables = .true.
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 <    return
549 <  end subroutine setSimVariables
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 +    
641 +     haveSIMvariables = .true.
642 +    
643 +     return
644 +   end subroutine setSimVariables
645 +
646    subroutine doReadyCheck(error)
647      integer, intent(out) :: error
648  
# Line 382 | Line 650 | contains
650  
651      error = 0
652  
653 <    if (.not. haveInteractionMap) then
653 >    if (.not. haveInteractionHash) then      
654 >       call createInteractionHash()      
655 >    endif
656  
657 <       myStatus = 0
657 >    if (.not. haveGtypeCutoffMap) then        
658 >       call createGtypeCutoffMap()      
659 >    endif
660  
389       call createInteractionMap(myStatus)
661  
662 <       if (myStatus .ne. 0) then
663 <          write(default_error, *) 'createInteractionMap failed in doForces!'
393 <          error = -1
394 <          return
395 <       endif
662 >    if (VisitCutoffsAfterComputing) then
663 >       call set_switch(GROUP_SWITCH, largestRcut, largestRcut)      
664      endif
665  
666 +
667      if (.not. haveSIMvariables) then
668         call setSimVariables()
669      endif
670  
671 <    if (.not. haveRlist) then
672 <       write(default_error, *) 'rList has not been set in doForces!'
673 <       error = -1
674 <       return
675 <    endif
671 >  !  if (.not. haveRlist) then
672 >  !     write(default_error, *) 'rList has not been set in doForces!'
673 >  !     error = -1
674 >  !     return
675 >  !  endif
676  
677      if (.not. haveNeighborList) then
678         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 428 | Line 697 | contains
697    end subroutine doReadyCheck
698  
699  
700 <  subroutine init_FF(use_RF_c, thisStat)
700 >  subroutine init_FF(thisStat)
701  
433    logical, intent(in) :: use_RF_c
434
702      integer, intent(out) :: thisStat  
703      integer :: my_status, nMatches
704      integer, pointer :: MatchList(:) => null()
438    real(kind=dp) :: rcut, rrf, rt, dielect
705  
706      !! assume things are copacetic, unless they aren't
707      thisStat = 0
442
443    !! Fortran's version of a cast:
444    FF_uses_RF = use_RF_c
708  
709      !! init_FF is called *after* all of the atom types have been
710      !! defined in atype_module using the new_atype subroutine.
# Line 450 | Line 713 | contains
713      !! interactions are used by the force field.    
714  
715      FF_uses_DirectionalAtoms = .false.
453    FF_uses_LennardJones = .false.
454    FF_uses_Electrostatics = .false.
455    FF_uses_Charges = .false.    
716      FF_uses_Dipoles = .false.
457    FF_uses_Sticky = .false.
458    FF_uses_StickyPower = .false.
717      FF_uses_GayBerne = .false.
718      FF_uses_EAM = .false.
461    FF_uses_Shapes = .false.
462    FF_uses_FLARB = .false.
719  
720      call getMatchingElementList(atypes, "is_Directional", .true., &
721           nMatches, MatchList)
722      if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
723  
468    call getMatchingElementList(atypes, "is_LennardJones", .true., &
469         nMatches, MatchList)
470    if (nMatches .gt. 0) FF_uses_LennardJones = .true.
471
472    call getMatchingElementList(atypes, "is_Electrostatic", .true., &
473         nMatches, MatchList)
474    if (nMatches .gt. 0) then
475       FF_uses_Electrostatics = .true.
476    endif
477
478    call getMatchingElementList(atypes, "is_Charge", .true., &
479         nMatches, MatchList)
480    if (nMatches .gt. 0) then
481       FF_uses_Charges = .true.  
482       FF_uses_Electrostatics = .true.
483    endif
484
724      call getMatchingElementList(atypes, "is_Dipole", .true., &
725           nMatches, MatchList)
726 <    if (nMatches .gt. 0) then
488 <       FF_uses_Dipoles = .true.
489 <       FF_uses_Electrostatics = .true.
490 <       FF_uses_DirectionalAtoms = .true.
491 <    endif
492 <
493 <    call getMatchingElementList(atypes, "is_Quadrupole", .true., &
494 <         nMatches, MatchList)
495 <    if (nMatches .gt. 0) then
496 <       FF_uses_Quadrupoles = .true.
497 <       FF_uses_Electrostatics = .true.
498 <       FF_uses_DirectionalAtoms = .true.
499 <    endif
500 <
501 <    call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
502 <         MatchList)
503 <    if (nMatches .gt. 0) then
504 <       FF_uses_Sticky = .true.
505 <       FF_uses_DirectionalAtoms = .true.
506 <    endif
507 <
508 <    call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
509 <         MatchList)
510 <    if (nMatches .gt. 0) then
511 <       FF_uses_StickyPower = .true.
512 <       FF_uses_DirectionalAtoms = .true.
513 <    endif
726 >    if (nMatches .gt. 0) FF_uses_Dipoles = .true.
727      
728      call getMatchingElementList(atypes, "is_GayBerne", .true., &
729           nMatches, MatchList)
730 <    if (nMatches .gt. 0) then
518 <       FF_uses_GayBerne = .true.
519 <       FF_uses_DirectionalAtoms = .true.
520 <    endif
730 >    if (nMatches .gt. 0) FF_uses_GayBerne = .true.
731  
732      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
733      if (nMatches .gt. 0) FF_uses_EAM = .true.
734  
525    call getMatchingElementList(atypes, "is_Shape", .true., &
526         nMatches, MatchList)
527    if (nMatches .gt. 0) then
528       FF_uses_Shapes = .true.
529       FF_uses_DirectionalAtoms = .true.
530    endif
735  
532    call getMatchingElementList(atypes, "is_FLARB", .true., &
533         nMatches, MatchList)
534    if (nMatches .gt. 0) FF_uses_FLARB = .true.
535
536    !! Assume sanity (for the sake of argument)
736      haveSaneForceField = .true.
538
539    !! check to make sure the FF_uses_RF setting makes sense
540
541    if (FF_uses_dipoles) then
542       if (FF_uses_RF) then
543          dielect = getDielect()
544          call initialize_rf(dielect)
545       endif
546    else
547       if (FF_uses_RF) then          
548          write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
549          thisStat = -1
550          haveSaneForceField = .false.
551          return
552       endif
553    endif
737  
555    !sticky module does not contain check_sticky_FF anymore
556    !if (FF_uses_sticky) then
557    !   call check_sticky_FF(my_status)
558    !   if (my_status /= 0) then
559    !      thisStat = -1
560    !      haveSaneForceField = .false.
561    !      return
562    !   end if
563    !endif
564
738      if (FF_uses_EAM) then
739         call init_EAM_FF(my_status)
740         if (my_status /= 0) then
# Line 572 | Line 745 | contains
745         end if
746      endif
747  
575    if (FF_uses_GayBerne) then
576       call check_gb_pair_FF(my_status)
577       if (my_status .ne. 0) then
578          thisStat = -1
579          haveSaneForceField = .false.
580          return
581       endif
582    endif
583
584    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
585    endif
586
748      if (.not. haveNeighborList) then
749         !! Create neighbor lists
750         call expandNeighborList(nLocal, my_status)
# Line 617 | Line 778 | contains
778  
779      !! Stress Tensor
780      real( kind = dp), dimension(9) :: tau  
781 <    real ( kind = dp ) :: pot
781 >    real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
782      logical ( kind = 2) :: do_pot_c, do_stress_c
783      logical :: do_pot
784      logical :: do_stress
785      logical :: in_switching_region
786   #ifdef IS_MPI
787 <    real( kind = DP ) :: pot_local
787 >    real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
788      integer :: nAtomsInRow
789      integer :: nAtomsInCol
790      integer :: nprocs
# Line 638 | Line 799 | contains
799      integer :: nlist
800      real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
801      real( kind = DP ) :: sw, dswdr, swderiv, mf
802 +    real( kind = DP ) :: rVal
803      real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
804      real(kind=dp) :: rfpot, mu_i, virial
805 +    real(kind=dp):: rCut
806      integer :: me_i, me_j, n_in_i, n_in_j
807      logical :: is_dp_i
808      integer :: neighborListSize
# Line 647 | Line 810 | contains
810      integer :: localError
811      integer :: propPack_i, propPack_j
812      integer :: loopStart, loopEnd, loop
813 <    integer :: iMap
814 <    real(kind=dp) :: listSkin = 1.0  
813 >    integer :: iHash
814 >    integer :: i1
815 >  
816  
817      !! initialize local variables  
818  
# Line 712 | Line 876 | contains
876         ! (but only on the first time through):
877         if (loop .eq. loopStart) then
878   #ifdef IS_MPI
879 <          call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
879 >          call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
880                 update_nlist)
881   #else
882 <          call checkNeighborList(nGroups, q_group, listSkin, &
882 >          call checkNeighborList(nGroups, q_group, skinThickness, &
883                 update_nlist)
884   #endif
885         endif
# Line 738 | Line 902 | contains
902         iend = nGroups - 1
903   #endif
904         outer: do i = istart, iend
741
742 #ifdef IS_MPI
743             me_i = atid_row(i)
744 #else
745             me_i = atid(i)
746 #endif
905  
906            if (update_nlist) point(i) = nlist + 1
907  
# Line 779 | Line 937 | contains
937               me_j = atid(j)
938               call get_interatomic_vector(q_group(:,i), &
939                    q_group(:,j), d_grp, rgrpsq)
940 < #endif
940 > #endif      
941  
942 <             if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
942 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
943                  if (update_nlist) then
944                     nlist = nlist + 1
945  
# Line 801 | Line 959 | contains
959  
960                     list(nlist) = j
961                  endif
962 +
963  
964 <                if (loop .eq. PAIR_LOOP) then
965 <                   vij = 0.0d0
807 <                   fij(1:3) = 0.0d0
808 <                endif
964 >                
965 >                if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
966  
967 <                call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
968 <                     in_switching_region)
969 <
970 <                n_in_j = groupStartCol(j+1) - groupStartCol(j)
971 <
972 <                do ia = groupStartRow(i), groupStartRow(i+1)-1
973 <
974 <                   atom1 = groupListRow(ia)
975 <
976 <                   inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
977 <
978 <                      atom2 = groupListCol(jb)
979 <
980 <                      if (skipThisPair(atom1, atom2)) cycle inner
981 <
982 <                      if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
983 <                         d_atm(1:3) = d_grp(1:3)
984 <                         ratmsq = rgrpsq
985 <                      else
986 < #ifdef IS_MPI
987 <                         call get_interatomic_vector(q_Row(:,atom1), &
988 <                              q_Col(:,atom2), d_atm, ratmsq)
989 < #else
990 <                         call get_interatomic_vector(q(:,atom1), &
991 <                              q(:,atom2), d_atm, ratmsq)
992 < #endif
993 <                      endif
994 <
995 <                      if (loop .eq. PREPAIR_LOOP) then
967 >                   rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
968 >                   if (loop .eq. PAIR_LOOP) then
969 >                      vij = 0.0d0
970 >                      fij(1:3) = 0.0d0
971 >                   endif
972 >                  
973 >                   call get_switch(rgrpsq, sw, dswdr, rgrp, &
974 >                        group_switch, in_switching_region)
975 >                  
976 >                   n_in_j = groupStartCol(j+1) - groupStartCol(j)
977 >                  
978 >                   do ia = groupStartRow(i), groupStartRow(i+1)-1
979 >                      
980 >                      atom1 = groupListRow(ia)
981 >                      
982 >                      inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
983 >                        
984 >                         atom2 = groupListCol(jb)
985 >                        
986 >                         if (skipThisPair(atom1, atom2))  cycle inner
987 >                        
988 >                         if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
989 >                            d_atm(1:3) = d_grp(1:3)
990 >                            ratmsq = rgrpsq
991 >                         else
992 > #ifdef IS_MPI
993 >                            call get_interatomic_vector(q_Row(:,atom1), &
994 >                                 q_Col(:,atom2), d_atm, ratmsq)
995 > #else
996 >                            call get_interatomic_vector(q(:,atom1), &
997 >                                 q(:,atom2), d_atm, ratmsq)
998 > #endif
999 >                         endif
1000 >                        
1001 >                         if (loop .eq. PREPAIR_LOOP) then
1002   #ifdef IS_MPI                      
1003 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1004 <                              rgrpsq, d_grp, do_pot, do_stress, &
1005 <                              eFrame, A, f, t, pot_local)
1003 >                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1004 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1005 >                                 eFrame, A, f, t, pot_local)
1006   #else
1007 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1008 <                              rgrpsq, d_grp, do_pot, do_stress, &
1009 <                              eFrame, A, f, t, pot)
1007 >                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1008 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1009 >                                 eFrame, A, f, t, pot)
1010   #endif                                              
1011 <                      else
1011 >                         else
1012   #ifdef IS_MPI                      
1013 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1014 <                              do_pot, &
1015 <                              eFrame, A, f, t, pot_local, vpair, fpair)
1013 >                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1014 >                                 do_pot, eFrame, A, f, t, pot_local, vpair, &
1015 >                                 fpair, d_grp, rgrp, rCut)
1016   #else
1017 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1018 <                              do_pot,  &
1019 <                              eFrame, A, f, t, pot, vpair, fpair)
1017 >                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1018 >                                 do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1019 >                                 d_grp, rgrp, rCut)
1020   #endif
1021 +                            vij = vij + vpair
1022 +                            fij(1:3) = fij(1:3) + fpair(1:3)
1023 +                         endif
1024 +                      enddo inner
1025 +                   enddo
1026  
1027 <                         vij = vij + vpair
1028 <                         fij(1:3) = fij(1:3) + fpair(1:3)
1029 <                      endif
1030 <                   enddo inner
1031 <                enddo
1032 <
1033 <                if (loop .eq. PAIR_LOOP) then
1034 <                   if (in_switching_region) then
1035 <                      swderiv = vij*dswdr/rgrp
1036 <                      fij(1) = fij(1) + swderiv*d_grp(1)
869 <                      fij(2) = fij(2) + swderiv*d_grp(2)
870 <                      fij(3) = fij(3) + swderiv*d_grp(3)
871 <
872 <                      do ia=groupStartRow(i), groupStartRow(i+1)-1
873 <                         atom1=groupListRow(ia)
874 <                         mf = mfactRow(atom1)
1027 >                   if (loop .eq. PAIR_LOOP) then
1028 >                      if (in_switching_region) then
1029 >                         swderiv = vij*dswdr/rgrp
1030 >                         fij(1) = fij(1) + swderiv*d_grp(1)
1031 >                         fij(2) = fij(2) + swderiv*d_grp(2)
1032 >                         fij(3) = fij(3) + swderiv*d_grp(3)
1033 >                        
1034 >                         do ia=groupStartRow(i), groupStartRow(i+1)-1
1035 >                            atom1=groupListRow(ia)
1036 >                            mf = mfactRow(atom1)
1037   #ifdef IS_MPI
1038 <                         f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1039 <                         f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1040 <                         f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1038 >                            f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1039 >                            f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1040 >                            f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1041   #else
1042 <                         f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1043 <                         f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1044 <                         f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1042 >                            f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1043 >                            f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1044 >                            f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1045   #endif
1046 <                      enddo
1047 <
1048 <                      do jb=groupStartCol(j), groupStartCol(j+1)-1
1049 <                         atom2=groupListCol(jb)
1050 <                         mf = mfactCol(atom2)
1046 >                         enddo
1047 >                        
1048 >                         do jb=groupStartCol(j), groupStartCol(j+1)-1
1049 >                            atom2=groupListCol(jb)
1050 >                            mf = mfactCol(atom2)
1051   #ifdef IS_MPI
1052 <                         f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1053 <                         f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1054 <                         f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1052 >                            f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1053 >                            f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1054 >                            f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1055   #else
1056 <                         f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1057 <                         f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1058 <                         f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1056 >                            f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1057 >                            f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1058 >                            f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1059   #endif
1060 <                      enddo
1061 <                   endif
1060 >                         enddo
1061 >                      endif
1062  
1063 <                   if (do_stress) call add_stress_tensor(d_grp, fij)
1063 >                      if (do_stress) call add_stress_tensor(d_grp, fij)
1064 >                   endif
1065                  endif
1066 <             end if
1066 >             endif
1067            enddo
1068 +          
1069         enddo outer
1070  
1071         if (update_nlist) then
# Line 961 | Line 1125 | contains
1125  
1126      if (do_pot) then
1127         ! scatter/gather pot_row into the members of my column
1128 <       call scatter(pot_Row, pot_Temp, plan_atom_row)
1129 <
1128 >       do i = 1,LR_POT_TYPES
1129 >          call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1130 >       end do
1131         ! scatter/gather pot_local into all other procs
1132         ! add resultant to get total pot
1133         do i = 1, nlocal
1134 <          pot_local = pot_local + pot_Temp(i)
1134 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1135 >               + pot_Temp(1:LR_POT_TYPES,i)
1136         enddo
1137  
1138         pot_Temp = 0.0_DP
1139 <
1140 <       call scatter(pot_Col, pot_Temp, plan_atom_col)
1139 >       do i = 1,LR_POT_TYPES
1140 >          call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1141 >       end do
1142         do i = 1, nlocal
1143 <          pot_local = pot_local + pot_Temp(i)
1143 >          pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1144 >               + pot_Temp(1:LR_POT_TYPES,i)
1145         enddo
1146  
1147      endif
1148   #endif
1149  
1150 <    if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1150 >    if (SIM_requires_postpair_calc) then
1151 >       do i = 1, nlocal            
1152 >          
1153 >          ! we loop only over the local atoms, so we don't need row and column
1154 >          ! lookups for the types
1155 >          
1156 >          me_i = atid(i)
1157 >          
1158 >          ! is the atom electrostatic?  See if it would have an
1159 >          ! electrostatic interaction with itself
1160 >          iHash = InteractionHash(me_i,me_i)
1161  
1162 <       if (FF_uses_RF .and. SIM_uses_RF) then
985 <
1162 >          if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1163   #ifdef IS_MPI
1164 <          call scatter(rf_Row,rf,plan_atom_row_3d)
1165 <          call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
989 <          do i = 1,nlocal
990 <             rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
991 <          end do
992 < #endif
993 <
994 <          do i = 1, nLocal
995 <
996 <             rfpot = 0.0_DP
997 < #ifdef IS_MPI
998 <             me_i = atid_row(i)
1164 >             call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1165 >                  t, do_pot)
1166   #else
1167 <             me_i = atid(i)
1167 >             call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1168 >                  t, do_pot)
1169   #endif
1170 <             iMap = InteractionMap(me_i, me_j)%InteractionHash
1170 >          endif
1171 >  
1172 >          
1173 >          if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1174              
1175 <             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1176 <
1177 <                mu_i = getDipoleMoment(me_i)
1178 <
1179 <                !! The reaction field needs to include a self contribution
1180 <                !! to the field:
1181 <                call accumulate_self_rf(i, mu_i, eFrame)
1182 <                !! Get the reaction field contribution to the
1183 <                !! potential and torques:
1184 <                call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1175 >             ! loop over the excludes to accumulate RF stuff we've
1176 >             ! left out of the normal pair loop
1177 >            
1178 >             do i1 = 1, nSkipsForAtom(i)
1179 >                j = skipsForAtom(i, i1)
1180 >                
1181 >                ! prevent overcounting of the skips
1182 >                if (i.lt.j) then
1183 >                   call get_interatomic_vector(q(:,i), &
1184 >                        q(:,j), d_atm, ratmsq)
1185 >                   rVal = dsqrt(ratmsq)
1186 >                   call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1187 >                        in_switching_region)
1188   #ifdef IS_MPI
1189 <                pot_local = pot_local + rfpot
1189 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1190 >                        vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1191   #else
1192 <                pot = pot + rfpot
1193 <
1192 >                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1193 >                        vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1194   #endif
1195 <             endif
1196 <          enddo
1197 <       endif
1195 >                endif
1196 >             enddo
1197 >          endif
1198 >       enddo
1199      endif
1200 <
1025 <
1200 >    
1201   #ifdef IS_MPI
1202 <
1202 >    
1203      if (do_pot) then
1204 <       pot = pot + pot_local
1205 <       !! we assume the c code will do the allreduce to get the total potential
1031 <       !! we could do it right here if we needed to...
1204 >       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1205 >            mpi_comm_world,mpi_err)            
1206      endif
1207 <
1207 >    
1208      if (do_stress) then
1209         call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1210              mpi_comm_world,mpi_err)
1211         call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1212              mpi_comm_world,mpi_err)
1213      endif
1214 <
1214 >    
1215   #else
1216 <
1216 >    
1217      if (do_stress) then
1218         tau = tau_Temp
1219         virial = virial_Temp
1220      endif
1221 <
1221 >    
1222   #endif
1223 <
1223 >    
1224    end subroutine do_force_loop
1225  
1226    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1227 <       eFrame, A, f, t, pot, vpair, fpair)
1227 >       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1228  
1229 <    real( kind = dp ) :: pot, vpair, sw
1229 >    real( kind = dp ) :: vpair, sw
1230 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1231      real( kind = dp ), dimension(3) :: fpair
1232      real( kind = dp ), dimension(nLocal)   :: mfact
1233      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 1063 | Line 1238 | contains
1238      logical, intent(inout) :: do_pot
1239      integer, intent(in) :: i, j
1240      real ( kind = dp ), intent(inout) :: rijsq
1241 <    real ( kind = dp )                :: r
1241 >    real ( kind = dp ), intent(inout) :: r_grp
1242      real ( kind = dp ), intent(inout) :: d(3)
1243 <    real ( kind = dp ) :: ebalance
1243 >    real ( kind = dp ), intent(inout) :: d_grp(3)
1244 >    real ( kind = dp ), intent(inout) :: rCut
1245 >    real ( kind = dp ) :: r
1246      integer :: me_i, me_j
1247  
1248 <    integer :: iMap
1248 >    integer :: iHash
1249  
1250      r = sqrt(rijsq)
1251      vpair = 0.0d0
# Line 1081 | Line 1258 | contains
1258      me_i = atid(i)
1259      me_j = atid(j)
1260   #endif
1084
1085    iMap = InteractionMap(me_i, me_j)%InteractionHash
1086
1087    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1088       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1089    endif
1090
1091    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1092       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1093            pot, eFrame, f, t, do_pot)
1261  
1262 <       if (FF_uses_RF .and. SIM_uses_RF) then
1263 <
1264 <          ! CHECK ME (RF needs to know about all electrostatic types)
1265 <          call accumulate_rf(i, j, r, eFrame, sw)
1266 <          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1100 <       endif
1101 <
1262 >    iHash = InteractionHash(me_i, me_j)
1263 >    
1264 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1265 >       call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1266 >            pot(VDW_POT), f, do_pot)
1267      endif
1268 <
1269 <    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1268 >    
1269 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1270 >       call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1271 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1272 >    endif
1273 >    
1274 >    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1275         call do_sticky_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 <
1279 <    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1278 >    
1279 >    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1280         call do_sticky_power_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 <
1284 <    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1283 >    
1284 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1285         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1286 <            pot, A, f, t, do_pot)
1286 >            pot(VDW_POT), A, f, t, do_pot)
1287      endif
1288      
1289 <    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1290 < !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1291 < !           pot, A, f, t, do_pot)
1289 >    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1290 >       call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1291 >            pot(VDW_POT), A, f, t, do_pot)
1292      endif
1293 <
1294 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1295 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1296 <            do_pot)
1293 >    
1294 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1295 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1296 >            pot(METALLIC_POT), f, do_pot)
1297      endif
1298 <
1299 <    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1298 >    
1299 >    if ( iand(iHash, SHAPE_PAIR).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(iMap, SHAPE_LJ).ne.0 ) then      
1303 >    
1304 >    if ( iand(iHash, SHAPE_LJ).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 +
1309 +    if ( iand(iHash, SC_PAIR).ne.0 ) then      
1310 +       call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1311 +            pot(METALLIC_POT), f, do_pot)
1312 +    endif
1313 +
1314      
1315 +    
1316    end subroutine do_pair
1317  
1318 <  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1318 >  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1319         do_pot, do_stress, eFrame, A, f, t, pot)
1320  
1321 <    real( kind = dp ) :: pot, sw
1321 >    real( kind = dp ) :: sw
1322 >    real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1323      real( kind = dp ), dimension(9,nLocal) :: eFrame
1324      real (kind=dp), dimension(9,nLocal) :: A
1325      real (kind=dp), dimension(3,nLocal) :: f
# Line 1149 | Line 1327 | contains
1327  
1328      logical, intent(inout) :: do_pot, do_stress
1329      integer, intent(in) :: i, j
1330 <    real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1330 >    real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1331      real ( kind = dp )                :: r, rc
1332      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1333  
1334 <    integer :: me_i, me_j, iMap
1334 >    integer :: me_i, me_j, iHash
1335  
1336 +    r = sqrt(rijsq)
1337 +
1338   #ifdef IS_MPI  
1339      me_i = atid_row(i)
1340      me_j = atid_col(j)  
# Line 1163 | Line 1343 | contains
1343      me_j = atid(j)  
1344   #endif
1345  
1346 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1346 >    iHash = InteractionHash(me_i, me_j)
1347  
1348 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1349 <            call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1348 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1349 >            call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1350      endif
1351 +
1352 +    if ( iand(iHash, SC_PAIR).ne.0 ) then      
1353 +            call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1354 +    endif
1355      
1356    end subroutine do_prepair
1357  
1358  
1359    subroutine do_preforce(nlocal,pot)
1360      integer :: nlocal
1361 <    real( kind = dp ) :: pot
1361 >    real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1362  
1363      if (FF_uses_EAM .and. SIM_uses_EAM) then
1364 <       call calc_EAM_preforce_Frho(nlocal,pot)
1364 >       call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1365      endif
1366 +    if (FF_uses_SC .and. SIM_uses_SC) then
1367 +       call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1368 +    endif
1369  
1370  
1371    end subroutine do_preforce
# Line 1263 | Line 1450 | contains
1450      pot_Col = 0.0_dp
1451      pot_Temp = 0.0_dp
1452  
1266    rf_Row = 0.0_dp
1267    rf_Col = 0.0_dp
1268    rf_Temp = 0.0_dp
1269
1453   #endif
1454  
1455      if (FF_uses_EAM .and. SIM_uses_EAM) then
1456         call clean_EAM()
1457      endif
1458  
1276    rf = 0.0_dp
1459      tau_Temp = 0.0_dp
1460      virial_Temp = 0.0_dp
1461    end subroutine zero_work_arrays
# Line 1362 | Line 1544 | contains
1544  
1545    function FF_UsesDirectionalAtoms() result(doesit)
1546      logical :: doesit
1547 <    doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1366 <         FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1367 <         FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1547 >    doesit = FF_uses_DirectionalAtoms
1548    end function FF_UsesDirectionalAtoms
1549  
1550    function FF_RequiresPrepairCalc() result(doesit)
1551      logical :: doesit
1552 <    doesit = FF_uses_EAM
1552 >    doesit = FF_uses_EAM .or. FF_uses_SC &
1553 >         .or. FF_uses_MEAM
1554    end function FF_RequiresPrepairCalc
1555  
1375  function FF_RequiresPostpairCalc() result(doesit)
1376    logical :: doesit
1377    doesit = FF_uses_RF
1378  end function FF_RequiresPostpairCalc
1379
1556   #ifdef PROFILE
1557    function getforcetime() result(totalforcetime)
1558      real(kind=dp) :: totalforcetime

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines