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 2270 by gezelter, Tue Aug 9 22:33:37 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.28 2005-08-09 22:33:37 gezelter Exp $, $Date: 2005-08-09 22:33:37 $, $Name: not supported by cvs2svn $, $Revision: 1.28 $
49  
50  
51   module doForces
# Line 73 | Line 73 | module doForces
73  
74   #define __FORTRAN90
75   #include "UseTheForce/fSwitchingFunction.h"
76 + #include "UseTheForce/DarkSide/fInteractionMap.h"
77  
78    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
79    INTEGER, PARAMETER:: PAIR_LOOP    = 2
80  
80 –  logical, save :: haveRlist = .false.
81    logical, save :: haveNeighborList = .false.
82    logical, save :: haveSIMvariables = .false.
83 –  logical, save :: havePropertyMap = .false.
83    logical, save :: haveSaneForceField = .false.
84 +  logical, save :: haveInteractionHash = .false.
85 +  logical, save :: haveGtypeCutoffMap = .false.
86  
87    logical, save :: FF_uses_DirectionalAtoms
87 –  logical, save :: FF_uses_LennardJones
88 –  logical, save :: FF_uses_Electrostatics
89 –  logical, save :: FF_uses_Charges
88    logical, save :: FF_uses_Dipoles
91 –  logical, save :: FF_uses_Quadrupoles
92 –  logical, save :: FF_uses_Sticky
93 –  logical, save :: FF_uses_StickyPower
89    logical, save :: FF_uses_GayBerne
90    logical, save :: FF_uses_EAM
96 –  logical, save :: FF_uses_Shapes
97 –  logical, save :: FF_uses_FLARB
91    logical, save :: FF_uses_RF
92  
93    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
94    logical, save :: SIM_uses_EAM
110 –  logical, save :: SIM_uses_Shapes
111 –  logical, save :: SIM_uses_FLARB
95    logical, save :: SIM_uses_RF
96    logical, save :: SIM_requires_postpair_calc
97    logical, save :: SIM_requires_prepair_calc
98    logical, save :: SIM_uses_PBC
116 –  logical, save :: SIM_uses_molecular_cutoffs
99  
118 –  real(kind=dp), save :: rlist, rlistsq
119 –
100    public :: init_FF
101    public :: do_force_loop
102 <  public :: setRlistDF
102 >  public :: createInteractionHash
103 >  public :: createGtypeCutoffMap
104  
105   #ifdef PROFILE
106    public :: getforcetime
# Line 127 | Line 108 | module doForces
108    real :: forceTimeInitial, forceTimeFinal
109    integer :: nLoops
110   #endif
111 <
112 <  type :: Properties
113 <     logical :: is_Directional   = .false.
114 <     logical :: is_LennardJones  = .false.
115 <     logical :: is_Electrostatic = .false.
116 <     logical :: is_Charge        = .false.
117 <     logical :: is_Dipole        = .false.
118 <     logical :: is_Quadrupole    = .false.
119 <     logical :: is_Sticky        = .false.
120 <     logical :: is_StickyPower   = .false.
121 <     logical :: is_GayBerne      = .false.
122 <     logical :: is_EAM           = .false.
123 <     logical :: is_Shape         = .false.
124 <     logical :: is_FLARB         = .false.
125 <  end type Properties
145 <
146 <  type(Properties), dimension(:),allocatable :: PropertyMap
147 <
111 >  
112 >  !! Variables for cutoff mapping and interaction mapping
113 >  ! Bit hash to determine pair-pair interactions.
114 >  integer, dimension(:,:), allocatable :: InteractionHash
115 >  real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
116 >  real(kind=dp), dimension(:), allocatable :: groupMaxCutoff
117 >  integer, dimension(:), allocatable :: groupToGtype
118 >  real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff
119 >  type ::gtypeCutoffs
120 >     real(kind=dp) :: rcut
121 >     real(kind=dp) :: rcutsq
122 >     real(kind=dp) :: rlistsq
123 >  end type gtypeCutoffs
124 >  type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
125 >  
126   contains
127  
128 <  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)
128 >  subroutine createInteractionHash(status)
129      integer :: nAtypes
130 <    integer :: status
130 >    integer, intent(out) :: status
131      integer :: i
132 <    logical :: thisProperty
133 <    real (kind=DP) :: thisDPproperty
132 >    integer :: j
133 >    integer :: iHash
134 >    !! Test Types
135 >    logical :: i_is_LJ
136 >    logical :: i_is_Elect
137 >    logical :: i_is_Sticky
138 >    logical :: i_is_StickyP
139 >    logical :: i_is_GB
140 >    logical :: i_is_EAM
141 >    logical :: i_is_Shape
142 >    logical :: j_is_LJ
143 >    logical :: j_is_Elect
144 >    logical :: j_is_Sticky
145 >    logical :: j_is_StickyP
146 >    logical :: j_is_GB
147 >    logical :: j_is_EAM
148 >    logical :: j_is_Shape
149 >    
150 >    status = 0  
151  
152 <    status = 0
153 <
152 >    if (.not. associated(atypes)) then
153 >       call handleError("atype", "atypes was not present before call of createInteractionHash!")
154 >       status = -1
155 >       return
156 >    endif
157 >    
158      nAtypes = getSize(atypes)
159 <
159 >    
160      if (nAtypes == 0) then
161         status = -1
162         return
163      end if
164  
165 <    if (.not. allocated(PropertyMap)) then
166 <       allocate(PropertyMap(nAtypes))
165 >    if (.not. allocated(InteractionHash)) then
166 >       allocate(InteractionHash(nAtypes,nAtypes))
167      endif
168  
169 +    if (.not. allocated(atypeMaxCutoff)) then
170 +       allocate(atypeMaxCutoff(nAtypes))
171 +    endif
172 +        
173      do i = 1, nAtypes
174 <       call getElementProperty(atypes, i, "is_Directional", thisProperty)
175 <       PropertyMap(i)%is_Directional = thisProperty
174 >       call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
175 >       call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
176 >       call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
177 >       call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
178 >       call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
179 >       call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
180 >       call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
181  
182 <       call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
183 <       PropertyMap(i)%is_LennardJones = thisProperty
182 >       if (i_is_LJ) then
183 >          thisCut = getDefaultLJCutoff(i)
184 >          if (thisCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisCut
185 >       endif
186  
187 <       call getElementProperty(atypes, i, "is_Electrostatic", thisProperty)
189 <       PropertyMap(i)%is_Electrostatic = thisProperty
187 >
188  
189 <       call getElementProperty(atypes, i, "is_Charge", thisProperty)
192 <       PropertyMap(i)%is_Charge = thisProperty
189 >       do j = i, nAtypes
190  
191 <       call getElementProperty(atypes, i, "is_Dipole", thisProperty)
192 <       PropertyMap(i)%is_Dipole = thisProperty
191 >          iHash = 0
192 >          myRcut = 0.0_dp
193  
194 <       call getElementProperty(atypes, i, "is_Quadrupole", thisProperty)
195 <       PropertyMap(i)%is_Quadrupole = thisProperty
194 >          call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ)
195 >          call getElementProperty(atypes, j, "is_Electrostatic", j_is_Elect)
196 >          call getElementProperty(atypes, j, "is_Sticky", j_is_Sticky)
197 >          call getElementProperty(atypes, j, "is_StickyPower", j_is_StickyP)
198 >          call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
199 >          call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
200 >          call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
201  
202 <       call getElementProperty(atypes, i, "is_Sticky", thisProperty)
203 <       PropertyMap(i)%is_Sticky = thisProperty
204 <      
205 <       call getElementProperty(atypes, i, "is_StickyPower", thisProperty)
206 <       PropertyMap(i)%is_StickyPower = thisProperty
202 >          if (i_is_LJ .and. j_is_LJ) then
203 >             iHash = ior(iHash, LJ_PAIR)            
204 >          endif
205 >          
206 >          if (i_is_Elect .and. j_is_Elect) then
207 >             iHash = ior(iHash, ELECTROSTATIC_PAIR)
208 >          endif
209 >          
210 >          if (i_is_Sticky .and. j_is_Sticky) then
211 >             iHash = ior(iHash, STICKY_PAIR)
212 >          endif
213  
214 <       call getElementProperty(atypes, i, "is_GayBerne", thisProperty)
215 <       PropertyMap(i)%is_GayBerne = thisProperty
214 >          if (i_is_StickyP .and. j_is_StickyP) then
215 >             iHash = ior(iHash, STICKYPOWER_PAIR)
216 >          endif
217  
218 <       call getElementProperty(atypes, i, "is_EAM", thisProperty)
219 <       PropertyMap(i)%is_EAM = thisProperty
218 >          if (i_is_EAM .and. j_is_EAM) then
219 >             iHash = ior(iHash, EAM_PAIR)
220 >          endif
221  
222 <       call getElementProperty(atypes, i, "is_Shape", thisProperty)
223 <       PropertyMap(i)%is_Shape = thisProperty
222 >          if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
223 >          if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
224 >          if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
225  
226 <       call getElementProperty(atypes, i, "is_FLARB", thisProperty)
227 <       PropertyMap(i)%is_FLARB = thisProperty
226 >          if (i_is_Shape .and. j_is_Shape) iHash = ior(iHash, SHAPE_PAIR)
227 >          if (i_is_Shape .and. j_is_LJ) iHash = ior(iHash, SHAPE_LJ)
228 >          if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
229 >
230 >
231 >          InteractionHash(i,j) = iHash
232 >          InteractionHash(j,i) = iHash
233 >
234 >       end do
235 >
236      end do
237  
238 <    havePropertyMap = .true.
238 >    haveInteractionHash = .true.
239 >  end subroutine createInteractionHash
240  
241 <  end subroutine createPropertyMap
241 >  subroutine createGtypeCutoffMap(defaultRcut, defaultSkinThickness, stat)
242  
243 +    real(kind=dp), intent(in), optional :: defaultRCut, defaultSkinThickness
244 +    integer, intent(out) :: stat
245 +
246 +    integer :: myStatus, nAtypes
247 +
248 +    stat = 0
249 +    if (.not. haveInteractionHash) then
250 +       call createInteractionHash(myStatus)      
251 +       if (myStatus .ne. 0) then
252 +          write(default_error, *) 'createInteractionHash failed in doForces!'
253 +          stat = -1
254 +          return
255 +       endif
256 +    endif
257 +
258 +    nAtypes = getSize(atypes)
259 +
260 +    do i = 1, nAtypes
261 +      
262 +       atypeMaxCutoff(i) =
263 +
264 +    
265 +
266 +
267 +
268 +     haveGtypeCutoffMap = .true.
269 +   end subroutine createGtypeCutoffMap
270 +
271    subroutine setSimVariables()
272      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()
273      SIM_uses_EAM = SimUsesEAM()
233 –    SIM_uses_Shapes = SimUsesShapes()
234 –    SIM_uses_FLARB = SimUsesFLARB()
274      SIM_uses_RF = SimUsesRF()
275      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
276      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
# Line 249 | Line 288 | contains
288  
289      error = 0
290  
291 <    if (.not. havePropertyMap) then
291 >    if (.not. haveInteractionHash) then      
292 >       myStatus = 0      
293 >       call createInteractionHash(myStatus)      
294 >       if (myStatus .ne. 0) then
295 >          write(default_error, *) 'createInteractionHash failed in doForces!'
296 >          error = -1
297 >          return
298 >       endif
299 >    endif
300  
301 <       myStatus = 0
302 <
303 <       call createPropertyMap(myStatus)
257 <
301 >    if (.not. haveGtypeCutoffMap) then        
302 >       myStatus = 0      
303 >       call createGtypeCutoffMap(myStatus)      
304         if (myStatus .ne. 0) then
305 <          write(default_error, *) 'createPropertyMap failed in doForces!'
305 >          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
306            error = -1
307            return
308         endif
# Line 317 | Line 363 | contains
363      !! interactions are used by the force field.    
364  
365      FF_uses_DirectionalAtoms = .false.
320 –    FF_uses_LennardJones = .false.
321 –    FF_uses_Electrostatics = .false.
322 –    FF_uses_Charges = .false.    
366      FF_uses_Dipoles = .false.
324 –    FF_uses_Sticky = .false.
325 –    FF_uses_StickyPower = .false.
367      FF_uses_GayBerne = .false.
368      FF_uses_EAM = .false.
328 –    FF_uses_Shapes = .false.
329 –    FF_uses_FLARB = .false.
369  
370      call getMatchingElementList(atypes, "is_Directional", .true., &
371           nMatches, MatchList)
372      if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
373  
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 –
374      call getMatchingElementList(atypes, "is_Dipole", .true., &
375           nMatches, MatchList)
376 <    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
376 >    if (nMatches .gt. 0) FF_uses_Dipoles = .true.
377      
378      call getMatchingElementList(atypes, "is_GayBerne", .true., &
379           nMatches, MatchList)
380 <    if (nMatches .gt. 0) then
385 <       FF_uses_GayBerne = .true.
386 <       FF_uses_DirectionalAtoms = .true.
387 <    endif
380 >    if (nMatches .gt. 0) FF_uses_GayBerne = .true.
381  
382      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
383      if (nMatches .gt. 0) FF_uses_EAM = .true.
384  
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
385  
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)
386      haveSaneForceField = .true.
387  
388      !! check to make sure the FF_uses_RF setting makes sense
389  
390 <    if (FF_uses_dipoles) then
390 >    if (FF_uses_Dipoles) then
391         if (FF_uses_RF) then
392            dielect = getDielect()
393            call initialize_rf(dielect)
# Line 419 | Line 401 | contains
401         endif
402      endif
403  
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 –
404      if (FF_uses_EAM) then
405         call init_EAM_FF(my_status)
406         if (my_status /= 0) then
# Line 446 | Line 418 | contains
418            haveSaneForceField = .false.
419            return
420         endif
449 –    endif
450 –
451 –    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
421      endif
422  
423      if (.not. haveNeighborList) then
# Line 514 | Line 483 | contains
483      integer :: localError
484      integer :: propPack_i, propPack_j
485      integer :: loopStart, loopEnd, loop
486 <
486 >    integer :: iHash
487      real(kind=dp) :: listSkin = 1.0  
488  
489      !! initialize local variables  
# Line 606 | Line 575 | contains
575   #endif
576         outer: do i = istart, iend
577  
578 + #ifdef IS_MPI
579 +             me_i = atid_row(i)
580 + #else
581 +             me_i = atid(i)
582 + #endif
583 +
584            if (update_nlist) point(i) = nlist + 1
585  
586            n_in_i = groupStartRow(i+1) - groupStartRow(i)
# Line 633 | Line 608 | contains
608               endif
609  
610   #ifdef IS_MPI
611 +             me_j = atid_col(j)
612               call get_interatomic_vector(q_group_Row(:,i), &
613                    q_group_Col(:,j), d_grp, rgrpsq)
614   #else
615 +             me_j = atid(j)
616               call get_interatomic_vector(q_group(:,i), &
617                    q_group(:,j), d_grp, rgrpsq)
618   #endif
619  
620 <             if (rgrpsq < rlistsq) then
620 >             if (rgrpsq < InteractionHash(me_i,me_j)%rListsq) then
621                  if (update_nlist) then
622                     nlist = nlist + 1
623  
# Line 858 | Line 835 | contains
835   #else
836               me_i = atid(i)
837   #endif
838 <
839 <             if (PropertyMap(me_i)%is_Dipole) then
838 >             iHash = InteractionHash(me_i,me_j)
839 >            
840 >             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
841  
842                  mu_i = getDipoleMoment(me_i)
843  
# Line 926 | Line 904 | contains
904      real ( kind = dp ) :: ebalance
905      integer :: me_i, me_j
906  
907 +    integer :: iHash
908 +
909      r = sqrt(rijsq)
910      vpair = 0.0d0
911      fpair(1:3) = 0.0d0
# Line 938 | Line 918 | contains
918      me_j = atid(j)
919   #endif
920  
921 <    !    write(*,*) i, j, me_i, me_j
921 >    iHash = InteractionHash(me_i, me_j)
922  
923 <    if (FF_uses_LennardJones .and. SIM_uses_LennardJones) then
924 <
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 <
923 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
924 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
925      endif
926  
927 <    if (FF_uses_Electrostatics .and. SIM_uses_Electrostatics) then
927 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
928 >       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
929 >            pot, eFrame, f, t, do_pot)
930  
931 <       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
931 >       if (FF_uses_RF .and. SIM_uses_RF) then
932  
933 <       if (FF_uses_dipoles .and. SIM_uses_dipoles) then      
934 <          if ( PropertyMap(me_i)%is_Dipole .and. &
935 <               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
933 >          ! CHECK ME (RF needs to know about all electrostatic types)
934 >          call accumulate_rf(i, j, r, eFrame, sw)
935 >          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
936         endif
937 +
938      endif
939  
940 +    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
941 +       call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
942 +            pot, A, f, t, do_pot)
943 +    endif
944  
945 <    if (FF_uses_Sticky .and. SIM_uses_sticky) then
946 <
947 <       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 <
945 >    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
946 >       call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
947 >            pot, A, f, t, do_pot)
948      endif
949  
950 <    if (FF_uses_StickyPower .and. SIM_uses_stickypower) then
951 <       if ( PropertyMap(me_i)%is_StickyPower .and. &
952 <            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
950 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
951 >       call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
952 >            pot, A, f, t, do_pot)
953      endif
954      
955 <    if (FF_uses_GayBerne .and. SIM_uses_GayBerne) then
956 <
957 <       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 <
955 >    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
956 > !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
957 > !           pot, A, f, t, do_pot)
958      endif
959  
960 <    if (FF_uses_EAM .and. SIM_uses_EAM) then
961 <
962 <       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 <
960 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
961 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
962 >            do_pot)
963      endif
964  
965 <
966 <    !    write(*,*) PropertyMap(me_i)%is_Shape,PropertyMap(me_j)%is_Shape
967 <
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
965 >    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
966 >       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
967 >            pot, A, f, t, do_pot)
968      endif
969  
970 +    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
971 +       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
972 +            pot, A, f, t, do_pot)
973 +    endif
974 +    
975    end subroutine do_pair
976  
977    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
# Line 1040 | Line 989 | contains
989      real ( kind = dp )                :: r, rc
990      real ( kind = dp ), intent(inout) :: d(3), dc(3)
991  
992 <    logical :: is_EAM_i, is_EAM_j
1044 <
1045 <    integer :: me_i, me_j
992 >    integer :: me_i, me_j, iHash
993  
1047 –
1048 –    r = sqrt(rijsq)
1049 –    if (SIM_uses_molecular_cutoffs) then
1050 –       rc = sqrt(rcijsq)
1051 –    else
1052 –       rc = r
1053 –    endif
1054 –
1055 –
994   #ifdef IS_MPI  
995      me_i = atid_row(i)
996      me_j = atid_col(j)  
# Line 1061 | Line 999 | contains
999      me_j = atid(j)  
1000   #endif
1001  
1002 <    if (FF_uses_EAM .and. SIM_uses_EAM) then
1065 <
1066 <       if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1067 <            call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1002 >    iHash = InteractionHash(me_i, me_j)
1003  
1004 +    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1005 +            call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1006      endif
1007 <
1007 >    
1008    end subroutine do_prepair
1009  
1010  
# Line 1261 | Line 1198 | contains
1198  
1199    function FF_UsesDirectionalAtoms() result(doesit)
1200      logical :: doesit
1201 <    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
1201 >    doesit = FF_uses_DirectionalAtoms
1202    end function FF_UsesDirectionalAtoms
1203  
1204    function FF_RequiresPrepairCalc() result(doesit)

Diff Legend

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