ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/UseTheForce/doForces.F90
(Generate patch)

Comparing trunk/OOPSE-2.0/src/UseTheForce/doForces.F90 (file contents):
Revision 2266 by chuckv, Thu Jul 28 22:12:45 2005 UTC vs.
Revision 2354 by chuckv, Tue Oct 11 22:00:30 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.24 2005-07-28 22:12:45 chuckv Exp $, $Date: 2005-07-28 22:12:45 $, $Name: not supported by cvs2svn $, $Revision: 1.24 $
48 > !! @version $Id: doForces.F90,v 1.53 2005-10-11 22:00:30 chuckv Exp $, $Date: 2005-10-11 22:00:30 $, $Name: not supported by cvs2svn $, $Revision: 1.53 $
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, 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 <    status = 0
171 <    
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 186 | 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 241 | 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
250 –  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,stat)
254 <    real(kind=dp), intent(in), optional :: defaultRList
255 <    integer :: iMap
256 <    integer :: map_i,map_j
257 <    real(kind=dp) :: thisRCut = 0.0_dp
258 <    real(kind=dp) :: actualCutoff = 0.0_dp
259 <    integer, intent(out) :: stat
260 <    integer :: nAtypes
261 <    integer :: myStatus
258 >    haveInteractionHash = .true.
259 >  end subroutine createInteractionHash
260  
261 <    stat = 0
264 <    if (.not. haveInteractionMap) then
261 >  subroutine createGtypeCutoffMap(stat)
262  
263 <       call createInteractionMap(myStatus)
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 +    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, *) 'createInteractionMap failed in doForces!'
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 +
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 +          
342 +          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
343 +             biggestAtypeCutoff = atypeMaxCutoff(i)
344 +          endif
345 +
346         endif
347 +    enddo
348 +  
349 +
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 <    nAtypes = getSize(atypes)
384 < ! If we pass a default rcut, set all atypes to that cutoff distance
385 <    if(present(defaultRList)) then
386 <       InteractionMap(:,:)%rList = defaultRList
387 <       InteractionMap(:,:)%rListSq = defaultRList*defaultRList
388 <       haveRlist = .true.
389 <       return
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 <    do map_i = 1,nAtypes
394 <       do map_j = map_i,nAtypes
395 <          iMap = InteractionMap(map_i, map_j)%InteractionHash
396 <          
397 <          if ( iand(iMap, LJ_PAIR).ne.0 ) then
398 < !            thisRCut = getLJCutOff(map_i,map_j)
399 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
400 <          endif
401 <          
402 <          if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
403 < !            thisRCut = getElectrostaticCutOff(map_i,map_j)
404 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
297 <          endif
298 <          
299 <          if ( iand(iMap, STICKY_PAIR).ne.0 ) then
300 < !             thisRCut = getStickyCutOff(map_i,map_j)
301 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
302 <           endif
303 <          
304 <           if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
305 < !              thisRCut = getStickyPowerCutOff(map_i,map_j)
306 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
307 <           endif
308 <          
309 <           if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
310 < !              thisRCut = getGayberneCutOff(map_i,map_j)
311 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
312 <           endif
313 <          
314 <           if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
315 < !              thisRCut = getGaybrneLJCutOff(map_i,map_j)
316 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
317 <           endif
318 <          
319 <           if ( iand(iMap, EAM_PAIR).ne.0 ) then      
320 < !              thisRCut = getEAMCutOff(map_i,map_j)
321 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
322 <           endif
323 <          
324 <           if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
325 < !              thisRCut = getShapeCutOff(map_i,map_j)
326 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
327 <           endif
328 <          
329 <           if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
330 < !              thisRCut = getShapeLJCutOff(map_i,map_j)
331 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
332 <           endif
333 <           InteractionMap(map_i, map_j)%rList = actualCutoff
334 <           InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
335 <        end do
336 <     end do
337 <          haveRlist = .true.
338 <  end subroutine createRcuts
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 < !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
410 < !!$  subroutine setRlistDF( this_rlist )
411 < !!$
412 < !!$   real(kind=dp) :: this_rlist
413 < !!$
414 < !!$    rlist = this_rlist
415 < !!$    rlistsq = rlist * rlist
416 < !!$
417 < !!$    haveRlist = .true.
418 < !!$
419 < !!$  end subroutine setRlistDF
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()
356 –    SIM_uses_LennardJones = SimUsesLennardJones()
357 –    SIM_uses_Electrostatics = SimUsesElectrostatics()
358 –    SIM_uses_Charges = SimUsesCharges()
359 –    SIM_uses_Dipoles = SimUsesDipoles()
360 –    SIM_uses_Sticky = SimUsesSticky()
361 –    SIM_uses_StickyPower = SimUsesStickyPower()
362 –    SIM_uses_GayBerne = SimUsesGayBerne()
572      SIM_uses_EAM = SimUsesEAM()
364 –    SIM_uses_Shapes = SimUsesShapes()
365 –    SIM_uses_FLARB = SimUsesFLARB()
366 –    SIM_uses_RF = SimUsesRF()
573      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
574      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
575      SIM_uses_PBC = SimUsesPBC()
# Line 380 | 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)
388 <
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 397 | 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 426 | 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
432 <
644 >    integer, intent(in) :: thisESM
645      integer, intent(out) :: thisStat  
646      integer :: my_status, nMatches
647      integer, pointer :: MatchList(:) => null()
# Line 438 | Line 650 | contains
650      !! assume things are copacetic, unless they aren't
651      thisStat = 0
652  
653 <    !! Fortran's version of a cast:
442 <    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 448 | Line 659 | contains
659      !! interactions are used by the force field.    
660  
661      FF_uses_DirectionalAtoms = .false.
451 –    FF_uses_LennardJones = .false.
452 –    FF_uses_Electrostatics = .false.
453 –    FF_uses_Charges = .false.    
662      FF_uses_Dipoles = .false.
455 –    FF_uses_Sticky = .false.
456 –    FF_uses_StickyPower = .false.
663      FF_uses_GayBerne = .false.
664      FF_uses_EAM = .false.
459 –    FF_uses_Shapes = .false.
460 –    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  
466 –    call getMatchingElementList(atypes, "is_LennardJones", .true., &
467 –         nMatches, MatchList)
468 –    if (nMatches .gt. 0) FF_uses_LennardJones = .true.
469 –
470 –    call getMatchingElementList(atypes, "is_Electrostatic", .true., &
471 –         nMatches, MatchList)
472 –    if (nMatches .gt. 0) then
473 –       FF_uses_Electrostatics = .true.
474 –    endif
475 –
476 –    call getMatchingElementList(atypes, "is_Charge", .true., &
477 –         nMatches, MatchList)
478 –    if (nMatches .gt. 0) then
479 –       FF_uses_Charges = .true.  
480 –       FF_uses_Electrostatics = .true.
481 –    endif
482 –
670      call getMatchingElementList(atypes, "is_Dipole", .true., &
671           nMatches, MatchList)
672 <    if (nMatches .gt. 0) then
486 <       FF_uses_Dipoles = .true.
487 <       FF_uses_Electrostatics = .true.
488 <       FF_uses_DirectionalAtoms = .true.
489 <    endif
490 <
491 <    call getMatchingElementList(atypes, "is_Quadrupole", .true., &
492 <         nMatches, MatchList)
493 <    if (nMatches .gt. 0) then
494 <       FF_uses_Quadrupoles = .true.
495 <       FF_uses_Electrostatics = .true.
496 <       FF_uses_DirectionalAtoms = .true.
497 <    endif
498 <
499 <    call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
500 <         MatchList)
501 <    if (nMatches .gt. 0) then
502 <       FF_uses_Sticky = .true.
503 <       FF_uses_DirectionalAtoms = .true.
504 <    endif
505 <
506 <    call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
507 <         MatchList)
508 <    if (nMatches .gt. 0) then
509 <       FF_uses_StickyPower = .true.
510 <       FF_uses_DirectionalAtoms = .true.
511 <    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
516 <       FF_uses_GayBerne = .true.
517 <       FF_uses_DirectionalAtoms = .true.
518 <    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  
523 –    call getMatchingElementList(atypes, "is_Shape", .true., &
524 –         nMatches, MatchList)
525 –    if (nMatches .gt. 0) then
526 –       FF_uses_Shapes = .true.
527 –       FF_uses_DirectionalAtoms = .true.
528 –    endif
681  
530 –    call getMatchingElementList(atypes, "is_FLARB", .true., &
531 –         nMatches, MatchList)
532 –    if (nMatches .gt. 0) FF_uses_FLARB = .true.
533 –
534 –    !! 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 550 | Line 697 | contains
697         endif
698      endif
699  
553 –    !sticky module does not contain check_sticky_FF anymore
554 –    !if (FF_uses_sticky) then
555 –    !   call check_sticky_FF(my_status)
556 –    !   if (my_status /= 0) then
557 –    !      thisStat = -1
558 –    !      haveSaneForceField = .false.
559 –    !      return
560 –    !   end if
561 –    !endif
562 –
700      if (FF_uses_EAM) then
701         call init_EAM_FF(my_status)
702         if (my_status /= 0) then
# Line 579 | Line 716 | contains
716         endif
717      endif
718  
582 –    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
583 –    endif
584 –
719      if (.not. haveNeighborList) then
720         !! Create neighbor lists
721         call expandNeighborList(nLocal, my_status)
# Line 645 | 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 743 | Line 877 | contains
877  
878            if (update_nlist) then
879   #ifdef IS_MPI
746 –             me_i = atid_row(i)
880               jstart = 1
881               jend = nGroupsInCol
882   #else
750 –             me_i = atid(i)
883               jstart = i+1
884               jend = nGroups
885   #endif
# Line 773 | Line 905 | contains
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 896 | Line 1028 | contains
1028                  endif
1029               end if
1030            enddo
1031 +
1032         enddo outer
1033  
1034         if (update_nlist) then
# Line 975 | 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 993 | 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 1059 | Line 1192 | contains
1192      real ( kind = dp ), intent(inout) :: rijsq
1193      real ( kind = dp )                :: r
1194      real ( kind = dp ), intent(inout) :: d(3)
1062 –    real ( kind = dp ) :: ebalance
1195      integer :: me_i, me_j
1196  
1197 <    integer :: iMap
1197 >    integer :: iHash
1198  
1199      r = sqrt(rijsq)
1200      vpair = 0.0d0
# Line 1076 | Line 1208 | contains
1208      me_j = atid(j)
1209   #endif
1210  
1211 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1211 >    iHash = InteractionHash(me_i, me_j)
1212  
1213 <    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1213 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1214         call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1215      endif
1216  
1217 <    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1217 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1218         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1219              pot, eFrame, f, t, do_pot)
1220  
1221 <       if (FF_uses_RF .and. SIM_uses_RF) then
1221 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1222  
1223            ! CHECK ME (RF needs to know about all electrostatic types)
1224            call accumulate_rf(i, j, r, eFrame, sw)
# Line 1095 | Line 1227 | contains
1227  
1228      endif
1229  
1230 <    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1230 >    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1231         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1232              pot, A, f, t, do_pot)
1233      endif
1234  
1235 <    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1235 >    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1236         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1237              pot, A, f, t, do_pot)
1238      endif
1239  
1240 <    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1240 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1241         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1242              pot, A, f, t, do_pot)
1243      endif
1244      
1245 <    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1245 >    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1246   !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1247   !           pot, A, f, t, do_pot)
1248      endif
1249  
1250 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1250 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1251         call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1252              do_pot)
1253      endif
1254  
1255 <    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1255 >    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1256         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1257              pot, A, f, t, do_pot)
1258      endif
1259  
1260 <    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1260 >    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1261         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1262              pot, A, f, t, do_pot)
1263      endif
# Line 1147 | Line 1279 | contains
1279      real ( kind = dp )                :: r, rc
1280      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1281  
1282 <    integer :: me_i, me_j, iMap
1282 >    integer :: me_i, me_j, iHash
1283  
1284 +    r = sqrt(rijsq)
1285 +
1286   #ifdef IS_MPI  
1287      me_i = atid_row(i)
1288      me_j = atid_col(j)  
# Line 1157 | Line 1291 | contains
1291      me_j = atid(j)  
1292   #endif
1293  
1294 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1294 >    iHash = InteractionHash(me_i, me_j)
1295  
1296 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1296 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1297              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1298      endif
1299      
# Line 1356 | Line 1490 | contains
1490  
1491    function FF_UsesDirectionalAtoms() result(doesit)
1492      logical :: doesit
1493 <    doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1360 <         FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1361 <         FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1493 >    doesit = FF_uses_DirectionalAtoms
1494    end function FF_UsesDirectionalAtoms
1495  
1496    function FF_RequiresPrepairCalc() result(doesit)
# Line 1368 | Line 1500 | contains
1500  
1501    function FF_RequiresPostpairCalc() result(doesit)
1502      logical :: doesit
1503 <    doesit = FF_uses_RF
1503 >    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1504    end function FF_RequiresPostpairCalc
1505  
1506   #ifdef PROFILE

Diff Legend

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