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 2266 by chuckv, Thu Jul 28 22:12:45 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.24 2005-07-28 22:12:45 chuckv Exp $, $Date: 2005-07-28 22:12:45 $, $Name: not supported by cvs2svn $, $Revision: 1.24 $
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
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_Quadrupole", thisProperty)
216 <       PropertyMap(i)%is_Quadrupole = 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_Sticky", thisProperty)
228 <       PropertyMap(i)%is_Sticky = thisProperty
229 <      
203 <       call getElementProperty(atypes, i, "is_StickyPower", thisProperty)
204 <       PropertyMap(i)%is_StickyPower = 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_GayBerne", thisProperty)
232 <       PropertyMap(i)%is_GayBerne = 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_EAM", thisProperty)
236 <       PropertyMap(i)%is_EAM = 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_Shape", thisProperty)
240 <       PropertyMap(i)%is_Shape = 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 <       call getElementProperty(atypes, i, "is_FLARB", thisProperty)
244 <       PropertyMap(i)%is_FLARB = thisProperty
243 >
244 >          InteractionMap(i,j)%InteractionHash = iHash
245 >          InteractionMap(j,i)%InteractionHash = iHash
246 >
247 >       end do
248 >
249      end do
250 +  end subroutine createInteractionMap
251  
252 <    havePropertyMap = .true.
252 > ! 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.
253 >  subroutine createRcuts(defaultRList,stat)
254 >    real(kind=dp), intent(in), optional :: defaultRList
255 >    integer :: iMap
256 >    integer :: map_i,map_j
257 >    real(kind=dp) :: thisRCut = 0.0_dp
258 >    real(kind=dp) :: actualCutoff = 0.0_dp
259 >    integer, intent(out) :: stat
260 >    integer :: nAtypes
261 >    integer :: myStatus
262  
263 <  end subroutine createPropertyMap
263 >    stat = 0
264 >    if (.not. haveInteractionMap) then
265  
266 +       call createInteractionMap(myStatus)
267 +
268 +       if (myStatus .ne. 0) then
269 +          write(default_error, *) 'createInteractionMap failed in doForces!'
270 +          stat = -1
271 +          return
272 +       endif
273 +    endif
274 +
275 +
276 +    nAtypes = getSize(atypes)
277 + ! If we pass a default rcut, set all atypes to that cutoff distance
278 +    if(present(defaultRList)) then
279 +       InteractionMap(:,:)%rList = defaultRList
280 +       InteractionMap(:,:)%rListSq = defaultRList*defaultRList
281 +       haveRlist = .true.
282 +       return
283 +    end if
284 +
285 +    do map_i = 1,nAtypes
286 +       do map_j = map_i,nAtypes
287 +          iMap = InteractionMap(map_i, map_j)%InteractionHash
288 +          
289 +          if ( iand(iMap, LJ_PAIR).ne.0 ) then
290 + !            thisRCut = getLJCutOff(map_i,map_j)
291 +             if (thisRcut > actualCutoff) actualCutoff = thisRcut
292 +          endif
293 +          
294 +          if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
295 + !            thisRCut = getElectrostaticCutOff(map_i,map_j)
296 +             if (thisRcut > actualCutoff) actualCutoff = thisRcut
297 +          endif
298 +          
299 +          if ( iand(iMap, STICKY_PAIR).ne.0 ) then
300 + !             thisRCut = getStickyCutOff(map_i,map_j)
301 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
302 +           endif
303 +          
304 +           if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
305 + !              thisRCut = getStickyPowerCutOff(map_i,map_j)
306 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
307 +           endif
308 +          
309 +           if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
310 + !              thisRCut = getGayberneCutOff(map_i,map_j)
311 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
312 +           endif
313 +          
314 +           if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
315 + !              thisRCut = getGaybrneLJCutOff(map_i,map_j)
316 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
317 +           endif
318 +          
319 +           if ( iand(iMap, EAM_PAIR).ne.0 ) then      
320 + !              thisRCut = getEAMCutOff(map_i,map_j)
321 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
322 +           endif
323 +          
324 +           if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
325 + !              thisRCut = getShapeCutOff(map_i,map_j)
326 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
327 +           endif
328 +          
329 +           if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
330 + !              thisRCut = getShapeLJCutOff(map_i,map_j)
331 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
332 +           endif
333 +           InteractionMap(map_i, map_j)%rList = actualCutoff
334 +           InteractionMap(map_i, map_j)%rListSq = actualCutoff * actualCutoff
335 +        end do
336 +     end do
337 +          haveRlist = .true.
338 +  end subroutine createRcuts
339 +
340 +
341 + !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
342 + !!$  subroutine setRlistDF( this_rlist )
343 + !!$
344 + !!$   real(kind=dp) :: this_rlist
345 + !!$
346 + !!$    rlist = this_rlist
347 + !!$    rlistsq = rlist * rlist
348 + !!$
349 + !!$    haveRlist = .true.
350 + !!$
351 + !!$  end subroutine setRlistDF
352 +
353 +
354    subroutine setSimVariables()
355      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
356      SIM_uses_LennardJones = SimUsesLennardJones()
# Line 249 | Line 380 | contains
380  
381      error = 0
382  
383 <    if (.not. havePropertyMap) then
383 >    if (.not. haveInteractionMap) then
384  
385         myStatus = 0
386  
387 <       call createPropertyMap(myStatus)
387 >       call createInteractionMap(myStatus)
388  
389         if (myStatus .ne. 0) then
390 <          write(default_error, *) 'createPropertyMap failed in doForces!'
390 >          write(default_error, *) 'createInteractionMap failed in doForces!'
391            error = -1
392            return
393         endif
# Line 514 | Line 645 | contains
645      integer :: localError
646      integer :: propPack_i, propPack_j
647      integer :: loopStart, loopEnd, loop
648 <
648 >    integer :: iMap
649      real(kind=dp) :: listSkin = 1.0  
650  
651      !! initialize local variables  
# Line 612 | Line 743 | contains
743  
744            if (update_nlist) then
745   #ifdef IS_MPI
746 +             me_i = atid_row(i)
747               jstart = 1
748               jend = nGroupsInCol
749   #else
750 +             me_i = atid(i)
751               jstart = i+1
752               jend = nGroups
753   #endif
# Line 633 | Line 766 | contains
766               endif
767  
768   #ifdef IS_MPI
769 +             me_j = atid_col(j)
770               call get_interatomic_vector(q_group_Row(:,i), &
771                    q_group_Col(:,j), d_grp, rgrpsq)
772   #else
773 +             me_j = atid(j)
774               call get_interatomic_vector(q_group(:,i), &
775                    q_group(:,j), d_grp, rgrpsq)
776   #endif
777  
778 <             if (rgrpsq < rlistsq) then
778 >             if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
779                  if (update_nlist) then
780                     nlist = nlist + 1
781  
# Line 858 | Line 993 | contains
993   #else
994               me_i = atid(i)
995   #endif
996 +             iMap = InteractionMap(me_i, me_j)%InteractionHash
997 +            
998 +             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
999  
862 –             if (PropertyMap(me_i)%is_Dipole) then
863 –
1000                  mu_i = getDipoleMoment(me_i)
1001  
1002                  !! The reaction field needs to include a self contribution
# Line 926 | Line 1062 | contains
1062      real ( kind = dp ) :: ebalance
1063      integer :: me_i, me_j
1064  
1065 +    integer :: iMap
1066 +
1067      r = sqrt(rijsq)
1068      vpair = 0.0d0
1069      fpair(1:3) = 0.0d0
# Line 938 | Line 1076 | contains
1076      me_j = atid(j)
1077   #endif
1078  
1079 <    !    write(*,*) i, j, me_i, me_j
1079 >    iMap = InteractionMap(me_i, me_j)%InteractionHash
1080  
1081 <    if (FF_uses_LennardJones .and. SIM_uses_LennardJones) then
1081 >    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1082 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1083 >    endif
1084  
1085 <       if ( PropertyMap(me_i)%is_LennardJones .and. &
1086 <            PropertyMap(me_j)%is_LennardJones ) then
1087 <          call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1085 >    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1086 >       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1087 >            pot, eFrame, f, t, do_pot)
1088 >
1089 >       if (FF_uses_RF .and. SIM_uses_RF) then
1090 >
1091 >          ! CHECK ME (RF needs to know about all electrostatic types)
1092 >          call accumulate_rf(i, j, r, eFrame, sw)
1093 >          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1094         endif
1095  
1096      endif
1097  
1098 <    if (FF_uses_Electrostatics .and. SIM_uses_Electrostatics) then
1098 >    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1099 >       call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1100 >            pot, A, f, t, do_pot)
1101 >    endif
1102  
1103 <       if (PropertyMap(me_i)%is_Electrostatic .and. &
1104 <            PropertyMap(me_j)%is_Electrostatic) then
1105 <          call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1106 <               pot, eFrame, f, t, do_pot, ebalance)
958 <       endif
1103 >    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1104 >       call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1105 >            pot, A, f, t, do_pot)
1106 >    endif
1107  
1108 <       if (FF_uses_dipoles .and. SIM_uses_dipoles) then      
1109 <          if ( PropertyMap(me_i)%is_Dipole .and. &
1110 <               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
1108 >    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1109 >       call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1110 >            pot, A, f, t, do_pot)
1111      endif
1112 <
1113 <
1114 <    if (FF_uses_Sticky .and. SIM_uses_sticky) then
1115 <
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 <
1112 >    
1113 >    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1114 > !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1115 > !           pot, A, f, t, do_pot)
1116      endif
1117  
1118 <    if (FF_uses_StickyPower .and. SIM_uses_stickypower) then
1119 <       if ( PropertyMap(me_i)%is_StickyPower .and. &
1120 <            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
1118 >    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1119 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1120 >            do_pot)
1121      endif
988 –    
989 –    if (FF_uses_GayBerne .and. SIM_uses_GayBerne) then
1122  
1123 <       if ( PropertyMap(me_i)%is_GayBerne .and. &
1124 <            PropertyMap(me_j)%is_GayBerne) then
1125 <          call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
994 <               pot, A, f, t, do_pot)
995 <       endif
996 <
1123 >    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1124 >       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1125 >            pot, A, f, t, do_pot)
1126      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
1127  
1128 +    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1129 +       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1130 +            pot, A, f, t, do_pot)
1131      endif
1132 <
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 <
1132 >    
1133    end subroutine do_pair
1134  
1135    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
# Line 1040 | Line 1147 | contains
1147      real ( kind = dp )                :: r, rc
1148      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1149  
1150 <    logical :: is_EAM_i, is_EAM_j
1044 <
1045 <    integer :: me_i, me_j
1150 >    integer :: me_i, me_j, iMap
1151  
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 –
1152   #ifdef IS_MPI  
1153      me_i = atid_row(i)
1154      me_j = atid_col(j)  
# Line 1061 | Line 1157 | contains
1157      me_j = atid(j)  
1158   #endif
1159  
1160 <    if (FF_uses_EAM .and. SIM_uses_EAM) then
1160 >    iMap = InteractionMap(me_i, me_j)%InteractionHash
1161  
1162 <       if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1162 >    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1163              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1068 –
1164      endif
1165 <
1165 >    
1166    end subroutine do_prepair
1167  
1168  

Diff Legend

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