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

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2259 by gezelter, Mon Jun 27 21:01:36 2005 UTC vs.
Revision 2277 by chrisfen, Fri Aug 26 21:30:41 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.20 2005-06-27 21:01:30 gezelter Exp $, $Date: 2005-06-27 21:01:30 $, $Name: not supported by cvs2svn $, $Revision: 1.20 $
48 > !! @version $Id: doForces.F90,v 1.32 2005-08-26 21:30:30 chrisfen Exp $, $Date: 2005-08-26 21:30:30 $, $Name: not supported by cvs2svn $, $Revision: 1.32 $
49  
50  
51   module doForces
# Line 73 | Line 73 | module doForces
73  
74   #define __FORTRAN90
75   #include "UseTheForce/fSwitchingFunction.h"
76 + #include "UseTheForce/fCutoffPolicy.h"
77   #include "UseTheForce/DarkSide/fInteractionMap.h"
78  
79 +
80    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
81    INTEGER, PARAMETER:: PAIR_LOOP    = 2
82  
81  logical, save :: haveRlist = .false.
83    logical, save :: haveNeighborList = .false.
84    logical, save :: haveSIMvariables = .false.
84  logical, save :: havePropertyMap = .false.
85    logical, save :: haveSaneForceField = .false.
86 +  logical, save :: haveInteractionHash = .false.
87 +  logical, save :: haveGtypeCutoffMap = .false.
88  
89    logical, save :: FF_uses_DirectionalAtoms
88  logical, save :: FF_uses_LennardJones
89  logical, save :: FF_uses_Electrostatics
90  logical, save :: FF_uses_Charges
90    logical, save :: FF_uses_Dipoles
92  logical, save :: FF_uses_Quadrupoles
93  logical, save :: FF_uses_Sticky
94  logical, save :: FF_uses_StickyPower
91    logical, save :: FF_uses_GayBerne
92    logical, save :: FF_uses_EAM
97  logical, save :: FF_uses_Shapes
98  logical, save :: FF_uses_FLARB
93    logical, save :: FF_uses_RF
94  
95    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
96    logical, save :: SIM_uses_EAM
111  logical, save :: SIM_uses_Shapes
112  logical, save :: SIM_uses_FLARB
97    logical, save :: SIM_uses_RF
98    logical, save :: SIM_requires_postpair_calc
99    logical, save :: SIM_requires_prepair_calc
100    logical, save :: SIM_uses_PBC
117  logical, save :: SIM_uses_molecular_cutoffs
101  
119  real(kind=dp), save :: rlist, rlistsq
120
102    public :: init_FF
103 +  public :: setDefaultCutoffs
104    public :: do_force_loop
105 <  public :: setRlistDF
105 >  public :: createInteractionHash
106 >  public :: createGtypeCutoffMap
107 >  public :: getStickyCut
108 >  public :: getStickyPowerCut
109 >  public :: getGayBerneCut
110 >  public :: getEAMCut
111 >  public :: getShapeCut
112  
113   #ifdef PROFILE
114    public :: getforcetime
# Line 128 | Line 116 | module doForces
116    real :: forceTimeInitial, forceTimeFinal
117    integer :: nLoops
118   #endif
119 +  
120 +  !! Variables for cutoff mapping and interaction mapping
121 +  ! Bit hash to determine pair-pair interactions.
122 +  integer, dimension(:,:), allocatable :: InteractionHash
123 +  real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
124 +  real(kind=dp), dimension(:), allocatable :: groupMaxCutoff
125 +  integer, dimension(:), allocatable :: groupToGtype
126 +  real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff
127 +  type ::gtypeCutoffs
128 +     real(kind=dp) :: rcut
129 +     real(kind=dp) :: rcutsq
130 +     real(kind=dp) :: rlistsq
131 +  end type gtypeCutoffs
132 +  type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
133  
134 <  type :: Properties
135 <     logical :: is_Directional   = .false.
136 <     logical :: is_LennardJones  = .false.
135 <     logical :: is_Electrostatic = .false.
136 <     logical :: is_Charge        = .false.
137 <     logical :: is_Dipole        = .false.
138 <     logical :: is_Quadrupole    = .false.
139 <     logical :: is_Sticky        = .false.
140 <     logical :: is_StickyPower   = .false.
141 <     logical :: is_GayBerne      = .false.
142 <     logical :: is_EAM           = .false.
143 <     logical :: is_Shape         = .false.
144 <     logical :: is_FLARB         = .false.
145 <  end type Properties
146 <
147 <  type(Properties), dimension(:),allocatable :: PropertyMap
148 <
134 >  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
135 >  real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
136 >  
137   contains
138  
139 <  subroutine setRlistDF( this_rlist )
152 <
153 <    real(kind=dp) :: this_rlist
154 <
155 <    rlist = this_rlist
156 <    rlistsq = rlist * rlist
157 <
158 <    haveRlist = .true.
159 <
160 <  end subroutine setRlistDF
161 <
162 <  subroutine createPropertyMap(status)
139 >  subroutine createInteractionHash(status)
140      integer :: nAtypes
141 <    integer :: status
141 >    integer, intent(out) :: status
142      integer :: i
143 <    logical :: thisProperty
144 <    real (kind=DP) :: thisDPproperty
143 >    integer :: j
144 >    integer :: iHash
145 >    !! Test Types
146 >    logical :: i_is_LJ
147 >    logical :: i_is_Elect
148 >    logical :: i_is_Sticky
149 >    logical :: i_is_StickyP
150 >    logical :: i_is_GB
151 >    logical :: i_is_EAM
152 >    logical :: i_is_Shape
153 >    logical :: j_is_LJ
154 >    logical :: j_is_Elect
155 >    logical :: j_is_Sticky
156 >    logical :: j_is_StickyP
157 >    logical :: j_is_GB
158 >    logical :: j_is_EAM
159 >    logical :: j_is_Shape
160 >    real(kind=dp) :: myRcut
161  
162 <    status = 0
162 >    status = 0  
163  
164 +    if (.not. associated(atypes)) then
165 +       call handleError("atype", "atypes was not present before call of createInteractionHash!")
166 +       status = -1
167 +       return
168 +    endif
169 +    
170      nAtypes = getSize(atypes)
171 <
171 >    
172      if (nAtypes == 0) then
173         status = -1
174         return
175      end if
176  
177 <    if (.not. allocated(PropertyMap)) then
178 <       allocate(PropertyMap(nAtypes))
177 >    if (.not. allocated(InteractionHash)) then
178 >       allocate(InteractionHash(nAtypes,nAtypes))
179      endif
180  
181 +    if (.not. allocated(atypeMaxCutoff)) then
182 +       allocate(atypeMaxCutoff(nAtypes))
183 +    endif
184 +        
185      do i = 1, nAtypes
186 <       call getElementProperty(atypes, i, "is_Directional", thisProperty)
187 <       PropertyMap(i)%is_Directional = thisProperty
188 <
189 <       call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
190 <       PropertyMap(i)%is_LennardJones = thisProperty
191 <
192 <       call getElementProperty(atypes, i, "is_Electrostatic", thisProperty)
190 <       PropertyMap(i)%is_Electrostatic = thisProperty
186 >       call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
187 >       call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
188 >       call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
189 >       call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
190 >       call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
191 >       call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
192 >       call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
193  
194 <       call getElementProperty(atypes, i, "is_Charge", thisProperty)
193 <       PropertyMap(i)%is_Charge = thisProperty
194 >       do j = i, nAtypes
195  
196 <       call getElementProperty(atypes, i, "is_Dipole", thisProperty)
197 <       PropertyMap(i)%is_Dipole = thisProperty
196 >          iHash = 0
197 >          myRcut = 0.0_dp
198  
199 <       call getElementProperty(atypes, i, "is_Quadrupole", thisProperty)
200 <       PropertyMap(i)%is_Quadrupole = thisProperty
199 >          call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ)
200 >          call getElementProperty(atypes, j, "is_Electrostatic", j_is_Elect)
201 >          call getElementProperty(atypes, j, "is_Sticky", j_is_Sticky)
202 >          call getElementProperty(atypes, j, "is_StickyPower", j_is_StickyP)
203 >          call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
204 >          call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
205 >          call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
206  
207 <       call getElementProperty(atypes, i, "is_Sticky", thisProperty)
208 <       PropertyMap(i)%is_Sticky = thisProperty
209 <      
210 <       call getElementProperty(atypes, i, "is_StickyPower", thisProperty)
211 <       PropertyMap(i)%is_StickyPower = thisProperty
207 >          if (i_is_LJ .and. j_is_LJ) then
208 >             iHash = ior(iHash, LJ_PAIR)            
209 >          endif
210 >          
211 >          if (i_is_Elect .and. j_is_Elect) then
212 >             iHash = ior(iHash, ELECTROSTATIC_PAIR)
213 >          endif
214 >          
215 >          if (i_is_Sticky .and. j_is_Sticky) then
216 >             iHash = ior(iHash, STICKY_PAIR)
217 >          endif
218  
219 <       call getElementProperty(atypes, i, "is_GayBerne", thisProperty)
220 <       PropertyMap(i)%is_GayBerne = thisProperty
219 >          if (i_is_StickyP .and. j_is_StickyP) then
220 >             iHash = ior(iHash, STICKYPOWER_PAIR)
221 >          endif
222  
223 <       call getElementProperty(atypes, i, "is_EAM", thisProperty)
224 <       PropertyMap(i)%is_EAM = thisProperty
223 >          if (i_is_EAM .and. j_is_EAM) then
224 >             iHash = ior(iHash, EAM_PAIR)
225 >          endif
226  
227 <       call getElementProperty(atypes, i, "is_Shape", thisProperty)
228 <       PropertyMap(i)%is_Shape = thisProperty
227 >          if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
228 >          if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
229 >          if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
230  
231 <       call getElementProperty(atypes, i, "is_FLARB", thisProperty)
232 <       PropertyMap(i)%is_FLARB = thisProperty
231 >          if (i_is_Shape .and. j_is_Shape) iHash = ior(iHash, SHAPE_PAIR)
232 >          if (i_is_Shape .and. j_is_LJ) iHash = ior(iHash, SHAPE_LJ)
233 >          if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
234 >
235 >
236 >          InteractionHash(i,j) = iHash
237 >          InteractionHash(j,i) = iHash
238 >
239 >       end do
240 >
241      end do
242  
243 <    havePropertyMap = .true.
243 >    haveInteractionHash = .true.
244 >  end subroutine createInteractionHash
245  
246 <  end subroutine createPropertyMap
246 >  subroutine createGtypeCutoffMap(stat)
247 >
248 >    integer, intent(out), optional :: stat
249 >    logical :: i_is_LJ
250 >    logical :: i_is_Elect
251 >    logical :: i_is_Sticky
252 >    logical :: i_is_StickyP
253 >    logical :: i_is_GB
254 >    logical :: i_is_EAM
255 >    logical :: i_is_Shape
256 >
257 >    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
258 >    integer :: n_in_i
259 >    real(kind=dp):: thisSigma, bigSigma, thisRcut
260 >    real(kind=dp) :: biggestAtypeCutoff
261 >
262 >    stat = 0
263 >    if (.not. haveInteractionHash) then
264 >       call createInteractionHash(myStatus)      
265 >       if (myStatus .ne. 0) then
266 >          write(default_error, *) 'createInteractionHash failed in doForces!'
267 >          stat = -1
268 >          return
269 >       endif
270 >    endif
271 >
272 >    nAtypes = getSize(atypes)
273 >    
274 >    do i = 1, nAtypes
275 >       if (SimHasAtype(i)) then          
276 >          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
277 >          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
278 >          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
279 >          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
280 >          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
281 >          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
282 >          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
283 >          
284 >          if (i_is_LJ) then
285 >             thisRcut = getSigma(i) * 2.5_dp
286 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
287 >          endif
288 >          if (i_is_Elect) then
289 >             thisRcut = defaultRcut
290 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
291 >          endif
292 >          if (i_is_Sticky) then
293 >             thisRcut = getStickyCut(i)
294 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
295 >          endif
296 >          if (i_is_StickyP) then
297 >             thisRcut = getStickyPowerCut(i)
298 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
299 >          endif
300 >          if (i_is_GB) then
301 >             thisRcut = getGayBerneCut(i)
302 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
303 >          endif
304 >          if (i_is_EAM) then
305 >             thisRcut = getEAMCut(i)
306 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
307 >          endif
308 >          if (i_is_Shape) then
309 >             thisRcut = getShapeCut(i)
310 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
311 >          endif
312 >          
313 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
314 >             biggestAtypeCutoff = atypeMaxCutoff(i)
315 >          endif
316 >       endif
317 >    enddo
318  
319 +    istart = 1
320 + #ifdef IS_MPI
321 +    iend = nGroupsInRow
322 + #else
323 +    iend = nGroups
324 + #endif
325 +    outer: do i = istart, iend
326 +      
327 +       n_in_i = groupStartRow(i+1) - groupStartRow(i)
328 +      
329 + #ifdef IS_MPI
330 +       jstart = 1
331 +       jend = nGroupsInCol
332 + #else
333 +       jstart = i+1
334 +       jend = nGroups
335 + #endif
336 +      
337 +      
338 +      
339 +      
340 +      
341 +      
342 +    enddo outer        
343 +    
344 +     haveGtypeCutoffMap = .true.
345 +   end subroutine createGtypeCutoffMap
346 +
347 +   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
348 +     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
349 +     integer, intent(in) :: cutPolicy
350 +
351 +     defaultRcut = defRcut
352 +     defaultRsw = defRsw
353 +     defaultRlist = defRlist
354 +     cutoffPolicy = cutPolicy
355 +   end subroutine setDefaultCutoffs
356 +
357 +   subroutine setCutoffPolicy(cutPolicy)
358 +
359 +     integer, intent(in) :: cutPolicy
360 +     cutoffPolicy = cutPolicy
361 +     call createGtypeCutoffMap()
362 +
363 +   end subroutine setCutoffPolicy
364 +    
365 +    
366    subroutine setSimVariables()
367      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
226    SIM_uses_LennardJones = SimUsesLennardJones()
227    SIM_uses_Electrostatics = SimUsesElectrostatics()
228    SIM_uses_Charges = SimUsesCharges()
229    SIM_uses_Dipoles = SimUsesDipoles()
230    SIM_uses_Sticky = SimUsesSticky()
231    SIM_uses_StickyPower = SimUsesStickyPower()
232    SIM_uses_GayBerne = SimUsesGayBerne()
368      SIM_uses_EAM = SimUsesEAM()
234    SIM_uses_Shapes = SimUsesShapes()
235    SIM_uses_FLARB = SimUsesFLARB()
369      SIM_uses_RF = SimUsesRF()
370      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
371      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
# Line 250 | Line 383 | contains
383  
384      error = 0
385  
386 <    if (.not. havePropertyMap) then
386 >    if (.not. haveInteractionHash) then      
387 >       myStatus = 0      
388 >       call createInteractionHash(myStatus)      
389 >       if (myStatus .ne. 0) then
390 >          write(default_error, *) 'createInteractionHash failed in doForces!'
391 >          error = -1
392 >          return
393 >       endif
394 >    endif
395  
396 <       myStatus = 0
397 <
398 <       call createPropertyMap(myStatus)
258 <
396 >    if (.not. haveGtypeCutoffMap) then        
397 >       myStatus = 0      
398 >       call createGtypeCutoffMap(myStatus)      
399         if (myStatus .ne. 0) then
400 <          write(default_error, *) 'createPropertyMap failed in doForces!'
400 >          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
401            error = -1
402            return
403         endif
# Line 318 | Line 458 | contains
458      !! interactions are used by the force field.    
459  
460      FF_uses_DirectionalAtoms = .false.
321    FF_uses_LennardJones = .false.
322    FF_uses_Electrostatics = .false.
323    FF_uses_Charges = .false.    
461      FF_uses_Dipoles = .false.
325    FF_uses_Sticky = .false.
326    FF_uses_StickyPower = .false.
462      FF_uses_GayBerne = .false.
463      FF_uses_EAM = .false.
329    FF_uses_Shapes = .false.
330    FF_uses_FLARB = .false.
464  
465      call getMatchingElementList(atypes, "is_Directional", .true., &
466           nMatches, MatchList)
467      if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
468  
336    call getMatchingElementList(atypes, "is_LennardJones", .true., &
337         nMatches, MatchList)
338    if (nMatches .gt. 0) FF_uses_LennardJones = .true.
339
340    call getMatchingElementList(atypes, "is_Electrostatic", .true., &
341         nMatches, MatchList)
342    if (nMatches .gt. 0) then
343       FF_uses_Electrostatics = .true.
344    endif
345
346    call getMatchingElementList(atypes, "is_Charge", .true., &
347         nMatches, MatchList)
348    if (nMatches .gt. 0) then
349       FF_uses_Charges = .true.  
350       FF_uses_Electrostatics = .true.
351    endif
352
469      call getMatchingElementList(atypes, "is_Dipole", .true., &
470           nMatches, MatchList)
471 <    if (nMatches .gt. 0) then
356 <       FF_uses_Dipoles = .true.
357 <       FF_uses_Electrostatics = .true.
358 <       FF_uses_DirectionalAtoms = .true.
359 <    endif
360 <
361 <    call getMatchingElementList(atypes, "is_Quadrupole", .true., &
362 <         nMatches, MatchList)
363 <    if (nMatches .gt. 0) then
364 <       FF_uses_Quadrupoles = .true.
365 <       FF_uses_Electrostatics = .true.
366 <       FF_uses_DirectionalAtoms = .true.
367 <    endif
368 <
369 <    call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
370 <         MatchList)
371 <    if (nMatches .gt. 0) then
372 <       FF_uses_Sticky = .true.
373 <       FF_uses_DirectionalAtoms = .true.
374 <    endif
375 <
376 <    call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
377 <         MatchList)
378 <    if (nMatches .gt. 0) then
379 <       FF_uses_StickyPower = .true.
380 <       FF_uses_DirectionalAtoms = .true.
381 <    endif
471 >    if (nMatches .gt. 0) FF_uses_Dipoles = .true.
472      
473      call getMatchingElementList(atypes, "is_GayBerne", .true., &
474           nMatches, MatchList)
475 <    if (nMatches .gt. 0) then
386 <       FF_uses_GayBerne = .true.
387 <       FF_uses_DirectionalAtoms = .true.
388 <    endif
475 >    if (nMatches .gt. 0) FF_uses_GayBerne = .true.
476  
477      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
478      if (nMatches .gt. 0) FF_uses_EAM = .true.
392
393    call getMatchingElementList(atypes, "is_Shape", .true., &
394         nMatches, MatchList)
395    if (nMatches .gt. 0) then
396       FF_uses_Shapes = .true.
397       FF_uses_DirectionalAtoms = .true.
398    endif
479  
400    call getMatchingElementList(atypes, "is_FLARB", .true., &
401         nMatches, MatchList)
402    if (nMatches .gt. 0) FF_uses_FLARB = .true.
480  
404    !! Assume sanity (for the sake of argument)
481      haveSaneForceField = .true.
482  
483      !! check to make sure the FF_uses_RF setting makes sense
484  
485 <    if (FF_uses_dipoles) then
485 >    if (FF_uses_Dipoles) then
486         if (FF_uses_RF) then
487            dielect = getDielect()
488            call initialize_rf(dielect)
# Line 420 | Line 496 | contains
496         endif
497      endif
498  
423    !sticky module does not contain check_sticky_FF anymore
424    !if (FF_uses_sticky) then
425    !   call check_sticky_FF(my_status)
426    !   if (my_status /= 0) then
427    !      thisStat = -1
428    !      haveSaneForceField = .false.
429    !      return
430    !   end if
431    !endif
432
499      if (FF_uses_EAM) then
500         call init_EAM_FF(my_status)
501         if (my_status /= 0) then
# Line 449 | Line 515 | contains
515         endif
516      endif
517  
452    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
453    endif
454
518      if (.not. haveNeighborList) then
519         !! Create neighbor lists
520         call expandNeighborList(nLocal, my_status)
# Line 515 | Line 578 | contains
578      integer :: localError
579      integer :: propPack_i, propPack_j
580      integer :: loopStart, loopEnd, loop
581 <
581 >    integer :: iHash
582      real(kind=dp) :: listSkin = 1.0  
583  
584      !! initialize local variables  
# Line 634 | Line 697 | contains
697               endif
698  
699   #ifdef IS_MPI
700 +             me_j = atid_col(j)
701               call get_interatomic_vector(q_group_Row(:,i), &
702                    q_group_Col(:,j), d_grp, rgrpsq)
703   #else
704 +             me_j = atid(j)
705               call get_interatomic_vector(q_group(:,i), &
706                    q_group(:,j), d_grp, rgrpsq)
707   #endif
708  
709 <             if (rgrpsq < rlistsq) then
709 >             if (rgrpsq < InteractionHash(me_i,me_j)%rListsq) then
710                  if (update_nlist) then
711                     nlist = nlist + 1
712  
# Line 859 | Line 924 | contains
924   #else
925               me_i = atid(i)
926   #endif
927 +             iHash = InteractionHash(me_i,me_j)
928 +            
929 +             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
930  
863             if (PropertyMap(me_i)%is_Dipole) then
864
931                  mu_i = getDipoleMoment(me_i)
932  
933                  !! The reaction field needs to include a self contribution
# Line 927 | Line 993 | contains
993      real ( kind = dp ) :: ebalance
994      integer :: me_i, me_j
995  
996 <    integer :: iMap
996 >    integer :: iHash
997  
998      r = sqrt(rijsq)
999      vpair = 0.0d0
# Line 941 | Line 1007 | contains
1007      me_j = atid(j)
1008   #endif
1009  
1010 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1010 >    iHash = InteractionHash(me_i, me_j)
1011  
1012 <    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1012 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
1013         call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1014      endif
1015  
1016 <    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1016 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1017         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1018              pot, eFrame, f, t, do_pot)
1019  
1020 <       if (FF_uses_dipoles .and. SIM_uses_dipoles) then                
1021 <          if ( PropertyMap(me_i)%is_Dipole .and. &
1022 <               PropertyMap(me_j)%is_Dipole) then
1023 <             if (FF_uses_RF .and. SIM_uses_RF) then
1024 <                call accumulate_rf(i, j, r, eFrame, sw)
959 <                call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
960 <             endif
961 <          endif
1020 >       if (FF_uses_RF .and. SIM_uses_RF) then
1021 >
1022 >          ! CHECK ME (RF needs to know about all electrostatic types)
1023 >          call accumulate_rf(i, j, r, eFrame, sw)
1024 >          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1025         endif
1026 +
1027      endif
1028  
1029 <    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1029 >    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1030         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1031              pot, A, f, t, do_pot)
1032      endif
1033  
1034 <    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1034 >    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1035         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1036              pot, A, f, t, do_pot)
1037      endif
1038  
1039 <    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1039 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1040         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1041              pot, A, f, t, do_pot)
1042      endif
1043      
1044 <    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1045 <       call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1046 <            pot, A, f, t, do_pot)
1044 >    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1045 > !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1046 > !           pot, A, f, t, do_pot)
1047      endif
1048  
1049 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1049 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1050         call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1051              do_pot)
1052      endif
1053  
1054 <    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1054 >    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1055         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1056              pot, A, f, t, do_pot)
1057      endif
1058  
1059 <    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1059 >    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1060         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1061              pot, A, f, t, do_pot)
1062      endif
# Line 1014 | Line 1078 | contains
1078      real ( kind = dp )                :: r, rc
1079      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1080  
1081 <    integer :: me_i, me_j, iMap
1081 >    integer :: me_i, me_j, iHash
1082  
1083   #ifdef IS_MPI  
1084      me_i = atid_row(i)
# Line 1024 | Line 1088 | contains
1088      me_j = atid(j)  
1089   #endif
1090  
1091 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1091 >    iHash = InteractionHash(me_i, me_j)
1092  
1093 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1093 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1094              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1095      endif
1096      
# Line 1223 | Line 1287 | contains
1287  
1288    function FF_UsesDirectionalAtoms() result(doesit)
1289      logical :: doesit
1290 <    doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1227 <         FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1228 <         FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1290 >    doesit = FF_uses_DirectionalAtoms
1291    end function FF_UsesDirectionalAtoms
1292  
1293    function FF_RequiresPrepairCalc() result(doesit)

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines