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

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2262 by chuckv, Sun Jul 3 20:53:43 2005 UTC vs.
Revision 2290 by gezelter, Wed Sep 14 20:28: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.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.41 2005-09-14 20:28:05 gezelter Exp $, $Date: 2005-09-14 20:28:05 $, $Name: not supported by cvs2svn $, $Revision: 1.41 $
49  
50  
51   module doForces
# 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    
140   contains
141  
142 <
151 <  subroutine createInteractionMap(status)
142 >  subroutine createInteractionHash(status)
143      integer :: nAtypes
144 <    integer :: status
144 >    integer, intent(out) :: status
145      integer :: i
146      integer :: j
147 <    integer :: ihash
148 <    real(kind=dp) :: myRcut
158 < ! Test Types
147 >    integer :: iHash
148 >    !! Test Types
149      logical :: i_is_LJ
150      logical :: i_is_Elect
151      logical :: i_is_Sticky
# Line 170 | Line 160 | contains
160      logical :: j_is_GB
161      logical :: j_is_EAM
162      logical :: j_is_Shape
163 <    
164 <    
163 >    real(kind=dp) :: myRcut
164 >
165 >    status = 0  
166 >
167      if (.not. associated(atypes)) then
168 <       call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
168 >       call handleError("atype", "atypes was not present before call of createInteractionHash!")
169         status = -1
170         return
171      endif
# Line 185 | Line 177 | contains
177         return
178      end if
179  
180 <    if (.not. allocated(InteractionMap)) then
181 <       allocate(InteractionMap(nAtypes,nAtypes))
180 >    if (.not. allocated(InteractionHash)) then
181 >       allocate(InteractionHash(nAtypes,nAtypes))
182      endif
183 +
184 +    if (.not. allocated(atypeMaxCutoff)) then
185 +       allocate(atypeMaxCutoff(nAtypes))
186 +    endif
187          
188      do i = 1, nAtypes
189         call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
# Line 240 | Line 236 | contains
236            if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
237  
238  
239 <          InteractionMap(i,j)%InteractionHash = iHash
240 <          InteractionMap(j,i)%InteractionHash = iHash
239 >          InteractionHash(i,j) = iHash
240 >          InteractionHash(j,i) = iHash
241  
242         end do
243  
244      end do
249 –  end subroutine createInteractionMap
245  
246 < ! 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.
247 <  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
246 >    haveInteractionHash = .true.
247 >  end subroutine createInteractionHash
248  
249 <    if(.not. allocated(InteractionMap)) return
249 >  subroutine createGtypeCutoffMap(stat)
250  
251 <    nAtypes = getSize(atypes)
252 < ! If we pass a default rcut, set all atypes to that cutoff distance
253 <    if(present(defaultRList)) then
254 <       InteractionMap(:,:)%rList = defaultRList
255 <       InteractionMap(:,:)%rListSq = defaultRList*defaultRList
256 <       haveRlist = .true.
257 <       return
258 <    end if
251 >    integer, intent(out), optional :: stat
252 >    logical :: i_is_LJ
253 >    logical :: i_is_Elect
254 >    logical :: i_is_Sticky
255 >    logical :: i_is_StickyP
256 >    logical :: i_is_GB
257 >    logical :: i_is_EAM
258 >    logical :: i_is_Shape
259 >    logical :: GtypeFound
260  
261 <    do map_i = 1,nAtypes
262 <       do map_j = map_i,nAtypes
263 <          iMap = InteractionMap(map_i, map_j)%InteractionHash
264 <          
265 <          if ( iand(iMap, LJ_PAIR).ne.0 ) then
266 < !            thisRCut = getLJCutOff(map_i,map_j)
267 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
268 <          endif
269 <          
270 <          if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
271 < !            thisRCut = getElectrostaticCutOff(map_i,map_j)
272 <             if (thisRcut > actualCutoff) actualCutoff = thisRcut
261 >    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
262 >    integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
263 >    integer :: nGroupsInRow
264 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
265 >    real(kind=dp) :: biggestAtypeCutoff
266 >
267 >    stat = 0
268 >    if (.not. haveInteractionHash) then
269 >       call createInteractionHash(myStatus)      
270 >       if (myStatus .ne. 0) then
271 >          write(default_error, *) 'createInteractionHash failed in doForces!'
272 >          stat = -1
273 >          return
274 >       endif
275 >    endif
276 > #ifdef IS_MPI
277 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
278 > #endif
279 >    nAtypes = getSize(atypes)
280 >    
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 >          atypeMaxCutoff(i) = 0.0_dp
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 +          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 <
326 <
327 < !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
328 < !!$  subroutine setRlistDF( this_rlist )
329 < !!$
330 < !!$   real(kind=dp) :: this_rlist
331 < !!$
332 < !!$    rlist = this_rlist
333 < !!$    rlistsq = rlist * rlist
334 < !!$
335 < !!$    haveRlist = .true.
336 < !!$
337 < !!$  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 >    endif
342 >    !! first we do a single loop over the cutoff groups to find the
343 >    !! largest cutoff for any atypes present in this group.  We also
344 >    !! create gtypes at this point.
345 >    
346 >    tol = 1.0d-6
347 >    
348 >    do i = istart, iend      
349 >       n_in_i = groupStartRow(i+1) - groupStartRow(i)
350 >       groupMaxCutoff(i) = 0.0_dp
351 >       do ia = groupStartRow(i), groupStartRow(i+1)-1
352 >          atom1 = groupListRow(ia)
353 > #ifdef IS_MPI
354 >          me_i = atid_row(atom1)
355 > #else
356 >          me_i = atid(atom1)
357 > #endif          
358 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
359 >             groupMaxCutoff(i)=atypeMaxCutoff(me_i)
360 >          endif          
361 >       enddo
362  
363 +       if (nGroupTypes.eq.0) then
364 +          nGroupTypes = nGroupTypes + 1
365 +          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
366 +          groupToGtype(i) = nGroupTypes
367 +       else
368 +          GtypeFound = .false.
369 +          do g = 1, nGroupTypes
370 +             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
371 +                groupToGtype(i) = g
372 +                GtypeFound = .true.
373 +             endif
374 +          enddo
375 +          if (.not.GtypeFound) then            
376 +             nGroupTypes = nGroupTypes + 1
377 +             gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
378 +             groupToGtype(i) = nGroupTypes
379 +          endif
380 +       endif
381 +    enddo    
382  
383 +    !! allocate the gtypeCutoffMap here.
384 +    allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
385 +    !! then we do a double loop over all the group TYPES to find the cutoff
386 +    !! map between groups of two types
387 +    
388 +    do i = 1, nGroupTypes
389 +       do j = 1, nGroupTypes
390 +      
391 +          select case(cutoffPolicy)
392 +          case(TRADITIONAL_CUTOFF_POLICY)
393 +             thisRcut = maxval(gtypeMaxCutoff)
394 +          case(MIX_CUTOFF_POLICY)
395 +             thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
396 +          case(MAX_CUTOFF_POLICY)
397 +             thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
398 +          case default
399 +             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
400 +             return
401 +          end select
402 +          gtypeCutoffMap(i,j)%rcut = thisRcut
403 +          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
404 +          skin = defaultRlist - defaultRcut
405 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
406 +
407 +       enddo
408 +    enddo
409 +    
410 +    haveGtypeCutoffMap = .true.
411 +    
412 +  end subroutine createGtypeCutoffMap
413 +  
414 +  subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
415 +    real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
416 +    integer, intent(in) :: cutPolicy
417 +    
418 +    defaultRcut = defRcut
419 +    defaultRsw = defRsw
420 +    defaultRlist = defRlist
421 +    cutoffPolicy = cutPolicy
422 +  end subroutine setDefaultCutoffs
423 +  
424 +  subroutine setCutoffPolicy(cutPolicy)
425 +
426 +     integer, intent(in) :: cutPolicy
427 +     cutoffPolicy = cutPolicy
428 +     call createGtypeCutoffMap()
429 +
430 +   end subroutine setCutoffPolicy
431 +    
432 +    
433    subroutine setSimVariables()
434      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()
435      SIM_uses_EAM = SimUsesEAM()
350 –    SIM_uses_Shapes = SimUsesShapes()
351 –    SIM_uses_FLARB = SimUsesFLARB()
436      SIM_uses_RF = SimUsesRF()
437      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
438      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
# Line 366 | 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)
374 <
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 383 | 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 412 | Line 503 | contains
503    end subroutine doReadyCheck
504  
505  
506 <  subroutine init_FF(use_RF_c, thisStat)
506 >  subroutine init_FF(use_RF, use_UW, use_DW, thisStat)
507  
508 <    logical, intent(in) :: use_RF_c
509 <
508 >    logical, intent(in) :: use_RF
509 >    logical, intent(in) :: use_UW
510 >    logical, intent(in) :: use_DW
511      integer, intent(out) :: thisStat  
512      integer :: my_status, nMatches
513 +    integer :: corrMethod
514      integer, pointer :: MatchList(:) => null()
515      real(kind=dp) :: rcut, rrf, rt, dielect
516  
# Line 425 | Line 518 | contains
518      thisStat = 0
519  
520      !! Fortran's version of a cast:
521 <    FF_uses_RF = use_RF_c
521 >    FF_uses_RF = use_RF
522  
523 +    !! set the electrostatic correction method
524 +    if (use_UW) then
525 +       corrMethod = 1
526 +    elseif (use_DW) then
527 +       corrMethod = 2
528 +    else
529 +       corrMethod = 0
530 +    endif
531 +    
532      !! init_FF is called *after* all of the atom types have been
533      !! defined in atype_module using the new_atype subroutine.
534      !!
# Line 434 | Line 536 | contains
536      !! interactions are used by the force field.    
537  
538      FF_uses_DirectionalAtoms = .false.
437 –    FF_uses_LennardJones = .false.
438 –    FF_uses_Electrostatics = .false.
439 –    FF_uses_Charges = .false.    
539      FF_uses_Dipoles = .false.
441 –    FF_uses_Sticky = .false.
442 –    FF_uses_StickyPower = .false.
540      FF_uses_GayBerne = .false.
541      FF_uses_EAM = .false.
445 –    FF_uses_Shapes = .false.
446 –    FF_uses_FLARB = .false.
542  
543      call getMatchingElementList(atypes, "is_Directional", .true., &
544           nMatches, MatchList)
545      if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
546  
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 –
547      call getMatchingElementList(atypes, "is_Dipole", .true., &
548           nMatches, MatchList)
549 <    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
549 >    if (nMatches .gt. 0) FF_uses_Dipoles = .true.
550      
551      call getMatchingElementList(atypes, "is_GayBerne", .true., &
552           nMatches, MatchList)
553 <    if (nMatches .gt. 0) then
502 <       FF_uses_GayBerne = .true.
503 <       FF_uses_DirectionalAtoms = .true.
504 <    endif
553 >    if (nMatches .gt. 0) FF_uses_GayBerne = .true.
554  
555      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
556      if (nMatches .gt. 0) FF_uses_EAM = .true.
557  
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
558  
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)
559      haveSaneForceField = .true.
560  
561      !! check to make sure the FF_uses_RF setting makes sense
562  
563 <    if (FF_uses_dipoles) then
563 >    if (FF_uses_Dipoles) then
564         if (FF_uses_RF) then
565            dielect = getDielect()
566            call initialize_rf(dielect)
# Line 536 | Line 574 | contains
574         endif
575      endif
576  
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 –
577      if (FF_uses_EAM) then
578         call init_EAM_FF(my_status)
579         if (my_status /= 0) then
# Line 565 | Line 593 | contains
593         endif
594      endif
595  
568 –    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
569 –    endif
570 –
596      if (.not. haveNeighborList) then
597         !! Create neighbor lists
598         call expandNeighborList(nLocal, my_status)
# Line 631 | Line 656 | contains
656      integer :: localError
657      integer :: propPack_i, propPack_j
658      integer :: loopStart, loopEnd, loop
659 <    integer :: iMap
659 >    integer :: iHash
660      real(kind=dp) :: listSkin = 1.0  
661  
662      !! initialize local variables  
# Line 750 | Line 775 | contains
775               endif
776  
777   #ifdef IS_MPI
778 +             me_j = atid_col(j)
779               call get_interatomic_vector(q_group_Row(:,i), &
780                    q_group_Col(:,j), d_grp, rgrpsq)
781   #else
782 +             me_j = atid(j)
783               call get_interatomic_vector(q_group(:,i), &
784                    q_group(:,j), d_grp, rgrpsq)
785   #endif
786  
787 <             if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
787 >             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
788                  if (update_nlist) then
789                     nlist = nlist + 1
790  
# Line 975 | Line 1002 | contains
1002   #else
1003               me_i = atid(i)
1004   #endif
1005 <             iMap = InteractionMap(me_i, me_j)%InteractionHash
1005 >             iHash = InteractionHash(me_i,me_j)
1006              
1007 <             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1007 >             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1008  
1009                  mu_i = getDipoleMoment(me_i)
1010  
# Line 1044 | Line 1071 | contains
1071      real ( kind = dp ) :: ebalance
1072      integer :: me_i, me_j
1073  
1074 <    integer :: iMap
1074 >    integer :: iHash
1075  
1076      r = sqrt(rijsq)
1077      vpair = 0.0d0
# Line 1058 | Line 1085 | contains
1085      me_j = atid(j)
1086   #endif
1087  
1088 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1088 >    iHash = InteractionHash(me_i, me_j)
1089  
1090 <    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1090 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1091         call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1092      endif
1093  
1094 <    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1094 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1095         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1096 <            pot, eFrame, f, t, do_pot)
1096 >            pot, eFrame, f, t, do_pot, corrMethod)
1097  
1098         if (FF_uses_RF .and. SIM_uses_RF) then
1099  
# Line 1077 | Line 1104 | contains
1104  
1105      endif
1106  
1107 <    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1107 >    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1108         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1109              pot, A, f, t, do_pot)
1110      endif
1111  
1112 <    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1112 >    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1113         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1114              pot, A, f, t, do_pot)
1115      endif
1116  
1117 <    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1117 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1118         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1119              pot, A, f, t, do_pot)
1120      endif
1121      
1122 <    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1122 >    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1123   !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1124   !           pot, A, f, t, do_pot)
1125      endif
1126  
1127 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1127 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1128         call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1129              do_pot)
1130      endif
1131  
1132 <    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1132 >    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1133         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1134              pot, A, f, t, do_pot)
1135      endif
1136  
1137 <    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1137 >    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1138         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1139              pot, A, f, t, do_pot)
1140      endif
# Line 1129 | Line 1156 | contains
1156      real ( kind = dp )                :: r, rc
1157      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1158  
1159 <    integer :: me_i, me_j, iMap
1159 >    integer :: me_i, me_j, iHash
1160  
1161 +    r = sqrt(rijsq)
1162 +
1163   #ifdef IS_MPI  
1164      me_i = atid_row(i)
1165      me_j = atid_col(j)  
# Line 1139 | Line 1168 | contains
1168      me_j = atid(j)  
1169   #endif
1170  
1171 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1171 >    iHash = InteractionHash(me_i, me_j)
1172  
1173 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1173 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1174              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1175      endif
1176      
# Line 1338 | Line 1367 | contains
1367  
1368    function FF_UsesDirectionalAtoms() result(doesit)
1369      logical :: doesit
1370 <    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
1370 >    doesit = FF_uses_DirectionalAtoms
1371    end function FF_UsesDirectionalAtoms
1372  
1373    function FF_RequiresPrepairCalc() result(doesit)

Diff Legend

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