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

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2267 by tim, Fri Jul 29 17:34:06 2005 UTC vs.
Revision 2306 by chrisfen, Fri Sep 16 20:37:05 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.25 2005-07-29 17:34:06 tim Exp $, $Date: 2005-07-29 17:34:06 $, $Name: not supported by cvs2svn $, $Revision: 1.25 $
48 > !! @version $Id: doForces.F90,v 1.45 2005-09-16 20:36:55 chrisfen Exp $, $Date: 2005-09-16 20:36:55 $, $Name: not supported by cvs2svn $, $Revision: 1.45 $
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  
79 +
80    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
81    INTEGER, PARAMETER:: PAIR_LOOP    = 2
82  
81  logical, save :: haveRlist = .false.
83    logical, save :: haveNeighborList = .false.
84    logical, save :: haveSIMvariables = .false.
85    logical, save :: haveSaneForceField = .false.
86 <  logical, save :: haveInteractionMap = .false.
86 >  logical, save :: haveInteractionHash = .false.
87 >  logical, save :: haveGtypeCutoffMap = .false.
88 >  logical, save :: haveRlist = .false.
89  
90    logical, save :: FF_uses_DirectionalAtoms
88  logical, save :: FF_uses_LennardJones
89  logical, save :: FF_uses_Electrostatics
90  logical, save :: FF_uses_Charges
91    logical, save :: FF_uses_Dipoles
92  logical, save :: FF_uses_Quadrupoles
93  logical, save :: FF_uses_Sticky
94  logical, save :: FF_uses_StickyPower
92    logical, save :: FF_uses_GayBerne
93    logical, save :: FF_uses_EAM
97  logical, save :: FF_uses_Shapes
98  logical, save :: FF_uses_FLARB
99  logical, save :: FF_uses_RF
94  
95    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
96    logical, save :: SIM_uses_EAM
111  logical, save :: SIM_uses_Shapes
112  logical, save :: SIM_uses_FLARB
113  logical, save :: SIM_uses_RF
97    logical, save :: SIM_requires_postpair_calc
98    logical, save :: SIM_requires_prepair_calc
99    logical, save :: SIM_uses_PBC
117  logical, save :: SIM_uses_molecular_cutoffs
100  
101 <  !!!GO AWAY---------
120 <  !!!!!real(kind=dp), save :: rlist, rlistsq
101 >  integer, save :: electrostaticSummationMethod
102  
103    public :: init_FF
104 +  public :: setDefaultCutoffs
105    public :: do_force_loop
106 < !  public :: setRlistDF
107 <  !public :: addInteraction
108 <  !public :: setInteractionHash
109 <  !public :: getInteractionHash
110 <  public :: createInteractionMap
111 <  public :: createRcuts
106 >  public :: createInteractionHash
107 >  public :: createGtypeCutoffMap
108 >  public :: getStickyCut
109 >  public :: getStickyPowerCut
110 >  public :: getGayBerneCut
111 >  public :: getEAMCut
112 >  public :: getShapeCut
113  
114   #ifdef PROFILE
115    public :: getforcetime
# Line 134 | Line 117 | module doForces
117    real :: forceTimeInitial, forceTimeFinal
118    integer :: nLoops
119   #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
120    
121 <  type(Interaction), dimension(:,:),allocatable :: InteractionMap
122 <  
121 >  !! Variables for cutoff mapping and interaction mapping
122 >  ! Bit hash to determine pair-pair interactions.
123 >  integer, dimension(:,:), allocatable :: InteractionHash
124 >  real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
125 >  real(kind=dp), dimension(:), allocatable :: groupMaxCutoff
126 >  integer, dimension(:), allocatable :: groupToGtype
127 >  real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff
128 >  type ::gtypeCutoffs
129 >     real(kind=dp) :: rcut
130 >     real(kind=dp) :: rcutsq
131 >     real(kind=dp) :: rlistsq
132 >  end type gtypeCutoffs
133 >  type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
134  
135 +  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
136 +  real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
137 +  real(kind=dp),save :: rcuti
138    
139   contains
140  
141 <
151 <  subroutine createInteractionMap(status)
141 >  subroutine createInteractionHash(status)
142      integer :: nAtypes
143      integer, intent(out) :: status
144      integer :: i
145      integer :: j
146 <    integer :: ihash
157 <    real(kind=dp) :: myRcut
146 >    integer :: iHash
147      !! Test Types
148      logical :: i_is_LJ
149      logical :: i_is_Elect
# Line 170 | Line 159 | contains
159      logical :: j_is_GB
160      logical :: j_is_EAM
161      logical :: j_is_Shape
162 <    
163 <    status = 0
164 <    
162 >    real(kind=dp) :: myRcut
163 >
164 >    status = 0  
165 >
166      if (.not. associated(atypes)) then
167 <       call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
167 >       call handleError("atype", "atypes was not present before call of createInteractionHash!")
168         status = -1
169         return
170      endif
# Line 186 | Line 176 | contains
176         return
177      end if
178  
179 <    if (.not. allocated(InteractionMap)) then
180 <       allocate(InteractionMap(nAtypes,nAtypes))
179 >    if (.not. allocated(InteractionHash)) then
180 >       allocate(InteractionHash(nAtypes,nAtypes))
181      endif
182 +
183 +    if (.not. allocated(atypeMaxCutoff)) then
184 +       allocate(atypeMaxCutoff(nAtypes))
185 +    endif
186          
187      do i = 1, nAtypes
188         call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
# Line 241 | Line 235 | contains
235            if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
236  
237  
238 <          InteractionMap(i,j)%InteractionHash = iHash
239 <          InteractionMap(j,i)%InteractionHash = iHash
238 >          InteractionHash(i,j) = iHash
239 >          InteractionHash(j,i) = iHash
240  
241         end do
242  
243      end do
244  
245 <    haveInteractionMap = .true.
246 <  end subroutine createInteractionMap
245 >    haveInteractionHash = .true.
246 >  end subroutine createInteractionHash
247  
248 < ! Query each potential and return the cutoff for that potential. We build the neighbor list based on the largest cutoff value for that atype. Each potential can decide whether to calculate the force for that atype based upon it's own cutoff.
255 <  subroutine createRcuts(defaultRList,stat)
256 <    real(kind=dp), intent(in), optional :: defaultRList
257 <    integer :: iMap
258 <    integer :: map_i,map_j
259 <    real(kind=dp) :: thisRCut = 0.0_dp
260 <    real(kind=dp) :: actualCutoff = 0.0_dp
261 <    integer, intent(out) :: stat
262 <    integer :: nAtypes
263 <    integer :: myStatus
248 >  subroutine createGtypeCutoffMap(stat)
249  
250 <    stat = 0
251 <    if (.not. haveInteractionMap) then
250 >    integer, intent(out), optional :: stat
251 >    logical :: i_is_LJ
252 >    logical :: i_is_Elect
253 >    logical :: i_is_Sticky
254 >    logical :: i_is_StickyP
255 >    logical :: i_is_GB
256 >    logical :: i_is_EAM
257 >    logical :: i_is_Shape
258 >    logical :: GtypeFound
259  
260 <       call createInteractionMap(myStatus)
260 >    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
261 >    integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
262 >    integer :: nGroupsInRow
263 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
264 >    real(kind=dp) :: biggestAtypeCutoff
265  
266 +    stat = 0
267 +    if (.not. haveInteractionHash) then
268 +       call createInteractionHash(myStatus)      
269         if (myStatus .ne. 0) then
270 <          write(default_error, *) 'createInteractionMap failed in doForces!'
270 >          write(default_error, *) 'createInteractionHash failed in doForces!'
271            stat = -1
272            return
273         endif
274      endif
275 <
276 <
275 > #ifdef IS_MPI
276 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
277 > #endif
278      nAtypes = getSize(atypes)
279 <    !! If we pass a default rcut, set all atypes to that cutoff distance
280 <    if(present(defaultRList)) then
281 <       InteractionMap(:,:)%rList = defaultRList
282 <       InteractionMap(:,:)%rListSq = defaultRList*defaultRList
283 <       haveRlist = .true.
284 <       return
285 <    end if
286 <
287 <    do map_i = 1,nAtypes
288 <       do map_j = map_i,nAtypes
289 <          iMap = InteractionMap(map_i, map_j)%InteractionHash
279 > ! Set all of the initial cutoffs to zero.
280 >    atypeMaxCutoff = 0.0_dp
281 >    do i = 1, nAtypes
282 >       if (SimHasAtype(i)) then    
283 >          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
284 >          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
285 >          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
286 >          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
287 >          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
288 >          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
289 >          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
290            
291 <          if ( iand(iMap, LJ_PAIR).ne.0 ) then
292 <             ! thisRCut = getLJCutOff(map_i,map_j)
293 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
291 >
292 >          if (i_is_LJ) then
293 >             thisRcut = getSigma(i) * 2.5_dp
294 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
295            endif
296 <          
297 <          if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
298 <             ! thisRCut = getElectrostaticCutOff(map_i,map_j)
298 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
296 >          if (i_is_Elect) then
297 >             thisRcut = defaultRcut
298 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
299            endif
300 +          if (i_is_Sticky) then
301 +             thisRcut = getStickyCut(i)
302 +             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
303 +          endif
304 +          if (i_is_StickyP) then
305 +             thisRcut = getStickyPowerCut(i)
306 +             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
307 +          endif
308 +          if (i_is_GB) then
309 +             thisRcut = getGayBerneCut(i)
310 +             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
311 +          endif
312 +          if (i_is_EAM) then
313 +             thisRcut = getEAMCut(i)
314 +             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
315 +          endif
316 +          if (i_is_Shape) then
317 +             thisRcut = getShapeCut(i)
318 +             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
319 +          endif
320            
321 <          if ( iand(iMap, STICKY_PAIR).ne.0 ) then
322 <             ! thisRCut = getStickyCutOff(map_i,map_j)
323 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
324 <           endif
325 <          
326 <           if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
327 <              ! thisRCut = getStickyPowerCutOff(map_i,map_j)
328 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
329 <           endif
330 <          
331 <           if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
332 <              ! thisRCut = getGayberneCutOff(map_i,map_j)
333 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
334 <           endif
335 <          
336 <           if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
337 < !              thisRCut = getGaybrneLJCutOff(map_i,map_j)
338 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
339 <           endif
340 <          
341 <           if ( iand(iMap, EAM_PAIR).ne.0 ) then      
342 < !              thisRCut = getEAMCutOff(map_i,map_j)
343 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
344 <           endif
345 <          
346 <           if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
347 < !              thisRCut = getShapeCutOff(map_i,map_j)
348 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
349 <           endif
350 <          
351 <           if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
352 < !              thisRCut = getShapeLJCutOff(map_i,map_j)
353 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
354 <           endif
355 <           InteractionMap(map_i, map_j)%rList = actualCutoff
356 <           InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
357 <        end do
358 <     end do
359 <     haveRlist = .true.
360 <  end subroutine createRcuts
361 <
362 <
363 < !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
344 < !!$  subroutine setRlistDF( this_rlist )
345 < !!$
346 < !!$   real(kind=dp) :: this_rlist
347 < !!$
348 < !!$    rlist = this_rlist
349 < !!$    rlistsq = rlist * rlist
350 < !!$
351 < !!$    haveRlist = .true.
352 < !!$
353 < !!$  end subroutine setRlistDF
321 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
322 >             biggestAtypeCutoff = atypeMaxCutoff(i)
323 >          endif
324 >       endif
325 >    enddo
326 >  
327 >    nGroupTypes = 0
328 >    
329 >    istart = 1
330 > #ifdef IS_MPI
331 >    iend = nGroupsInRow
332 > #else
333 >    iend = nGroups
334 > #endif
335 >    
336 >    !! allocate the groupToGtype and gtypeMaxCutoff here.
337 >    if(.not.allocated(groupToGtype)) then
338 >       allocate(groupToGtype(iend))
339 >       allocate(groupMaxCutoff(iend))
340 >       allocate(gtypeMaxCutoff(iend))
341 >       groupMaxCutoff = 0.0_dp
342 >       gtypeMaxCutoff = 0.0_dp
343 >    endif
344 >    !! first we do a single loop over the cutoff groups to find the
345 >    !! largest cutoff for any atypes present in this group.  We also
346 >    !! create gtypes at this point.
347 >    
348 >    tol = 1.0d-6
349 >    
350 >    do i = istart, iend      
351 >       n_in_i = groupStartRow(i+1) - groupStartRow(i)
352 >       groupMaxCutoff(i) = 0.0_dp
353 >       do ia = groupStartRow(i), groupStartRow(i+1)-1
354 >          atom1 = groupListRow(ia)
355 > #ifdef IS_MPI
356 >          me_i = atid_row(atom1)
357 > #else
358 >          me_i = atid(atom1)
359 > #endif          
360 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
361 >             groupMaxCutoff(i)=atypeMaxCutoff(me_i)
362 >          endif          
363 >       enddo
364  
365 +       if (nGroupTypes.eq.0) then
366 +          nGroupTypes = nGroupTypes + 1
367 +          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
368 +          groupToGtype(i) = nGroupTypes
369 +       else
370 +          GtypeFound = .false.
371 +          do g = 1, nGroupTypes
372 +             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
373 +                groupToGtype(i) = g
374 +                GtypeFound = .true.
375 +             endif
376 +          enddo
377 +          if (.not.GtypeFound) then            
378 +             nGroupTypes = nGroupTypes + 1
379 +             gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
380 +             groupToGtype(i) = nGroupTypes
381 +          endif
382 +       endif
383 +    enddo    
384  
385 +    !! allocate the gtypeCutoffMap here.
386 +    allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
387 +    !! then we do a double loop over all the group TYPES to find the cutoff
388 +    !! map between groups of two types
389 +    
390 +    do i = 1, nGroupTypes
391 +       do j = 1, nGroupTypes
392 +      
393 +          select case(cutoffPolicy)
394 +          case(TRADITIONAL_CUTOFF_POLICY)
395 +             thisRcut = maxval(gtypeMaxCutoff)
396 +          case(MIX_CUTOFF_POLICY)
397 +             thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
398 +          case(MAX_CUTOFF_POLICY)
399 +             thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
400 +          case default
401 +             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
402 +             return
403 +          end select
404 +          gtypeCutoffMap(i,j)%rcut = thisRcut
405 +          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
406 +          skin = defaultRlist - defaultRcut
407 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
408 +
409 +       enddo
410 +    enddo
411 +    
412 +    haveGtypeCutoffMap = .true.
413 +   end subroutine createGtypeCutoffMap
414 +
415 +   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
416 +     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
417 +     integer, intent(in) :: cutPolicy
418 +
419 +     defaultRcut = defRcut
420 +     defaultRsw = defRsw
421 +     defaultRlist = defRlist
422 +     cutoffPolicy = cutPolicy
423 +     rcuti = 1.0_dp / defaultRcut
424 +   end subroutine setDefaultCutoffs
425 +
426 +   subroutine setCutoffPolicy(cutPolicy)
427 +
428 +     integer, intent(in) :: cutPolicy
429 +     cutoffPolicy = cutPolicy
430 +     call createGtypeCutoffMap()
431 +   end subroutine setCutoffPolicy
432 +    
433 +    
434    subroutine setSimVariables()
435      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
358    SIM_uses_LennardJones = SimUsesLennardJones()
359    SIM_uses_Electrostatics = SimUsesElectrostatics()
360    SIM_uses_Charges = SimUsesCharges()
361    SIM_uses_Dipoles = SimUsesDipoles()
362    SIM_uses_Sticky = SimUsesSticky()
363    SIM_uses_StickyPower = SimUsesStickyPower()
364    SIM_uses_GayBerne = SimUsesGayBerne()
436      SIM_uses_EAM = SimUsesEAM()
366    SIM_uses_Shapes = SimUsesShapes()
367    SIM_uses_FLARB = SimUsesFLARB()
368    SIM_uses_RF = SimUsesRF()
437      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
438      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
439      SIM_uses_PBC = SimUsesPBC()
# Line 382 | Line 450 | contains
450  
451      error = 0
452  
453 <    if (.not. haveInteractionMap) then
453 >    if (.not. haveInteractionHash) then      
454 >       myStatus = 0      
455 >       call createInteractionHash(myStatus)      
456 >       if (myStatus .ne. 0) then
457 >          write(default_error, *) 'createInteractionHash failed in doForces!'
458 >          error = -1
459 >          return
460 >       endif
461 >    endif
462  
463 <       myStatus = 0
464 <
465 <       call createInteractionMap(myStatus)
390 <
463 >    if (.not. haveGtypeCutoffMap) then        
464 >       myStatus = 0      
465 >       call createGtypeCutoffMap(myStatus)      
466         if (myStatus .ne. 0) then
467 <          write(default_error, *) 'createInteractionMap failed in doForces!'
467 >          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
468            error = -1
469            return
470         endif
# Line 399 | Line 474 | contains
474         call setSimVariables()
475      endif
476  
477 <    if (.not. haveRlist) then
478 <       write(default_error, *) 'rList has not been set in doForces!'
479 <       error = -1
480 <       return
481 <    endif
477 >  !  if (.not. haveRlist) then
478 >  !     write(default_error, *) 'rList has not been set in doForces!'
479 >  !     error = -1
480 >  !     return
481 >  !  endif
482  
483      if (.not. haveNeighborList) then
484         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 428 | Line 503 | contains
503    end subroutine doReadyCheck
504  
505  
506 <  subroutine init_FF(use_RF_c, thisStat)
506 >  subroutine init_FF(thisESM, dampingAlpha, thisStat)
507  
508 <    logical, intent(in) :: use_RF_c
509 <
508 >    integer, intent(in) :: thisESM
509 >    real(kind=dp), intent(in) :: dampingAlpha
510      integer, intent(out) :: thisStat  
511      integer :: my_status, nMatches
512      integer, pointer :: MatchList(:) => null()
# Line 440 | Line 515 | contains
515      !! assume things are copacetic, unless they aren't
516      thisStat = 0
517  
518 <    !! Fortran's version of a cast:
444 <    FF_uses_RF = use_RF_c
518 >    electrostaticSummationMethod = thisESM
519  
520      !! init_FF is called *after* all of the atom types have been
521      !! defined in atype_module using the new_atype subroutine.
# Line 450 | Line 524 | contains
524      !! interactions are used by the force field.    
525  
526      FF_uses_DirectionalAtoms = .false.
453    FF_uses_LennardJones = .false.
454    FF_uses_Electrostatics = .false.
455    FF_uses_Charges = .false.    
527      FF_uses_Dipoles = .false.
457    FF_uses_Sticky = .false.
458    FF_uses_StickyPower = .false.
528      FF_uses_GayBerne = .false.
529      FF_uses_EAM = .false.
461    FF_uses_Shapes = .false.
462    FF_uses_FLARB = .false.
530  
531      call getMatchingElementList(atypes, "is_Directional", .true., &
532           nMatches, MatchList)
533      if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
534  
468    call getMatchingElementList(atypes, "is_LennardJones", .true., &
469         nMatches, MatchList)
470    if (nMatches .gt. 0) FF_uses_LennardJones = .true.
471
472    call getMatchingElementList(atypes, "is_Electrostatic", .true., &
473         nMatches, MatchList)
474    if (nMatches .gt. 0) then
475       FF_uses_Electrostatics = .true.
476    endif
477
478    call getMatchingElementList(atypes, "is_Charge", .true., &
479         nMatches, MatchList)
480    if (nMatches .gt. 0) then
481       FF_uses_Charges = .true.  
482       FF_uses_Electrostatics = .true.
483    endif
484
535      call getMatchingElementList(atypes, "is_Dipole", .true., &
536           nMatches, MatchList)
537 <    if (nMatches .gt. 0) then
488 <       FF_uses_Dipoles = .true.
489 <       FF_uses_Electrostatics = .true.
490 <       FF_uses_DirectionalAtoms = .true.
491 <    endif
492 <
493 <    call getMatchingElementList(atypes, "is_Quadrupole", .true., &
494 <         nMatches, MatchList)
495 <    if (nMatches .gt. 0) then
496 <       FF_uses_Quadrupoles = .true.
497 <       FF_uses_Electrostatics = .true.
498 <       FF_uses_DirectionalAtoms = .true.
499 <    endif
500 <
501 <    call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
502 <         MatchList)
503 <    if (nMatches .gt. 0) then
504 <       FF_uses_Sticky = .true.
505 <       FF_uses_DirectionalAtoms = .true.
506 <    endif
507 <
508 <    call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
509 <         MatchList)
510 <    if (nMatches .gt. 0) then
511 <       FF_uses_StickyPower = .true.
512 <       FF_uses_DirectionalAtoms = .true.
513 <    endif
537 >    if (nMatches .gt. 0) FF_uses_Dipoles = .true.
538      
539      call getMatchingElementList(atypes, "is_GayBerne", .true., &
540           nMatches, MatchList)
541 <    if (nMatches .gt. 0) then
518 <       FF_uses_GayBerne = .true.
519 <       FF_uses_DirectionalAtoms = .true.
520 <    endif
541 >    if (nMatches .gt. 0) FF_uses_GayBerne = .true.
542  
543      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
544      if (nMatches .gt. 0) FF_uses_EAM = .true.
545  
525    call getMatchingElementList(atypes, "is_Shape", .true., &
526         nMatches, MatchList)
527    if (nMatches .gt. 0) then
528       FF_uses_Shapes = .true.
529       FF_uses_DirectionalAtoms = .true.
530    endif
546  
532    call getMatchingElementList(atypes, "is_FLARB", .true., &
533         nMatches, MatchList)
534    if (nMatches .gt. 0) FF_uses_FLARB = .true.
535
536    !! Assume sanity (for the sake of argument)
547      haveSaneForceField = .true.
548  
549 <    !! check to make sure the FF_uses_RF setting makes sense
549 >    !! check to make sure the reaction field setting makes sense
550  
551 <    if (FF_uses_dipoles) then
552 <       if (FF_uses_RF) then
551 >    if (FF_uses_Dipoles) then
552 >       if (electrostaticSummationMethod == 3) then
553            dielect = getDielect()
554            call initialize_rf(dielect)
555         endif
556      else
557 <       if (FF_uses_RF) then          
557 >       if (electrostaticSummationMethod == 3) then
558            write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
559            thisStat = -1
560            haveSaneForceField = .false.
# Line 552 | Line 562 | contains
562         endif
563      endif
564  
555    !sticky module does not contain check_sticky_FF anymore
556    !if (FF_uses_sticky) then
557    !   call check_sticky_FF(my_status)
558    !   if (my_status /= 0) then
559    !      thisStat = -1
560    !      haveSaneForceField = .false.
561    !      return
562    !   end if
563    !endif
564
565      if (FF_uses_EAM) then
566         call init_EAM_FF(my_status)
567         if (my_status /= 0) then
# Line 581 | Line 581 | contains
581         endif
582      endif
583  
584    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
585    endif
586
584      if (.not. haveNeighborList) then
585         !! Create neighbor lists
586         call expandNeighborList(nLocal, my_status)
# Line 647 | Line 644 | contains
644      integer :: localError
645      integer :: propPack_i, propPack_j
646      integer :: loopStart, loopEnd, loop
647 <    integer :: iMap
647 >    integer :: iHash
648      real(kind=dp) :: listSkin = 1.0  
649  
650      !! initialize local variables  
# Line 739 | Line 736 | contains
736   #endif
737         outer: do i = istart, iend
738  
742 #ifdef IS_MPI
743             me_i = atid_row(i)
744 #else
745             me_i = atid(i)
746 #endif
747
739            if (update_nlist) point(i) = nlist + 1
740  
741            n_in_i = groupStartRow(i+1) - groupStartRow(i)
# Line 781 | Line 772 | contains
772                    q_group(:,j), d_grp, rgrpsq)
773   #endif
774  
775 <             if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
775 >             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
776                  if (update_nlist) then
777                     nlist = nlist + 1
778  
# Line 981 | Line 972 | contains
972  
973      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
974  
975 <       if (FF_uses_RF .and. SIM_uses_RF) then
975 >       if (electrostaticSummationMethod == 3) then
976  
977   #ifdef IS_MPI
978            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 999 | Line 990 | contains
990   #else
991               me_i = atid(i)
992   #endif
993 <             iMap = InteractionMap(me_i, me_j)%InteractionHash
993 >             iHash = InteractionHash(me_i,me_j)
994              
995 <             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
995 >             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
996  
997                  mu_i = getDipoleMoment(me_i)
998  
# Line 1065 | Line 1056 | contains
1056      real ( kind = dp ), intent(inout) :: rijsq
1057      real ( kind = dp )                :: r
1058      real ( kind = dp ), intent(inout) :: d(3)
1068    real ( kind = dp ) :: ebalance
1059      integer :: me_i, me_j
1060  
1061 <    integer :: iMap
1061 >    integer :: iHash
1062  
1063      r = sqrt(rijsq)
1064      vpair = 0.0d0
# Line 1082 | Line 1072 | contains
1072      me_j = atid(j)
1073   #endif
1074  
1075 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1075 >    iHash = InteractionHash(me_i, me_j)
1076  
1077 <    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1077 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1078         call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1079      endif
1080  
1081 <    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1081 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1082         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1083              pot, eFrame, f, t, do_pot)
1084  
1085 <       if (FF_uses_RF .and. SIM_uses_RF) then
1085 >       if (electrostaticSummationMethod == 3) then
1086  
1087            ! CHECK ME (RF needs to know about all electrostatic types)
1088            call accumulate_rf(i, j, r, eFrame, sw)
# Line 1101 | Line 1091 | contains
1091  
1092      endif
1093  
1094 <    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1094 >    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1095         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1096              pot, A, f, t, do_pot)
1097      endif
1098  
1099 <    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1099 >    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1100         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1101              pot, A, f, t, do_pot)
1102      endif
1103  
1104 <    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1104 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1105         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1106              pot, A, f, t, do_pot)
1107      endif
1108      
1109 <    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1109 >    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1110   !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1111   !           pot, A, f, t, do_pot)
1112      endif
1113  
1114 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1114 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1115         call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1116              do_pot)
1117      endif
1118  
1119 <    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1119 >    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1120         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1121              pot, A, f, t, do_pot)
1122      endif
1123  
1124 <    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1124 >    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1125         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1126              pot, A, f, t, do_pot)
1127      endif
# Line 1153 | Line 1143 | contains
1143      real ( kind = dp )                :: r, rc
1144      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1145  
1146 <    integer :: me_i, me_j, iMap
1146 >    integer :: me_i, me_j, iHash
1147  
1148 +    r = sqrt(rijsq)
1149 +
1150   #ifdef IS_MPI  
1151      me_i = atid_row(i)
1152      me_j = atid_col(j)  
# Line 1163 | Line 1155 | contains
1155      me_j = atid(j)  
1156   #endif
1157  
1158 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1158 >    iHash = InteractionHash(me_i, me_j)
1159  
1160 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1160 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1161              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1162      endif
1163      
# Line 1362 | Line 1354 | contains
1354  
1355    function FF_UsesDirectionalAtoms() result(doesit)
1356      logical :: doesit
1357 <    doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1366 <         FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1367 <         FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1357 >    doesit = FF_uses_DirectionalAtoms
1358    end function FF_UsesDirectionalAtoms
1359  
1360    function FF_RequiresPrepairCalc() result(doesit)
# Line 1374 | Line 1364 | contains
1364  
1365    function FF_RequiresPostpairCalc() result(doesit)
1366      logical :: doesit
1367 <    doesit = FF_uses_RF
1367 >    if (electrostaticSummationMethod == 3) doesit = .true.
1368    end function FF_RequiresPostpairCalc
1369  
1370   #ifdef PROFILE

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines