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

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2229 by chrisfen, Tue May 17 22:35:01 2005 UTC vs.
Revision 2260 by chuckv, Mon Jun 27 22:21: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.17 2005-05-17 22:35:01 chrisfen Exp $, $Date: 2005-05-17 22:35:01 $, $Name: not supported by cvs2svn $, $Revision: 1.17 $
48 > !! @version $Id: doForces.F90,v 1.21 2005-06-27 22:21:37 chuckv Exp $, $Date: 2005-06-27 22:21:37 $, $Name: not supported by cvs2svn $, $Revision: 1.21 $
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 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  
126   #ifdef PROFILE
127    public :: getforcetime
# Line 145 | Line 147 | contains
147  
148    type(Properties), dimension(:),allocatable :: PropertyMap
149  
148 – contains
150  
151 <  subroutine setRlistDF( this_rlist )
151 >  
152 >  type, public :: Interaction
153 >     integer :: InteractionHash
154 >     real(kind=dp) :: rCut
155 >  end type Interaction
156 >  
157 >  type(Interaction), public, dimension(:,:), allocatable :: InteractionMap
158 >  
159 >  !public :: addInteraction
160 >  !public :: setInteractionHash
161 >  !public :: getInteractionHash
162 >  public :: createInteractionMap
163  
164 <    real(kind=dp) :: this_rlist
164 > contains
165  
166 <    rlist = this_rlist
167 <    rlistsq = rlist * rlist
166 >
167 >  subroutine createInteractionMap(status)
168 >    integer :: nAtypes
169 >    integer :: status
170 >    integer :: i
171 >    integer :: j
172 >    integer :: ihash
173 >    real(kind=dp) :: myRcut
174 > ! Test Types
175 >    logical :: i_is_LJ
176 >    logical :: i_is_Elect
177 >    logical :: i_is_Sticky
178 >    logical :: i_is_StickyP
179 >    logical :: i_is_GB
180 >    logical :: i_is_EAM
181 >    logical :: i_is_Shape
182 >    logical :: j_is_LJ
183 >    logical :: j_is_Elect
184 >    logical :: j_is_Sticky
185 >    logical :: j_is_StickyP
186 >    logical :: j_is_GB
187 >    logical :: j_is_EAM
188 >    logical :: j_is_Shape
189 >    
190 >    
191 >    if (.not. associated(atypes)) then
192 >       call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
193 >       status = -1
194 >       return
195 >    endif
196 >    
197 >    nAtypes = getSize(atypes)
198 >    
199 >    if (nAtypes == 0) then
200 >       status = -1
201 >       return
202 >    end if
203  
204 <    haveRlist = .true.
204 >    if (.not. allocated(InteractionMap)) then
205 >       allocate(InteractionMap(nAtypes,nAtypes))
206 >    endif
207 >        
208 >    do i = 1, nAtypes
209 >       call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
210 >       call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
211 >       call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
212 >       call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
213 >       call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
214 >       call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
215 >       call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
216  
217 <  end subroutine setRlistDF
217 >       do j = i, nAtypes
218 >
219 >          iHash = 0
220 >          myRcut = 0.0_dp
221 >
222 >          call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ)
223 >          call getElementProperty(atypes, j, "is_Electrostatic", j_is_Elect)
224 >          call getElementProperty(atypes, j, "is_Sticky", j_is_Sticky)
225 >          call getElementProperty(atypes, j, "is_StickyPower", j_is_StickyP)
226 >          call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
227 >          call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
228 >          call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
229 >
230 >          if (i_is_LJ .and. j_is_LJ) then
231 >             iHash = ior(iHash, LJ_PAIR)
232 >            
233 >
234 >
235 >          endif
236 >
237 >
238 >
239 >          if (i_is_Elect .and. j_is_Elect) iHash = ior(iHash, ELECTROSTATIC_PAIR)
240 >          if (i_is_Sticky .and. j_is_Sticky) iHash = ior(iHash, STICKY_PAIR)
241 >          if (i_is_StickyP .and. j_is_StickyP) iHash = ior(iHash, STICKYPOWER_PAIR)
242 >
243 >          if (i_is_EAM .and. j_is_EAM) iHash = ior(iHash, EAM_PAIR)
244 >
245 >          if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
246 >          if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
247 >          if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
248  
249 +          if (i_is_Shape .and. j_is_Shape) iHash = ior(iHash, SHAPE_PAIR)
250 +          if (i_is_Shape .and. j_is_LJ) iHash = ior(iHash, SHAPE_LJ)
251 +          if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
252 +
253 +
254 +          InteractionMap(i,j)%InteractionHash = iHash
255 +          InteractionMap(j,i)%InteractionHash = iHash
256 +
257 +       end do
258 +
259 +    end do
260 +  end subroutine createInteractionMap
261 +
262 +
263 +
264 + !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
265 + !!$  subroutine setRlistDF( this_rlist )
266 + !!$
267 + !!$   real(kind=dp) :: this_rlist
268 + !!$
269 + !!$    rlist = this_rlist
270 + !!$    rlistsq = rlist * rlist
271 + !!$
272 + !!$    haveRlist = .true.
273 + !!$
274 + !!$  end subroutine setRlistDF
275 +
276    subroutine createPropertyMap(status)
277      integer :: nAtypes
278      integer :: status
# Line 923 | Line 1038 | contains
1038      real ( kind = dp ), intent(inout) :: rijsq
1039      real ( kind = dp )                :: r
1040      real ( kind = dp ), intent(inout) :: d(3)
1041 +    real ( kind = dp ) :: ebalance
1042      integer :: me_i, me_j
1043  
1044 +    integer :: iMap
1045 +
1046      r = sqrt(rijsq)
1047      vpair = 0.0d0
1048      fpair(1:3) = 0.0d0
# Line 937 | Line 1055 | contains
1055      me_j = atid(j)
1056   #endif
1057  
1058 <    !    write(*,*) i, j, me_i, me_j
1058 >    iMap = InteractionMap(me_i, me_j)%InteractionHash
1059  
1060 <    if (FF_uses_LennardJones .and. SIM_uses_LennardJones) then
1061 <
944 <       if ( PropertyMap(me_i)%is_LennardJones .and. &
945 <            PropertyMap(me_j)%is_LennardJones ) then
946 <          call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
947 <       endif
948 <
1060 >    if ( iand(iMap, LJ_PAIR).ne.0 ) then
1061 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1062      endif
1063  
1064 <    if (FF_uses_Electrostatics .and. SIM_uses_Electrostatics) then
1064 >    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1065 >       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1066 >            pot, eFrame, f, t, do_pot)
1067  
1068 <       if (PropertyMap(me_i)%is_Electrostatic .and. &
954 <            PropertyMap(me_j)%is_Electrostatic) then
955 <          call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
956 <               pot, eFrame, f, t, do_pot)
957 <       endif
958 <
959 <       if (FF_uses_dipoles .and. SIM_uses_dipoles) then      
1068 >       if (FF_uses_dipoles .and. SIM_uses_dipoles) then                
1069            if ( PropertyMap(me_i)%is_Dipole .and. &
1070                 PropertyMap(me_j)%is_Dipole) then
1071               if (FF_uses_RF .and. SIM_uses_RF) then
# Line 967 | Line 1076 | contains
1076         endif
1077      endif
1078  
1079 +    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1080 +       call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1081 +            pot, A, f, t, do_pot)
1082 +    endif
1083  
1084 <    if (FF_uses_Sticky .and. SIM_uses_sticky) then
1085 <
1086 <       if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
974 <          call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
975 <               pot, A, f, t, do_pot)
976 <       endif
977 <
1084 >    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1085 >       call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1086 >            pot, A, f, t, do_pot)
1087      endif
1088  
1089 <    if (FF_uses_StickyPower .and. SIM_uses_stickypower) then
1090 <       if ( PropertyMap(me_i)%is_StickyPower .and. &
1091 <            PropertyMap(me_j)%is_StickyPower) then
983 <          call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
984 <               pot, A, f, t, do_pot)
985 <       endif
1089 >    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1090 >       call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1091 >            pot, A, f, t, do_pot)
1092      endif
1093      
1094 <    if (FF_uses_GayBerne .and. SIM_uses_GayBerne) then
1095 <
1096 <       if ( PropertyMap(me_i)%is_GayBerne .and. &
991 <            PropertyMap(me_j)%is_GayBerne) then
992 <          call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
993 <               pot, A, f, t, do_pot)
994 <       endif
995 <
1094 >    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1095 >       call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1096 >            pot, A, f, t, do_pot)
1097      endif
1098  
1099 <    if (FF_uses_EAM .and. SIM_uses_EAM) then
1100 <
1101 <       if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
1001 <          call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1002 <               do_pot)
1003 <       endif
1004 <
1099 >    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1100 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1101 >            do_pot)
1102      endif
1103  
1104 <
1105 <    !    write(*,*) PropertyMap(me_i)%is_Shape,PropertyMap(me_j)%is_Shape
1106 <
1010 <    if (FF_uses_Shapes .and. SIM_uses_Shapes) then
1011 <       if ( PropertyMap(me_i)%is_Shape .and. &
1012 <            PropertyMap(me_j)%is_Shape ) then
1013 <          call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1014 <               pot, A, f, t, do_pot)
1015 <       endif
1016 <       if ( (PropertyMap(me_i)%is_Shape .and. &
1017 <            PropertyMap(me_j)%is_LennardJones) .or. &
1018 <            (PropertyMap(me_i)%is_LennardJones .and. &
1019 <            PropertyMap(me_j)%is_Shape) ) then
1020 <          call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1021 <               pot, A, f, t, do_pot)
1022 <       endif
1104 >    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1105 >       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1106 >            pot, A, f, t, do_pot)
1107      endif
1108  
1109 +    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1110 +       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1111 +            pot, A, f, t, do_pot)
1112 +    endif
1113 +    
1114    end subroutine do_pair
1115  
1116    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
# Line 1039 | Line 1128 | contains
1128      real ( kind = dp )                :: r, rc
1129      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1130  
1131 <    logical :: is_EAM_i, is_EAM_j
1131 >    integer :: me_i, me_j, iMap
1132  
1044 –    integer :: me_i, me_j
1045 –
1046 –
1047 –    r = sqrt(rijsq)
1048 –    if (SIM_uses_molecular_cutoffs) then
1049 –       rc = sqrt(rcijsq)
1050 –    else
1051 –       rc = r
1052 –    endif
1053 –
1054 –
1133   #ifdef IS_MPI  
1134      me_i = atid_row(i)
1135      me_j = atid_col(j)  
# Line 1060 | Line 1138 | contains
1138      me_j = atid(j)  
1139   #endif
1140  
1141 <    if (FF_uses_EAM .and. SIM_uses_EAM) then
1141 >    iMap = InteractionMap(me_i, me_j)%InteractionHash
1142  
1143 <       if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1143 >    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1144              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1067 –
1145      endif
1146 <
1146 >    
1147    end subroutine do_prepair
1148  
1149  

Diff Legend

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