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

Comparing trunk/OOPSE-3.0/src/UseTheForce/doForces.F90 (file contents):
Revision 2262 by chuckv, Sun Jul 3 20:53:43 2005 UTC vs.
Revision 2350 by chuckv, Mon Oct 10 21:20:46 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.51 2005-10-10 21:20:46 chuckv Exp $, $Date: 2005-10-10 21:20:46 $, $Name: not supported by cvs2svn $, $Revision: 1.51 $
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      endif
189 +
190 +    if (.not. allocated(atypeMaxCutoff)) then
191 +       allocate(atypeMaxCutoff(nAtypes))
192 +    endif
193          
194      do i = 1, nAtypes
195         call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
# Line 240 | Line 242 | contains
242            if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
243  
244  
245 <          InteractionMap(i,j)%InteractionHash = iHash
246 <          InteractionMap(j,i)%InteractionHash = iHash
245 >          InteractionHash(i,j) = iHash
246 >          InteractionHash(j,i) = iHash
247  
248         end do
249  
250      end do
249 –  end subroutine createInteractionMap
251  
252 < ! 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.
253 <  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
252 >    haveInteractionHash = .true.
253 >  end subroutine createInteractionHash
254  
255 <    if(.not. allocated(InteractionMap)) return
255 >  subroutine createGtypeCutoffMap(stat)
256  
257 <    nAtypes = getSize(atypes)
258 < ! If we pass a default rcut, set all atypes to that cutoff distance
259 <    if(present(defaultRList)) then
260 <       InteractionMap(:,:)%rList = defaultRList
261 <       InteractionMap(:,:)%rListSq = defaultRList*defaultRList
262 <       haveRlist = .true.
263 <       return
264 <    end if
257 >    integer, intent(out), optional :: stat
258 >    logical :: i_is_LJ
259 >    logical :: i_is_Elect
260 >    logical :: i_is_Sticky
261 >    logical :: i_is_StickyP
262 >    logical :: i_is_GB
263 >    logical :: i_is_EAM
264 >    logical :: i_is_Shape
265 >    logical :: GtypeFound
266  
267 <    do map_i = 1,nAtypes
268 <       do map_j = map_i,nAtypes
269 <          iMap = InteractionMap(map_i, map_j)%InteractionHash
270 <          
271 <          if ( iand(iMap, LJ_PAIR).ne.0 ) then
272 < !            thisRCut = getLJCutOff(map_i,map_j)
273 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
274 <          endif
275 <          
276 <          if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
277 < !            thisRCut = getElectrostaticCutOff(map_i,map_j)
278 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
267 >    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
268 >    integer :: n_in_i, me_i, ia, g, atom1
269 >    integer :: nGroupsInRow
270 >    integer :: nGroupsInCol
271 >    integer :: nGroupTypesRow,nGroupTypesCol
272 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin
273 >    real(kind=dp) :: biggestAtypeCutoff
274 >
275 >    stat = 0
276 >    if (.not. haveInteractionHash) then
277 >       call createInteractionHash(myStatus)      
278 >       if (myStatus .ne. 0) then
279 >          write(default_error, *) 'createInteractionHash failed in doForces!'
280 >          stat = -1
281 >          return
282 >       endif
283 >    endif
284 > #ifdef IS_MPI
285 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
286 >    nGroupsInCol = getNgroupsInCol(plan_group_col)
287 > #endif
288 >    nAtypes = getSize(atypes)
289 > ! Set all of the initial cutoffs to zero.
290 >    atypeMaxCutoff = 0.0_dp
291 >    do i = 1, nAtypes
292 >       if (SimHasAtype(i)) then    
293 >          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
294 >          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
295 >          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
296 >          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
297 >          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
298 >          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
299 >          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
300 >          
301 >
302 >          if (haveDefaultCutoffs) then
303 >             atypeMaxCutoff(i) = defaultRcut
304 >          else
305 >             if (i_is_LJ) then          
306 >                thisRcut = getSigma(i) * 2.5_dp
307 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
308 >             endif
309 >             if (i_is_Elect) then
310 >                thisRcut = defaultRcut
311 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
312 >             endif
313 >             if (i_is_Sticky) then
314 >                thisRcut = getStickyCut(i)
315 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
316 >             endif
317 >             if (i_is_StickyP) then
318 >                thisRcut = getStickyPowerCut(i)
319 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
320 >             endif
321 >             if (i_is_GB) then
322 >                thisRcut = getGayBerneCut(i)
323 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
324 >             endif
325 >             if (i_is_EAM) then
326 >                thisRcut = getEAMCut(i)
327 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
328 >             endif
329 >             if (i_is_Shape) then
330 >                thisRcut = getShapeCut(i)
331 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
332 >             endif
333            endif
334            
335 <          if ( iand(iMap, STICKY_PAIR).ne.0 ) then
336 < !             thisRCut = getStickyCutOff(map_i,map_j)
337 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
338 <           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
335 >          
336 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
337 >             biggestAtypeCutoff = atypeMaxCutoff(i)
338 >          endif
339  
340 +       endif
341 +    enddo
342 +  
343  
344 < !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
345 < !!$  subroutine setRlistDF( this_rlist )
346 < !!$
347 < !!$   real(kind=dp) :: this_rlist
348 < !!$
349 < !!$    rlist = this_rlist
350 < !!$    rlistsq = rlist * rlist
351 < !!$
352 < !!$    haveRlist = .true.
353 < !!$
354 < !!$  end subroutine setRlistDF
344 >    
345 >    istart = 1
346 >    jstart = 1
347 > #ifdef IS_MPI
348 >    iend = nGroupsInRow
349 >    jend = nGroupsInCol
350 > #else
351 >    iend = nGroups
352 >    jend = nGroups
353 > #endif
354 >    
355 >    !! allocate the groupToGtype and gtypeMaxCutoff here.
356 >    if(.not.allocated(groupToGtypeRow)) then
357 >     !  allocate(groupToGtype(iend))
358 >       allocate(groupToGtypeRow(iend))
359 >    else
360 >       deallocate(groupToGtypeRow)
361 >       allocate(groupToGtypeRow(iend))
362 >    endif
363 >    if(.not.allocated(groupMaxCutoffRow)) then
364 >       allocate(groupMaxCutoffRow(iend))
365 >    else
366 >       deallocate(groupMaxCutoffRow)
367 >       allocate(groupMaxCutoffRow(iend))
368 >    end if
369 >
370 >    if(.not.allocated(gtypeMaxCutoffRow)) then
371 >       allocate(gtypeMaxCutoffRow(iend))
372 >    else
373 >       deallocate(gtypeMaxCutoffRow)
374 >       allocate(gtypeMaxCutoffRow(iend))
375 >    endif
376 >
377 >
378 > #ifdef IS_MPI
379 >       ! We only allocate new storage if we are in MPI because Ncol /= Nrow
380 >    if(.not.allocated(groupToGtypeCol)) then
381 >       allocate(groupToGtypeCol(jend))
382 >    else
383 >       deallocate(groupToGtypeCol)
384 >       allocate(groupToGtypeCol(jend))
385 >    end if
386 >
387 >    if(.not.allocated(groupToGtypeCol)) then
388 >       allocate(groupToGtypeCol(jend))
389 >    else
390 >       deallocate(groupToGtypeCol)
391 >       allocate(groupToGtypeCol(jend))
392 >    end if
393 >    if(.not.allocated(gtypeMaxCutoffCol)) then
394 >       allocate(gtypeMaxCutoffCol(jend))
395 >    else
396 >       deallocate(gtypeMaxCutoffCol)      
397 >       allocate(gtypeMaxCutoffCol(jend))
398 >    end if
399 >
400 >       groupMaxCutoffCol = 0.0_dp
401 >       gtypeMaxCutoffCol = 0.0_dp
402 >
403 > #endif
404 >       groupMaxCutoffRow = 0.0_dp
405 >       gtypeMaxCutoffRow = 0.0_dp
406 >
407 >
408 >    !! first we do a single loop over the cutoff groups to find the
409 >    !! largest cutoff for any atypes present in this group.  We also
410 >    !! create gtypes at this point.
411 >    
412 >    tol = 1.0d-6
413 >    nGroupTypesRow = 0
414 >
415 >    do i = istart, iend      
416 >       n_in_i = groupStartRow(i+1) - groupStartRow(i)
417 >       groupMaxCutoffRow(i) = 0.0_dp
418 >       do ia = groupStartRow(i), groupStartRow(i+1)-1
419 >          atom1 = groupListRow(ia)
420 > #ifdef IS_MPI
421 >          me_i = atid_row(atom1)
422 > #else
423 >          me_i = atid(atom1)
424 > #endif          
425 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
426 >             groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
427 >          endif          
428 >       enddo
429  
430 +       if (nGroupTypesRow.eq.0) then
431 +          nGroupTypesRow = nGroupTypesRow + 1
432 +          gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
433 +          groupToGtypeRow(i) = nGroupTypesRow
434 +       else
435 +          GtypeFound = .false.
436 +          do g = 1, nGroupTypesRow
437 +             if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
438 +                groupToGtypeRow(i) = g
439 +                GtypeFound = .true.
440 +             endif
441 +          enddo
442 +          if (.not.GtypeFound) then            
443 +             nGroupTypesRow = nGroupTypesRow + 1
444 +             gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
445 +             groupToGtypeRow(i) = nGroupTypesRow
446 +          endif
447 +       endif
448 +    enddo    
449  
450 + #ifdef IS_MPI
451 +    do j = jstart, jend      
452 +       n_in_j = groupStartCol(j+1) - groupStartCol(j)
453 +       groupMaxCutoffCol(j) = 0.0_dp
454 +       do ja = groupStartCol(j), groupStartCol(j+1)-1
455 +          atom1 = groupListCol(ja)
456 +
457 +          me_j = atid_col(atom1)
458 +
459 +          if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
460 +             groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
461 +          endif          
462 +       enddo
463 +
464 +       if (nGroupTypesCol.eq.0) then
465 +          nGroupTypesCol = nGroupTypesCol + 1
466 +          gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
467 +          groupToGtypeCol(j) = nGroupTypesCol
468 +       else
469 +          GtypeFound = .false.
470 +          do g = 1, nGroupTypesCol
471 +             if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
472 +                groupToGtypeCol(j) = g
473 +                GtypeFound = .true.
474 +             endif
475 +          enddo
476 +          if (.not.GtypeFound) then            
477 +             nGroupTypesCol = nGroupTypesCol + 1
478 +             gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
479 +             groupToGtypeCol(j) = nGroupTypesCol
480 +          endif
481 +       endif
482 +    enddo    
483 +
484 + #else
485 + ! Set pointers to information we just found
486 +    nGroupTypesCol = nGroupTypesRow
487 +    groupToGtypeCol => groupToGtypeRow
488 +    gtypeMaxCutoffCol => gtypeMaxCutoffRow
489 +    groupMaxCutoffCol => groupMaxCutoffRow
490 + #endif
491 +
492 +
493 +
494 +
495 +
496 +    !! allocate the gtypeCutoffMap here.
497 +    allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
498 +    !! then we do a double loop over all the group TYPES to find the cutoff
499 +    !! map between groups of two types
500 +    tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
501 +
502 +    do i = 1, nGroupTypesRow
503 +       do j = 1, nGroupTypesCol
504 +      
505 +          select case(cutoffPolicy)
506 +          case(TRADITIONAL_CUTOFF_POLICY)
507 +             thisRcut = tradRcut
508 +          case(MIX_CUTOFF_POLICY)
509 +             thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
510 +          case(MAX_CUTOFF_POLICY)
511 +             thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
512 +          case default
513 +             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
514 +             return
515 +          end select
516 +          gtypeCutoffMap(i,j)%rcut = thisRcut
517 +          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
518 +          skin = defaultRlist - defaultRcut
519 +          listSkin = skin ! set neighbor list skin thickness
520 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
521 +
522 +          ! sanity check
523 +
524 +          if (haveDefaultCutoffs) then
525 +             if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
526 +                call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
527 +             endif
528 +          endif
529 +       enddo
530 +    enddo
531 +    if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
532 +    if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
533 +    if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
534 + #ifdef IS_MPI
535 +    if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
536 +    if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
537 + #endif
538 +    groupMaxCutoffCol => null()
539 +    gtypeMaxCutoffCol => null()
540 +    
541 +    haveGtypeCutoffMap = .true.
542 +   end subroutine createGtypeCutoffMap
543 +
544 +   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
545 +     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
546 +     integer, intent(in) :: cutPolicy
547 +
548 +     defaultRcut = defRcut
549 +     defaultRsw = defRsw
550 +     defaultRlist = defRlist
551 +     cutoffPolicy = cutPolicy
552 +
553 +     haveDefaultCutoffs = .true.
554 +   end subroutine setDefaultCutoffs
555 +
556 +   subroutine setCutoffPolicy(cutPolicy)
557 +
558 +     integer, intent(in) :: cutPolicy
559 +     cutoffPolicy = cutPolicy
560 +     call createGtypeCutoffMap()
561 +   end subroutine setCutoffPolicy
562 +    
563 +    
564    subroutine setSimVariables()
565      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()
566      SIM_uses_EAM = SimUsesEAM()
350 –    SIM_uses_Shapes = SimUsesShapes()
351 –    SIM_uses_FLARB = SimUsesFLARB()
352 –    SIM_uses_RF = SimUsesRF()
567      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
568      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
569      SIM_uses_PBC = SimUsesPBC()
# Line 366 | Line 580 | contains
580  
581      error = 0
582  
583 <    if (.not. haveInteractionMap) then
583 >    if (.not. haveInteractionHash) then      
584 >       myStatus = 0      
585 >       call createInteractionHash(myStatus)      
586 >       if (myStatus .ne. 0) then
587 >          write(default_error, *) 'createInteractionHash failed in doForces!'
588 >          error = -1
589 >          return
590 >       endif
591 >    endif
592  
593 <       myStatus = 0
594 <
595 <       call createInteractionMap(myStatus)
374 <
593 >    if (.not. haveGtypeCutoffMap) then        
594 >       myStatus = 0      
595 >       call createGtypeCutoffMap(myStatus)      
596         if (myStatus .ne. 0) then
597 <          write(default_error, *) 'createInteractionMap failed in doForces!'
597 >          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
598            error = -1
599            return
600         endif
# Line 383 | Line 604 | contains
604         call setSimVariables()
605      endif
606  
607 <    if (.not. haveRlist) then
608 <       write(default_error, *) 'rList has not been set in doForces!'
609 <       error = -1
610 <       return
611 <    endif
607 >  !  if (.not. haveRlist) then
608 >  !     write(default_error, *) 'rList has not been set in doForces!'
609 >  !     error = -1
610 >  !     return
611 >  !  endif
612  
613      if (.not. haveNeighborList) then
614         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 412 | Line 633 | contains
633    end subroutine doReadyCheck
634  
635  
636 <  subroutine init_FF(use_RF_c, thisStat)
636 >  subroutine init_FF(thisESM, thisStat)
637  
638 <    logical, intent(in) :: use_RF_c
418 <
638 >    integer, intent(in) :: thisESM
639      integer, intent(out) :: thisStat  
640      integer :: my_status, nMatches
641      integer, pointer :: MatchList(:) => null()
# Line 424 | Line 644 | contains
644      !! assume things are copacetic, unless they aren't
645      thisStat = 0
646  
647 <    !! Fortran's version of a cast:
428 <    FF_uses_RF = use_RF_c
647 >    electrostaticSummationMethod = thisESM
648  
649      !! init_FF is called *after* all of the atom types have been
650      !! defined in atype_module using the new_atype subroutine.
# Line 434 | Line 653 | contains
653      !! interactions are used by the force field.    
654  
655      FF_uses_DirectionalAtoms = .false.
437 –    FF_uses_LennardJones = .false.
438 –    FF_uses_Electrostatics = .false.
439 –    FF_uses_Charges = .false.    
656      FF_uses_Dipoles = .false.
441 –    FF_uses_Sticky = .false.
442 –    FF_uses_StickyPower = .false.
657      FF_uses_GayBerne = .false.
658      FF_uses_EAM = .false.
445 –    FF_uses_Shapes = .false.
446 –    FF_uses_FLARB = .false.
659  
660      call getMatchingElementList(atypes, "is_Directional", .true., &
661           nMatches, MatchList)
662      if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
663  
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 –
664      call getMatchingElementList(atypes, "is_Dipole", .true., &
665           nMatches, MatchList)
666 <    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
666 >    if (nMatches .gt. 0) FF_uses_Dipoles = .true.
667      
668      call getMatchingElementList(atypes, "is_GayBerne", .true., &
669           nMatches, MatchList)
670 <    if (nMatches .gt. 0) then
502 <       FF_uses_GayBerne = .true.
503 <       FF_uses_DirectionalAtoms = .true.
504 <    endif
670 >    if (nMatches .gt. 0) FF_uses_GayBerne = .true.
671  
672      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
673      if (nMatches .gt. 0) FF_uses_EAM = .true.
674  
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
675  
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)
676      haveSaneForceField = .true.
677  
678 <    !! check to make sure the FF_uses_RF setting makes sense
678 >    !! check to make sure the reaction field setting makes sense
679  
680 <    if (FF_uses_dipoles) then
681 <       if (FF_uses_RF) then
680 >    if (FF_uses_Dipoles) then
681 >       if (electrostaticSummationMethod == REACTION_FIELD) then
682            dielect = getDielect()
683            call initialize_rf(dielect)
684         endif
685      else
686 <       if (FF_uses_RF) then          
686 >       if (electrostaticSummationMethod == REACTION_FIELD) then
687            write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
688            thisStat = -1
689            haveSaneForceField = .false.
# Line 536 | Line 691 | contains
691         endif
692      endif
693  
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 –
694      if (FF_uses_EAM) then
695         call init_EAM_FF(my_status)
696         if (my_status /= 0) then
# Line 565 | Line 710 | contains
710         endif
711      endif
712  
568 –    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
569 –    endif
570 –
713      if (.not. haveNeighborList) then
714         !! Create neighbor lists
715         call expandNeighborList(nLocal, my_status)
# Line 631 | Line 773 | contains
773      integer :: localError
774      integer :: propPack_i, propPack_j
775      integer :: loopStart, loopEnd, loop
776 <    integer :: iMap
777 <    real(kind=dp) :: listSkin = 1.0  
776 >    integer :: iHash
777 >  
778  
779      !! initialize local variables  
780  
# Line 750 | Line 892 | contains
892               endif
893  
894   #ifdef IS_MPI
895 +             me_j = atid_col(j)
896               call get_interatomic_vector(q_group_Row(:,i), &
897                    q_group_Col(:,j), d_grp, rgrpsq)
898   #else
899 +             me_j = atid(j)
900               call get_interatomic_vector(q_group(:,i), &
901                    q_group(:,j), d_grp, rgrpsq)
902 < #endif
902 > #endif      
903  
904 <             if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
904 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
905                  if (update_nlist) then
906                     nlist = nlist + 1
907  
# Line 878 | Line 1022 | contains
1022                  endif
1023               end if
1024            enddo
1025 +
1026         enddo outer
1027  
1028         if (update_nlist) then
# Line 957 | Line 1102 | contains
1102  
1103      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1104  
1105 <       if (FF_uses_RF .and. SIM_uses_RF) then
1105 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1106  
1107   #ifdef IS_MPI
1108            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 975 | Line 1120 | contains
1120   #else
1121               me_i = atid(i)
1122   #endif
1123 <             iMap = InteractionMap(me_i, me_j)%InteractionHash
1123 >             iHash = InteractionHash(me_i,me_j)
1124              
1125 <             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1125 >             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1126  
1127                  mu_i = getDipoleMoment(me_i)
1128  
# Line 1041 | Line 1186 | contains
1186      real ( kind = dp ), intent(inout) :: rijsq
1187      real ( kind = dp )                :: r
1188      real ( kind = dp ), intent(inout) :: d(3)
1044 –    real ( kind = dp ) :: ebalance
1189      integer :: me_i, me_j
1190  
1191 <    integer :: iMap
1191 >    integer :: iHash
1192  
1193      r = sqrt(rijsq)
1194      vpair = 0.0d0
# Line 1058 | Line 1202 | contains
1202      me_j = atid(j)
1203   #endif
1204  
1205 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1205 >    iHash = InteractionHash(me_i, me_j)
1206  
1207 <    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1207 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1208         call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1209      endif
1210  
1211 <    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1211 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1212         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1213              pot, eFrame, f, t, do_pot)
1214  
1215 <       if (FF_uses_RF .and. SIM_uses_RF) then
1215 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1216  
1217            ! CHECK ME (RF needs to know about all electrostatic types)
1218            call accumulate_rf(i, j, r, eFrame, sw)
# Line 1077 | Line 1221 | contains
1221  
1222      endif
1223  
1224 <    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1224 >    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1225         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1226              pot, A, f, t, do_pot)
1227      endif
1228  
1229 <    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1229 >    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1230         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1231              pot, A, f, t, do_pot)
1232      endif
1233  
1234 <    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1234 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1235         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1236              pot, A, f, t, do_pot)
1237      endif
1238      
1239 <    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1239 >    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1240   !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1241   !           pot, A, f, t, do_pot)
1242      endif
1243  
1244 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1244 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1245         call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1246              do_pot)
1247      endif
1248  
1249 <    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1249 >    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1250         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1251              pot, A, f, t, do_pot)
1252      endif
1253  
1254 <    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1254 >    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1255         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1256              pot, A, f, t, do_pot)
1257      endif
# Line 1129 | Line 1273 | contains
1273      real ( kind = dp )                :: r, rc
1274      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1275  
1276 <    integer :: me_i, me_j, iMap
1276 >    integer :: me_i, me_j, iHash
1277  
1278 +    r = sqrt(rijsq)
1279 +
1280   #ifdef IS_MPI  
1281      me_i = atid_row(i)
1282      me_j = atid_col(j)  
# Line 1139 | Line 1285 | contains
1285      me_j = atid(j)  
1286   #endif
1287  
1288 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1288 >    iHash = InteractionHash(me_i, me_j)
1289  
1290 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1290 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1291              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1292      endif
1293      
# Line 1338 | Line 1484 | contains
1484  
1485    function FF_UsesDirectionalAtoms() result(doesit)
1486      logical :: doesit
1487 <    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
1487 >    doesit = FF_uses_DirectionalAtoms
1488    end function FF_UsesDirectionalAtoms
1489  
1490    function FF_RequiresPrepairCalc() result(doesit)
# Line 1350 | Line 1494 | contains
1494  
1495    function FF_RequiresPostpairCalc() result(doesit)
1496      logical :: doesit
1497 <    doesit = FF_uses_RF
1497 >    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1498    end function FF_RequiresPostpairCalc
1499  
1500   #ifdef PROFILE

Diff Legend

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