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

Comparing trunk/OOPSE-3.0/src/UseTheForce/doForces.F90 (file contents):
Revision 2260 by chuckv, Mon Jun 27 22:21:37 2005 UTC vs.
Revision 2350 by chuckv, Mon Oct 10 21:20:46 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.21 2005-06-27 22:21:37 chuckv Exp $, $Date: 2005-06-27 22:21:37 $, $Name: not supported by cvs2svn $, $Revision: 1.21 $
48 > !! @version $Id: doForces.F90,v 1.51 2005-10-10 21:20:46 chuckv Exp $, $Date: 2005-10-10 21:20:46 $, $Name: not supported by cvs2svn $, $Revision: 1.51 $
49  
50  
51   module doForces
# Line 58 | Line 58 | module doForces
58    use lj
59    use sticky
60    use electrostatic_module
61 <  use reaction_field
61 >  use reaction_field_module
62    use gb_pair
63    use shapes
64    use vector_class
# Line 73 | Line 73 | module doForces
73  
74   #define __FORTRAN90
75   #include "UseTheForce/fSwitchingFunction.h"
76 + #include "UseTheForce/fCutoffPolicy.h"
77   #include "UseTheForce/DarkSide/fInteractionMap.h"
78 + #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
79  
80 +
81    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
82    INTEGER, PARAMETER:: PAIR_LOOP    = 2
83  
81 –  logical, save :: haveRlist = .false.
84    logical, save :: haveNeighborList = .false.
85    logical, save :: haveSIMvariables = .false.
84 –  logical, save :: havePropertyMap = .false.
86    logical, save :: haveSaneForceField = .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
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 129 | Line 119 | module doForces
119    real :: forceTimeInitial, forceTimeFinal
120    integer :: nLoops
121   #endif
122 +  
123 +  !! Variables for cutoff mapping and interaction mapping
124 +  ! Bit hash to determine pair-pair interactions.
125 +  integer, dimension(:,:), allocatable :: InteractionHash
126 +  real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
127 +  real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
128 +  real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
129  
130 <  type :: Properties
131 <     logical :: is_Directional   = .false.
135 <     logical :: is_LennardJones  = .false.
136 <     logical :: is_Electrostatic = .false.
137 <     logical :: is_Charge        = .false.
138 <     logical :: is_Dipole        = .false.
139 <     logical :: is_Quadrupole    = .false.
140 <     logical :: is_Sticky        = .false.
141 <     logical :: is_StickyPower   = .false.
142 <     logical :: is_GayBerne      = .false.
143 <     logical :: is_EAM           = .false.
144 <     logical :: is_Shape         = .false.
145 <     logical :: is_FLARB         = .false.
146 <  end type Properties
130 >  integer, dimension(:), allocatable, target :: groupToGtypeRow
131 >  integer, dimension(:), pointer :: groupToGtypeCol => null()
132  
133 <  type(Properties), dimension(:),allocatable :: PropertyMap
133 >  real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
134 >  real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
135 >  type ::gtypeCutoffs
136 >     real(kind=dp) :: rcut
137 >     real(kind=dp) :: rcutsq
138 >     real(kind=dp) :: rlistsq
139 >  end type gtypeCutoffs
140 >  type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
141  
142 <
142 >  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
143 >  real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
144 >  real(kind=dp),save :: listSkin
145    
152 –  type, public :: Interaction
153 –     integer :: InteractionHash
154 –     real(kind=dp) :: rCut
155 –  end type Interaction
156 –  
157 –  type(Interaction), public, dimension(:,:), allocatable :: InteractionMap
158 –  
159 –  !public :: addInteraction
160 –  !public :: setInteractionHash
161 –  !public :: getInteractionHash
162 –  public :: createInteractionMap
163 –
146   contains
147  
148 <
167 <  subroutine createInteractionMap(status)
148 >  subroutine createInteractionHash(status)
149      integer :: nAtypes
150 <    integer :: status
150 >    integer, intent(out) :: status
151      integer :: i
152      integer :: j
153 <    integer :: ihash
154 <    real(kind=dp) :: myRcut
174 < ! Test Types
153 >    integer :: iHash
154 >    !! Test Types
155      logical :: i_is_LJ
156      logical :: i_is_Elect
157      logical :: i_is_Sticky
# Line 186 | Line 166 | contains
166      logical :: j_is_GB
167      logical :: j_is_EAM
168      logical :: j_is_Shape
169 <    
170 <    
169 >    real(kind=dp) :: myRcut
170 >
171 >    status = 0  
172 >
173      if (.not. associated(atypes)) then
174 <       call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
174 >       call handleError("atype", "atypes was not present before call of createInteractionHash!")
175         status = -1
176         return
177      endif
# Line 201 | Line 183 | contains
183         return
184      end if
185  
186 <    if (.not. allocated(InteractionMap)) then
187 <       allocate(InteractionMap(nAtypes,nAtypes))
186 >    if (.not. allocated(InteractionHash)) then
187 >       allocate(InteractionHash(nAtypes,nAtypes))
188      endif
189 +
190 +    if (.not. allocated(atypeMaxCutoff)) then
191 +       allocate(atypeMaxCutoff(nAtypes))
192 +    endif
193          
194      do i = 1, nAtypes
195         call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
# Line 228 | Line 214 | contains
214            call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
215  
216            if (i_is_LJ .and. j_is_LJ) then
217 <             iHash = ior(iHash, LJ_PAIR)
218 <            
217 >             iHash = ior(iHash, LJ_PAIR)            
218 >          endif
219 >          
220 >          if (i_is_Elect .and. j_is_Elect) then
221 >             iHash = ior(iHash, ELECTROSTATIC_PAIR)
222 >          endif
223 >          
224 >          if (i_is_Sticky .and. j_is_Sticky) then
225 >             iHash = ior(iHash, STICKY_PAIR)
226 >          endif
227  
228 +          if (i_is_StickyP .and. j_is_StickyP) then
229 +             iHash = ior(iHash, STICKYPOWER_PAIR)
230 +          endif
231  
232 +          if (i_is_EAM .and. j_is_EAM) then
233 +             iHash = ior(iHash, EAM_PAIR)
234            endif
235  
237 –
238 –
239 –          if (i_is_Elect .and. j_is_Elect) iHash = ior(iHash, ELECTROSTATIC_PAIR)
240 –          if (i_is_Sticky .and. j_is_Sticky) iHash = ior(iHash, STICKY_PAIR)
241 –          if (i_is_StickyP .and. j_is_StickyP) iHash = ior(iHash, STICKYPOWER_PAIR)
242 –
243 –          if (i_is_EAM .and. j_is_EAM) iHash = ior(iHash, EAM_PAIR)
244 –
236            if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
237            if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
238            if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
# Line 251 | Line 242 | contains
242            if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
243  
244  
245 <          InteractionMap(i,j)%InteractionHash = iHash
246 <          InteractionMap(j,i)%InteractionHash = iHash
245 >          InteractionHash(i,j) = iHash
246 >          InteractionHash(j,i) = iHash
247  
248         end do
249  
250      end do
260 –  end subroutine createInteractionMap
251  
252 +    haveInteractionHash = .true.
253 +  end subroutine createInteractionHash
254  
255 +  subroutine createGtypeCutoffMap(stat)
256  
257 < !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
258 < !!$  subroutine setRlistDF( this_rlist )
259 < !!$
260 < !!$   real(kind=dp) :: this_rlist
261 < !!$
262 < !!$    rlist = this_rlist
263 < !!$    rlistsq = rlist * rlist
264 < !!$
265 < !!$    haveRlist = .true.
273 < !!$
274 < !!$  end subroutine setRlistDF
257 >    integer, intent(out), optional :: stat
258 >    logical :: i_is_LJ
259 >    logical :: i_is_Elect
260 >    logical :: i_is_Sticky
261 >    logical :: i_is_StickyP
262 >    logical :: i_is_GB
263 >    logical :: i_is_EAM
264 >    logical :: i_is_Shape
265 >    logical :: GtypeFound
266  
267 <  subroutine createPropertyMap(status)
268 <    integer :: nAtypes
269 <    integer :: status
270 <    integer :: i
271 <    logical :: thisProperty
272 <    real (kind=DP) :: thisDPproperty
273 <
283 <    status = 0
267 >    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
268 >    integer :: n_in_i, me_i, ia, g, atom1
269 >    integer :: nGroupsInRow
270 >    integer :: nGroupsInCol
271 >    integer :: nGroupTypesRow,nGroupTypesCol
272 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin
273 >    real(kind=dp) :: biggestAtypeCutoff
274  
275 +    stat = 0
276 +    if (.not. haveInteractionHash) then
277 +       call createInteractionHash(myStatus)      
278 +       if (myStatus .ne. 0) then
279 +          write(default_error, *) 'createInteractionHash failed in doForces!'
280 +          stat = -1
281 +          return
282 +       endif
283 +    endif
284 + #ifdef IS_MPI
285 +    nGroupsInRow = getNgroupsInRow(plan_group_row)
286 +    nGroupsInCol = getNgroupsInCol(plan_group_col)
287 + #endif
288      nAtypes = getSize(atypes)
289 + ! Set all of the initial cutoffs to zero.
290 +    atypeMaxCutoff = 0.0_dp
291 +    do i = 1, nAtypes
292 +       if (SimHasAtype(i)) then    
293 +          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
294 +          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
295 +          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
296 +          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
297 +          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
298 +          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
299 +          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
300 +          
301 +
302 +          if (haveDefaultCutoffs) then
303 +             atypeMaxCutoff(i) = defaultRcut
304 +          else
305 +             if (i_is_LJ) then          
306 +                thisRcut = getSigma(i) * 2.5_dp
307 +                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
308 +             endif
309 +             if (i_is_Elect) then
310 +                thisRcut = defaultRcut
311 +                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
312 +             endif
313 +             if (i_is_Sticky) then
314 +                thisRcut = getStickyCut(i)
315 +                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
316 +             endif
317 +             if (i_is_StickyP) then
318 +                thisRcut = getStickyPowerCut(i)
319 +                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
320 +             endif
321 +             if (i_is_GB) then
322 +                thisRcut = getGayBerneCut(i)
323 +                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
324 +             endif
325 +             if (i_is_EAM) then
326 +                thisRcut = getEAMCut(i)
327 +                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
328 +             endif
329 +             if (i_is_Shape) then
330 +                thisRcut = getShapeCut(i)
331 +                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
332 +             endif
333 +          endif
334 +          
335 +          
336 +          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
337 +             biggestAtypeCutoff = atypeMaxCutoff(i)
338 +          endif
339  
340 <    if (nAtypes == 0) then
341 <       status = -1
342 <       return
340 >       endif
341 >    enddo
342 >  
343 >
344 >    
345 >    istart = 1
346 >    jstart = 1
347 > #ifdef IS_MPI
348 >    iend = nGroupsInRow
349 >    jend = nGroupsInCol
350 > #else
351 >    iend = nGroups
352 >    jend = nGroups
353 > #endif
354 >    
355 >    !! allocate the groupToGtype and gtypeMaxCutoff here.
356 >    if(.not.allocated(groupToGtypeRow)) then
357 >     !  allocate(groupToGtype(iend))
358 >       allocate(groupToGtypeRow(iend))
359 >    else
360 >       deallocate(groupToGtypeRow)
361 >       allocate(groupToGtypeRow(iend))
362 >    endif
363 >    if(.not.allocated(groupMaxCutoffRow)) then
364 >       allocate(groupMaxCutoffRow(iend))
365 >    else
366 >       deallocate(groupMaxCutoffRow)
367 >       allocate(groupMaxCutoffRow(iend))
368      end if
369  
370 <    if (.not. allocated(PropertyMap)) then
371 <       allocate(PropertyMap(nAtypes))
370 >    if(.not.allocated(gtypeMaxCutoffRow)) then
371 >       allocate(gtypeMaxCutoffRow(iend))
372 >    else
373 >       deallocate(gtypeMaxCutoffRow)
374 >       allocate(gtypeMaxCutoffRow(iend))
375      endif
376  
296 –    do i = 1, nAtypes
297 –       call getElementProperty(atypes, i, "is_Directional", thisProperty)
298 –       PropertyMap(i)%is_Directional = thisProperty
377  
378 <       call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
379 <       PropertyMap(i)%is_LennardJones = thisProperty
378 > #ifdef IS_MPI
379 >       ! We only allocate new storage if we are in MPI because Ncol /= Nrow
380 >    if(.not.allocated(groupToGtypeCol)) then
381 >       allocate(groupToGtypeCol(jend))
382 >    else
383 >       deallocate(groupToGtypeCol)
384 >       allocate(groupToGtypeCol(jend))
385 >    end if
386  
387 <       call getElementProperty(atypes, i, "is_Electrostatic", thisProperty)
388 <       PropertyMap(i)%is_Electrostatic = thisProperty
387 >    if(.not.allocated(groupToGtypeCol)) then
388 >       allocate(groupToGtypeCol(jend))
389 >    else
390 >       deallocate(groupToGtypeCol)
391 >       allocate(groupToGtypeCol(jend))
392 >    end if
393 >    if(.not.allocated(gtypeMaxCutoffCol)) then
394 >       allocate(gtypeMaxCutoffCol(jend))
395 >    else
396 >       deallocate(gtypeMaxCutoffCol)      
397 >       allocate(gtypeMaxCutoffCol(jend))
398 >    end if
399  
400 <       call getElementProperty(atypes, i, "is_Charge", thisProperty)
401 <       PropertyMap(i)%is_Charge = thisProperty
400 >       groupMaxCutoffCol = 0.0_dp
401 >       gtypeMaxCutoffCol = 0.0_dp
402  
403 <       call getElementProperty(atypes, i, "is_Dipole", thisProperty)
404 <       PropertyMap(i)%is_Dipole = thisProperty
403 > #endif
404 >       groupMaxCutoffRow = 0.0_dp
405 >       gtypeMaxCutoffRow = 0.0_dp
406  
312 –       call getElementProperty(atypes, i, "is_Quadrupole", thisProperty)
313 –       PropertyMap(i)%is_Quadrupole = thisProperty
407  
408 <       call getElementProperty(atypes, i, "is_Sticky", thisProperty)
409 <       PropertyMap(i)%is_Sticky = thisProperty
410 <      
411 <       call getElementProperty(atypes, i, "is_StickyPower", thisProperty)
412 <       PropertyMap(i)%is_StickyPower = thisProperty
408 >    !! first we do a single loop over the cutoff groups to find the
409 >    !! largest cutoff for any atypes present in this group.  We also
410 >    !! create gtypes at this point.
411 >    
412 >    tol = 1.0d-6
413 >    nGroupTypesRow = 0
414  
415 <       call getElementProperty(atypes, i, "is_GayBerne", thisProperty)
416 <       PropertyMap(i)%is_GayBerne = thisProperty
415 >    do i = istart, iend      
416 >       n_in_i = groupStartRow(i+1) - groupStartRow(i)
417 >       groupMaxCutoffRow(i) = 0.0_dp
418 >       do ia = groupStartRow(i), groupStartRow(i+1)-1
419 >          atom1 = groupListRow(ia)
420 > #ifdef IS_MPI
421 >          me_i = atid_row(atom1)
422 > #else
423 >          me_i = atid(atom1)
424 > #endif          
425 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
426 >             groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
427 >          endif          
428 >       enddo
429  
430 <       call getElementProperty(atypes, i, "is_EAM", thisProperty)
431 <       PropertyMap(i)%is_EAM = thisProperty
432 <
433 <       call getElementProperty(atypes, i, "is_Shape", thisProperty)
434 <       PropertyMap(i)%is_Shape = thisProperty
430 >       if (nGroupTypesRow.eq.0) then
431 >          nGroupTypesRow = nGroupTypesRow + 1
432 >          gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
433 >          groupToGtypeRow(i) = nGroupTypesRow
434 >       else
435 >          GtypeFound = .false.
436 >          do g = 1, nGroupTypesRow
437 >             if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
438 >                groupToGtypeRow(i) = g
439 >                GtypeFound = .true.
440 >             endif
441 >          enddo
442 >          if (.not.GtypeFound) then            
443 >             nGroupTypesRow = nGroupTypesRow + 1
444 >             gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
445 >             groupToGtypeRow(i) = nGroupTypesRow
446 >          endif
447 >       endif
448 >    enddo    
449  
450 <       call getElementProperty(atypes, i, "is_FLARB", thisProperty)
451 <       PropertyMap(i)%is_FLARB = thisProperty
452 <    end do
450 > #ifdef IS_MPI
451 >    do j = jstart, jend      
452 >       n_in_j = groupStartCol(j+1) - groupStartCol(j)
453 >       groupMaxCutoffCol(j) = 0.0_dp
454 >       do ja = groupStartCol(j), groupStartCol(j+1)-1
455 >          atom1 = groupListCol(ja)
456  
457 <    havePropertyMap = .true.
457 >          me_j = atid_col(atom1)
458  
459 <  end subroutine createPropertyMap
459 >          if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
460 >             groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
461 >          endif          
462 >       enddo
463  
464 +       if (nGroupTypesCol.eq.0) then
465 +          nGroupTypesCol = nGroupTypesCol + 1
466 +          gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
467 +          groupToGtypeCol(j) = nGroupTypesCol
468 +       else
469 +          GtypeFound = .false.
470 +          do g = 1, nGroupTypesCol
471 +             if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
472 +                groupToGtypeCol(j) = g
473 +                GtypeFound = .true.
474 +             endif
475 +          enddo
476 +          if (.not.GtypeFound) then            
477 +             nGroupTypesCol = nGroupTypesCol + 1
478 +             gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
479 +             groupToGtypeCol(j) = nGroupTypesCol
480 +          endif
481 +       endif
482 +    enddo    
483 +
484 + #else
485 + ! Set pointers to information we just found
486 +    nGroupTypesCol = nGroupTypesRow
487 +    groupToGtypeCol => groupToGtypeRow
488 +    gtypeMaxCutoffCol => gtypeMaxCutoffRow
489 +    groupMaxCutoffCol => groupMaxCutoffRow
490 + #endif
491 +
492 +
493 +
494 +
495 +
496 +    !! allocate the gtypeCutoffMap here.
497 +    allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
498 +    !! then we do a double loop over all the group TYPES to find the cutoff
499 +    !! map between groups of two types
500 +    tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
501 +
502 +    do i = 1, nGroupTypesRow
503 +       do j = 1, nGroupTypesCol
504 +      
505 +          select case(cutoffPolicy)
506 +          case(TRADITIONAL_CUTOFF_POLICY)
507 +             thisRcut = tradRcut
508 +          case(MIX_CUTOFF_POLICY)
509 +             thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
510 +          case(MAX_CUTOFF_POLICY)
511 +             thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
512 +          case default
513 +             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
514 +             return
515 +          end select
516 +          gtypeCutoffMap(i,j)%rcut = thisRcut
517 +          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
518 +          skin = defaultRlist - defaultRcut
519 +          listSkin = skin ! set neighbor list skin thickness
520 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
521 +
522 +          ! sanity check
523 +
524 +          if (haveDefaultCutoffs) then
525 +             if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
526 +                call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
527 +             endif
528 +          endif
529 +       enddo
530 +    enddo
531 +    if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
532 +    if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
533 +    if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
534 + #ifdef IS_MPI
535 +    if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
536 +    if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
537 + #endif
538 +    groupMaxCutoffCol => null()
539 +    gtypeMaxCutoffCol => null()
540 +    
541 +    haveGtypeCutoffMap = .true.
542 +   end subroutine createGtypeCutoffMap
543 +
544 +   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
545 +     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
546 +     integer, intent(in) :: cutPolicy
547 +
548 +     defaultRcut = defRcut
549 +     defaultRsw = defRsw
550 +     defaultRlist = defRlist
551 +     cutoffPolicy = cutPolicy
552 +
553 +     haveDefaultCutoffs = .true.
554 +   end subroutine setDefaultCutoffs
555 +
556 +   subroutine setCutoffPolicy(cutPolicy)
557 +
558 +     integer, intent(in) :: cutPolicy
559 +     cutoffPolicy = cutPolicy
560 +     call createGtypeCutoffMap()
561 +   end subroutine setCutoffPolicy
562 +    
563 +    
564    subroutine setSimVariables()
565      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
340 –    SIM_uses_LennardJones = SimUsesLennardJones()
341 –    SIM_uses_Electrostatics = SimUsesElectrostatics()
342 –    SIM_uses_Charges = SimUsesCharges()
343 –    SIM_uses_Dipoles = SimUsesDipoles()
344 –    SIM_uses_Sticky = SimUsesSticky()
345 –    SIM_uses_StickyPower = SimUsesStickyPower()
346 –    SIM_uses_GayBerne = SimUsesGayBerne()
566      SIM_uses_EAM = SimUsesEAM()
348 –    SIM_uses_Shapes = SimUsesShapes()
349 –    SIM_uses_FLARB = SimUsesFLARB()
350 –    SIM_uses_RF = SimUsesRF()
567      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
568      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
569      SIM_uses_PBC = SimUsesPBC()
# Line 364 | Line 580 | contains
580  
581      error = 0
582  
583 <    if (.not. havePropertyMap) then
583 >    if (.not. haveInteractionHash) then      
584 >       myStatus = 0      
585 >       call createInteractionHash(myStatus)      
586 >       if (myStatus .ne. 0) then
587 >          write(default_error, *) 'createInteractionHash failed in doForces!'
588 >          error = -1
589 >          return
590 >       endif
591 >    endif
592  
593 <       myStatus = 0
594 <
595 <       call createPropertyMap(myStatus)
372 <
593 >    if (.not. haveGtypeCutoffMap) then        
594 >       myStatus = 0      
595 >       call createGtypeCutoffMap(myStatus)      
596         if (myStatus .ne. 0) then
597 <          write(default_error, *) 'createPropertyMap failed in doForces!'
597 >          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
598            error = -1
599            return
600         endif
# Line 381 | Line 604 | contains
604         call setSimVariables()
605      endif
606  
607 <    if (.not. haveRlist) then
608 <       write(default_error, *) 'rList has not been set in doForces!'
609 <       error = -1
610 <       return
611 <    endif
607 >  !  if (.not. haveRlist) then
608 >  !     write(default_error, *) 'rList has not been set in doForces!'
609 >  !     error = -1
610 >  !     return
611 >  !  endif
612  
613      if (.not. haveNeighborList) then
614         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 410 | Line 633 | contains
633    end subroutine doReadyCheck
634  
635  
636 <  subroutine init_FF(use_RF_c, thisStat)
414 <
415 <    logical, intent(in) :: use_RF_c
636 >  subroutine init_FF(thisESM, thisStat)
637  
638 +    integer, intent(in) :: thisESM
639      integer, intent(out) :: thisStat  
640      integer :: my_status, nMatches
641      integer, pointer :: MatchList(:) => null()
# Line 422 | Line 644 | contains
644      !! assume things are copacetic, unless they aren't
645      thisStat = 0
646  
647 <    !! Fortran's version of a cast:
426 <    FF_uses_RF = use_RF_c
647 >    electrostaticSummationMethod = thisESM
648  
649      !! init_FF is called *after* all of the atom types have been
650      !! defined in atype_module using the new_atype subroutine.
# Line 432 | Line 653 | contains
653      !! interactions are used by the force field.    
654  
655      FF_uses_DirectionalAtoms = .false.
435 –    FF_uses_LennardJones = .false.
436 –    FF_uses_Electrostatics = .false.
437 –    FF_uses_Charges = .false.    
656      FF_uses_Dipoles = .false.
439 –    FF_uses_Sticky = .false.
440 –    FF_uses_StickyPower = .false.
657      FF_uses_GayBerne = .false.
658      FF_uses_EAM = .false.
443 –    FF_uses_Shapes = .false.
444 –    FF_uses_FLARB = .false.
659  
660      call getMatchingElementList(atypes, "is_Directional", .true., &
661           nMatches, MatchList)
662      if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
663  
450 –    call getMatchingElementList(atypes, "is_LennardJones", .true., &
451 –         nMatches, MatchList)
452 –    if (nMatches .gt. 0) FF_uses_LennardJones = .true.
453 –
454 –    call getMatchingElementList(atypes, "is_Electrostatic", .true., &
455 –         nMatches, MatchList)
456 –    if (nMatches .gt. 0) then
457 –       FF_uses_Electrostatics = .true.
458 –    endif
459 –
460 –    call getMatchingElementList(atypes, "is_Charge", .true., &
461 –         nMatches, MatchList)
462 –    if (nMatches .gt. 0) then
463 –       FF_uses_Charges = .true.  
464 –       FF_uses_Electrostatics = .true.
465 –    endif
466 –
664      call getMatchingElementList(atypes, "is_Dipole", .true., &
665           nMatches, MatchList)
666 <    if (nMatches .gt. 0) then
470 <       FF_uses_Dipoles = .true.
471 <       FF_uses_Electrostatics = .true.
472 <       FF_uses_DirectionalAtoms = .true.
473 <    endif
474 <
475 <    call getMatchingElementList(atypes, "is_Quadrupole", .true., &
476 <         nMatches, MatchList)
477 <    if (nMatches .gt. 0) then
478 <       FF_uses_Quadrupoles = .true.
479 <       FF_uses_Electrostatics = .true.
480 <       FF_uses_DirectionalAtoms = .true.
481 <    endif
482 <
483 <    call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
484 <         MatchList)
485 <    if (nMatches .gt. 0) then
486 <       FF_uses_Sticky = .true.
487 <       FF_uses_DirectionalAtoms = .true.
488 <    endif
489 <
490 <    call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
491 <         MatchList)
492 <    if (nMatches .gt. 0) then
493 <       FF_uses_StickyPower = .true.
494 <       FF_uses_DirectionalAtoms = .true.
495 <    endif
666 >    if (nMatches .gt. 0) FF_uses_Dipoles = .true.
667      
668      call getMatchingElementList(atypes, "is_GayBerne", .true., &
669           nMatches, MatchList)
670 <    if (nMatches .gt. 0) then
500 <       FF_uses_GayBerne = .true.
501 <       FF_uses_DirectionalAtoms = .true.
502 <    endif
670 >    if (nMatches .gt. 0) FF_uses_GayBerne = .true.
671  
672      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
673      if (nMatches .gt. 0) FF_uses_EAM = .true.
674  
507 –    call getMatchingElementList(atypes, "is_Shape", .true., &
508 –         nMatches, MatchList)
509 –    if (nMatches .gt. 0) then
510 –       FF_uses_Shapes = .true.
511 –       FF_uses_DirectionalAtoms = .true.
512 –    endif
675  
514 –    call getMatchingElementList(atypes, "is_FLARB", .true., &
515 –         nMatches, MatchList)
516 –    if (nMatches .gt. 0) FF_uses_FLARB = .true.
517 –
518 –    !! Assume sanity (for the sake of argument)
676      haveSaneForceField = .true.
677  
678 <    !! check to make sure the FF_uses_RF setting makes sense
678 >    !! check to make sure the reaction field setting makes sense
679  
680 <    if (FF_uses_dipoles) then
681 <       if (FF_uses_RF) then
680 >    if (FF_uses_Dipoles) then
681 >       if (electrostaticSummationMethod == REACTION_FIELD) then
682            dielect = getDielect()
683            call initialize_rf(dielect)
684         endif
685      else
686 <       if (FF_uses_RF) then          
686 >       if (electrostaticSummationMethod == REACTION_FIELD) then
687            write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
688            thisStat = -1
689            haveSaneForceField = .false.
690            return
691         endif
692      endif
536 –
537 –    !sticky module does not contain check_sticky_FF anymore
538 –    !if (FF_uses_sticky) then
539 –    !   call check_sticky_FF(my_status)
540 –    !   if (my_status /= 0) then
541 –    !      thisStat = -1
542 –    !      haveSaneForceField = .false.
543 –    !      return
544 –    !   end if
545 –    !endif
693  
694      if (FF_uses_EAM) then
695         call init_EAM_FF(my_status)
# Line 563 | Line 710 | contains
710         endif
711      endif
712  
566 –    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
567 –    endif
568 –
713      if (.not. haveNeighborList) then
714         !! Create neighbor lists
715         call expandNeighborList(nLocal, my_status)
# Line 629 | Line 773 | contains
773      integer :: localError
774      integer :: propPack_i, propPack_j
775      integer :: loopStart, loopEnd, loop
776 +    integer :: iHash
777 +  
778  
633 –    real(kind=dp) :: listSkin = 1.0  
634 –
779      !! initialize local variables  
780  
781   #ifdef IS_MPI
# Line 748 | Line 892 | contains
892               endif
893  
894   #ifdef IS_MPI
895 +             me_j = atid_col(j)
896               call get_interatomic_vector(q_group_Row(:,i), &
897                    q_group_Col(:,j), d_grp, rgrpsq)
898   #else
899 +             me_j = atid(j)
900               call get_interatomic_vector(q_group(:,i), &
901                    q_group(:,j), d_grp, rgrpsq)
902 < #endif
902 > #endif      
903  
904 <             if (rgrpsq < rlistsq) then
904 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
905                  if (update_nlist) then
906                     nlist = nlist + 1
907  
# Line 876 | Line 1022 | contains
1022                  endif
1023               end if
1024            enddo
1025 +
1026         enddo outer
1027  
1028         if (update_nlist) then
# Line 955 | Line 1102 | contains
1102  
1103      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1104  
1105 <       if (FF_uses_RF .and. SIM_uses_RF) then
1105 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1106  
1107   #ifdef IS_MPI
1108            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 973 | Line 1120 | contains
1120   #else
1121               me_i = atid(i)
1122   #endif
1123 <
1124 <             if (PropertyMap(me_i)%is_Dipole) then
1123 >             iHash = InteractionHash(me_i,me_j)
1124 >            
1125 >             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1126  
1127                  mu_i = getDipoleMoment(me_i)
1128  
# Line 1038 | Line 1186 | contains
1186      real ( kind = dp ), intent(inout) :: rijsq
1187      real ( kind = dp )                :: r
1188      real ( kind = dp ), intent(inout) :: d(3)
1041 –    real ( kind = dp ) :: ebalance
1189      integer :: me_i, me_j
1190  
1191 <    integer :: iMap
1191 >    integer :: iHash
1192  
1193      r = sqrt(rijsq)
1194      vpair = 0.0d0
# Line 1055 | Line 1202 | contains
1202      me_j = atid(j)
1203   #endif
1204  
1205 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1205 >    iHash = InteractionHash(me_i, me_j)
1206  
1207 <    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1207 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1208         call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1209      endif
1210  
1211 <    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1211 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1212         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1213              pot, eFrame, f, t, do_pot)
1214  
1215 <       if (FF_uses_dipoles .and. SIM_uses_dipoles) then                
1216 <          if ( PropertyMap(me_i)%is_Dipole .and. &
1217 <               PropertyMap(me_j)%is_Dipole) then
1218 <             if (FF_uses_RF .and. SIM_uses_RF) then
1219 <                call accumulate_rf(i, j, r, eFrame, sw)
1073 <                call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1074 <             endif
1075 <          endif
1215 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1216 >
1217 >          ! CHECK ME (RF needs to know about all electrostatic types)
1218 >          call accumulate_rf(i, j, r, eFrame, sw)
1219 >          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1220         endif
1221 +
1222      endif
1223  
1224 <    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1224 >    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1225         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1226              pot, A, f, t, do_pot)
1227      endif
1228  
1229 <    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1229 >    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1230         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1231              pot, A, f, t, do_pot)
1232      endif
1233  
1234 <    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1234 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1235         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1236              pot, A, f, t, do_pot)
1237      endif
1238      
1239 <    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1240 <       call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1241 <            pot, A, f, t, do_pot)
1239 >    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1240 > !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1241 > !           pot, A, f, t, do_pot)
1242      endif
1243  
1244 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1244 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1245         call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1246              do_pot)
1247      endif
1248  
1249 <    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1249 >    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1250         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1251              pot, A, f, t, do_pot)
1252      endif
1253  
1254 <    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1254 >    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1255         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1256              pot, A, f, t, do_pot)
1257      endif
# Line 1128 | Line 1273 | contains
1273      real ( kind = dp )                :: r, rc
1274      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1275  
1276 <    integer :: me_i, me_j, iMap
1276 >    integer :: me_i, me_j, iHash
1277  
1278 +    r = sqrt(rijsq)
1279 +
1280   #ifdef IS_MPI  
1281      me_i = atid_row(i)
1282      me_j = atid_col(j)  
# Line 1138 | Line 1285 | contains
1285      me_j = atid(j)  
1286   #endif
1287  
1288 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1288 >    iHash = InteractionHash(me_i, me_j)
1289  
1290 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1290 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1291              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1292      endif
1293      
# Line 1337 | Line 1484 | contains
1484  
1485    function FF_UsesDirectionalAtoms() result(doesit)
1486      logical :: doesit
1487 <    doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1341 <         FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1342 <         FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1487 >    doesit = FF_uses_DirectionalAtoms
1488    end function FF_UsesDirectionalAtoms
1489  
1490    function FF_RequiresPrepairCalc() result(doesit)
# Line 1349 | Line 1494 | contains
1494  
1495    function FF_RequiresPostpairCalc() result(doesit)
1496      logical :: doesit
1497 <    doesit = FF_uses_RF
1497 >    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1498    end function FF_RequiresPostpairCalc
1499  
1500   #ifdef PROFILE

Diff Legend

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