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 2267 by tim, Fri Jul 29 17:34:06 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.25 2005-07-29 17:34:06 tim Exp $, $Date: 2005-07-29 17:34:06 $, $Name: not supported by cvs2svn $, $Revision: 1.25 $
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.
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
138 >  type, public :: Interaction
139 >     integer :: InteractionHash
140 >     real(kind=dp) :: rList = 0.0_dp
141 >     real(kind=dp) :: rListSq = 0.0_dp
142 >  end type Interaction
143 >  
144 >  type(Interaction), dimension(:,:),allocatable :: InteractionMap
145 >  
146  
147 <  type(Properties), dimension(:),allocatable :: PropertyMap
147 <
147 >  
148   contains
149  
150 <  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)
150 >
151 >  subroutine createInteractionMap(status)
152      integer :: nAtypes
153 <    integer :: status
153 >    integer, intent(out) :: status
154      integer :: i
155 <    logical :: thisProperty
156 <    real (kind=DP) :: thisDPproperty
157 <
155 >    integer :: j
156 >    integer :: ihash
157 >    real(kind=dp) :: myRcut
158 >    !! Test Types
159 >    logical :: i_is_LJ
160 >    logical :: i_is_Elect
161 >    logical :: i_is_Sticky
162 >    logical :: i_is_StickyP
163 >    logical :: i_is_GB
164 >    logical :: i_is_EAM
165 >    logical :: i_is_Shape
166 >    logical :: j_is_LJ
167 >    logical :: j_is_Elect
168 >    logical :: j_is_Sticky
169 >    logical :: j_is_StickyP
170 >    logical :: j_is_GB
171 >    logical :: j_is_EAM
172 >    logical :: j_is_Shape
173 >    
174      status = 0
175 <
175 >    
176 >    if (.not. associated(atypes)) then
177 >       call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
178 >       status = -1
179 >       return
180 >    endif
181 >    
182      nAtypes = getSize(atypes)
183 <
183 >    
184      if (nAtypes == 0) then
185         status = -1
186         return
187      end if
188  
189 <    if (.not. allocated(PropertyMap)) then
190 <       allocate(PropertyMap(nAtypes))
189 >    if (.not. allocated(InteractionMap)) then
190 >       allocate(InteractionMap(nAtypes,nAtypes))
191      endif
192 <
192 >        
193      do i = 1, nAtypes
194 <       call getElementProperty(atypes, i, "is_Directional", thisProperty)
195 <       PropertyMap(i)%is_Directional = thisProperty
194 >       call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
195 >       call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
196 >       call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
197 >       call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
198 >       call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
199 >       call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
200 >       call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
201  
202 <       call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
186 <       PropertyMap(i)%is_LennardJones = thisProperty
202 >       do j = i, nAtypes
203  
204 <       call getElementProperty(atypes, i, "is_Electrostatic", thisProperty)
205 <       PropertyMap(i)%is_Electrostatic = thisProperty
204 >          iHash = 0
205 >          myRcut = 0.0_dp
206  
207 <       call getElementProperty(atypes, i, "is_Charge", thisProperty)
208 <       PropertyMap(i)%is_Charge = thisProperty
209 <
210 <       call getElementProperty(atypes, i, "is_Dipole", thisProperty)
211 <       PropertyMap(i)%is_Dipole = thisProperty
212 <
213 <       call getElementProperty(atypes, i, "is_Quadrupole", thisProperty)
198 <       PropertyMap(i)%is_Quadrupole = thisProperty
207 >          call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ)
208 >          call getElementProperty(atypes, j, "is_Electrostatic", j_is_Elect)
209 >          call getElementProperty(atypes, j, "is_Sticky", j_is_Sticky)
210 >          call getElementProperty(atypes, j, "is_StickyPower", j_is_StickyP)
211 >          call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
212 >          call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
213 >          call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
214  
215 <       call getElementProperty(atypes, i, "is_Sticky", thisProperty)
216 <       PropertyMap(i)%is_Sticky = thisProperty
217 <      
218 <       call getElementProperty(atypes, i, "is_StickyPower", thisProperty)
219 <       PropertyMap(i)%is_StickyPower = thisProperty
215 >          if (i_is_LJ .and. j_is_LJ) then
216 >             iHash = ior(iHash, LJ_PAIR)            
217 >          endif
218 >          
219 >          if (i_is_Elect .and. j_is_Elect) then
220 >             iHash = ior(iHash, ELECTROSTATIC_PAIR)
221 >          endif
222 >          
223 >          if (i_is_Sticky .and. j_is_Sticky) then
224 >             iHash = ior(iHash, STICKY_PAIR)
225 >          endif
226  
227 <       call getElementProperty(atypes, i, "is_GayBerne", thisProperty)
228 <       PropertyMap(i)%is_GayBerne = thisProperty
227 >          if (i_is_StickyP .and. j_is_StickyP) then
228 >             iHash = ior(iHash, STICKYPOWER_PAIR)
229 >          endif
230  
231 <       call getElementProperty(atypes, i, "is_EAM", thisProperty)
232 <       PropertyMap(i)%is_EAM = thisProperty
231 >          if (i_is_EAM .and. j_is_EAM) then
232 >             iHash = ior(iHash, EAM_PAIR)
233 >          endif
234  
235 <       call getElementProperty(atypes, i, "is_Shape", thisProperty)
236 <       PropertyMap(i)%is_Shape = thisProperty
235 >          if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
236 >          if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
237 >          if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
238  
239 <       call getElementProperty(atypes, i, "is_FLARB", thisProperty)
240 <       PropertyMap(i)%is_FLARB = thisProperty
239 >          if (i_is_Shape .and. j_is_Shape) iHash = ior(iHash, SHAPE_PAIR)
240 >          if (i_is_Shape .and. j_is_LJ) iHash = ior(iHash, SHAPE_LJ)
241 >          if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
242 >
243 >
244 >          InteractionMap(i,j)%InteractionHash = iHash
245 >          InteractionMap(j,i)%InteractionHash = iHash
246 >
247 >       end do
248 >
249      end do
250  
251 <    havePropertyMap = .true.
251 >    haveInteractionMap = .true.
252 >  end subroutine createInteractionMap
253  
254 <  end subroutine createPropertyMap
254 > ! Query each potential and return the cutoff for that potential. We build the neighbor list based on the largest cutoff value for that atype. Each potential can decide whether to calculate the force for that atype based upon it's own cutoff.
255 >  subroutine createRcuts(defaultRList,stat)
256 >    real(kind=dp), intent(in), optional :: defaultRList
257 >    integer :: iMap
258 >    integer :: map_i,map_j
259 >    real(kind=dp) :: thisRCut = 0.0_dp
260 >    real(kind=dp) :: actualCutoff = 0.0_dp
261 >    integer, intent(out) :: stat
262 >    integer :: nAtypes
263 >    integer :: myStatus
264  
265 +    stat = 0
266 +    if (.not. haveInteractionMap) then
267 +
268 +       call createInteractionMap(myStatus)
269 +
270 +       if (myStatus .ne. 0) then
271 +          write(default_error, *) 'createInteractionMap failed in doForces!'
272 +          stat = -1
273 +          return
274 +       endif
275 +    endif
276 +
277 +
278 +    nAtypes = getSize(atypes)
279 +    !! If we pass a default rcut, set all atypes to that cutoff distance
280 +    if(present(defaultRList)) then
281 +       InteractionMap(:,:)%rList = defaultRList
282 +       InteractionMap(:,:)%rListSq = defaultRList*defaultRList
283 +       haveRlist = .true.
284 +       return
285 +    end if
286 +
287 +    do map_i = 1,nAtypes
288 +       do map_j = map_i,nAtypes
289 +          iMap = InteractionMap(map_i, map_j)%InteractionHash
290 +          
291 +          if ( iand(iMap, LJ_PAIR).ne.0 ) then
292 +             ! thisRCut = getLJCutOff(map_i,map_j)
293 +             if (thisRcut > actualCutoff) actualCutoff = thisRcut
294 +          endif
295 +          
296 +          if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
297 +             ! thisRCut = getElectrostaticCutOff(map_i,map_j)
298 +             if (thisRcut > actualCutoff) actualCutoff = thisRcut
299 +          endif
300 +          
301 +          if ( iand(iMap, STICKY_PAIR).ne.0 ) then
302 +             ! thisRCut = getStickyCutOff(map_i,map_j)
303 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
304 +           endif
305 +          
306 +           if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
307 +              ! thisRCut = getStickyPowerCutOff(map_i,map_j)
308 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
309 +           endif
310 +          
311 +           if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
312 +              ! thisRCut = getGayberneCutOff(map_i,map_j)
313 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
314 +           endif
315 +          
316 +           if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
317 + !              thisRCut = getGaybrneLJCutOff(map_i,map_j)
318 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
319 +           endif
320 +          
321 +           if ( iand(iMap, EAM_PAIR).ne.0 ) then      
322 + !              thisRCut = getEAMCutOff(map_i,map_j)
323 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
324 +           endif
325 +          
326 +           if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
327 + !              thisRCut = getShapeCutOff(map_i,map_j)
328 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
329 +           endif
330 +          
331 +           if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
332 + !              thisRCut = getShapeLJCutOff(map_i,map_j)
333 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
334 +           endif
335 +           InteractionMap(map_i, map_j)%rList = actualCutoff
336 +           InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
337 +        end do
338 +     end do
339 +     haveRlist = .true.
340 +  end subroutine createRcuts
341 +
342 +
343 + !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
344 + !!$  subroutine setRlistDF( this_rlist )
345 + !!$
346 + !!$   real(kind=dp) :: this_rlist
347 + !!$
348 + !!$    rlist = this_rlist
349 + !!$    rlistsq = rlist * rlist
350 + !!$
351 + !!$    haveRlist = .true.
352 + !!$
353 + !!$  end subroutine setRlistDF
354 +
355 +
356    subroutine setSimVariables()
357      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
358      SIM_uses_LennardJones = SimUsesLennardJones()
# Line 249 | Line 382 | contains
382  
383      error = 0
384  
385 <    if (.not. havePropertyMap) then
385 >    if (.not. haveInteractionMap) then
386  
387         myStatus = 0
388  
389 <       call createPropertyMap(myStatus)
389 >       call createInteractionMap(myStatus)
390  
391         if (myStatus .ne. 0) then
392 <          write(default_error, *) 'createPropertyMap failed in doForces!'
392 >          write(default_error, *) 'createInteractionMap failed in doForces!'
393            error = -1
394            return
395         endif
# Line 514 | Line 647 | contains
647      integer :: localError
648      integer :: propPack_i, propPack_j
649      integer :: loopStart, loopEnd, loop
650 <
650 >    integer :: iMap
651      real(kind=dp) :: listSkin = 1.0  
652  
653      !! initialize local variables  
# Line 606 | Line 739 | contains
739   #endif
740         outer: do i = istart, iend
741  
742 + #ifdef IS_MPI
743 +             me_i = atid_row(i)
744 + #else
745 +             me_i = atid(i)
746 + #endif
747 +
748            if (update_nlist) point(i) = nlist + 1
749  
750            n_in_i = groupStartRow(i+1) - groupStartRow(i)
# Line 633 | Line 772 | contains
772               endif
773  
774   #ifdef IS_MPI
775 +             me_j = atid_col(j)
776               call get_interatomic_vector(q_group_Row(:,i), &
777                    q_group_Col(:,j), d_grp, rgrpsq)
778   #else
779 +             me_j = atid(j)
780               call get_interatomic_vector(q_group(:,i), &
781                    q_group(:,j), d_grp, rgrpsq)
782   #endif
783  
784 <             if (rgrpsq < rlistsq) then
784 >             if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
785                  if (update_nlist) then
786                     nlist = nlist + 1
787  
# Line 858 | Line 999 | contains
999   #else
1000               me_i = atid(i)
1001   #endif
1002 <
1003 <             if (PropertyMap(me_i)%is_Dipole) then
1002 >             iMap = InteractionMap(me_i, me_j)%InteractionHash
1003 >            
1004 >             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1005  
1006                  mu_i = getDipoleMoment(me_i)
1007  
# Line 926 | Line 1068 | contains
1068      real ( kind = dp ) :: ebalance
1069      integer :: me_i, me_j
1070  
1071 +    integer :: iMap
1072 +
1073      r = sqrt(rijsq)
1074      vpair = 0.0d0
1075      fpair(1:3) = 0.0d0
# Line 938 | Line 1082 | contains
1082      me_j = atid(j)
1083   #endif
1084  
1085 <    !    write(*,*) i, j, me_i, me_j
1085 >    iMap = InteractionMap(me_i, me_j)%InteractionHash
1086  
1087 <    if (FF_uses_LennardJones .and. SIM_uses_LennardJones) then
1087 >    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1088 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1089 >    endif
1090  
1091 <       if ( PropertyMap(me_i)%is_LennardJones .and. &
1092 <            PropertyMap(me_j)%is_LennardJones ) then
1093 <          call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1091 >    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1092 >       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1093 >            pot, eFrame, f, t, do_pot)
1094 >
1095 >       if (FF_uses_RF .and. SIM_uses_RF) then
1096 >
1097 >          ! CHECK ME (RF needs to know about all electrostatic types)
1098 >          call accumulate_rf(i, j, r, eFrame, sw)
1099 >          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1100         endif
1101  
1102      endif
1103  
1104 <    if (FF_uses_Electrostatics .and. SIM_uses_Electrostatics) then
1104 >    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1105 >       call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1106 >            pot, A, f, t, do_pot)
1107 >    endif
1108  
1109 <       if (PropertyMap(me_i)%is_Electrostatic .and. &
1110 <            PropertyMap(me_j)%is_Electrostatic) then
1111 <          call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1112 <               pot, eFrame, f, t, do_pot, ebalance)
958 <       endif
1109 >    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1110 >       call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1111 >            pot, A, f, t, do_pot)
1112 >    endif
1113  
1114 <       if (FF_uses_dipoles .and. SIM_uses_dipoles) then      
1115 <          if ( PropertyMap(me_i)%is_Dipole .and. &
1116 <               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
1114 >    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1115 >       call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1116 >            pot, A, f, t, do_pot)
1117      endif
1118 <
1119 <
1120 <    if (FF_uses_Sticky .and. SIM_uses_sticky) then
1121 <
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 <
1118 >    
1119 >    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1120 > !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1121 > !           pot, A, f, t, do_pot)
1122      endif
1123  
1124 <    if (FF_uses_StickyPower .and. SIM_uses_stickypower) then
1125 <       if ( PropertyMap(me_i)%is_StickyPower .and. &
1126 <            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
1124 >    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1125 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1126 >            do_pot)
1127      endif
988 –    
989 –    if (FF_uses_GayBerne .and. SIM_uses_GayBerne) then
1128  
1129 <       if ( PropertyMap(me_i)%is_GayBerne .and. &
1130 <            PropertyMap(me_j)%is_GayBerne) then
1131 <          call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
994 <               pot, A, f, t, do_pot)
995 <       endif
996 <
1129 >    if ( iand(iMap, SHAPE_PAIR).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
998 –
999 –    if (FF_uses_EAM .and. SIM_uses_EAM) then
1000 –
1001 –       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
1133  
1134 +    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1135 +       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1136 +            pot, A, f, t, do_pot)
1137      endif
1138 <
1008 <
1009 <    !    write(*,*) PropertyMap(me_i)%is_Shape,PropertyMap(me_j)%is_Shape
1010 <
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
1024 <    endif
1025 <
1138 >    
1139    end subroutine do_pair
1140  
1141    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
# Line 1040 | Line 1153 | contains
1153      real ( kind = dp )                :: r, rc
1154      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1155  
1156 <    logical :: is_EAM_i, is_EAM_j
1044 <
1045 <    integer :: me_i, me_j
1156 >    integer :: me_i, me_j, iMap
1157  
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 –
1158   #ifdef IS_MPI  
1159      me_i = atid_row(i)
1160      me_j = atid_col(j)  
# Line 1061 | Line 1163 | contains
1163      me_j = atid(j)  
1164   #endif
1165  
1166 <    if (FF_uses_EAM .and. SIM_uses_EAM) then
1166 >    iMap = InteractionMap(me_i, me_j)%InteractionHash
1167  
1168 <       if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1168 >    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1169              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1068 –
1170      endif
1171 <
1171 >    
1172    end subroutine do_prepair
1173  
1174  

Diff Legend

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