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 2211 by chrisfen, Thu Apr 21 14:12:19 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.14 2005-04-21 14:12:19 chrisfen Exp $, $Date: 2005-04-21 14:12:19 $, $Name: not supported by cvs2svn $, $Revision: 1.14 $
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 89 | Line 90 | module doForces
90    logical, save :: FF_uses_Charges
91    logical, save :: FF_uses_Dipoles
92    logical, save :: FF_uses_Quadrupoles
93 <  logical, save :: FF_uses_sticky
93 >  logical, save :: FF_uses_Sticky
94 >  logical, save :: FF_uses_StickyPower
95    logical, save :: FF_uses_GayBerne
96    logical, save :: FF_uses_EAM
97    logical, save :: FF_uses_Shapes
# Line 103 | Line 105 | module doForces
105    logical, save :: SIM_uses_Dipoles
106    logical, save :: SIM_uses_Quadrupoles
107    logical, save :: SIM_uses_Sticky
108 +  logical, save :: SIM_uses_StickyPower
109    logical, save :: SIM_uses_GayBerne
110    logical, save :: SIM_uses_EAM
111    logical, save :: SIM_uses_Shapes
# Line 113 | 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 134 | Line 138 | module doForces
138       logical :: is_Dipole        = .false.
139       logical :: is_Quadrupole    = .false.
140       logical :: is_Sticky        = .false.
141 +     logical :: is_StickyPower   = .false.
142       logical :: is_GayBerne      = .false.
143       logical :: is_EAM           = .false.
144       logical :: is_Shape         = .false.
# Line 142 | Line 147 | contains
147  
148    type(Properties), dimension(:),allocatable :: PropertyMap
149  
145 – 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 196 | Line 314 | contains
314  
315         call getElementProperty(atypes, i, "is_Sticky", thisProperty)
316         PropertyMap(i)%is_Sticky = thisProperty
317 +      
318 +       call getElementProperty(atypes, i, "is_StickyPower", thisProperty)
319 +       PropertyMap(i)%is_StickyPower = thisProperty
320  
321         call getElementProperty(atypes, i, "is_GayBerne", thisProperty)
322         PropertyMap(i)%is_GayBerne = thisProperty
# Line 221 | Line 342 | contains
342      SIM_uses_Charges = SimUsesCharges()
343      SIM_uses_Dipoles = SimUsesDipoles()
344      SIM_uses_Sticky = SimUsesSticky()
345 +    SIM_uses_StickyPower = SimUsesStickyPower()
346      SIM_uses_GayBerne = SimUsesGayBerne()
347      SIM_uses_EAM = SimUsesEAM()
348      SIM_uses_Shapes = SimUsesShapes()
# Line 315 | Line 437 | contains
437      FF_uses_Charges = .false.    
438      FF_uses_Dipoles = .false.
439      FF_uses_Sticky = .false.
440 +    FF_uses_StickyPower = .false.
441      FF_uses_GayBerne = .false.
442      FF_uses_EAM = .false.
443      FF_uses_Shapes = .false.
# Line 364 | Line 487 | contains
487         FF_uses_DirectionalAtoms = .true.
488      endif
489  
490 +    call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
491 +         MatchList)
492 +    if (nMatches .gt. 0) then
493 +       FF_uses_StickyPower = .true.
494 +       FF_uses_DirectionalAtoms = .true.
495 +    endif
496 +    
497      call getMatchingElementList(atypes, "is_GayBerne", .true., &
498           nMatches, MatchList)
499      if (nMatches .gt. 0) then
# Line 908 | 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 922 | 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 <
929 <       if ( PropertyMap(me_i)%is_LennardJones .and. &
930 <            PropertyMap(me_j)%is_LennardJones ) then
931 <          call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
932 <       endif
933 <
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
1065 <
1066 <       if (PropertyMap(me_i)%is_Electrostatic .and. &
939 <            PropertyMap(me_j)%is_Electrostatic) then
940 <          call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
941 <               pot, eFrame, f, t, do_pot)
942 <       endif
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 (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 952 | Line 1076 | contains
1076         endif
1077      endif
1078  
1079 <
1080 <    if (FF_uses_Sticky .and. SIM_uses_sticky) then
1081 <
958 <       if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
959 <          call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
960 <               pot, A, f, t, do_pot)
961 <       endif
962 <
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 <
1085 <    if (FF_uses_GayBerne .and. SIM_uses_GayBerne) then
1086 <
968 <       if ( PropertyMap(me_i)%is_GayBerne .and. &
969 <            PropertyMap(me_j)%is_GayBerne) then
970 <          call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
971 <               pot, A, f, t, do_pot)
972 <       endif
973 <
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_EAM .and. SIM_uses_EAM) then
1090 <
1091 <       if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
1092 <          call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1093 <               do_pot)
1094 <       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 ( 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 ( 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 <
988 <    if (FF_uses_Shapes .and. SIM_uses_Shapes) then
989 <       if ( PropertyMap(me_i)%is_Shape .and. &
990 <            PropertyMap(me_j)%is_Shape ) then
991 <          call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
992 <               pot, A, f, t, do_pot)
993 <       endif
994 <       if ( (PropertyMap(me_i)%is_Shape .and. &
995 <            PropertyMap(me_j)%is_LennardJones) .or. &
996 <            (PropertyMap(me_i)%is_LennardJones .and. &
997 <            PropertyMap(me_j)%is_Shape) ) then
998 <          call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
999 <               pot, A, f, t, do_pot)
1000 <       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 1017 | 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
1021 <
1022 <    integer :: me_i, me_j
1023 <
1024 <
1025 <    r = sqrt(rijsq)
1026 <    if (SIM_uses_molecular_cutoffs) then
1027 <       rc = sqrt(rcijsq)
1028 <    else
1029 <       rc = r
1030 <    endif
1031 <
1131 >    integer :: me_i, me_j, iMap
1132  
1133   #ifdef IS_MPI  
1134      me_i = atid_row(i)
# Line 1038 | 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 )
1045 –
1145      endif
1146 <
1146 >    
1147    end subroutine do_prepair
1148  
1149  
# Line 1240 | Line 1339 | contains
1339      logical :: doesit
1340      doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1341           FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1342 <         FF_uses_GayBerne .or. FF_uses_Shapes
1342 >         FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1343    end function FF_UsesDirectionalAtoms
1344  
1345    function FF_RequiresPrepairCalc() result(doesit)

Diff Legend

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