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 2301 by gezelter, Thu Sep 15 22:05:21 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.44 2005-09-15 22:05:17 gezelter Exp $, $Date: 2005-09-15 22:05:17 $, $Name: not supported by cvs2svn $, $Revision: 1.44 $
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
94    logical, save :: FF_uses_RF
95  
96    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
97    logical, save :: SIM_uses_EAM
111  logical, save :: SIM_uses_Shapes
112  logical, save :: SIM_uses_FLARB
98    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 :: corrMethod
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 :: groupMaxCutoff
128 >  integer, dimension(:), allocatable :: groupToGtype
129 >  real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff
130 >  type ::gtypeCutoffs
131 >     real(kind=dp) :: rcut
132 >     real(kind=dp) :: rcutsq
133 >     real(kind=dp) :: rlistsq
134 >  end type gtypeCutoffs
135 >  type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
136  
137 +  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
138 +  real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
139 +  real(kind=dp),save :: rcuti
140    
141   contains
142  
143 <
151 <  subroutine createInteractionMap(status)
143 >  subroutine createInteractionHash(status)
144      integer :: nAtypes
145      integer, intent(out) :: status
146      integer :: i
147      integer :: j
148 <    integer :: ihash
157 <    real(kind=dp) :: myRcut
148 >    integer :: iHash
149      !! Test Types
150      logical :: i_is_LJ
151      logical :: i_is_Elect
# Line 170 | Line 161 | contains
161      logical :: j_is_GB
162      logical :: j_is_EAM
163      logical :: j_is_Shape
164 <    
165 <    status = 0
166 <    
164 >    real(kind=dp) :: myRcut
165 >
166 >    status = 0  
167 >
168      if (.not. associated(atypes)) then
169 <       call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
169 >       call handleError("atype", "atypes was not present before call of createInteractionHash!")
170         status = -1
171         return
172      endif
# Line 186 | Line 178 | contains
178         return
179      end if
180  
181 <    if (.not. allocated(InteractionMap)) then
182 <       allocate(InteractionMap(nAtypes,nAtypes))
181 >    if (.not. allocated(InteractionHash)) then
182 >       allocate(InteractionHash(nAtypes,nAtypes))
183      endif
184 +
185 +    if (.not. allocated(atypeMaxCutoff)) then
186 +       allocate(atypeMaxCutoff(nAtypes))
187 +    endif
188          
189      do i = 1, nAtypes
190         call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
# Line 241 | Line 237 | contains
237            if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
238  
239  
240 <          InteractionMap(i,j)%InteractionHash = iHash
241 <          InteractionMap(j,i)%InteractionHash = iHash
240 >          InteractionHash(i,j) = iHash
241 >          InteractionHash(j,i) = iHash
242  
243         end do
244  
245      end do
246  
247 <    haveInteractionMap = .true.
248 <  end subroutine createInteractionMap
247 >    haveInteractionHash = .true.
248 >  end subroutine createInteractionHash
249  
250 < ! 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
250 >  subroutine createGtypeCutoffMap(stat)
251  
252 <    stat = 0
253 <    if (.not. haveInteractionMap) then
252 >    integer, intent(out), optional :: stat
253 >    logical :: i_is_LJ
254 >    logical :: i_is_Elect
255 >    logical :: i_is_Sticky
256 >    logical :: i_is_StickyP
257 >    logical :: i_is_GB
258 >    logical :: i_is_EAM
259 >    logical :: i_is_Shape
260 >    logical :: GtypeFound
261  
262 <       call createInteractionMap(myStatus)
262 >    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
263 >    integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
264 >    integer :: nGroupsInRow
265 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
266 >    real(kind=dp) :: biggestAtypeCutoff
267  
268 +    stat = 0
269 +    if (.not. haveInteractionHash) then
270 +       call createInteractionHash(myStatus)      
271         if (myStatus .ne. 0) then
272 <          write(default_error, *) 'createInteractionMap failed in doForces!'
272 >          write(default_error, *) 'createInteractionHash failed in doForces!'
273            stat = -1
274            return
275         endif
276      endif
277 <
278 <
277 > #ifdef IS_MPI
278 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
279 > #endif
280      nAtypes = getSize(atypes)
281 <    !! If we pass a default rcut, set all atypes to that cutoff distance
282 <    if(present(defaultRList)) then
283 <       InteractionMap(:,:)%rList = defaultRList
284 <       InteractionMap(:,:)%rListSq = defaultRList*defaultRList
285 <       haveRlist = .true.
286 <       return
287 <    end if
288 <
289 <    do map_i = 1,nAtypes
290 <       do map_j = map_i,nAtypes
291 <          iMap = InteractionMap(map_i, map_j)%InteractionHash
281 > ! Set all of the initial cutoffs to zero.
282 >    atypeMaxCutoff = 0.0_dp
283 >    do i = 1, nAtypes
284 >       if (SimHasAtype(i)) then    
285 >          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
286 >          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
287 >          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
288 >          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
289 >          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
290 >          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
291 >          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
292            
293 <          if ( iand(iMap, LJ_PAIR).ne.0 ) then
294 <             ! thisRCut = getLJCutOff(map_i,map_j)
295 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
293 >
294 >          if (i_is_LJ) then
295 >             thisRcut = getSigma(i) * 2.5_dp
296 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
297            endif
298 <          
299 <          if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
300 <             ! thisRCut = getElectrostaticCutOff(map_i,map_j)
298 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
298 >          if (i_is_Elect) then
299 >             thisRcut = defaultRcut
300 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
301            endif
302 +          if (i_is_Sticky) then
303 +             thisRcut = getStickyCut(i)
304 +             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
305 +          endif
306 +          if (i_is_StickyP) then
307 +             thisRcut = getStickyPowerCut(i)
308 +             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
309 +          endif
310 +          if (i_is_GB) then
311 +             thisRcut = getGayBerneCut(i)
312 +             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
313 +          endif
314 +          if (i_is_EAM) then
315 +             thisRcut = getEAMCut(i)
316 +             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
317 +          endif
318 +          if (i_is_Shape) then
319 +             thisRcut = getShapeCut(i)
320 +             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
321 +          endif
322            
323 <          if ( iand(iMap, STICKY_PAIR).ne.0 ) then
324 <             ! thisRCut = getStickyCutOff(map_i,map_j)
325 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
326 <           endif
327 <          
328 <           if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
329 <              ! thisRCut = getStickyPowerCutOff(map_i,map_j)
330 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
331 <           endif
332 <          
333 <           if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
334 <              ! thisRCut = getGayberneCutOff(map_i,map_j)
335 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
336 <           endif
337 <          
338 <           if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
339 < !              thisRCut = getGaybrneLJCutOff(map_i,map_j)
340 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
341 <           endif
342 <          
343 <           if ( iand(iMap, EAM_PAIR).ne.0 ) then      
344 < !              thisRCut = getEAMCutOff(map_i,map_j)
345 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
346 <           endif
347 <          
348 <           if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
349 < !              thisRCut = getShapeCutOff(map_i,map_j)
350 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
351 <           endif
352 <          
353 <           if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
354 < !              thisRCut = getShapeLJCutOff(map_i,map_j)
355 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
356 <           endif
357 <           InteractionMap(map_i, map_j)%rList = actualCutoff
358 <           InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
359 <        end do
360 <     end do
361 <     haveRlist = .true.
362 <  end subroutine createRcuts
363 <
364 <
365 < !!! 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
323 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
324 >             biggestAtypeCutoff = atypeMaxCutoff(i)
325 >          endif
326 >       endif
327 >    enddo
328 >  
329 >    nGroupTypes = 0
330 >    
331 >    istart = 1
332 > #ifdef IS_MPI
333 >    iend = nGroupsInRow
334 > #else
335 >    iend = nGroups
336 > #endif
337 >    
338 >    !! allocate the groupToGtype and gtypeMaxCutoff here.
339 >    if(.not.allocated(groupToGtype)) then
340 >       allocate(groupToGtype(iend))
341 >       allocate(groupMaxCutoff(iend))
342 >       allocate(gtypeMaxCutoff(iend))
343 >       groupMaxCutoff = 0.0_dp
344 >       gtypeMaxCutoff = 0.0_dp
345 >    endif
346 >    !! first we do a single loop over the cutoff groups to find the
347 >    !! largest cutoff for any atypes present in this group.  We also
348 >    !! create gtypes at this point.
349 >    
350 >    tol = 1.0d-6
351 >    
352 >    do i = istart, iend      
353 >       n_in_i = groupStartRow(i+1) - groupStartRow(i)
354 >       groupMaxCutoff(i) = 0.0_dp
355 >       do ia = groupStartRow(i), groupStartRow(i+1)-1
356 >          atom1 = groupListRow(ia)
357 > #ifdef IS_MPI
358 >          me_i = atid_row(atom1)
359 > #else
360 >          me_i = atid(atom1)
361 > #endif          
362 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
363 >             groupMaxCutoff(i)=atypeMaxCutoff(me_i)
364 >          endif          
365 >       enddo
366  
367 +       if (nGroupTypes.eq.0) then
368 +          nGroupTypes = nGroupTypes + 1
369 +          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
370 +          groupToGtype(i) = nGroupTypes
371 +       else
372 +          GtypeFound = .false.
373 +          do g = 1, nGroupTypes
374 +             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
375 +                groupToGtype(i) = g
376 +                GtypeFound = .true.
377 +             endif
378 +          enddo
379 +          if (.not.GtypeFound) then            
380 +             nGroupTypes = nGroupTypes + 1
381 +             gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
382 +             groupToGtype(i) = nGroupTypes
383 +          endif
384 +       endif
385 +    enddo    
386  
387 +    !! allocate the gtypeCutoffMap here.
388 +    allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
389 +    !! then we do a double loop over all the group TYPES to find the cutoff
390 +    !! map between groups of two types
391 +    
392 +    do i = 1, nGroupTypes
393 +       do j = 1, nGroupTypes
394 +      
395 +          select case(cutoffPolicy)
396 +          case(TRADITIONAL_CUTOFF_POLICY)
397 +             thisRcut = maxval(gtypeMaxCutoff)
398 +          case(MIX_CUTOFF_POLICY)
399 +             thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
400 +          case(MAX_CUTOFF_POLICY)
401 +             thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
402 +          case default
403 +             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
404 +             return
405 +          end select
406 +          gtypeCutoffMap(i,j)%rcut = thisRcut
407 +          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
408 +          skin = defaultRlist - defaultRcut
409 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
410 +
411 +       enddo
412 +    enddo
413 +    
414 +    haveGtypeCutoffMap = .true.
415 +   end subroutine createGtypeCutoffMap
416 +
417 +   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
418 +     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
419 +     integer, intent(in) :: cutPolicy
420 +
421 +     defaultRcut = defRcut
422 +     defaultRsw = defRsw
423 +     defaultRlist = defRlist
424 +     cutoffPolicy = cutPolicy
425 +     rcuti = 1.0_dp / defaultRcut
426 +   end subroutine setDefaultCutoffs
427 +
428 +   subroutine setCutoffPolicy(cutPolicy)
429 +
430 +     integer, intent(in) :: cutPolicy
431 +     cutoffPolicy = cutPolicy
432 +     call createGtypeCutoffMap()
433 +   end subroutine setCutoffPolicy
434 +    
435 +    
436    subroutine setSimVariables()
437      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()
438      SIM_uses_EAM = SimUsesEAM()
366    SIM_uses_Shapes = SimUsesShapes()
367    SIM_uses_FLARB = SimUsesFLARB()
368    SIM_uses_RF = SimUsesRF()
439      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
440      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
441      SIM_uses_PBC = SimUsesPBC()
442 +    SIM_uses_RF = SimUsesRF()
443  
444      haveSIMvariables = .true.
445  
# Line 382 | Line 453 | contains
453  
454      error = 0
455  
456 <    if (.not. haveInteractionMap) then
456 >    if (.not. haveInteractionHash) then      
457 >       myStatus = 0      
458 >       call createInteractionHash(myStatus)      
459 >       if (myStatus .ne. 0) then
460 >          write(default_error, *) 'createInteractionHash failed in doForces!'
461 >          error = -1
462 >          return
463 >       endif
464 >    endif
465  
466 <       myStatus = 0
467 <
468 <       call createInteractionMap(myStatus)
390 <
466 >    if (.not. haveGtypeCutoffMap) then        
467 >       myStatus = 0      
468 >       call createGtypeCutoffMap(myStatus)      
469         if (myStatus .ne. 0) then
470 <          write(default_error, *) 'createInteractionMap failed in doForces!'
470 >          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
471            error = -1
472            return
473         endif
# Line 399 | Line 477 | contains
477         call setSimVariables()
478      endif
479  
480 <    if (.not. haveRlist) then
481 <       write(default_error, *) 'rList has not been set in doForces!'
482 <       error = -1
483 <       return
484 <    endif
480 >  !  if (.not. haveRlist) then
481 >  !     write(default_error, *) 'rList has not been set in doForces!'
482 >  !     error = -1
483 >  !     return
484 >  !  endif
485  
486      if (.not. haveNeighborList) then
487         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 428 | Line 506 | contains
506    end subroutine doReadyCheck
507  
508  
509 <  subroutine init_FF(use_RF_c, thisStat)
509 >  subroutine init_FF(use_RF, correctionMethod, dampingAlpha, thisStat)
510  
511 <    logical, intent(in) :: use_RF_c
512 <
511 >    logical, intent(in) :: use_RF
512 >    integer, intent(in) :: correctionMethod
513 >    real(kind=dp), intent(in) :: dampingAlpha
514      integer, intent(out) :: thisStat  
515      integer :: my_status, nMatches
516      integer, pointer :: MatchList(:) => null()
# Line 441 | Line 520 | contains
520      thisStat = 0
521  
522      !! Fortran's version of a cast:
523 <    FF_uses_RF = use_RF_c
523 >    FF_uses_RF = use_RF
524  
525 +        
526      !! init_FF is called *after* all of the atom types have been
527      !! defined in atype_module using the new_atype subroutine.
528      !!
# Line 450 | Line 530 | contains
530      !! interactions are used by the force field.    
531  
532      FF_uses_DirectionalAtoms = .false.
453    FF_uses_LennardJones = .false.
454    FF_uses_Electrostatics = .false.
455    FF_uses_Charges = .false.    
533      FF_uses_Dipoles = .false.
457    FF_uses_Sticky = .false.
458    FF_uses_StickyPower = .false.
534      FF_uses_GayBerne = .false.
535      FF_uses_EAM = .false.
461    FF_uses_Shapes = .false.
462    FF_uses_FLARB = .false.
536  
537      call getMatchingElementList(atypes, "is_Directional", .true., &
538           nMatches, MatchList)
539      if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
467
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
540  
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
541      call getMatchingElementList(atypes, "is_Dipole", .true., &
542           nMatches, MatchList)
543 <    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
543 >    if (nMatches .gt. 0) FF_uses_Dipoles = .true.
544      
545      call getMatchingElementList(atypes, "is_GayBerne", .true., &
546           nMatches, MatchList)
547 <    if (nMatches .gt. 0) then
518 <       FF_uses_GayBerne = .true.
519 <       FF_uses_DirectionalAtoms = .true.
520 <    endif
547 >    if (nMatches .gt. 0) FF_uses_GayBerne = .true.
548  
549      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
550      if (nMatches .gt. 0) FF_uses_EAM = .true.
551  
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
552  
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)
553      haveSaneForceField = .true.
554  
555      !! check to make sure the FF_uses_RF setting makes sense
556  
557 <    if (FF_uses_dipoles) then
557 >    if (FF_uses_Dipoles) then
558         if (FF_uses_RF) then
559            dielect = getDielect()
560            call initialize_rf(dielect)
561         endif
562      else
563 <       if (FF_uses_RF) then          
563 >       if ((corrMethod == 3) .or. FF_uses_RF) then
564            write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
565            thisStat = -1
566            haveSaneForceField = .false.
567            return
568         endif
569      endif
554
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
570  
571      if (FF_uses_EAM) then
572         call init_EAM_FF(my_status)
# Line 581 | Line 587 | contains
587         endif
588      endif
589  
584    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
585    endif
586
590      if (.not. haveNeighborList) then
591         !! Create neighbor lists
592         call expandNeighborList(nLocal, my_status)
# Line 647 | Line 650 | contains
650      integer :: localError
651      integer :: propPack_i, propPack_j
652      integer :: loopStart, loopEnd, loop
653 <    integer :: iMap
653 >    integer :: iHash
654      real(kind=dp) :: listSkin = 1.0  
655  
656      !! initialize local variables  
# Line 739 | Line 742 | contains
742   #endif
743         outer: do i = istart, iend
744  
742 #ifdef IS_MPI
743             me_i = atid_row(i)
744 #else
745             me_i = atid(i)
746 #endif
747
745            if (update_nlist) point(i) = nlist + 1
746  
747            n_in_i = groupStartRow(i+1) - groupStartRow(i)
# Line 781 | Line 778 | contains
778                    q_group(:,j), d_grp, rgrpsq)
779   #endif
780  
781 <             if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
781 >             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
782                  if (update_nlist) then
783                     nlist = nlist + 1
784  
# Line 981 | Line 978 | contains
978  
979      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
980  
981 <       if (FF_uses_RF .and. SIM_uses_RF) then
981 >       if ((FF_uses_RF .and. SIM_uses_RF) .or. (corrMethod == 3)) then
982  
983   #ifdef IS_MPI
984            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 999 | Line 996 | contains
996   #else
997               me_i = atid(i)
998   #endif
999 <             iMap = InteractionMap(me_i, me_j)%InteractionHash
999 >             iHash = InteractionHash(me_i,me_j)
1000              
1001 <             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1001 >             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1002  
1003                  mu_i = getDipoleMoment(me_i)
1004  
# Line 1065 | Line 1062 | contains
1062      real ( kind = dp ), intent(inout) :: rijsq
1063      real ( kind = dp )                :: r
1064      real ( kind = dp ), intent(inout) :: d(3)
1068    real ( kind = dp ) :: ebalance
1065      integer :: me_i, me_j
1066  
1067 <    integer :: iMap
1067 >    integer :: iHash
1068  
1069      r = sqrt(rijsq)
1070      vpair = 0.0d0
# Line 1082 | Line 1078 | contains
1078      me_j = atid(j)
1079   #endif
1080  
1081 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1081 >    iHash = InteractionHash(me_i, me_j)
1082  
1083 <    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1083 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1084         call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1085      endif
1086  
1087 <    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1087 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1088         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1089 <            pot, eFrame, f, t, do_pot)
1089 >            pot, eFrame, f, t, do_pot, corrMethod, rcuti)
1090  
1091 <       if (FF_uses_RF .and. SIM_uses_RF) then
1091 >       if ((FF_uses_RF .and. SIM_uses_RF) .or. (corrMethod == 3)) then
1092  
1093            ! CHECK ME (RF needs to know about all electrostatic types)
1094            call accumulate_rf(i, j, r, eFrame, sw)
# Line 1101 | Line 1097 | contains
1097  
1098      endif
1099  
1100 <    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1100 >    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1101         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1102              pot, A, f, t, do_pot)
1103      endif
1104  
1105 <    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1105 >    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1106         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1107              pot, A, f, t, do_pot)
1108      endif
1109  
1110 <    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1110 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1111         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1112              pot, A, f, t, do_pot)
1113      endif
1114      
1115 <    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1115 >    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1116   !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1117   !           pot, A, f, t, do_pot)
1118      endif
1119  
1120 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1120 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1121         call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1122              do_pot)
1123      endif
1124  
1125 <    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1125 >    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1126         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1127              pot, A, f, t, do_pot)
1128      endif
1129  
1130 <    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1130 >    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1131         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1132              pot, A, f, t, do_pot)
1133      endif
# Line 1153 | Line 1149 | contains
1149      real ( kind = dp )                :: r, rc
1150      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1151  
1152 <    integer :: me_i, me_j, iMap
1152 >    integer :: me_i, me_j, iHash
1153  
1154 +    r = sqrt(rijsq)
1155 +
1156   #ifdef IS_MPI  
1157      me_i = atid_row(i)
1158      me_j = atid_col(j)  
# Line 1163 | Line 1161 | contains
1161      me_j = atid(j)  
1162   #endif
1163  
1164 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1164 >    iHash = InteractionHash(me_i, me_j)
1165  
1166 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1166 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1167              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1168      endif
1169      
# Line 1362 | Line 1360 | contains
1360  
1361    function FF_UsesDirectionalAtoms() result(doesit)
1362      logical :: doesit
1363 <    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
1363 >    doesit = FF_uses_DirectionalAtoms
1364    end function FF_UsesDirectionalAtoms
1365  
1366    function FF_RequiresPrepairCalc() result(doesit)
# Line 1375 | Line 1371 | contains
1371    function FF_RequiresPostpairCalc() result(doesit)
1372      logical :: doesit
1373      doesit = FF_uses_RF
1374 +    if (corrMethod == 3) doesit = .true.
1375    end function FF_RequiresPostpairCalc
1376  
1377   #ifdef PROFILE

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines