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 2284 by gezelter, Wed Sep 7 19:18:54 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.37 2005-09-07 19:18:54 gezelter Exp $, $Date: 2005-09-07 19:18:54 $, $Name: not supported by cvs2svn $, $Revision: 1.37 $
49  
50  
51   module doForces
# Line 73 | Line 73 | module doForces
73  
74   #define __FORTRAN90
75   #include "UseTheForce/fSwitchingFunction.h"
76 + #include "UseTheForce/fCutoffPolicy.h"
77 + #include "UseTheForce/DarkSide/fInteractionMap.h"
78  
79 +
80    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
81    INTEGER, PARAMETER:: PAIR_LOOP    = 2
82  
80  logical, save :: haveRlist = .false.
83    logical, save :: haveNeighborList = .false.
84    logical, save :: haveSIMvariables = .false.
83  logical, save :: havePropertyMap = .false.
85    logical, save :: haveSaneForceField = .false.
86 +  logical, save :: haveInteractionHash = .false.
87 +  logical, save :: haveGtypeCutoffMap = .false.
88 +  logical, save :: haveRlist = .false.
89  
90    logical, save :: FF_uses_DirectionalAtoms
87  logical, save :: FF_uses_LennardJones
88  logical, save :: FF_uses_Electrostatics
89  logical, save :: FF_uses_Charges
91    logical, save :: FF_uses_Dipoles
91  logical, save :: FF_uses_Quadrupoles
92  logical, save :: FF_uses_Sticky
93  logical, save :: FF_uses_StickyPower
92    logical, save :: FF_uses_GayBerne
93    logical, save :: FF_uses_EAM
96  logical, save :: FF_uses_Shapes
97  logical, save :: FF_uses_FLARB
94    logical, save :: FF_uses_RF
95  
96    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
97    logical, save :: SIM_uses_EAM
110  logical, save :: SIM_uses_Shapes
111  logical, save :: SIM_uses_FLARB
98    logical, save :: SIM_uses_RF
99    logical, save :: SIM_requires_postpair_calc
100    logical, save :: SIM_requires_prepair_calc
101    logical, save :: SIM_uses_PBC
116  logical, save :: SIM_uses_molecular_cutoffs
102  
103 <  real(kind=dp), save :: rlist, rlistsq
103 >  integer, save :: corrMethod
104  
105    public :: init_FF
106 +  public :: setDefaultCutoffs
107    public :: do_force_loop
108 <  public :: setRlistDF
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.
134 <     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 >  
140   contains
141  
142 <  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)
142 >  subroutine createInteractionHash(status)
143      integer :: nAtypes
144 <    integer :: status
144 >    integer, intent(out) :: status
145      integer :: i
146 <    logical :: thisProperty
147 <    real (kind=DP) :: thisDPproperty
146 >    integer :: j
147 >    integer :: iHash
148 >    !! Test Types
149 >    logical :: i_is_LJ
150 >    logical :: i_is_Elect
151 >    logical :: i_is_Sticky
152 >    logical :: i_is_StickyP
153 >    logical :: i_is_GB
154 >    logical :: i_is_EAM
155 >    logical :: i_is_Shape
156 >    logical :: j_is_LJ
157 >    logical :: j_is_Elect
158 >    logical :: j_is_Sticky
159 >    logical :: j_is_StickyP
160 >    logical :: j_is_GB
161 >    logical :: j_is_EAM
162 >    logical :: j_is_Shape
163 >    real(kind=dp) :: myRcut
164  
165 <    status = 0
165 >    status = 0  
166  
167 +    if (.not. associated(atypes)) then
168 +       call handleError("atype", "atypes was not present before call of createInteractionHash!")
169 +       status = -1
170 +       return
171 +    endif
172 +    
173      nAtypes = getSize(atypes)
174 <
174 >    
175      if (nAtypes == 0) then
176         status = -1
177         return
178      end if
179  
180 <    if (.not. allocated(PropertyMap)) then
181 <       allocate(PropertyMap(nAtypes))
180 >    if (.not. allocated(InteractionHash)) then
181 >       allocate(InteractionHash(nAtypes,nAtypes))
182      endif
183  
184 +    if (.not. allocated(atypeMaxCutoff)) then
185 +       allocate(atypeMaxCutoff(nAtypes))
186 +    endif
187 +        
188      do i = 1, nAtypes
189 <       call getElementProperty(atypes, i, "is_Directional", thisProperty)
190 <       PropertyMap(i)%is_Directional = thisProperty
189 >       call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
190 >       call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
191 >       call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
192 >       call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
193 >       call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
194 >       call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
195 >       call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
196  
197 <       call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
186 <       PropertyMap(i)%is_LennardJones = thisProperty
197 >       do j = i, nAtypes
198  
199 <       call getElementProperty(atypes, i, "is_Electrostatic", thisProperty)
200 <       PropertyMap(i)%is_Electrostatic = thisProperty
199 >          iHash = 0
200 >          myRcut = 0.0_dp
201  
202 <       call getElementProperty(atypes, i, "is_Charge", thisProperty)
203 <       PropertyMap(i)%is_Charge = thisProperty
202 >          call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ)
203 >          call getElementProperty(atypes, j, "is_Electrostatic", j_is_Elect)
204 >          call getElementProperty(atypes, j, "is_Sticky", j_is_Sticky)
205 >          call getElementProperty(atypes, j, "is_StickyPower", j_is_StickyP)
206 >          call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
207 >          call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
208 >          call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
209  
210 <       call getElementProperty(atypes, i, "is_Dipole", thisProperty)
211 <       PropertyMap(i)%is_Dipole = thisProperty
210 >          if (i_is_LJ .and. j_is_LJ) then
211 >             iHash = ior(iHash, LJ_PAIR)            
212 >          endif
213 >          
214 >          if (i_is_Elect .and. j_is_Elect) then
215 >             iHash = ior(iHash, ELECTROSTATIC_PAIR)
216 >          endif
217 >          
218 >          if (i_is_Sticky .and. j_is_Sticky) then
219 >             iHash = ior(iHash, STICKY_PAIR)
220 >          endif
221  
222 <       call getElementProperty(atypes, i, "is_Quadrupole", thisProperty)
223 <       PropertyMap(i)%is_Quadrupole = thisProperty
222 >          if (i_is_StickyP .and. j_is_StickyP) then
223 >             iHash = ior(iHash, STICKYPOWER_PAIR)
224 >          endif
225  
226 <       call getElementProperty(atypes, i, "is_Sticky", thisProperty)
227 <       PropertyMap(i)%is_Sticky = thisProperty
228 <      
203 <       call getElementProperty(atypes, i, "is_StickyPower", thisProperty)
204 <       PropertyMap(i)%is_StickyPower = thisProperty
226 >          if (i_is_EAM .and. j_is_EAM) then
227 >             iHash = ior(iHash, EAM_PAIR)
228 >          endif
229  
230 <       call getElementProperty(atypes, i, "is_GayBerne", thisProperty)
231 <       PropertyMap(i)%is_GayBerne = thisProperty
230 >          if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
231 >          if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
232 >          if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
233  
234 <       call getElementProperty(atypes, i, "is_EAM", thisProperty)
235 <       PropertyMap(i)%is_EAM = thisProperty
234 >          if (i_is_Shape .and. j_is_Shape) iHash = ior(iHash, SHAPE_PAIR)
235 >          if (i_is_Shape .and. j_is_LJ) iHash = ior(iHash, SHAPE_LJ)
236 >          if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
237  
212       call getElementProperty(atypes, i, "is_Shape", thisProperty)
213       PropertyMap(i)%is_Shape = thisProperty
238  
239 <       call getElementProperty(atypes, i, "is_FLARB", thisProperty)
240 <       PropertyMap(i)%is_FLARB = thisProperty
239 >          InteractionHash(i,j) = iHash
240 >          InteractionHash(j,i) = iHash
241 >
242 >       end do
243 >
244      end do
245  
246 <    havePropertyMap = .true.
246 >    haveInteractionHash = .true.
247 >  end subroutine createInteractionHash
248  
249 <  end subroutine createPropertyMap
249 >  subroutine createGtypeCutoffMap(stat)
250 >
251 >    integer, intent(out), optional :: stat
252 >    logical :: i_is_LJ
253 >    logical :: i_is_Elect
254 >    logical :: i_is_Sticky
255 >    logical :: i_is_StickyP
256 >    logical :: i_is_GB
257 >    logical :: i_is_EAM
258 >    logical :: i_is_Shape
259 >
260 >    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
261 >    integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
262 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
263 >    real(kind=dp) :: biggestAtypeCutoff
264 >
265 >    stat = 0
266 >    if (.not. haveInteractionHash) then
267 >       call createInteractionHash(myStatus)      
268 >       if (myStatus .ne. 0) then
269 >          write(default_error, *) 'createInteractionHash failed in doForces!'
270 >          stat = -1
271 >          return
272 >       endif
273 >    endif
274 >
275 >    nAtypes = getSize(atypes)
276 >    
277 >    do i = 1, nAtypes
278 >       if (SimHasAtype(i)) then    
279 >          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
280 >          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
281 >          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
282 >          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
283 >          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
284 >          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
285 >          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
286 >          
287 >          atypeMaxCutoff(i) = 0.0_dp
288 >          if (i_is_LJ) then
289 >             thisRcut = getSigma(i) * 2.5_dp
290 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
291 >          endif
292 >          if (i_is_Elect) then
293 >             thisRcut = defaultRcut
294 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
295 >          endif
296 >          if (i_is_Sticky) then
297 >             thisRcut = getStickyCut(i)
298 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
299 >          endif
300 >          if (i_is_StickyP) then
301 >             thisRcut = getStickyPowerCut(i)
302 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
303 >          endif
304 >          if (i_is_GB) then
305 >             thisRcut = getGayBerneCut(i)
306 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
307 >          endif
308 >          if (i_is_EAM) then
309 >             thisRcut = getEAMCut(i)
310 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
311 >          endif
312 >          if (i_is_Shape) then
313 >             thisRcut = getShapeCut(i)
314 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
315 >          endif
316 >          
317 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
318 >             biggestAtypeCutoff = atypeMaxCutoff(i)
319 >          endif
320 >       endif
321 >    enddo
322 >  
323 >    nGroupTypes = 0
324 >    
325 >    istart = 1
326 > #ifdef IS_MPI
327 >    iend = nGroupsInRow
328 > #else
329 >    iend = nGroups
330 > #endif
331 >    
332 >    !! allocate the groupToGtype and gtypeMaxCutoff here.
333 >    if(.not.allocated(groupToGtype)) then
334 >       allocate(groupToGtype(iend))
335 >       allocate(groupMaxCutoff(iend))
336 >       allocate(gtypeMaxCutoff(iend))
337 >    endif
338 >    !! first we do a single loop over the cutoff groups to find the
339 >    !! largest cutoff for any atypes present in this group.  We also
340 >    !! create gtypes at this point.
341 >    
342 >    tol = 1.0d-6
343 >    
344 >    do i = istart, iend      
345 >       n_in_i = groupStartRow(i+1) - groupStartRow(i)
346 >       groupMaxCutoff(i) = 0.0_dp
347 >       do ia = groupStartRow(i), groupStartRow(i+1)-1
348 >          atom1 = groupListRow(ia)
349 > #ifdef IS_MPI
350 >          me_i = atid_row(atom1)
351 > #else
352 >          me_i = atid(atom1)
353 > #endif          
354 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
355 >             groupMaxCutoff(i)=atypeMaxCutoff(me_i)
356 >          endif
357 >       enddo
358 >       if (nGroupTypes.eq.0) then
359 >          nGroupTypes = nGroupTypes + 1
360 >          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
361 >          groupToGtype(i) = nGroupTypes
362 >       else
363 >          do g = 1, nGroupTypes
364 >             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).gt.tol) then
365 >                nGroupTypes = nGroupTypes + 1
366 >                gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
367 >                groupToGtype(i) = nGroupTypes
368 >             else
369 >                groupToGtype(i) = g
370 >             endif
371 >          enddo
372 >       endif
373 >    enddo
374 >    
375 >    !! allocate the gtypeCutoffMap here.
376 >    allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
377 >    !! then we do a double loop over all the group TYPES to find the cutoff
378 >    !! map between groups of two types
379 >    
380 >    do i = 1, nGroupTypes
381 >       do j = 1, nGroupTypes
382 >      
383 >          select case(cutoffPolicy)
384 >          case(TRADITIONAL_CUTOFF_POLICY)
385 >             thisRcut = maxval(gtypeMaxCutoff)
386 >          case(MIX_CUTOFF_POLICY)
387 >             thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
388 >          case(MAX_CUTOFF_POLICY)
389 >             thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
390 >          case default
391 >             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
392 >             return
393 >          end select
394 >          gtypeCutoffMap(i,j)%rcut = thisRcut
395 >          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
396 >          skin = defaultRlist - defaultRcut
397 >          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
398  
399 +       enddo
400 +    enddo
401 +    
402 +    haveGtypeCutoffMap = .true.
403 +    
404 +  end subroutine createGtypeCutoffMap
405 +  
406 +  subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
407 +    real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
408 +    integer, intent(in) :: cutPolicy
409 +    
410 +    defaultRcut = defRcut
411 +    defaultRsw = defRsw
412 +    defaultRlist = defRlist
413 +    cutoffPolicy = cutPolicy
414 +  end subroutine setDefaultCutoffs
415 +  
416 +  subroutine setCutoffPolicy(cutPolicy)
417 +
418 +     integer, intent(in) :: cutPolicy
419 +     cutoffPolicy = cutPolicy
420 +     call createGtypeCutoffMap()
421 +
422 +   end subroutine setCutoffPolicy
423 +    
424 +    
425    subroutine setSimVariables()
426      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()
427      SIM_uses_EAM = SimUsesEAM()
233    SIM_uses_Shapes = SimUsesShapes()
234    SIM_uses_FLARB = SimUsesFLARB()
428      SIM_uses_RF = SimUsesRF()
429      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
430      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
# Line 249 | Line 442 | contains
442  
443      error = 0
444  
445 <    if (.not. havePropertyMap) then
445 >    if (.not. haveInteractionHash) then      
446 >       myStatus = 0      
447 >       call createInteractionHash(myStatus)      
448 >       if (myStatus .ne. 0) then
449 >          write(default_error, *) 'createInteractionHash failed in doForces!'
450 >          error = -1
451 >          return
452 >       endif
453 >    endif
454  
455 <       myStatus = 0
456 <
457 <       call createPropertyMap(myStatus)
257 <
455 >    if (.not. haveGtypeCutoffMap) then        
456 >       myStatus = 0      
457 >       call createGtypeCutoffMap(myStatus)      
458         if (myStatus .ne. 0) then
459 <          write(default_error, *) 'createPropertyMap failed in doForces!'
459 >          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
460            error = -1
461            return
462         endif
# Line 266 | Line 466 | contains
466         call setSimVariables()
467      endif
468  
469 <    if (.not. haveRlist) then
470 <       write(default_error, *) 'rList has not been set in doForces!'
471 <       error = -1
472 <       return
473 <    endif
469 >  !  if (.not. haveRlist) then
470 >  !     write(default_error, *) 'rList has not been set in doForces!'
471 >  !     error = -1
472 >  !     return
473 >  !  endif
474  
475      if (.not. haveNeighborList) then
476         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 295 | Line 495 | contains
495    end subroutine doReadyCheck
496  
497  
498 <  subroutine init_FF(use_RF_c, thisStat)
498 >  subroutine init_FF(use_RF, use_UW, use_DW, thisStat)
499  
500 <    logical, intent(in) :: use_RF_c
501 <
500 >    logical, intent(in) :: use_RF
501 >    logical, intent(in) :: use_UW
502 >    logical, intent(in) :: use_DW
503      integer, intent(out) :: thisStat  
504      integer :: my_status, nMatches
505 +    integer :: corrMethod
506      integer, pointer :: MatchList(:) => null()
507      real(kind=dp) :: rcut, rrf, rt, dielect
508  
# Line 308 | Line 510 | contains
510      thisStat = 0
511  
512      !! Fortran's version of a cast:
513 <    FF_uses_RF = use_RF_c
513 >    FF_uses_RF = use_RF
514  
515 +    !! set the electrostatic correction method
516 +    if (use_UW) then
517 +       corrMethod = 1
518 +    elseif (use_DW) then
519 +       corrMethod = 2
520 +    else
521 +       corrMethod = 0
522 +    endif
523 +    
524      !! init_FF is called *after* all of the atom types have been
525      !! defined in atype_module using the new_atype subroutine.
526      !!
# Line 317 | Line 528 | contains
528      !! interactions are used by the force field.    
529  
530      FF_uses_DirectionalAtoms = .false.
320    FF_uses_LennardJones = .false.
321    FF_uses_Electrostatics = .false.
322    FF_uses_Charges = .false.    
531      FF_uses_Dipoles = .false.
324    FF_uses_Sticky = .false.
325    FF_uses_StickyPower = .false.
532      FF_uses_GayBerne = .false.
533      FF_uses_EAM = .false.
328    FF_uses_Shapes = .false.
329    FF_uses_FLARB = .false.
534  
535      call getMatchingElementList(atypes, "is_Directional", .true., &
536           nMatches, MatchList)
537      if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
538  
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
539      call getMatchingElementList(atypes, "is_Dipole", .true., &
540           nMatches, MatchList)
541 <    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
541 >    if (nMatches .gt. 0) FF_uses_Dipoles = .true.
542      
543      call getMatchingElementList(atypes, "is_GayBerne", .true., &
544           nMatches, MatchList)
545 <    if (nMatches .gt. 0) then
385 <       FF_uses_GayBerne = .true.
386 <       FF_uses_DirectionalAtoms = .true.
387 <    endif
545 >    if (nMatches .gt. 0) FF_uses_GayBerne = .true.
546  
547      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
548      if (nMatches .gt. 0) FF_uses_EAM = .true.
549  
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
550  
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)
551      haveSaneForceField = .true.
552  
553      !! check to make sure the FF_uses_RF setting makes sense
554  
555 <    if (FF_uses_dipoles) then
555 >    if (FF_uses_Dipoles) then
556         if (FF_uses_RF) then
557            dielect = getDielect()
558            call initialize_rf(dielect)
# Line 418 | Line 565 | contains
565            return
566         endif
567      endif
421
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
568  
569      if (FF_uses_EAM) then
570         call init_EAM_FF(my_status)
# Line 448 | Line 585 | contains
585         endif
586      endif
587  
451    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
452    endif
453
588      if (.not. haveNeighborList) then
589         !! Create neighbor lists
590         call expandNeighborList(nLocal, my_status)
# Line 514 | Line 648 | contains
648      integer :: localError
649      integer :: propPack_i, propPack_j
650      integer :: loopStart, loopEnd, loop
651 <
651 >    integer :: iHash
652      real(kind=dp) :: listSkin = 1.0  
653  
654      !! initialize local variables  
# Line 633 | Line 767 | contains
767               endif
768  
769   #ifdef IS_MPI
770 +             me_j = atid_col(j)
771               call get_interatomic_vector(q_group_Row(:,i), &
772                    q_group_Col(:,j), d_grp, rgrpsq)
773   #else
774 +             me_j = atid(j)
775               call get_interatomic_vector(q_group(:,i), &
776                    q_group(:,j), d_grp, rgrpsq)
777   #endif
778  
779 <             if (rgrpsq < rlistsq) then
779 >             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
780                  if (update_nlist) then
781                     nlist = nlist + 1
782  
# Line 858 | Line 994 | contains
994   #else
995               me_i = atid(i)
996   #endif
997 <
998 <             if (PropertyMap(me_i)%is_Dipole) then
997 >             iHash = InteractionHash(me_i,me_j)
998 >            
999 >             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1000  
1001                  mu_i = getDipoleMoment(me_i)
1002  
# Line 926 | Line 1063 | contains
1063      real ( kind = dp ) :: ebalance
1064      integer :: me_i, me_j
1065  
1066 +    integer :: iHash
1067 +
1068      r = sqrt(rijsq)
1069      vpair = 0.0d0
1070      fpair(1:3) = 0.0d0
# Line 938 | Line 1077 | contains
1077      me_j = atid(j)
1078   #endif
1079  
1080 <    !    write(*,*) i, j, me_i, me_j
1080 >    iHash = InteractionHash(me_i, me_j)
1081  
1082 <    if (FF_uses_LennardJones .and. SIM_uses_LennardJones) then
1083 <
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 <
1082 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1083 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1084      endif
1085  
1086 <    if (FF_uses_Electrostatics .and. SIM_uses_Electrostatics) then
1086 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1087 >       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1088 >            pot, eFrame, f, t, do_pot, corrMethod)
1089  
1090 <       if (PropertyMap(me_i)%is_Electrostatic .and. &
1091 <            PropertyMap(me_j)%is_Electrostatic) then
1092 <          call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1093 <               pot, eFrame, f, t, do_pot, ebalance)
1090 >       if (FF_uses_RF .and. SIM_uses_RF) then
1091 >
1092 >          ! CHECK ME (RF needs to know about all electrostatic types)
1093 >          call accumulate_rf(i, j, r, eFrame, sw)
1094 >          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1095         endif
1096  
960       if (FF_uses_dipoles .and. SIM_uses_dipoles) then      
961          if ( PropertyMap(me_i)%is_Dipole .and. &
962               PropertyMap(me_j)%is_Dipole) then
963             if (FF_uses_RF .and. SIM_uses_RF) then
964                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
968       endif
1097      endif
1098  
1099 <
1100 <    if (FF_uses_Sticky .and. SIM_uses_sticky) then
1101 <
974 <       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 <
1099 >    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1100 >       call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1101 >            pot, A, f, t, do_pot)
1102      endif
1103  
1104 <    if (FF_uses_StickyPower .and. SIM_uses_stickypower) then
1105 <       if ( PropertyMap(me_i)%is_StickyPower .and. &
1106 <            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
1104 >    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1105 >       call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1106 >            pot, A, f, t, do_pot)
1107      endif
988    
989    if (FF_uses_GayBerne .and. SIM_uses_GayBerne) then
1108  
1109 <       if ( PropertyMap(me_i)%is_GayBerne .and. &
1110 <            PropertyMap(me_j)%is_GayBerne) then
1111 <          call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
994 <               pot, A, f, t, do_pot)
995 <       endif
996 <
1109 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1110 >       call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1111 >            pot, A, f, t, do_pot)
1112      endif
1113 +    
1114 +    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1115 + !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1116 + !           pot, A, f, t, do_pot)
1117 +    endif
1118  
1119 <    if (FF_uses_EAM .and. SIM_uses_EAM) then
1120 <
1121 <       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 <
1119 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1120 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1121 >            do_pot)
1122      endif
1123  
1124 <
1125 <    !    write(*,*) PropertyMap(me_i)%is_Shape,PropertyMap(me_j)%is_Shape
1126 <
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
1124 >    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1125 >       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1126 >            pot, A, f, t, do_pot)
1127      endif
1128  
1129 +    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1130 +       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1131 +            pot, A, f, t, do_pot)
1132 +    endif
1133 +    
1134    end subroutine do_pair
1135  
1136    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
# Line 1039 | Line 1147 | contains
1147      real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1148      real ( kind = dp )                :: r, rc
1149      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1042
1043    logical :: is_EAM_i, is_EAM_j
1044
1045    integer :: me_i, me_j
1046
1047
1048    r = sqrt(rijsq)
1049    if (SIM_uses_molecular_cutoffs) then
1050       rc = sqrt(rcijsq)
1051    else
1052       rc = r
1053    endif
1150  
1151 +    integer :: me_i, me_j, iHash
1152  
1153   #ifdef IS_MPI  
1154      me_i = atid_row(i)
# Line 1061 | Line 1158 | contains
1158      me_j = atid(j)  
1159   #endif
1160  
1161 <    if (FF_uses_EAM .and. SIM_uses_EAM) then
1161 >    iHash = InteractionHash(me_i, me_j)
1162  
1163 <       if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1163 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1164              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1068
1165      endif
1166 <
1166 >    
1167    end subroutine do_prepair
1168  
1169  
# Line 1261 | Line 1357 | contains
1357  
1358    function FF_UsesDirectionalAtoms() result(doesit)
1359      logical :: doesit
1360 <    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
1360 >    doesit = FF_uses_DirectionalAtoms
1361    end function FF_UsesDirectionalAtoms
1362  
1363    function FF_RequiresPrepairCalc() result(doesit)

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines