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 2268 by gezelter, Fri Jul 29 19:38:27 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.26 2005-07-29 19:38:27 gezelter Exp $, $Date: 2005-07-29 19:38:27 $, $Name: not supported by cvs2svn $, $Revision: 1.26 $
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
# Line 80 | Line 81 | module doForces
81    logical, save :: haveRlist = .false.
82    logical, save :: haveNeighborList = .false.
83    logical, save :: haveSIMvariables = .false.
83 –  logical, save :: havePropertyMap = .false.
84    logical, save :: haveSaneForceField = .false.
85 +  logical, save :: haveInteractionMap = .false.
86  
87    logical, save :: FF_uses_DirectionalAtoms
88    logical, save :: FF_uses_LennardJones
# Line 115 | Line 116 | module doForces
116    logical, save :: SIM_uses_PBC
117    logical, save :: SIM_uses_molecular_cutoffs
118  
119 <  real(kind=dp), save :: rlist, rlistsq
119 >  !!!GO AWAY---------
120 >  !!!!!real(kind=dp), save :: rlist, rlistsq
121  
122    public :: init_FF
123    public :: do_force_loop
124 <  public :: setRlistDF
124 > !  public :: setRlistDF
125 >  !public :: addInteraction
126 >  !public :: setInteractionHash
127 >  !public :: getInteractionHash
128 >  public :: createInteractionMap
129 >  public :: createRcuts
130  
131   #ifdef PROFILE
132    public :: getforcetime
# Line 128 | Line 135 | module doForces
135    integer :: nLoops
136   #endif
137  
138 <  type :: Properties
139 <     logical :: is_Directional   = .false.
140 <     logical :: is_LennardJones  = .false.
141 <     logical :: is_Electrostatic = .false.
142 <     logical :: is_Charge        = .false.
143 <     logical :: is_Dipole        = .false.
144 <     logical :: is_Quadrupole    = .false.
145 <     logical :: is_Sticky        = .false.
146 <     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
138 >  type, public :: Interaction
139 >     integer :: InteractionHash
140 >     real(kind=dp) :: rCut = 0.0_dp
141 >     real(kind=dp) :: rCutSq = 0.0_dp    
142 >     real(kind=dp) :: rListSq = 0.0_dp
143 >  end type Interaction
144 >  
145 >  type(Interaction), dimension(:,:),allocatable :: InteractionMap
146 >  
147  
148 <  type(Properties), dimension(:),allocatable :: PropertyMap
147 <
148 >  
149   contains
150  
151 <  subroutine setRlistDF( this_rlist )
152 <
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)
151 >
152 >  subroutine createInteractionMap(status)
153      integer :: nAtypes
154 <    integer :: status
154 >    integer, intent(out) :: status
155      integer :: i
156 <    logical :: thisProperty
157 <    real (kind=DP) :: thisDPproperty
156 >    integer :: j
157 >    integer :: ihash
158 >    real(kind=dp) :: myRcut
159 >    !! Test Types
160 >    logical :: i_is_LJ
161 >    logical :: i_is_Elect
162 >    logical :: i_is_Sticky
163 >    logical :: i_is_StickyP
164 >    logical :: i_is_GB
165 >    logical :: i_is_EAM
166 >    logical :: i_is_Shape
167 >    logical :: j_is_LJ
168 >    logical :: j_is_Elect
169 >    logical :: j_is_Sticky
170 >    logical :: j_is_StickyP
171 >    logical :: j_is_GB
172 >    logical :: j_is_EAM
173 >    logical :: j_is_Shape
174 >    
175 >    status = 0  
176  
177 <    status = 0
178 <
177 >    if (.not. associated(atypes)) then
178 >       call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
179 >       status = -1
180 >       return
181 >    endif
182 >    
183      nAtypes = getSize(atypes)
184 <
184 >    
185      if (nAtypes == 0) then
186         status = -1
187         return
188      end if
189  
190 <    if (.not. allocated(PropertyMap)) then
191 <       allocate(PropertyMap(nAtypes))
190 >    if (.not. allocated(InteractionMap)) then
191 >       allocate(InteractionMap(nAtypes,nAtypes))
192      endif
193 <
193 >        
194      do i = 1, nAtypes
195 <       call getElementProperty(atypes, i, "is_Directional", thisProperty)
196 <       PropertyMap(i)%is_Directional = thisProperty
195 >       call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
196 >       call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
197 >       call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
198 >       call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
199 >       call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
200 >       call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
201 >       call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
202  
203 <       call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
186 <       PropertyMap(i)%is_LennardJones = thisProperty
203 >       do j = i, nAtypes
204  
205 <       call getElementProperty(atypes, i, "is_Electrostatic", thisProperty)
206 <       PropertyMap(i)%is_Electrostatic = thisProperty
205 >          iHash = 0
206 >          myRcut = 0.0_dp
207  
208 <       call getElementProperty(atypes, i, "is_Charge", thisProperty)
209 <       PropertyMap(i)%is_Charge = thisProperty
208 >          call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ)
209 >          call getElementProperty(atypes, j, "is_Electrostatic", j_is_Elect)
210 >          call getElementProperty(atypes, j, "is_Sticky", j_is_Sticky)
211 >          call getElementProperty(atypes, j, "is_StickyPower", j_is_StickyP)
212 >          call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
213 >          call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
214 >          call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
215  
216 <       call getElementProperty(atypes, i, "is_Dipole", thisProperty)
217 <       PropertyMap(i)%is_Dipole = thisProperty
216 >          if (i_is_LJ .and. j_is_LJ) then
217 >             iHash = ior(iHash, LJ_PAIR)            
218 >          endif
219 >          
220 >          if (i_is_Elect .and. j_is_Elect) then
221 >             iHash = ior(iHash, ELECTROSTATIC_PAIR)
222 >          endif
223 >          
224 >          if (i_is_Sticky .and. j_is_Sticky) then
225 >             iHash = ior(iHash, STICKY_PAIR)
226 >          endif
227  
228 <       call getElementProperty(atypes, i, "is_Quadrupole", thisProperty)
229 <       PropertyMap(i)%is_Quadrupole = thisProperty
228 >          if (i_is_StickyP .and. j_is_StickyP) then
229 >             iHash = ior(iHash, STICKYPOWER_PAIR)
230 >          endif
231  
232 <       call getElementProperty(atypes, i, "is_Sticky", thisProperty)
233 <       PropertyMap(i)%is_Sticky = thisProperty
234 <      
203 <       call getElementProperty(atypes, i, "is_StickyPower", thisProperty)
204 <       PropertyMap(i)%is_StickyPower = thisProperty
232 >          if (i_is_EAM .and. j_is_EAM) then
233 >             iHash = ior(iHash, EAM_PAIR)
234 >          endif
235  
236 <       call getElementProperty(atypes, i, "is_GayBerne", thisProperty)
237 <       PropertyMap(i)%is_GayBerne = thisProperty
236 >          if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
237 >          if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
238 >          if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
239  
240 <       call getElementProperty(atypes, i, "is_EAM", thisProperty)
241 <       PropertyMap(i)%is_EAM = thisProperty
240 >          if (i_is_Shape .and. j_is_Shape) iHash = ior(iHash, SHAPE_PAIR)
241 >          if (i_is_Shape .and. j_is_LJ) iHash = ior(iHash, SHAPE_LJ)
242 >          if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
243  
212 –       call getElementProperty(atypes, i, "is_Shape", thisProperty)
213 –       PropertyMap(i)%is_Shape = thisProperty
244  
245 <       call getElementProperty(atypes, i, "is_FLARB", thisProperty)
246 <       PropertyMap(i)%is_FLARB = thisProperty
245 >          InteractionMap(i,j)%InteractionHash = iHash
246 >          InteractionMap(j,i)%InteractionHash = iHash
247 >
248 >       end do
249 >
250      end do
251  
252 <    havePropertyMap = .true.
252 >    haveInteractionMap = .true.
253 >  end subroutine createInteractionMap
254  
255 <  end subroutine createPropertyMap
255 >  ! Query each potential and return the cutoff for that potential. We
256 >  ! build the neighbor list based on the largest cutoff value for that
257 >  ! atype. Each potential can decide whether to calculate the force for
258 >  ! that atype based upon it's own cutoff.
259 >  
260 >  subroutine createRcuts(defaultRcut, defaultSkinThickness, stat)
261  
262 +    real(kind=dp), intent(in), optional :: defaultRCut, defaultSkinThickness
263 +    integer :: iMap
264 +    integer :: map_i,map_j
265 +    real(kind=dp) :: thisRCut = 0.0_dp
266 +    real(kind=dp) :: actualCutoff = 0.0_dp
267 +    integer, intent(out) :: stat
268 +    integer :: nAtypes
269 +    integer :: myStatus
270 +
271 +    stat = 0
272 +    if (.not. haveInteractionMap) then
273 +
274 +       call createInteractionMap(myStatus)
275 +
276 +       if (myStatus .ne. 0) then
277 +          write(default_error, *) 'createInteractionMap failed in doForces!'
278 +          stat = -1
279 +          return
280 +       endif
281 +    endif
282 +
283 +    nAtypes = getSize(atypes)
284 +    !! If we pass a default rcut, set all atypes to that cutoff distance
285 +    if(present(defaultRList)) then
286 +       InteractionMap(:,:)%rCut = defaultRCut
287 +       InteractionMap(:,:)%rCutSq = defaultRCut*defaultRCut
288 +       InteractionMap(:,:)%rListSq = (defaultRCut+defaultSkinThickness)**2
289 +       haveRlist = .true.
290 +       return
291 +    end if
292 +
293 +    do map_i = 1,nAtypes
294 +       do map_j = map_i,nAtypes
295 +          iMap = InteractionMap(map_i, map_j)%InteractionHash
296 +          
297 +          if ( iand(iMap, LJ_PAIR).ne.0 ) then
298 +             ! thisRCut = getLJCutOff(map_i,map_j)
299 +             if (thisRcut > actualCutoff) actualCutoff = thisRcut
300 +          endif
301 +          
302 +          if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
303 +             ! thisRCut = getElectrostaticCutOff(map_i,map_j)
304 +             if (thisRcut > actualCutoff) actualCutoff = thisRcut
305 +          endif
306 +          
307 +          if ( iand(iMap, STICKY_PAIR).ne.0 ) then
308 +             ! thisRCut = getStickyCutOff(map_i,map_j)
309 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
310 +           endif
311 +          
312 +           if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
313 +              ! thisRCut = getStickyPowerCutOff(map_i,map_j)
314 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
315 +           endif
316 +          
317 +           if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
318 +              ! thisRCut = getGayberneCutOff(map_i,map_j)
319 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
320 +           endif
321 +          
322 +           if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
323 + !              thisRCut = getGaybrneLJCutOff(map_i,map_j)
324 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
325 +           endif
326 +          
327 +           if ( iand(iMap, EAM_PAIR).ne.0 ) then      
328 + !              thisRCut = getEAMCutOff(map_i,map_j)
329 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
330 +           endif
331 +          
332 +           if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
333 + !              thisRCut = getShapeCutOff(map_i,map_j)
334 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
335 +           endif
336 +          
337 +           if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
338 + !              thisRCut = getShapeLJCutOff(map_i,map_j)
339 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
340 +           endif
341 +           InteractionMap(map_i, map_j)%rCut = actualCutoff
342 +           InteractionMap(map_i, map_j)%rCutSq = actualCutoff * actualCutoff
343 +           InteractionMap(map_i, map_j)%rListSq = (actualCutoff + skinThickness)**2
344 +
345 +           InteractionMap(map_j, map_i)%rCut = InteractionMap(map_i, map_j)%rCut
346 +           InteractionMap(map_j, map_i)%rCutSq = InteractionMap(map_i, map_j)%rCutSq
347 +           InteractionMap(map_j, map_i)%rListSq = InteractionMap(map_i, map_j)%rListSq
348 +        end do
349 +     end do
350 +     ! now the groups
351 +
352 +
353 +
354 +     haveRlist = .true.
355 +  end subroutine createRcuts
356 +
357 +
358 + !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
359 + !!$  subroutine setRlistDF( this_rlist )
360 + !!$
361 + !!$   real(kind=dp) :: this_rlist
362 + !!$
363 + !!$    rlist = this_rlist
364 + !!$    rlistsq = rlist * rlist
365 + !!$
366 + !!$    haveRlist = .true.
367 + !!$
368 + !!$  end subroutine setRlistDF
369 +
370 +
371    subroutine setSimVariables()
372      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
373      SIM_uses_LennardJones = SimUsesLennardJones()
# Line 249 | Line 397 | contains
397  
398      error = 0
399  
400 <    if (.not. havePropertyMap) then
401 <
402 <       myStatus = 0
403 <
404 <       call createPropertyMap(myStatus)
257 <
400 >    if (.not. haveInteractionMap) then
401 >      
402 >       myStatus = 0      
403 >       call createInteractionMap(myStatus)
404 >      
405         if (myStatus .ne. 0) then
406 <          write(default_error, *) 'createPropertyMap failed in doForces!'
406 >          write(default_error, *) 'createInteractionMap failed in doForces!'
407            error = -1
408            return
409         endif
# Line 514 | Line 661 | contains
661      integer :: localError
662      integer :: propPack_i, propPack_j
663      integer :: loopStart, loopEnd, loop
664 <
664 >    integer :: iMap
665      real(kind=dp) :: listSkin = 1.0  
666  
667      !! initialize local variables  
# Line 605 | Line 752 | contains
752         iend = nGroups - 1
753   #endif
754         outer: do i = istart, iend
755 +
756 + #ifdef IS_MPI
757 +             me_i = atid_row(i)
758 + #else
759 +             me_i = atid(i)
760 + #endif
761  
762            if (update_nlist) point(i) = nlist + 1
763  
# Line 633 | Line 786 | contains
786               endif
787  
788   #ifdef IS_MPI
789 +             me_j = atid_col(j)
790               call get_interatomic_vector(q_group_Row(:,i), &
791                    q_group_Col(:,j), d_grp, rgrpsq)
792   #else
793 +             me_j = atid(j)
794               call get_interatomic_vector(q_group(:,i), &
795                    q_group(:,j), d_grp, rgrpsq)
796   #endif
797  
798 <             if (rgrpsq < rlistsq) then
798 >             if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
799                  if (update_nlist) then
800                     nlist = nlist + 1
801  
# Line 858 | Line 1013 | contains
1013   #else
1014               me_i = atid(i)
1015   #endif
1016 <
1017 <             if (PropertyMap(me_i)%is_Dipole) then
1016 >             iMap = InteractionMap(me_i, me_j)%InteractionHash
1017 >            
1018 >             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1019  
1020                  mu_i = getDipoleMoment(me_i)
1021  
# Line 925 | Line 1081 | contains
1081      real ( kind = dp ), intent(inout) :: d(3)
1082      real ( kind = dp ) :: ebalance
1083      integer :: me_i, me_j
1084 +
1085 +    integer :: iMap
1086  
1087      r = sqrt(rijsq)
1088      vpair = 0.0d0
# Line 938 | Line 1096 | contains
1096      me_j = atid(j)
1097   #endif
1098  
1099 <    !    write(*,*) i, j, me_i, me_j
1099 >    iMap = InteractionMap(me_i, me_j)%InteractionHash
1100  
1101 <    if (FF_uses_LennardJones .and. SIM_uses_LennardJones) then
1102 <
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 <
1101 >    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1102 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1103      endif
1104  
1105 <    if (FF_uses_Electrostatics .and. SIM_uses_Electrostatics) then
1105 >    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1106 >       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1107 >            pot, eFrame, f, t, do_pot)
1108  
1109 <       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
1109 >       if (FF_uses_RF .and. SIM_uses_RF) then
1110  
1111 <       if (FF_uses_dipoles .and. SIM_uses_dipoles) then      
1112 <          if ( PropertyMap(me_i)%is_Dipole .and. &
1113 <               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
1111 >          ! CHECK ME (RF needs to know about all electrostatic types)
1112 >          call accumulate_rf(i, j, r, eFrame, sw)
1113 >          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1114         endif
1115 +
1116      endif
1117  
1118 +    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1119 +       call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1120 +            pot, A, f, t, do_pot)
1121 +    endif
1122  
1123 <    if (FF_uses_Sticky .and. SIM_uses_sticky) then
1124 <
1125 <       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 <
1123 >    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1124 >       call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1125 >            pot, A, f, t, do_pot)
1126      endif
1127  
1128 <    if (FF_uses_StickyPower .and. SIM_uses_stickypower) then
1129 <       if ( PropertyMap(me_i)%is_StickyPower .and. &
1130 <            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
1128 >    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1129 >       call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1130 >            pot, A, f, t, do_pot)
1131      endif
1132      
1133 <    if (FF_uses_GayBerne .and. SIM_uses_GayBerne) then
1134 <
1135 <       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 <
1133 >    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1134 > !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1135 > !           pot, A, f, t, do_pot)
1136      endif
1137  
1138 <    if (FF_uses_EAM .and. SIM_uses_EAM) then
1139 <
1140 <       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 <
1138 >    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1139 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1140 >            do_pot)
1141      endif
1142  
1143 <
1144 <    !    write(*,*) PropertyMap(me_i)%is_Shape,PropertyMap(me_j)%is_Shape
1145 <
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
1143 >    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1144 >       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1145 >            pot, A, f, t, do_pot)
1146      endif
1147  
1148 +    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1149 +       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1150 +            pot, A, f, t, do_pot)
1151 +    endif
1152 +    
1153    end subroutine do_pair
1154  
1155    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
# Line 1040 | Line 1167 | contains
1167      real ( kind = dp )                :: r, rc
1168      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1169  
1170 <    logical :: is_EAM_i, is_EAM_j
1044 <
1045 <    integer :: me_i, me_j
1170 >    integer :: me_i, me_j, iMap
1171  
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 –
1172   #ifdef IS_MPI  
1173      me_i = atid_row(i)
1174      me_j = atid_col(j)  
# Line 1061 | Line 1177 | contains
1177      me_j = atid(j)  
1178   #endif
1179  
1180 <    if (FF_uses_EAM .and. SIM_uses_EAM) then
1180 >    iMap = InteractionMap(me_i, me_j)%InteractionHash
1181  
1182 <       if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1182 >    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1183              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1068 –
1184      endif
1185 <
1185 >    
1186    end subroutine do_prepair
1187  
1188  

Diff Legend

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