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 2262 by chuckv, Sun Jul 3 20:53:43 2005 UTC vs.
Revision 2355 by chuckv, Wed Oct 12 18:59: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.23 2005-07-03 20:53:43 chuckv Exp $, $Date: 2005-07-03 20:53:43 $, $Name: not supported by cvs2svn $, $Revision: 1.23 $
48 > !! @version $Id: doForces.F90,v 1.54 2005-10-12 18:59:16 chuckv Exp $, $Date: 2005-10-12 18:59:16 $, $Name: not supported by cvs2svn $, $Revision: 1.54 $
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
61 >  use reaction_field_module
62    use gb_pair
63    use shapes
64    use vector_class
# 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 :: haveRlist = .false.
91  
92    logical, save :: FF_uses_DirectionalAtoms
88 –  logical, save :: FF_uses_LennardJones
89 –  logical, save :: FF_uses_Electrostatics
90 –  logical, save :: FF_uses_Charges
93    logical, save :: FF_uses_Dipoles
92 –  logical, save :: FF_uses_Quadrupoles
93 –  logical, save :: FF_uses_Sticky
94 –  logical, save :: FF_uses_StickyPower
94    logical, save :: FF_uses_GayBerne
95    logical, save :: FF_uses_EAM
97 –  logical, save :: FF_uses_Shapes
98 –  logical, save :: FF_uses_FLARB
99 –  logical, save :: FF_uses_RF
96  
97    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
98    logical, save :: SIM_uses_EAM
111 –  logical, save :: SIM_uses_Shapes
112 –  logical, save :: SIM_uses_FLARB
113 –  logical, save :: SIM_uses_RF
99    logical, save :: SIM_requires_postpair_calc
100    logical, save :: SIM_requires_prepair_calc
101    logical, save :: SIM_uses_PBC
117 –  logical, save :: SIM_uses_molecular_cutoffs
102  
103 <  !!!GO AWAY---------
120 <  !!!!!real(kind=dp), save :: rlist, rlistsq
103 >  integer, save :: electrostaticSummationMethod
104  
105    public :: init_FF
106 +  public :: setDefaultCutoffs
107    public :: do_force_loop
108 < !  public :: setRlistDF
109 <  !public :: addInteraction
110 <  !public :: setInteractionHash
111 <  !public :: getInteractionHash
112 <  public :: createInteractionMap
113 <  public :: createRcuts
108 >  public :: createInteractionHash
109 >  public :: createGtypeCutoffMap
110 >  public :: getStickyCut
111 >  public :: getStickyPowerCut
112 >  public :: getGayBerneCut
113 >  public :: getEAMCut
114 >  public :: getShapeCut
115  
116   #ifdef PROFILE
117    public :: getforcetime
# Line 134 | Line 119 | module doForces
119    real :: forceTimeInitial, forceTimeFinal
120    integer :: nLoops
121   #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
122    
123 <  type(Interaction), dimension(:,:),allocatable :: InteractionMap
124 <  
123 >  !! Variables for cutoff mapping and interaction mapping
124 >  ! Bit hash to determine pair-pair interactions.
125 >  integer, dimension(:,:), allocatable :: InteractionHash
126 >  real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
127 >  real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
128 >  real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
129  
130 +  integer, dimension(:), allocatable, target :: groupToGtypeRow
131 +  integer, dimension(:), pointer :: groupToGtypeCol => null()
132 +
133 +  real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
134 +  real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
135 +  type ::gtypeCutoffs
136 +     real(kind=dp) :: rcut
137 +     real(kind=dp) :: rcutsq
138 +     real(kind=dp) :: rlistsq
139 +  end type gtypeCutoffs
140 +  type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
141 +
142 +  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
143 +  real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
144 +  real(kind=dp),save :: listSkin
145    
146   contains
147  
148 <
151 <  subroutine createInteractionMap(status)
148 >  subroutine createInteractionHash(status)
149      integer :: nAtypes
150 <    integer :: status
150 >    integer, intent(out) :: status
151      integer :: i
152      integer :: j
153 <    integer :: ihash
154 <    real(kind=dp) :: myRcut
158 < ! Test Types
153 >    integer :: iHash
154 >    !! Test Types
155      logical :: i_is_LJ
156      logical :: i_is_Elect
157      logical :: i_is_Sticky
# Line 170 | Line 166 | contains
166      logical :: j_is_GB
167      logical :: j_is_EAM
168      logical :: j_is_Shape
169 <    
170 <    
169 >    real(kind=dp) :: myRcut
170 >
171 >    status = 0  
172 >
173      if (.not. associated(atypes)) then
174 <       call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
174 >       call handleError("atype", "atypes was not present before call of createInteractionHash!")
175         status = -1
176         return
177      endif
# Line 185 | Line 183 | contains
183         return
184      end if
185  
186 <    if (.not. allocated(InteractionMap)) then
187 <       allocate(InteractionMap(nAtypes,nAtypes))
186 >    if (.not. allocated(InteractionHash)) then
187 >       allocate(InteractionHash(nAtypes,nAtypes))
188 >    else
189 >       deallocate(InteractionHash)
190 >       allocate(InteractionHash(nAtypes,nAtypes))
191      endif
192 +
193 +    if (.not. allocated(atypeMaxCutoff)) then
194 +       allocate(atypeMaxCutoff(nAtypes))
195 +    else
196 +       deallocate(atypeMaxCutoff)
197 +       allocate(atypeMaxCutoff(nAtypes))
198 +    endif
199          
200      do i = 1, nAtypes
201         call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
# Line 240 | Line 248 | contains
248            if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
249  
250  
251 <          InteractionMap(i,j)%InteractionHash = iHash
252 <          InteractionMap(j,i)%InteractionHash = iHash
251 >          InteractionHash(i,j) = iHash
252 >          InteractionHash(j,i) = iHash
253  
254         end do
255  
256      end do
249 –  end subroutine createInteractionMap
257  
258 < ! 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.
259 <  subroutine createRcuts(defaultRList)
253 <    real(kind=dp), intent(in), optional :: defaultRList
254 <    integer :: iMap
255 <    integer :: map_i,map_j
256 <    real(kind=dp) :: thisRCut = 0.0_dp
257 <    real(kind=dp) :: actualCutoff = 0.0_dp
258 <    integer :: nAtypes
258 >    haveInteractionHash = .true.
259 >  end subroutine createInteractionHash
260  
261 <    if(.not. allocated(InteractionMap)) return
261 >  subroutine createGtypeCutoffMap(stat)
262  
263 <    nAtypes = getSize(atypes)
264 < ! If we pass a default rcut, set all atypes to that cutoff distance
265 <    if(present(defaultRList)) then
266 <       InteractionMap(:,:)%rList = defaultRList
267 <       InteractionMap(:,:)%rListSq = defaultRList*defaultRList
268 <       haveRlist = .true.
269 <       return
270 <    end if
263 >    integer, intent(out), optional :: stat
264 >    logical :: i_is_LJ
265 >    logical :: i_is_Elect
266 >    logical :: i_is_Sticky
267 >    logical :: i_is_StickyP
268 >    logical :: i_is_GB
269 >    logical :: i_is_EAM
270 >    logical :: i_is_Shape
271 >    logical :: GtypeFound
272  
273 <    do map_i = 1,nAtypes
274 <       do map_j = map_i,nAtypes
275 <          iMap = InteractionMap(map_i, map_j)%InteractionHash
273 >    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
274 >    integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
275 >    integer :: nGroupsInRow
276 >    integer :: nGroupsInCol
277 >    integer :: nGroupTypesRow,nGroupTypesCol
278 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin
279 >    real(kind=dp) :: biggestAtypeCutoff
280 >
281 >    stat = 0
282 >    if (.not. haveInteractionHash) then
283 >       call createInteractionHash(myStatus)      
284 >       if (myStatus .ne. 0) then
285 >          write(default_error, *) 'createInteractionHash failed in doForces!'
286 >          stat = -1
287 >          return
288 >       endif
289 >    endif
290 > #ifdef IS_MPI
291 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
292 >    nGroupsInCol = getNgroupsInCol(plan_group_col)
293 > #endif
294 >    nAtypes = getSize(atypes)
295 > ! Set all of the initial cutoffs to zero.
296 >    atypeMaxCutoff = 0.0_dp
297 >    do i = 1, nAtypes
298 >       if (SimHasAtype(i)) then    
299 >          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
300 >          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
301 >          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
302 >          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
303 >          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
304 >          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
305 >          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
306            
307 <          if ( iand(iMap, LJ_PAIR).ne.0 ) then
308 < !            thisRCut = getLJCutOff(map_i,map_j)
309 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
310 <          endif
311 <          
312 <          if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
313 < !            thisRCut = getElectrostaticCutOff(map_i,map_j)
314 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
307 >
308 >          if (haveDefaultCutoffs) then
309 >             atypeMaxCutoff(i) = defaultRcut
310 >          else
311 >             if (i_is_LJ) then          
312 >                thisRcut = getSigma(i) * 2.5_dp
313 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
314 >             endif
315 >             if (i_is_Elect) then
316 >                thisRcut = defaultRcut
317 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
318 >             endif
319 >             if (i_is_Sticky) then
320 >                thisRcut = getStickyCut(i)
321 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
322 >             endif
323 >             if (i_is_StickyP) then
324 >                thisRcut = getStickyPowerCut(i)
325 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
326 >             endif
327 >             if (i_is_GB) then
328 >                thisRcut = getGayBerneCut(i)
329 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
330 >             endif
331 >             if (i_is_EAM) then
332 >                thisRcut = getEAMCut(i)
333 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
334 >             endif
335 >             if (i_is_Shape) then
336 >                thisRcut = getShapeCut(i)
337 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
338 >             endif
339            endif
340            
341 <          if ( iand(iMap, STICKY_PAIR).ne.0 ) then
342 < !             thisRCut = getStickyCutOff(map_i,map_j)
343 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
344 <           endif
289 <          
290 <           if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
291 < !              thisRCut = getStickyPowerCutOff(map_i,map_j)
292 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
293 <           endif
294 <          
295 <           if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
296 < !              thisRCut = getGayberneCutOff(map_i,map_j)
297 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
298 <           endif
299 <          
300 <           if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
301 < !              thisRCut = getGaybrneLJCutOff(map_i,map_j)
302 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
303 <           endif
304 <          
305 <           if ( iand(iMap, EAM_PAIR).ne.0 ) then      
306 < !              thisRCut = getEAMCutOff(map_i,map_j)
307 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
308 <           endif
309 <          
310 <           if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
311 < !              thisRCut = getShapeCutOff(map_i,map_j)
312 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
313 <           endif
314 <          
315 <           if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
316 < !              thisRCut = getShapeLJCutOff(map_i,map_j)
317 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
318 <           endif
319 <           InteractionMap(map_i, map_j)%rList = actualCutoff
320 <           InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
321 <        end do
322 <     end do
323 <          haveRlist = .true.
324 <  end subroutine createRcuts
341 >          
342 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
343 >             biggestAtypeCutoff = atypeMaxCutoff(i)
344 >          endif
345  
346 +       endif
347 +    enddo
348 +  
349  
350 < !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
351 < !!$  subroutine setRlistDF( this_rlist )
352 < !!$
353 < !!$   real(kind=dp) :: this_rlist
354 < !!$
355 < !!$    rlist = this_rlist
356 < !!$    rlistsq = rlist * rlist
357 < !!$
358 < !!$    haveRlist = .true.
359 < !!$
360 < !!$  end subroutine setRlistDF
350 >    
351 >    istart = 1
352 >    jstart = 1
353 > #ifdef IS_MPI
354 >    iend = nGroupsInRow
355 >    jend = nGroupsInCol
356 > #else
357 >    iend = nGroups
358 >    jend = nGroups
359 > #endif
360 >    
361 >    !! allocate the groupToGtype and gtypeMaxCutoff here.
362 >    if(.not.allocated(groupToGtypeRow)) then
363 >     !  allocate(groupToGtype(iend))
364 >       allocate(groupToGtypeRow(iend))
365 >    else
366 >       deallocate(groupToGtypeRow)
367 >       allocate(groupToGtypeRow(iend))
368 >    endif
369 >    if(.not.allocated(groupMaxCutoffRow)) then
370 >       allocate(groupMaxCutoffRow(iend))
371 >    else
372 >       deallocate(groupMaxCutoffRow)
373 >       allocate(groupMaxCutoffRow(iend))
374 >    end if
375  
376 +    if(.not.allocated(gtypeMaxCutoffRow)) then
377 +       allocate(gtypeMaxCutoffRow(iend))
378 +    else
379 +       deallocate(gtypeMaxCutoffRow)
380 +       allocate(gtypeMaxCutoffRow(iend))
381 +    endif
382  
383 +
384 + #ifdef IS_MPI
385 +       ! We only allocate new storage if we are in MPI because Ncol /= Nrow
386 +    if(.not.associated(groupToGtypeCol)) then
387 +       allocate(groupToGtypeCol(jend))
388 +    else
389 +       deallocate(groupToGtypeCol)
390 +       allocate(groupToGtypeCol(jend))
391 +    end if
392 +
393 +    if(.not.associated(groupToGtypeCol)) then
394 +       allocate(groupToGtypeCol(jend))
395 +    else
396 +       deallocate(groupToGtypeCol)
397 +       allocate(groupToGtypeCol(jend))
398 +    end if
399 +    if(.not.associated(gtypeMaxCutoffCol)) then
400 +       allocate(gtypeMaxCutoffCol(jend))
401 +    else
402 +       deallocate(gtypeMaxCutoffCol)      
403 +       allocate(gtypeMaxCutoffCol(jend))
404 +    end if
405 +
406 +       groupMaxCutoffCol = 0.0_dp
407 +       gtypeMaxCutoffCol = 0.0_dp
408 +
409 + #endif
410 +       groupMaxCutoffRow = 0.0_dp
411 +       gtypeMaxCutoffRow = 0.0_dp
412 +
413 +
414 +    !! first we do a single loop over the cutoff groups to find the
415 +    !! largest cutoff for any atypes present in this group.  We also
416 +    !! create gtypes at this point.
417 +    
418 +    tol = 1.0d-6
419 +    nGroupTypesRow = 0
420 +
421 +    do i = istart, iend      
422 +       n_in_i = groupStartRow(i+1) - groupStartRow(i)
423 +       groupMaxCutoffRow(i) = 0.0_dp
424 +       do ia = groupStartRow(i), groupStartRow(i+1)-1
425 +          atom1 = groupListRow(ia)
426 + #ifdef IS_MPI
427 +          me_i = atid_row(atom1)
428 + #else
429 +          me_i = atid(atom1)
430 + #endif          
431 +          if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
432 +             groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
433 +          endif          
434 +       enddo
435 +
436 +       if (nGroupTypesRow.eq.0) then
437 +          nGroupTypesRow = nGroupTypesRow + 1
438 +          gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
439 +          groupToGtypeRow(i) = nGroupTypesRow
440 +       else
441 +          GtypeFound = .false.
442 +          do g = 1, nGroupTypesRow
443 +             if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
444 +                groupToGtypeRow(i) = g
445 +                GtypeFound = .true.
446 +             endif
447 +          enddo
448 +          if (.not.GtypeFound) then            
449 +             nGroupTypesRow = nGroupTypesRow + 1
450 +             gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
451 +             groupToGtypeRow(i) = nGroupTypesRow
452 +          endif
453 +       endif
454 +    enddo    
455 +
456 + #ifdef IS_MPI
457 +    do j = jstart, jend      
458 +       n_in_j = groupStartCol(j+1) - groupStartCol(j)
459 +       groupMaxCutoffCol(j) = 0.0_dp
460 +       do ja = groupStartCol(j), groupStartCol(j+1)-1
461 +          atom1 = groupListCol(ja)
462 +
463 +          me_j = atid_col(atom1)
464 +
465 +          if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
466 +             groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
467 +          endif          
468 +       enddo
469 +
470 +       if (nGroupTypesCol.eq.0) then
471 +          nGroupTypesCol = nGroupTypesCol + 1
472 +          gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
473 +          groupToGtypeCol(j) = nGroupTypesCol
474 +       else
475 +          GtypeFound = .false.
476 +          do g = 1, nGroupTypesCol
477 +             if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
478 +                groupToGtypeCol(j) = g
479 +                GtypeFound = .true.
480 +             endif
481 +          enddo
482 +          if (.not.GtypeFound) then            
483 +             nGroupTypesCol = nGroupTypesCol + 1
484 +             gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
485 +             groupToGtypeCol(j) = nGroupTypesCol
486 +          endif
487 +       endif
488 +    enddo    
489 +
490 + #else
491 + ! Set pointers to information we just found
492 +    nGroupTypesCol = nGroupTypesRow
493 +    groupToGtypeCol => groupToGtypeRow
494 +    gtypeMaxCutoffCol => gtypeMaxCutoffRow
495 +    groupMaxCutoffCol => groupMaxCutoffRow
496 + #endif
497 +
498 +
499 +
500 +
501 +
502 +    !! allocate the gtypeCutoffMap here.
503 +    allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
504 +    !! then we do a double loop over all the group TYPES to find the cutoff
505 +    !! map between groups of two types
506 +    tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
507 +
508 +    do i = 1, nGroupTypesRow
509 +       do j = 1, nGroupTypesCol
510 +      
511 +          select case(cutoffPolicy)
512 +          case(TRADITIONAL_CUTOFF_POLICY)
513 +             thisRcut = tradRcut
514 +          case(MIX_CUTOFF_POLICY)
515 +             thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
516 +          case(MAX_CUTOFF_POLICY)
517 +             thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
518 +          case default
519 +             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
520 +             return
521 +          end select
522 +          gtypeCutoffMap(i,j)%rcut = thisRcut
523 +          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
524 +          skin = defaultRlist - defaultRcut
525 +          listSkin = skin ! set neighbor list skin thickness
526 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
527 +
528 +          ! sanity check
529 +
530 +          if (haveDefaultCutoffs) then
531 +             if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
532 +                call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
533 +             endif
534 +          endif
535 +       enddo
536 +    enddo
537 +    if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
538 +    if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
539 +    if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
540 + #ifdef IS_MPI
541 +    if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
542 +    if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
543 + #endif
544 +    groupMaxCutoffCol => null()
545 +    gtypeMaxCutoffCol => null()
546 +    
547 +    haveGtypeCutoffMap = .true.
548 +   end subroutine createGtypeCutoffMap
549 +
550 +   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
551 +     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
552 +     integer, intent(in) :: cutPolicy
553 +
554 +     defaultRcut = defRcut
555 +     defaultRsw = defRsw
556 +     defaultRlist = defRlist
557 +     cutoffPolicy = cutPolicy
558 +
559 +     haveDefaultCutoffs = .true.
560 +   end subroutine setDefaultCutoffs
561 +
562 +   subroutine setCutoffPolicy(cutPolicy)
563 +
564 +     integer, intent(in) :: cutPolicy
565 +     cutoffPolicy = cutPolicy
566 +     call createGtypeCutoffMap()
567 +   end subroutine setCutoffPolicy
568 +    
569 +    
570    subroutine setSimVariables()
571      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
342 –    SIM_uses_LennardJones = SimUsesLennardJones()
343 –    SIM_uses_Electrostatics = SimUsesElectrostatics()
344 –    SIM_uses_Charges = SimUsesCharges()
345 –    SIM_uses_Dipoles = SimUsesDipoles()
346 –    SIM_uses_Sticky = SimUsesSticky()
347 –    SIM_uses_StickyPower = SimUsesStickyPower()
348 –    SIM_uses_GayBerne = SimUsesGayBerne()
572      SIM_uses_EAM = SimUsesEAM()
350 –    SIM_uses_Shapes = SimUsesShapes()
351 –    SIM_uses_FLARB = SimUsesFLARB()
352 –    SIM_uses_RF = SimUsesRF()
573      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
574      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
575      SIM_uses_PBC = SimUsesPBC()
# Line 366 | Line 586 | contains
586  
587      error = 0
588  
589 <    if (.not. haveInteractionMap) then
589 >    if (.not. haveInteractionHash) then      
590 >       myStatus = 0      
591 >       call createInteractionHash(myStatus)      
592 >       if (myStatus .ne. 0) then
593 >          write(default_error, *) 'createInteractionHash failed in doForces!'
594 >          error = -1
595 >          return
596 >       endif
597 >    endif
598  
599 <       myStatus = 0
600 <
601 <       call createInteractionMap(myStatus)
374 <
599 >    if (.not. haveGtypeCutoffMap) then        
600 >       myStatus = 0      
601 >       call createGtypeCutoffMap(myStatus)      
602         if (myStatus .ne. 0) then
603 <          write(default_error, *) 'createInteractionMap failed in doForces!'
603 >          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
604            error = -1
605            return
606         endif
# Line 383 | Line 610 | contains
610         call setSimVariables()
611      endif
612  
613 <    if (.not. haveRlist) then
614 <       write(default_error, *) 'rList has not been set in doForces!'
615 <       error = -1
616 <       return
617 <    endif
613 >  !  if (.not. haveRlist) then
614 >  !     write(default_error, *) 'rList has not been set in doForces!'
615 >  !     error = -1
616 >  !     return
617 >  !  endif
618  
619      if (.not. haveNeighborList) then
620         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 412 | Line 639 | contains
639    end subroutine doReadyCheck
640  
641  
642 <  subroutine init_FF(use_RF_c, thisStat)
642 >  subroutine init_FF(thisESM, thisStat)
643  
644 <    logical, intent(in) :: use_RF_c
418 <
644 >    integer, intent(in) :: thisESM
645      integer, intent(out) :: thisStat  
646      integer :: my_status, nMatches
647      integer, pointer :: MatchList(:) => null()
# Line 424 | Line 650 | contains
650      !! assume things are copacetic, unless they aren't
651      thisStat = 0
652  
653 <    !! Fortran's version of a cast:
428 <    FF_uses_RF = use_RF_c
653 >    electrostaticSummationMethod = thisESM
654  
655      !! init_FF is called *after* all of the atom types have been
656      !! defined in atype_module using the new_atype subroutine.
# Line 434 | Line 659 | contains
659      !! interactions are used by the force field.    
660  
661      FF_uses_DirectionalAtoms = .false.
437 –    FF_uses_LennardJones = .false.
438 –    FF_uses_Electrostatics = .false.
439 –    FF_uses_Charges = .false.    
662      FF_uses_Dipoles = .false.
441 –    FF_uses_Sticky = .false.
442 –    FF_uses_StickyPower = .false.
663      FF_uses_GayBerne = .false.
664      FF_uses_EAM = .false.
445 –    FF_uses_Shapes = .false.
446 –    FF_uses_FLARB = .false.
665  
666      call getMatchingElementList(atypes, "is_Directional", .true., &
667           nMatches, MatchList)
668      if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
669  
452 –    call getMatchingElementList(atypes, "is_LennardJones", .true., &
453 –         nMatches, MatchList)
454 –    if (nMatches .gt. 0) FF_uses_LennardJones = .true.
455 –
456 –    call getMatchingElementList(atypes, "is_Electrostatic", .true., &
457 –         nMatches, MatchList)
458 –    if (nMatches .gt. 0) then
459 –       FF_uses_Electrostatics = .true.
460 –    endif
461 –
462 –    call getMatchingElementList(atypes, "is_Charge", .true., &
463 –         nMatches, MatchList)
464 –    if (nMatches .gt. 0) then
465 –       FF_uses_Charges = .true.  
466 –       FF_uses_Electrostatics = .true.
467 –    endif
468 –
670      call getMatchingElementList(atypes, "is_Dipole", .true., &
671           nMatches, MatchList)
672 <    if (nMatches .gt. 0) then
472 <       FF_uses_Dipoles = .true.
473 <       FF_uses_Electrostatics = .true.
474 <       FF_uses_DirectionalAtoms = .true.
475 <    endif
476 <
477 <    call getMatchingElementList(atypes, "is_Quadrupole", .true., &
478 <         nMatches, MatchList)
479 <    if (nMatches .gt. 0) then
480 <       FF_uses_Quadrupoles = .true.
481 <       FF_uses_Electrostatics = .true.
482 <       FF_uses_DirectionalAtoms = .true.
483 <    endif
484 <
485 <    call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
486 <         MatchList)
487 <    if (nMatches .gt. 0) then
488 <       FF_uses_Sticky = .true.
489 <       FF_uses_DirectionalAtoms = .true.
490 <    endif
491 <
492 <    call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
493 <         MatchList)
494 <    if (nMatches .gt. 0) then
495 <       FF_uses_StickyPower = .true.
496 <       FF_uses_DirectionalAtoms = .true.
497 <    endif
672 >    if (nMatches .gt. 0) FF_uses_Dipoles = .true.
673      
674      call getMatchingElementList(atypes, "is_GayBerne", .true., &
675           nMatches, MatchList)
676 <    if (nMatches .gt. 0) then
502 <       FF_uses_GayBerne = .true.
503 <       FF_uses_DirectionalAtoms = .true.
504 <    endif
676 >    if (nMatches .gt. 0) FF_uses_GayBerne = .true.
677  
678      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
679      if (nMatches .gt. 0) FF_uses_EAM = .true.
680  
509 –    call getMatchingElementList(atypes, "is_Shape", .true., &
510 –         nMatches, MatchList)
511 –    if (nMatches .gt. 0) then
512 –       FF_uses_Shapes = .true.
513 –       FF_uses_DirectionalAtoms = .true.
514 –    endif
681  
516 –    call getMatchingElementList(atypes, "is_FLARB", .true., &
517 –         nMatches, MatchList)
518 –    if (nMatches .gt. 0) FF_uses_FLARB = .true.
519 –
520 –    !! Assume sanity (for the sake of argument)
682      haveSaneForceField = .true.
683  
684 <    !! check to make sure the FF_uses_RF setting makes sense
684 >    !! check to make sure the reaction field setting makes sense
685  
686 <    if (FF_uses_dipoles) then
687 <       if (FF_uses_RF) then
686 >    if (FF_uses_Dipoles) then
687 >       if (electrostaticSummationMethod == REACTION_FIELD) then
688            dielect = getDielect()
689            call initialize_rf(dielect)
690         endif
691      else
692 <       if (FF_uses_RF) then          
692 >       if (electrostaticSummationMethod == REACTION_FIELD) then
693            write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
694            thisStat = -1
695            haveSaneForceField = .false.
# Line 536 | Line 697 | contains
697         endif
698      endif
699  
539 –    !sticky module does not contain check_sticky_FF anymore
540 –    !if (FF_uses_sticky) then
541 –    !   call check_sticky_FF(my_status)
542 –    !   if (my_status /= 0) then
543 –    !      thisStat = -1
544 –    !      haveSaneForceField = .false.
545 –    !      return
546 –    !   end if
547 –    !endif
548 –
700      if (FF_uses_EAM) then
701         call init_EAM_FF(my_status)
702         if (my_status /= 0) then
# Line 565 | Line 716 | contains
716         endif
717      endif
718  
568 –    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
569 –    endif
570 –
719      if (.not. haveNeighborList) then
720         !! Create neighbor lists
721         call expandNeighborList(nLocal, my_status)
# Line 601 | Line 749 | contains
749  
750      !! Stress Tensor
751      real( kind = dp), dimension(9) :: tau  
752 <    real ( kind = dp ) :: pot
752 >    real ( kind = dp ),dimension(POT_ARRAY_SIZE) :: pot
753      logical ( kind = 2) :: do_pot_c, do_stress_c
754      logical :: do_pot
755      logical :: do_stress
756      logical :: in_switching_region
757   #ifdef IS_MPI
758 <    real( kind = DP ) :: pot_local
758 >    real( kind = DP ), dimension(POT_ARRAY_SIZE) :: pot_local
759      integer :: nAtomsInRow
760      integer :: nAtomsInCol
761      integer :: nprocs
# Line 631 | Line 779 | contains
779      integer :: localError
780      integer :: propPack_i, propPack_j
781      integer :: loopStart, loopEnd, loop
782 <    integer :: iMap
783 <    real(kind=dp) :: listSkin = 1.0  
782 >    integer :: iHash
783 >  
784  
785      !! initialize local variables  
786  
# Line 750 | Line 898 | contains
898               endif
899  
900   #ifdef IS_MPI
901 +             me_j = atid_col(j)
902               call get_interatomic_vector(q_group_Row(:,i), &
903                    q_group_Col(:,j), d_grp, rgrpsq)
904   #else
905 +             me_j = atid(j)
906               call get_interatomic_vector(q_group(:,i), &
907                    q_group(:,j), d_grp, rgrpsq)
908 < #endif
908 > #endif      
909  
910 <             if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
910 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
911                  if (update_nlist) then
912                     nlist = nlist + 1
913  
# Line 878 | Line 1028 | contains
1028                  endif
1029               end if
1030            enddo
1031 +
1032         enddo outer
1033  
1034         if (update_nlist) then
# Line 957 | Line 1108 | contains
1108  
1109      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1110  
1111 <       if (FF_uses_RF .and. SIM_uses_RF) then
1111 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1112  
1113   #ifdef IS_MPI
1114            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 975 | Line 1126 | contains
1126   #else
1127               me_i = atid(i)
1128   #endif
1129 <             iMap = InteractionMap(me_i, me_j)%InteractionHash
1129 >             iHash = InteractionHash(me_i,me_j)
1130              
1131 <             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1131 >             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1132  
1133                  mu_i = getDipoleMoment(me_i)
1134  
# Line 988 | Line 1139 | contains
1139                  !! potential and torques:
1140                  call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1141   #ifdef IS_MPI
1142 <                pot_local = pot_local + rfpot
1142 >                pot_local(RF_POT) = pot_local(RF_POT) + rfpot
1143   #else
1144 <                pot = pot + rfpot
1144 >                pot(RF_POT) = pot(RF_POT) + rfpot
1145  
1146   #endif
1147               endif
# Line 1002 | Line 1153 | contains
1153   #ifdef IS_MPI
1154  
1155      if (do_pot) then
1156 <       pot = pot + pot_local
1156 >       pot(1:SIZE_POT_ARRAY) = pot(1:SIZE_POT_ARRAY) &
1157 >            + pot_local(1:SIZE_POT_ARRAY)
1158         !! we assume the c code will do the allreduce to get the total potential
1159         !! we could do it right here if we needed to...
1160      endif
# Line 1028 | Line 1180 | contains
1180    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1181         eFrame, A, f, t, pot, vpair, fpair)
1182  
1183 <    real( kind = dp ) :: pot, vpair, sw
1183 >    real( kind = dp ) :: vpair, sw
1184 >    real( kind = dp ), dimension(POT_ARRAY_SIZE) :: pot
1185      real( kind = dp ), dimension(3) :: fpair
1186      real( kind = dp ), dimension(nLocal)   :: mfact
1187      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 1041 | Line 1194 | contains
1194      real ( kind = dp ), intent(inout) :: rijsq
1195      real ( kind = dp )                :: r
1196      real ( kind = dp ), intent(inout) :: d(3)
1044 –    real ( kind = dp ) :: ebalance
1197      integer :: me_i, me_j
1198  
1199 <    integer :: iMap
1199 >    integer :: iHash
1200  
1201      r = sqrt(rijsq)
1202      vpair = 0.0d0
# Line 1058 | Line 1210 | contains
1210      me_j = atid(j)
1211   #endif
1212  
1213 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1213 >    iHash = InteractionHash(me_i, me_j)
1214  
1215 <    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1216 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1215 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1216 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot(LJ_POT), f, do_pot)
1217      endif
1218  
1219 <    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1219 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1220         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1221 <            pot, eFrame, f, t, do_pot)
1221 >            pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1222  
1223 <       if (FF_uses_RF .and. SIM_uses_RF) then
1223 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1224  
1225            ! CHECK ME (RF needs to know about all electrostatic types)
1226            call accumulate_rf(i, j, r, eFrame, sw)
# Line 1077 | Line 1229 | contains
1229  
1230      endif
1231  
1232 <    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1232 >    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1233         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1234 <            pot, A, f, t, do_pot)
1234 >            pot(STICKY_POT), A, f, t, do_pot)
1235      endif
1236  
1237 <    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1237 >    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1238         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1239 <            pot, A, f, t, do_pot)
1239 >            pot(STICKYPOWER_POT), A, f, t, do_pot)
1240      endif
1241  
1242 <    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1242 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1243         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1244 <            pot, A, f, t, do_pot)
1244 >            pot(GAYBERNE_POT), A, f, t, do_pot)
1245      endif
1246      
1247 <    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1247 >    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1248   !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1249 < !           pot, A, f, t, do_pot)
1249 > !           pot(GAYBERNE_LJ_POT), A, f, t, do_pot)
1250      endif
1251  
1252 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1253 <       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1252 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1253 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot(EAM_POT), f, &
1254              do_pot)
1255      endif
1256  
1257 <    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1257 >    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1258         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1259 <            pot, A, f, t, do_pot)
1259 >            pot(SHAPE_POT), A, f, t, do_pot)
1260      endif
1261  
1262 <    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1262 >    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1263         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1264 <            pot, A, f, t, do_pot)
1264 >            pot(SHAPE_LJ_POT), A, f, t, do_pot)
1265      endif
1266      
1267    end subroutine do_pair
# Line 1117 | Line 1269 | contains
1269    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1270         do_pot, do_stress, eFrame, A, f, t, pot)
1271  
1272 <    real( kind = dp ) :: pot, sw
1272 >    real( kind = dp ) :: sw
1273 >    real( kind = dp ), dimension(POT_ARRAY_SIZE) :: pot
1274      real( kind = dp ), dimension(9,nLocal) :: eFrame
1275      real (kind=dp), dimension(9,nLocal) :: A
1276      real (kind=dp), dimension(3,nLocal) :: f
# Line 1129 | Line 1282 | contains
1282      real ( kind = dp )                :: r, rc
1283      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1284  
1285 <    integer :: me_i, me_j, iMap
1285 >    integer :: me_i, me_j, iHash
1286  
1287 +    r = sqrt(rijsq)
1288 +
1289   #ifdef IS_MPI  
1290      me_i = atid_row(i)
1291      me_j = atid_col(j)  
# Line 1139 | Line 1294 | contains
1294      me_j = atid(j)  
1295   #endif
1296  
1297 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1297 >    iHash = InteractionHash(me_i, me_j)
1298  
1299 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1299 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1300              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1301      endif
1302      
# Line 1150 | Line 1305 | contains
1305  
1306    subroutine do_preforce(nlocal,pot)
1307      integer :: nlocal
1308 <    real( kind = dp ) :: pot
1308 >    real( kind = dp ),dimension(POT_ARRAY_SIZE) :: pot
1309  
1310      if (FF_uses_EAM .and. SIM_uses_EAM) then
1311 <       call calc_EAM_preforce_Frho(nlocal,pot)
1311 >       call calc_EAM_preforce_Frho(nlocal,pot(EAM_POT))
1312      endif
1313  
1314  
# Line 1338 | Line 1493 | contains
1493  
1494    function FF_UsesDirectionalAtoms() result(doesit)
1495      logical :: doesit
1496 <    doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1342 <         FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1343 <         FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1496 >    doesit = FF_uses_DirectionalAtoms
1497    end function FF_UsesDirectionalAtoms
1498  
1499    function FF_RequiresPrepairCalc() result(doesit)
# Line 1350 | Line 1503 | contains
1503  
1504    function FF_RequiresPostpairCalc() result(doesit)
1505      logical :: doesit
1506 <    doesit = FF_uses_RF
1506 >    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1507    end function FF_RequiresPostpairCalc
1508  
1509   #ifdef PROFILE

Diff Legend

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