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 2231 by chrisfen, Wed May 18 18:31:40 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.18 2005-05-18 18:31:40 chrisfen Exp $, $Date: 2005-05-18 18:31:40 $, $Name: not supported by cvs2svn $, $Revision: 1.18 $
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  
80 –  logical, save :: haveRlist = .false.
84    logical, save :: haveNeighborList = .false.
85    logical, save :: haveSIMvariables = .false.
83 –  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
87 –  logical, save :: FF_uses_LennardJones
88 –  logical, save :: FF_uses_Electrostatics
89 –  logical, save :: FF_uses_Charges
93    logical, save :: FF_uses_Dipoles
91 –  logical, save :: FF_uses_Quadrupoles
92 –  logical, save :: FF_uses_Sticky
93 –  logical, save :: FF_uses_StickyPower
94    logical, save :: FF_uses_GayBerne
95    logical, save :: FF_uses_EAM
96 –  logical, save :: FF_uses_Shapes
97 –  logical, save :: FF_uses_FLARB
98 –  logical, save :: FF_uses_RF
96  
97    logical, save :: SIM_uses_DirectionalAtoms
101 –  logical, save :: SIM_uses_LennardJones
102 –  logical, save :: SIM_uses_Electrostatics
103 –  logical, save :: SIM_uses_Charges
104 –  logical, save :: SIM_uses_Dipoles
105 –  logical, save :: SIM_uses_Quadrupoles
106 –  logical, save :: SIM_uses_Sticky
107 –  logical, save :: SIM_uses_StickyPower
108 –  logical, save :: SIM_uses_GayBerne
98    logical, save :: SIM_uses_EAM
110 –  logical, save :: SIM_uses_Shapes
111 –  logical, save :: SIM_uses_FLARB
112 –  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
116 –  logical, save :: SIM_uses_molecular_cutoffs
102  
103 <  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 127 | 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 :: 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 <  type :: Properties
138 <     logical :: is_Directional   = .false.
139 <     logical :: is_LennardJones  = .false.
140 <     logical :: is_Electrostatic = .false.
135 <     logical :: is_Charge        = .false.
136 <     logical :: is_Dipole        = .false.
137 <     logical :: is_Quadrupole    = .false.
138 <     logical :: is_Sticky        = .false.
139 <     logical :: is_StickyPower   = .false.
140 <     logical :: is_GayBerne      = .false.
141 <     logical :: is_EAM           = .false.
142 <     logical :: is_Shape         = .false.
143 <     logical :: is_FLARB         = .false.
144 <  end type Properties
145 <
146 <  type(Properties), dimension(:),allocatable :: PropertyMap
147 <
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 <  subroutine setRlistDF( this_rlist )
151 <
152 <    real(kind=dp) :: this_rlist
153 <
154 <    rlist = this_rlist
155 <    rlistsq = rlist * rlist
156 <
157 <    haveRlist = .true.
158 <
159 <  end subroutine setRlistDF
160 <
161 <  subroutine createPropertyMap(status)
143 >  subroutine createInteractionHash(status)
144      integer :: nAtypes
145 <    integer :: status
145 >    integer, intent(out) :: status
146      integer :: i
147 <    logical :: thisProperty
148 <    real (kind=DP) :: thisDPproperty
147 >    integer :: j
148 >    integer :: iHash
149 >    !! Test Types
150 >    logical :: i_is_LJ
151 >    logical :: i_is_Elect
152 >    logical :: i_is_Sticky
153 >    logical :: i_is_StickyP
154 >    logical :: i_is_GB
155 >    logical :: i_is_EAM
156 >    logical :: i_is_Shape
157 >    logical :: j_is_LJ
158 >    logical :: j_is_Elect
159 >    logical :: j_is_Sticky
160 >    logical :: j_is_StickyP
161 >    logical :: j_is_GB
162 >    logical :: j_is_EAM
163 >    logical :: j_is_Shape
164 >    real(kind=dp) :: myRcut
165  
166 <    status = 0
166 >    status = 0  
167  
168 +    if (.not. associated(atypes)) then
169 +       call handleError("atype", "atypes was not present before call of createInteractionHash!")
170 +       status = -1
171 +       return
172 +    endif
173 +    
174      nAtypes = getSize(atypes)
175 <
175 >    
176      if (nAtypes == 0) then
177         status = -1
178         return
179      end if
180  
181 <    if (.not. allocated(PropertyMap)) then
182 <       allocate(PropertyMap(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_Directional", thisProperty)
191 <       PropertyMap(i)%is_Directional = thisProperty
190 >       call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
191 >       call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
192 >       call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
193 >       call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
194 >       call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
195 >       call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
196 >       call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
197  
198 <       call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
186 <       PropertyMap(i)%is_LennardJones = thisProperty
198 >       do j = i, nAtypes
199  
200 <       call getElementProperty(atypes, i, "is_Electrostatic", thisProperty)
201 <       PropertyMap(i)%is_Electrostatic = thisProperty
200 >          iHash = 0
201 >          myRcut = 0.0_dp
202  
203 <       call getElementProperty(atypes, i, "is_Charge", thisProperty)
204 <       PropertyMap(i)%is_Charge = thisProperty
203 >          call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ)
204 >          call getElementProperty(atypes, j, "is_Electrostatic", j_is_Elect)
205 >          call getElementProperty(atypes, j, "is_Sticky", j_is_Sticky)
206 >          call getElementProperty(atypes, j, "is_StickyPower", j_is_StickyP)
207 >          call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
208 >          call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
209 >          call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
210  
211 <       call getElementProperty(atypes, i, "is_Dipole", thisProperty)
212 <       PropertyMap(i)%is_Dipole = thisProperty
211 >          if (i_is_LJ .and. j_is_LJ) then
212 >             iHash = ior(iHash, LJ_PAIR)            
213 >          endif
214 >          
215 >          if (i_is_Elect .and. j_is_Elect) then
216 >             iHash = ior(iHash, ELECTROSTATIC_PAIR)
217 >          endif
218 >          
219 >          if (i_is_Sticky .and. j_is_Sticky) then
220 >             iHash = ior(iHash, STICKY_PAIR)
221 >          endif
222  
223 <       call getElementProperty(atypes, i, "is_Quadrupole", thisProperty)
224 <       PropertyMap(i)%is_Quadrupole = thisProperty
223 >          if (i_is_StickyP .and. j_is_StickyP) then
224 >             iHash = ior(iHash, STICKYPOWER_PAIR)
225 >          endif
226  
227 <       call getElementProperty(atypes, i, "is_Sticky", thisProperty)
228 <       PropertyMap(i)%is_Sticky = thisProperty
229 <      
203 <       call getElementProperty(atypes, i, "is_StickyPower", thisProperty)
204 <       PropertyMap(i)%is_StickyPower = thisProperty
205 <
206 <       call getElementProperty(atypes, i, "is_GayBerne", thisProperty)
207 <       PropertyMap(i)%is_GayBerne = thisProperty
227 >          if (i_is_EAM .and. j_is_EAM) then
228 >             iHash = ior(iHash, EAM_PAIR)
229 >          endif
230  
231 <       call getElementProperty(atypes, i, "is_EAM", thisProperty)
232 <       PropertyMap(i)%is_EAM = thisProperty
231 >          if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
232 >          if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
233 >          if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
234  
235 <       call getElementProperty(atypes, i, "is_Shape", thisProperty)
236 <       PropertyMap(i)%is_Shape = thisProperty
235 >          if (i_is_Shape .and. j_is_Shape) iHash = ior(iHash, SHAPE_PAIR)
236 >          if (i_is_Shape .and. j_is_LJ) iHash = ior(iHash, SHAPE_LJ)
237 >          if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
238  
239 <       call getElementProperty(atypes, i, "is_FLARB", thisProperty)
240 <       PropertyMap(i)%is_FLARB = thisProperty
239 >
240 >          InteractionHash(i,j) = iHash
241 >          InteractionHash(j,i) = iHash
242 >
243 >       end do
244 >
245      end do
246  
247 <    havePropertyMap = .true.
247 >    haveInteractionHash = .true.
248 >  end subroutine createInteractionHash
249  
250 <  end subroutine createPropertyMap
250 >  subroutine createGtypeCutoffMap(stat)
251  
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 +    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, *) 'createInteractionHash failed in doForces!'
273 +          stat = -1
274 +          return
275 +       endif
276 +    endif
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 +
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 +          
327 +          
328 +          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
329 +             biggestAtypeCutoff = atypeMaxCutoff(i)
330 +          endif
331 +
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()
225 –    SIM_uses_LennardJones = SimUsesLennardJones()
226 –    SIM_uses_Electrostatics = SimUsesElectrostatics()
227 –    SIM_uses_Charges = SimUsesCharges()
228 –    SIM_uses_Dipoles = SimUsesDipoles()
229 –    SIM_uses_Sticky = SimUsesSticky()
230 –    SIM_uses_StickyPower = SimUsesStickyPower()
231 –    SIM_uses_GayBerne = SimUsesGayBerne()
453      SIM_uses_EAM = SimUsesEAM()
233 –    SIM_uses_Shapes = SimUsesShapes()
234 –    SIM_uses_FLARB = SimUsesFLARB()
235 –    SIM_uses_RF = SimUsesRF()
454      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
455      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
456      SIM_uses_PBC = SimUsesPBC()
# Line 249 | Line 467 | contains
467  
468      error = 0
469  
470 <    if (.not. havePropertyMap) 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 createPropertyMap(myStatus)
257 <
480 >    if (.not. haveGtypeCutoffMap) then        
481 >       myStatus = 0      
482 >       call createGtypeCutoffMap(myStatus)      
483         if (myStatus .ne. 0) then
484 <          write(default_error, *) 'createPropertyMap failed in doForces!'
484 >          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
485            error = -1
486            return
487         endif
# Line 266 | 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 295 | 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
301 <
525 >    integer, intent(in) :: thisESM
526      integer, intent(out) :: thisStat  
527      integer :: my_status, nMatches
528      integer, pointer :: MatchList(:) => null()
# Line 307 | Line 531 | contains
531      !! assume things are copacetic, unless they aren't
532      thisStat = 0
533  
534 <    !! Fortran's version of a cast:
311 <    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 317 | Line 540 | contains
540      !! interactions are used by the force field.    
541  
542      FF_uses_DirectionalAtoms = .false.
320 –    FF_uses_LennardJones = .false.
321 –    FF_uses_Electrostatics = .false.
322 –    FF_uses_Charges = .false.    
543      FF_uses_Dipoles = .false.
324 –    FF_uses_Sticky = .false.
325 –    FF_uses_StickyPower = .false.
544      FF_uses_GayBerne = .false.
545      FF_uses_EAM = .false.
328 –    FF_uses_Shapes = .false.
329 –    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  
335 –    call getMatchingElementList(atypes, "is_LennardJones", .true., &
336 –         nMatches, MatchList)
337 –    if (nMatches .gt. 0) FF_uses_LennardJones = .true.
338 –
339 –    call getMatchingElementList(atypes, "is_Electrostatic", .true., &
340 –         nMatches, MatchList)
341 –    if (nMatches .gt. 0) then
342 –       FF_uses_Electrostatics = .true.
343 –    endif
344 –
345 –    call getMatchingElementList(atypes, "is_Charge", .true., &
346 –         nMatches, MatchList)
347 –    if (nMatches .gt. 0) then
348 –       FF_uses_Charges = .true.  
349 –       FF_uses_Electrostatics = .true.
350 –    endif
351 –
551      call getMatchingElementList(atypes, "is_Dipole", .true., &
552           nMatches, MatchList)
553 <    if (nMatches .gt. 0) then
355 <       FF_uses_Dipoles = .true.
356 <       FF_uses_Electrostatics = .true.
357 <       FF_uses_DirectionalAtoms = .true.
358 <    endif
359 <
360 <    call getMatchingElementList(atypes, "is_Quadrupole", .true., &
361 <         nMatches, MatchList)
362 <    if (nMatches .gt. 0) then
363 <       FF_uses_Quadrupoles = .true.
364 <       FF_uses_Electrostatics = .true.
365 <       FF_uses_DirectionalAtoms = .true.
366 <    endif
367 <
368 <    call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
369 <         MatchList)
370 <    if (nMatches .gt. 0) then
371 <       FF_uses_Sticky = .true.
372 <       FF_uses_DirectionalAtoms = .true.
373 <    endif
374 <
375 <    call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
376 <         MatchList)
377 <    if (nMatches .gt. 0) then
378 <       FF_uses_StickyPower = .true.
379 <       FF_uses_DirectionalAtoms = .true.
380 <    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
385 <       FF_uses_GayBerne = .true.
386 <       FF_uses_DirectionalAtoms = .true.
387 <    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  
392 –    call getMatchingElementList(atypes, "is_Shape", .true., &
393 –         nMatches, MatchList)
394 –    if (nMatches .gt. 0) then
395 –       FF_uses_Shapes = .true.
396 –       FF_uses_DirectionalAtoms = .true.
397 –    endif
562  
399 –    call getMatchingElementList(atypes, "is_FLARB", .true., &
400 –         nMatches, MatchList)
401 –    if (nMatches .gt. 0) FF_uses_FLARB = .true.
402 –
403 –    !! 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.
# Line 419 | Line 578 | contains
578         endif
579      endif
580  
422 –    !sticky module does not contain check_sticky_FF anymore
423 –    !if (FF_uses_sticky) then
424 –    !   call check_sticky_FF(my_status)
425 –    !   if (my_status /= 0) then
426 –    !      thisStat = -1
427 –    !      haveSaneForceField = .false.
428 –    !      return
429 –    !   end if
430 –    !endif
431 –
581      if (FF_uses_EAM) then
582         call init_EAM_FF(my_status)
583         if (my_status /= 0) then
# Line 446 | Line 595 | contains
595            haveSaneForceField = .false.
596            return
597         endif
449 –    endif
450 –
451 –    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
598      endif
599  
600      if (.not. haveNeighborList) then
# Line 514 | Line 660 | contains
660      integer :: localError
661      integer :: propPack_i, propPack_j
662      integer :: loopStart, loopEnd, loop
663 +    integer :: iHash
664 +  
665  
518 –    real(kind=dp) :: listSkin = 1.0  
519 –
666      !! initialize local variables  
667  
668   #ifdef IS_MPI
# Line 633 | Line 779 | contains
779               endif
780  
781   #ifdef IS_MPI
782 +             me_j = atid_col(j)
783               call get_interatomic_vector(q_group_Row(:,i), &
784                    q_group_Col(:,j), d_grp, rgrpsq)
785   #else
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 < rlistsq) then
791 >             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
792                  if (update_nlist) then
793                     nlist = nlist + 1
794  
# Line 840 | 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 858 | Line 1006 | contains
1006   #else
1007               me_i = atid(i)
1008   #endif
1009 +             iHash = InteractionHash(me_i,me_j)
1010 +            
1011 +             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1012  
862 –             if (PropertyMap(me_i)%is_Dipole) then
863 –
1013                  mu_i = getDipoleMoment(me_i)
1014  
1015                  !! The reaction field needs to include a self contribution
# Line 923 | Line 1072 | contains
1072      real ( kind = dp ), intent(inout) :: rijsq
1073      real ( kind = dp )                :: r
1074      real ( kind = dp ), intent(inout) :: d(3)
926 –    real ( kind = dp ) :: ebalance
1075      integer :: me_i, me_j
1076 +
1077 +    integer :: iHash
1078  
1079      r = sqrt(rijsq)
1080      vpair = 0.0d0
# Line 938 | Line 1088 | contains
1088      me_j = atid(j)
1089   #endif
1090  
1091 <    !    write(*,*) i, j, me_i, me_j
1091 >    iHash = InteractionHash(me_i, me_j)
1092  
1093 <    if (FF_uses_LennardJones .and. SIM_uses_LennardJones) then
1094 <
945 <       if ( PropertyMap(me_i)%is_LennardJones .and. &
946 <            PropertyMap(me_j)%is_LennardJones ) then
947 <          call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
948 <       endif
949 <
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 (FF_uses_Electrostatics .and. SIM_uses_Electrostatics) then
1098 <
1099 <       if (PropertyMap(me_i)%is_Electrostatic .and. &
955 <            PropertyMap(me_j)%is_Electrostatic) then
956 <          call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
957 <               pot, eFrame, f, t, do_pot, ebalance)
958 <       endif
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_dipoles .and. SIM_uses_dipoles) then      
1102 <          if ( PropertyMap(me_i)%is_Dipole .and. &
1103 <               PropertyMap(me_j)%is_Dipole) then
1104 <             if (FF_uses_RF .and. SIM_uses_RF) then
1105 <                call accumulate_rf(i, j, r, eFrame, sw)
965 <                call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
966 <             endif
967 <          endif
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)
1105 >          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1106         endif
1107 +
1108      endif
1109  
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 (FF_uses_Sticky .and. SIM_uses_sticky) then
1116 <
1117 <       if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
975 <          call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
976 <               pot, A, f, t, do_pot)
977 <       endif
978 <
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 (FF_uses_StickyPower .and. SIM_uses_stickypower) then
1121 <       if ( PropertyMap(me_i)%is_StickyPower .and. &
1122 <            PropertyMap(me_j)%is_StickyPower) then
984 <          call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
985 <               pot, A, f, t, do_pot, ebalance)
986 <       endif
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 (FF_uses_GayBerne .and. SIM_uses_GayBerne) then
1126 <
1127 <       if ( PropertyMap(me_i)%is_GayBerne .and. &
992 <            PropertyMap(me_j)%is_GayBerne) then
993 <          call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
994 <               pot, A, f, t, do_pot)
995 <       endif
996 <
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 (FF_uses_EAM .and. SIM_uses_EAM) then
1131 <
1132 <       if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
1002 <          call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1003 <               do_pot)
1004 <       endif
1005 <
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 <
1136 <    !    write(*,*) PropertyMap(me_i)%is_Shape,PropertyMap(me_j)%is_Shape
1137 <
1011 <    if (FF_uses_Shapes .and. SIM_uses_Shapes) then
1012 <       if ( PropertyMap(me_i)%is_Shape .and. &
1013 <            PropertyMap(me_j)%is_Shape ) then
1014 <          call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1015 <               pot, A, f, t, do_pot)
1016 <       endif
1017 <       if ( (PropertyMap(me_i)%is_Shape .and. &
1018 <            PropertyMap(me_j)%is_LennardJones) .or. &
1019 <            (PropertyMap(me_i)%is_LennardJones .and. &
1020 <            PropertyMap(me_j)%is_Shape) ) then
1021 <          call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1022 <               pot, A, f, t, do_pot)
1023 <       endif
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(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
1144 +    
1145    end subroutine do_pair
1146  
1147    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
# Line 1040 | Line 1159 | contains
1159      real ( kind = dp )                :: r, rc
1160      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1161  
1162 <    logical :: is_EAM_i, is_EAM_j
1162 >    integer :: me_i, me_j, iHash
1163  
1045 –    integer :: me_i, me_j
1046 –
1047 –
1164      r = sqrt(rijsq)
1049 –    if (SIM_uses_molecular_cutoffs) then
1050 –       rc = sqrt(rcijsq)
1051 –    else
1052 –       rc = r
1053 –    endif
1165  
1055 –
1166   #ifdef IS_MPI  
1167      me_i = atid_row(i)
1168      me_j = atid_col(j)  
# Line 1061 | Line 1171 | contains
1171      me_j = atid(j)  
1172   #endif
1173  
1174 <    if (FF_uses_EAM .and. SIM_uses_EAM) then
1174 >    iHash = InteractionHash(me_i, me_j)
1175  
1176 <       if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1176 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1177              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1068 –
1178      endif
1179 <
1179 >    
1180    end subroutine do_prepair
1181  
1182  
# Line 1261 | Line 1370 | contains
1370  
1371    function FF_UsesDirectionalAtoms() result(doesit)
1372      logical :: doesit
1373 <    doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1265 <         FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1266 <         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 1273 | 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