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 2323 by chuckv, Fri Sep 23 20:31:02 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.49 2005-09-23 20:31:02 chuckv Exp $, $Date: 2005-09-23 20:31:02 $, $Name: not supported by cvs2svn $, $Revision: 1.49 $
49  
50  
51   module doForces
# Line 58 | Line 58 | module doForces
58    use lj
59    use sticky
60    use electrostatic_module
61 <  use reaction_field
61 >  use reaction_field_module
62    use gb_pair
63    use shapes
64    use vector_class
# Line 73 | Line 73 | module doForces
73  
74   #define __FORTRAN90
75   #include "UseTheForce/fSwitchingFunction.h"
76 + #include "UseTheForce/fCutoffPolicy.h"
77   #include "UseTheForce/DarkSide/fInteractionMap.h"
78 + #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
79  
80 +
81    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
82    INTEGER, PARAMETER:: PAIR_LOOP    = 2
83  
81  logical, save :: haveRlist = .false.
84    logical, save :: haveNeighborList = .false.
85    logical, save :: haveSIMvariables = .false.
86    logical, save :: haveSaneForceField = .false.
87 <  logical, save :: haveInteractionMap = .false.
87 >  logical, save :: haveInteractionHash = .false.
88 >  logical, save :: haveGtypeCutoffMap = .false.
89 >  logical, save :: haveDefaultCutoffs = .false.
90 >  logical, save :: haveRlist = .false.
91  
92    logical, save :: FF_uses_DirectionalAtoms
88  logical, save :: FF_uses_LennardJones
89  logical, save :: FF_uses_Electrostatics
90  logical, save :: FF_uses_Charges
93    logical, save :: FF_uses_Dipoles
92  logical, save :: FF_uses_Quadrupoles
93  logical, save :: FF_uses_Sticky
94  logical, save :: FF_uses_StickyPower
94    logical, save :: FF_uses_GayBerne
95    logical, save :: FF_uses_EAM
97  logical, save :: FF_uses_Shapes
98  logical, save :: FF_uses_FLARB
99  logical, save :: FF_uses_RF
96  
97    logical, save :: SIM_uses_DirectionalAtoms
102  logical, save :: SIM_uses_LennardJones
103  logical, save :: SIM_uses_Electrostatics
104  logical, save :: SIM_uses_Charges
105  logical, save :: SIM_uses_Dipoles
106  logical, save :: SIM_uses_Quadrupoles
107  logical, save :: SIM_uses_Sticky
108  logical, save :: SIM_uses_StickyPower
109  logical, save :: SIM_uses_GayBerne
98    logical, save :: SIM_uses_EAM
111  logical, save :: SIM_uses_Shapes
112  logical, save :: SIM_uses_FLARB
113  logical, save :: SIM_uses_RF
99    logical, save :: SIM_requires_postpair_calc
100    logical, save :: SIM_requires_prepair_calc
101    logical, save :: SIM_uses_PBC
117  logical, save :: SIM_uses_molecular_cutoffs
102  
103 <  !!!GO AWAY---------
120 <  !!!!!real(kind=dp), save :: rlist, rlistsq
103 >  integer, save :: electrostaticSummationMethod
104  
105    public :: init_FF
106 +  public :: setDefaultCutoffs
107    public :: do_force_loop
108 < !  public :: setRlistDF
109 <  !public :: addInteraction
110 <  !public :: setInteractionHash
111 <  !public :: getInteractionHash
112 <  public :: createInteractionMap
113 <  public :: createRcuts
108 >  public :: createInteractionHash
109 >  public :: createGtypeCutoffMap
110 >  public :: getStickyCut
111 >  public :: getStickyPowerCut
112 >  public :: getGayBerneCut
113 >  public :: getEAMCut
114 >  public :: getShapeCut
115  
116   #ifdef PROFILE
117    public :: getforcetime
# Line 134 | Line 119 | module doForces
119    real :: forceTimeInitial, forceTimeFinal
120    integer :: nLoops
121   #endif
137
138  type, public :: Interaction
139     integer :: InteractionHash
140     real(kind=dp) :: rList = 0.0_dp
141     real(kind=dp) :: rListSq = 0.0_dp
142  end type Interaction
122    
123 <  type(Interaction), dimension(:,:),allocatable :: InteractionMap
124 <  
123 >  !! Variables for cutoff mapping and interaction mapping
124 >  ! Bit hash to determine pair-pair interactions.
125 >  integer, dimension(:,:), allocatable :: InteractionHash
126 >  real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
127 >  real(kind=dp), dimension(:), allocatable :: 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 :: listSkin
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 <
279 <    nAtypes = getSize(atypes)
280 <    !! If we pass a default rcut, set all atypes to that cutoff distance
281 <    if(present(defaultRList)) then
282 <       InteractionMap(:,:)%rList = defaultRList
283 <       InteractionMap(:,:)%rListSq = defaultRList*defaultRList
284 <       haveRlist = .true.
285 <       return
286 <    end if
287 <
288 <    do map_i = 1,nAtypes
289 <       do map_j = map_i,nAtypes
290 <          iMap = InteractionMap(map_i, map_j)%InteractionHash
277 > #ifdef IS_MPI
278 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
279 > #endif
280 >    nAtypes = getSize(atypes)
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 (haveDefaultCutoffs) then
295 >             atypeMaxCutoff(i) = defaultRcut
296 >          else
297 >             if (i_is_LJ) then          
298 >                thisRcut = getSigma(i) * 2.5_dp
299 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
300 >             endif
301 >             if (i_is_Elect) then
302 >                thisRcut = defaultRcut
303 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
304 >             endif
305 >             if (i_is_Sticky) then
306 >                thisRcut = getStickyCut(i)
307 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
308 >             endif
309 >             if (i_is_StickyP) then
310 >                thisRcut = getStickyPowerCut(i)
311 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
312 >             endif
313 >             if (i_is_GB) then
314 >                thisRcut = getGayBerneCut(i)
315 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
316 >             endif
317 >             if (i_is_EAM) then
318 >                thisRcut = getEAMCut(i)
319 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
320 >             endif
321 >             if (i_is_Shape) then
322 >                thisRcut = getShapeCut(i)
323 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
324 >             endif
325            endif
326            
296          if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
297             ! thisRCut = getElectrostaticCutOff(map_i,map_j)
298             if (thisRcut > actualCutoff) actualCutoff = thisRcut
299          endif
327            
328 <          if ( iand(iMap, STICKY_PAIR).ne.0 ) then
329 <             ! thisRCut = getStickyCutOff(map_i,map_j)
330 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
304 <           endif
305 <          
306 <           if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
307 <              ! thisRCut = getStickyPowerCutOff(map_i,map_j)
308 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
309 <           endif
310 <          
311 <           if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
312 <              ! thisRCut = getGayberneCutOff(map_i,map_j)
313 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
314 <           endif
315 <          
316 <           if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
317 < !              thisRCut = getGaybrneLJCutOff(map_i,map_j)
318 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
319 <           endif
320 <          
321 <           if ( iand(iMap, EAM_PAIR).ne.0 ) then      
322 < !              thisRCut = getEAMCutOff(map_i,map_j)
323 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
324 <           endif
325 <          
326 <           if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
327 < !              thisRCut = getShapeCutOff(map_i,map_j)
328 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
329 <           endif
330 <          
331 <           if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
332 < !              thisRCut = getShapeLJCutOff(map_i,map_j)
333 <              if (thisRcut > actualCutoff) actualCutoff = thisRcut
334 <           endif
335 <           InteractionMap(map_i, map_j)%rList = actualCutoff
336 <           InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
337 <        end do
338 <     end do
339 <     haveRlist = .true.
340 <  end subroutine createRcuts
328 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
329 >             biggestAtypeCutoff = atypeMaxCutoff(i)
330 >          endif
331  
332 <
333 < !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
334 < !!$  subroutine setRlistDF( this_rlist )
335 < !!$
336 < !!$   real(kind=dp) :: this_rlist
337 < !!$
338 < !!$    rlist = this_rlist
339 < !!$    rlistsq = rlist * rlist
340 < !!$
341 < !!$    haveRlist = .true.
342 < !!$
343 < !!$  end subroutine setRlistDF
332 >       endif
333 >    enddo
334 >  
335 >    nGroupTypes = 0
336 >    
337 >    istart = 1
338 > #ifdef IS_MPI
339 >    iend = nGroupsInRow
340 > #else
341 >    iend = nGroups
342 > #endif
343 >    
344 >    !! allocate the groupToGtype and gtypeMaxCutoff here.
345 >    if(.not.allocated(groupToGtype)) then
346 >       allocate(groupToGtype(iend))
347 >       allocate(groupMaxCutoff(iend))
348 >       allocate(gtypeMaxCutoff(iend))
349 >       groupMaxCutoff = 0.0_dp
350 >       gtypeMaxCutoff = 0.0_dp
351 >    endif
352 >    !! first we do a single loop over the cutoff groups to find the
353 >    !! largest cutoff for any atypes present in this group.  We also
354 >    !! create gtypes at this point.
355 >    
356 >    tol = 1.0d-6
357 >    
358 >    do i = istart, iend      
359 >       n_in_i = groupStartRow(i+1) - groupStartRow(i)
360 >       groupMaxCutoff(i) = 0.0_dp
361 >       do ia = groupStartRow(i), groupStartRow(i+1)-1
362 >          atom1 = groupListRow(ia)
363 > #ifdef IS_MPI
364 >          me_i = atid_row(atom1)
365 > #else
366 >          me_i = atid(atom1)
367 > #endif          
368 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
369 >             groupMaxCutoff(i)=atypeMaxCutoff(me_i)
370 >          endif          
371 >       enddo
372  
373 +       if (nGroupTypes.eq.0) then
374 +          nGroupTypes = nGroupTypes + 1
375 +          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
376 +          groupToGtype(i) = nGroupTypes
377 +       else
378 +          GtypeFound = .false.
379 +          do g = 1, nGroupTypes
380 +             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
381 +                groupToGtype(i) = g
382 +                GtypeFound = .true.
383 +             endif
384 +          enddo
385 +          if (.not.GtypeFound) then            
386 +             nGroupTypes = nGroupTypes + 1
387 +             gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
388 +             groupToGtype(i) = nGroupTypes
389 +          endif
390 +       endif
391 +    enddo    
392  
393 +    !! allocate the gtypeCutoffMap here.
394 +    allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
395 +    !! then we do a double loop over all the group TYPES to find the cutoff
396 +    !! map between groups of two types
397 +    
398 +    do i = 1, nGroupTypes
399 +       do j = 1, nGroupTypes
400 +      
401 +          select case(cutoffPolicy)
402 +          case(TRADITIONAL_CUTOFF_POLICY)
403 +             thisRcut = maxval(gtypeMaxCutoff)
404 +          case(MIX_CUTOFF_POLICY)
405 +             thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
406 +          case(MAX_CUTOFF_POLICY)
407 +             thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
408 +          case default
409 +             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
410 +             return
411 +          end select
412 +          gtypeCutoffMap(i,j)%rcut = thisRcut
413 +          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
414 +          skin = defaultRlist - defaultRcut
415 +          listSkin = skin ! set neighbor list skin thickness
416 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
417 +
418 +          ! sanity check
419 +
420 +          if (haveDefaultCutoffs) then
421 +             if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
422 +                call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
423 +             endif
424 +          endif
425 +       enddo
426 +    enddo
427 +
428 +    haveGtypeCutoffMap = .true.
429 +   end subroutine createGtypeCutoffMap
430 +
431 +   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
432 +     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
433 +     integer, intent(in) :: cutPolicy
434 +
435 +     defaultRcut = defRcut
436 +     defaultRsw = defRsw
437 +     defaultRlist = defRlist
438 +     cutoffPolicy = cutPolicy
439 +
440 +     haveDefaultCutoffs = .true.
441 +   end subroutine setDefaultCutoffs
442 +
443 +   subroutine setCutoffPolicy(cutPolicy)
444 +
445 +     integer, intent(in) :: cutPolicy
446 +     cutoffPolicy = cutPolicy
447 +     call createGtypeCutoffMap()
448 +   end subroutine setCutoffPolicy
449 +    
450 +    
451    subroutine setSimVariables()
452      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()
453      SIM_uses_EAM = SimUsesEAM()
366    SIM_uses_Shapes = SimUsesShapes()
367    SIM_uses_FLARB = SimUsesFLARB()
368    SIM_uses_RF = SimUsesRF()
454      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
455      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
456      SIM_uses_PBC = SimUsesPBC()
# Line 382 | Line 467 | contains
467  
468      error = 0
469  
470 <    if (.not. haveInteractionMap) then
470 >    if (.not. haveInteractionHash) then      
471 >       myStatus = 0      
472 >       call createInteractionHash(myStatus)      
473 >       if (myStatus .ne. 0) then
474 >          write(default_error, *) 'createInteractionHash failed in doForces!'
475 >          error = -1
476 >          return
477 >       endif
478 >    endif
479  
480 <       myStatus = 0
481 <
482 <       call createInteractionMap(myStatus)
390 <
480 >    if (.not. haveGtypeCutoffMap) then        
481 >       myStatus = 0      
482 >       call createGtypeCutoffMap(myStatus)      
483         if (myStatus .ne. 0) then
484 <          write(default_error, *) 'createInteractionMap failed in doForces!'
484 >          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
485            error = -1
486            return
487         endif
# Line 399 | Line 491 | contains
491         call setSimVariables()
492      endif
493  
494 <    if (.not. haveRlist) then
495 <       write(default_error, *) 'rList has not been set in doForces!'
496 <       error = -1
497 <       return
498 <    endif
494 >  !  if (.not. haveRlist) then
495 >  !     write(default_error, *) 'rList has not been set in doForces!'
496 >  !     error = -1
497 >  !     return
498 >  !  endif
499  
500      if (.not. haveNeighborList) then
501         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 428 | Line 520 | contains
520    end subroutine doReadyCheck
521  
522  
523 <  subroutine init_FF(use_RF_c, thisStat)
523 >  subroutine init_FF(thisESM, thisStat)
524  
525 <    logical, intent(in) :: use_RF_c
434 <
525 >    integer, intent(in) :: thisESM
526      integer, intent(out) :: thisStat  
527      integer :: my_status, nMatches
528      integer, pointer :: MatchList(:) => null()
# Line 440 | Line 531 | contains
531      !! assume things are copacetic, unless they aren't
532      thisStat = 0
533  
534 <    !! Fortran's version of a cast:
444 <    FF_uses_RF = use_RF_c
534 >    electrostaticSummationMethod = thisESM
535  
536      !! init_FF is called *after* all of the atom types have been
537      !! defined in atype_module using the new_atype subroutine.
# Line 450 | Line 540 | contains
540      !! interactions are used by the force field.    
541  
542      FF_uses_DirectionalAtoms = .false.
453    FF_uses_LennardJones = .false.
454    FF_uses_Electrostatics = .false.
455    FF_uses_Charges = .false.    
543      FF_uses_Dipoles = .false.
457    FF_uses_Sticky = .false.
458    FF_uses_StickyPower = .false.
544      FF_uses_GayBerne = .false.
545      FF_uses_EAM = .false.
461    FF_uses_Shapes = .false.
462    FF_uses_FLARB = .false.
546  
547      call getMatchingElementList(atypes, "is_Directional", .true., &
548           nMatches, MatchList)
549      if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
550  
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
551      call getMatchingElementList(atypes, "is_Dipole", .true., &
552           nMatches, MatchList)
553 <    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
553 >    if (nMatches .gt. 0) FF_uses_Dipoles = .true.
554      
555      call getMatchingElementList(atypes, "is_GayBerne", .true., &
556           nMatches, MatchList)
557 <    if (nMatches .gt. 0) then
518 <       FF_uses_GayBerne = .true.
519 <       FF_uses_DirectionalAtoms = .true.
520 <    endif
557 >    if (nMatches .gt. 0) FF_uses_GayBerne = .true.
558  
559      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
560      if (nMatches .gt. 0) FF_uses_EAM = .true.
561  
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
562  
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)
563      haveSaneForceField = .true.
564  
565 <    !! check to make sure the FF_uses_RF setting makes sense
565 >    !! check to make sure the reaction field setting makes sense
566  
567 <    if (FF_uses_dipoles) then
568 <       if (FF_uses_RF) then
567 >    if (FF_uses_Dipoles) then
568 >       if (electrostaticSummationMethod == REACTION_FIELD) then
569            dielect = getDielect()
570            call initialize_rf(dielect)
571         endif
572      else
573 <       if (FF_uses_RF) then          
573 >       if (electrostaticSummationMethod == REACTION_FIELD) then
574            write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
575            thisStat = -1
576            haveSaneForceField = .false.
577            return
578         endif
579      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
580  
581      if (FF_uses_EAM) then
582         call init_EAM_FF(my_status)
# Line 581 | Line 597 | contains
597         endif
598      endif
599  
584    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
585    endif
586
600      if (.not. haveNeighborList) then
601         !! Create neighbor lists
602         call expandNeighborList(nLocal, my_status)
# Line 647 | Line 660 | contains
660      integer :: localError
661      integer :: propPack_i, propPack_j
662      integer :: loopStart, loopEnd, loop
663 <    integer :: iMap
664 <    real(kind=dp) :: listSkin = 1.0  
663 >    integer :: iHash
664 >  
665  
666      !! initialize local variables  
667  
# Line 738 | Line 751 | contains
751         iend = nGroups - 1
752   #endif
753         outer: do i = istart, iend
741
742 #ifdef IS_MPI
743             me_i = atid_row(i)
744 #else
745             me_i = atid(i)
746 #endif
754  
755            if (update_nlist) point(i) = nlist + 1
756  
# Line 779 | Line 786 | contains
786               me_j = atid(j)
787               call get_interatomic_vector(q_group(:,i), &
788                    q_group(:,j), d_grp, rgrpsq)
789 < #endif
789 > #endif      
790  
791 <             if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
791 >             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
792                  if (update_nlist) then
793                     nlist = nlist + 1
794  
# Line 981 | Line 988 | contains
988  
989      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
990  
991 <       if (FF_uses_RF .and. SIM_uses_RF) then
991 >       if (electrostaticSummationMethod == REACTION_FIELD) then
992  
993   #ifdef IS_MPI
994            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 999 | Line 1006 | contains
1006   #else
1007               me_i = atid(i)
1008   #endif
1009 <             iMap = InteractionMap(me_i, me_j)%InteractionHash
1009 >             iHash = InteractionHash(me_i,me_j)
1010              
1011 <             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1011 >             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1012  
1013                  mu_i = getDipoleMoment(me_i)
1014  
# Line 1065 | Line 1072 | contains
1072      real ( kind = dp ), intent(inout) :: rijsq
1073      real ( kind = dp )                :: r
1074      real ( kind = dp ), intent(inout) :: d(3)
1068    real ( kind = dp ) :: ebalance
1075      integer :: me_i, me_j
1076  
1077 <    integer :: iMap
1077 >    integer :: iHash
1078  
1079      r = sqrt(rijsq)
1080      vpair = 0.0d0
# Line 1082 | Line 1088 | contains
1088      me_j = atid(j)
1089   #endif
1090  
1091 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1091 >    iHash = InteractionHash(me_i, me_j)
1092  
1093 <    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1093 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1094         call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1095      endif
1096  
1097 <    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1097 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1098         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1099              pot, eFrame, f, t, do_pot)
1100  
1101 <       if (FF_uses_RF .and. SIM_uses_RF) then
1101 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1102  
1103            ! CHECK ME (RF needs to know about all electrostatic types)
1104            call accumulate_rf(i, j, r, eFrame, sw)
# Line 1101 | Line 1107 | contains
1107  
1108      endif
1109  
1110 <    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1110 >    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1111         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1112              pot, A, f, t, do_pot)
1113      endif
1114  
1115 <    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1115 >    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1116         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1117              pot, A, f, t, do_pot)
1118      endif
1119  
1120 <    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1120 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1121         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1122              pot, A, f, t, do_pot)
1123      endif
1124      
1125 <    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1125 >    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1126   !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1127   !           pot, A, f, t, do_pot)
1128      endif
1129  
1130 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1130 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1131         call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1132              do_pot)
1133      endif
1134  
1135 <    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1135 >    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1136         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1137              pot, A, f, t, do_pot)
1138      endif
1139  
1140 <    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1140 >    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1141         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1142              pot, A, f, t, do_pot)
1143      endif
# Line 1153 | Line 1159 | contains
1159      real ( kind = dp )                :: r, rc
1160      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1161  
1162 <    integer :: me_i, me_j, iMap
1162 >    integer :: me_i, me_j, iHash
1163  
1164 +    r = sqrt(rijsq)
1165 +
1166   #ifdef IS_MPI  
1167      me_i = atid_row(i)
1168      me_j = atid_col(j)  
# Line 1163 | Line 1171 | contains
1171      me_j = atid(j)  
1172   #endif
1173  
1174 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1174 >    iHash = InteractionHash(me_i, me_j)
1175  
1176 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1176 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1177              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1178      endif
1179      
# Line 1362 | Line 1370 | contains
1370  
1371    function FF_UsesDirectionalAtoms() result(doesit)
1372      logical :: doesit
1373 <    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
1373 >    doesit = FF_uses_DirectionalAtoms
1374    end function FF_UsesDirectionalAtoms
1375  
1376    function FF_RequiresPrepairCalc() result(doesit)
# Line 1374 | Line 1380 | contains
1380  
1381    function FF_RequiresPostpairCalc() result(doesit)
1382      logical :: doesit
1383 <    doesit = FF_uses_RF
1383 >    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1384    end function FF_RequiresPostpairCalc
1385  
1386   #ifdef PROFILE

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines