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 2261 by gezelter, Tue Jun 28 13:58: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.22 2005-06-28 13:58:45 gezelter Exp $, $Date: 2005-06-28 13:58:45 $, $Name: not supported by cvs2svn $, $Revision: 1.22 $
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  
86    logical, save :: FF_uses_DirectionalAtoms
# Line 115 | Line 115 | module doForces
115    logical, save :: SIM_uses_PBC
116    logical, save :: SIM_uses_molecular_cutoffs
117  
118 <  real(kind=dp), save :: rlist, rlistsq
118 >  !!!GO AWAY---------
119 >  !!!!!real(kind=dp), save :: rlist, rlistsq
120  
121    public :: init_FF
122    public :: do_force_loop
123 <  public :: setRlistDF
123 > !  public :: setRlistDF
124  
125   #ifdef PROFILE
126    public :: getforcetime
# Line 128 | Line 129 | module doForces
129    integer :: nLoops
130   #endif
131  
132 <  type :: Properties
133 <     logical :: is_Directional   = .false.
134 <     logical :: is_LennardJones  = .false.
135 <     logical :: is_Electrostatic = .false.
136 <     logical :: is_Charge        = .false.
137 <     logical :: is_Dipole        = .false.
138 <     logical :: is_Quadrupole    = .false.
139 <     logical :: is_Sticky        = .false.
140 <     logical :: is_StickyPower   = .false.
141 <     logical :: is_GayBerne      = .false.
142 <     logical :: is_EAM           = .false.
143 <     logical :: is_Shape         = .false.
143 <     logical :: is_FLARB         = .false.
144 <  end type Properties
145 <
146 <  type(Properties), dimension(:),allocatable :: PropertyMap
147 <
132 >  type, public :: Interaction
133 >     integer :: InteractionHash
134 >     real(kind=dp) :: rCut
135 >  end type Interaction
136 >  
137 >  type(Interaction), dimension(:,:),allocatable :: InteractionMap
138 >  
139 >  !public :: addInteraction
140 >  !public :: setInteractionHash
141 >  !public :: getInteractionHash
142 >  public :: createInteractionMap
143 >  
144   contains
145  
146 <  subroutine setRlistDF( this_rlist )
147 <
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)
146 >
147 >  subroutine createInteractionMap(status)
148      integer :: nAtypes
149      integer :: status
150      integer :: i
151 <    logical :: thisProperty
152 <    real (kind=DP) :: thisDPproperty
153 <
154 <    status = 0
155 <
151 >    integer :: j
152 >    integer :: ihash
153 >    real(kind=dp) :: myRcut
154 > ! Test Types
155 >    logical :: i_is_LJ
156 >    logical :: i_is_Elect
157 >    logical :: i_is_Sticky
158 >    logical :: i_is_StickyP
159 >    logical :: i_is_GB
160 >    logical :: i_is_EAM
161 >    logical :: i_is_Shape
162 >    logical :: j_is_LJ
163 >    logical :: j_is_Elect
164 >    logical :: j_is_Sticky
165 >    logical :: j_is_StickyP
166 >    logical :: j_is_GB
167 >    logical :: j_is_EAM
168 >    logical :: j_is_Shape
169 >    
170 >    
171 >    if (.not. associated(atypes)) then
172 >       call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
173 >       status = -1
174 >       return
175 >    endif
176 >    
177      nAtypes = getSize(atypes)
178 <
178 >    
179      if (nAtypes == 0) then
180         status = -1
181         return
182      end if
183  
184 <    if (.not. allocated(PropertyMap)) then
185 <       allocate(PropertyMap(nAtypes))
184 >    if (.not. allocated(InteractionMap)) then
185 >       allocate(InteractionMap(nAtypes,nAtypes))
186      endif
187 <
187 >        
188      do i = 1, nAtypes
189 <       call getElementProperty(atypes, i, "is_Directional", thisProperty)
190 <       PropertyMap(i)%is_Directional = thisProperty
189 >       call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
190 >       call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
191 >       call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
192 >       call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
193 >       call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
194 >       call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
195 >       call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
196  
197 <       call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
186 <       PropertyMap(i)%is_LennardJones = thisProperty
197 >       do j = i, nAtypes
198  
199 <       call getElementProperty(atypes, i, "is_Electrostatic", thisProperty)
200 <       PropertyMap(i)%is_Electrostatic = thisProperty
199 >          iHash = 0
200 >          myRcut = 0.0_dp
201  
202 <       call getElementProperty(atypes, i, "is_Charge", thisProperty)
203 <       PropertyMap(i)%is_Charge = thisProperty
202 >          call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ)
203 >          call getElementProperty(atypes, j, "is_Electrostatic", j_is_Elect)
204 >          call getElementProperty(atypes, j, "is_Sticky", j_is_Sticky)
205 >          call getElementProperty(atypes, j, "is_StickyPower", j_is_StickyP)
206 >          call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
207 >          call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
208 >          call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
209  
210 <       call getElementProperty(atypes, i, "is_Dipole", thisProperty)
211 <       PropertyMap(i)%is_Dipole = thisProperty
210 >          if (i_is_LJ .and. j_is_LJ) then
211 >             iHash = ior(iHash, LJ_PAIR)            
212 >          endif
213 >          
214 >          if (i_is_Elect .and. j_is_Elect) then
215 >             iHash = ior(iHash, ELECTROSTATIC_PAIR)
216 >          endif
217 >          
218 >          if (i_is_Sticky .and. j_is_Sticky) then
219 >             iHash = ior(iHash, STICKY_PAIR)
220 >          endif
221  
222 <       call getElementProperty(atypes, i, "is_Quadrupole", thisProperty)
223 <       PropertyMap(i)%is_Quadrupole = thisProperty
222 >          if (i_is_StickyP .and. j_is_StickyP) then
223 >             iHash = ior(iHash, STICKYPOWER_PAIR)
224 >          endif
225  
226 <       call getElementProperty(atypes, i, "is_Sticky", thisProperty)
227 <       PropertyMap(i)%is_Sticky = thisProperty
228 <      
203 <       call getElementProperty(atypes, i, "is_StickyPower", thisProperty)
204 <       PropertyMap(i)%is_StickyPower = thisProperty
226 >          if (i_is_EAM .and. j_is_EAM) then
227 >             iHash = ior(iHash, EAM_PAIR)
228 >          endif
229  
230 <       call getElementProperty(atypes, i, "is_GayBerne", thisProperty)
231 <       PropertyMap(i)%is_GayBerne = thisProperty
230 >          if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
231 >          if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
232 >          if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
233  
234 <       call getElementProperty(atypes, i, "is_EAM", thisProperty)
235 <       PropertyMap(i)%is_EAM = thisProperty
234 >          if (i_is_Shape .and. j_is_Shape) iHash = ior(iHash, SHAPE_PAIR)
235 >          if (i_is_Shape .and. j_is_LJ) iHash = ior(iHash, SHAPE_LJ)
236 >          if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
237  
212 –       call getElementProperty(atypes, i, "is_Shape", thisProperty)
213 –       PropertyMap(i)%is_Shape = thisProperty
238  
239 <       call getElementProperty(atypes, i, "is_FLARB", thisProperty)
240 <       PropertyMap(i)%is_FLARB = thisProperty
239 >          InteractionMap(i,j)%InteractionHash = iHash
240 >          InteractionMap(j,i)%InteractionHash = iHash
241 >
242 >       end do
243 >
244      end do
245 +  end subroutine createInteractionMap
246  
219 –    havePropertyMap = .true.
247  
221 –  end subroutine createPropertyMap
248  
249 + !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
250 + !!$  subroutine setRlistDF( this_rlist )
251 + !!$
252 + !!$   real(kind=dp) :: this_rlist
253 + !!$
254 + !!$    rlist = this_rlist
255 + !!$    rlistsq = rlist * rlist
256 + !!$
257 + !!$    haveRlist = .true.
258 + !!$
259 + !!$  end subroutine setRlistDF
260 +
261 +
262    subroutine setSimVariables()
263      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
264      SIM_uses_LennardJones = SimUsesLennardJones()
# Line 249 | Line 288 | contains
288  
289      error = 0
290  
291 <    if (.not. havePropertyMap) then
291 >    if (.not. haveInteractionMap) then
292  
293         myStatus = 0
294  
295 <       call createPropertyMap(myStatus)
295 >       call createInteractionMap(myStatus)
296  
297         if (myStatus .ne. 0) then
298 <          write(default_error, *) 'createPropertyMap failed in doForces!'
298 >          write(default_error, *) 'createInteractionMap failed in doForces!'
299            error = -1
300            return
301         endif
# Line 858 | Line 897 | contains
897   #else
898               me_i = atid(i)
899   #endif
900 <
901 <             if (PropertyMap(me_i)%is_Dipole) then
900 >             iMap = InteractionMap(me_i, me_j)%InteractionHash
901 >            
902 >             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
903  
904                  mu_i = getDipoleMoment(me_i)
905  
# Line 926 | Line 966 | contains
966      real ( kind = dp ) :: ebalance
967      integer :: me_i, me_j
968  
969 +    integer :: iMap
970 +
971      r = sqrt(rijsq)
972      vpair = 0.0d0
973      fpair(1:3) = 0.0d0
# Line 938 | Line 980 | contains
980      me_j = atid(j)
981   #endif
982  
983 <    !    write(*,*) i, j, me_i, me_j
942 <
943 <    if (FF_uses_LennardJones .and. SIM_uses_LennardJones) then
944 <
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 <
950 <    endif
951 <
952 <    if (FF_uses_Electrostatics .and. SIM_uses_Electrostatics) then
953 <
954 <       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
983 >    iMap = InteractionMap(me_i, me_j)%InteractionHash
984  
985 <       if (FF_uses_dipoles .and. SIM_uses_dipoles) then      
986 <          if ( PropertyMap(me_i)%is_Dipole .and. &
962 <               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
985 >    if ( iand(iMap, LJ_PAIR).ne.0 ) then
986 >       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
987      endif
988  
989 +    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
990 +       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
991 +            pot, eFrame, f, t, do_pot)
992  
993 <    if (FF_uses_Sticky .and. SIM_uses_sticky) then
993 >       if (FF_uses_RF .and. SIM_uses_RF) then
994  
995 <       if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
996 <          call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
997 <               pot, A, f, t, do_pot)
995 >          ! CHECK ME (RF needs to know about all electrostatic types)
996 >          call accumulate_rf(i, j, r, eFrame, sw)
997 >          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
998         endif
999  
1000      endif
1001  
1002 <    if (FF_uses_StickyPower .and. SIM_uses_stickypower) then
1003 <       if ( PropertyMap(me_i)%is_StickyPower .and. &
1004 <            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
1002 >    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
1003 >       call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1004 >            pot, A, f, t, do_pot)
1005      endif
988 –    
989 –    if (FF_uses_GayBerne .and. SIM_uses_GayBerne) then
1006  
1007 <       if ( PropertyMap(me_i)%is_GayBerne .and. &
1008 <            PropertyMap(me_j)%is_GayBerne) then
1009 <          call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
994 <               pot, A, f, t, do_pot)
995 <       endif
996 <
1007 >    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
1008 >       call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1009 >            pot, A, f, t, do_pot)
1010      endif
1011  
1012 <    if (FF_uses_EAM .and. SIM_uses_EAM) then
1013 <
1014 <       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 <
1012 >    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
1013 >       call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1014 >            pot, A, f, t, do_pot)
1015      endif
1016 +    
1017 +    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1018 +       call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1019 +            pot, A, f, t, do_pot)
1020 +    endif
1021  
1022 <
1023 <    !    write(*,*) PropertyMap(me_i)%is_Shape,PropertyMap(me_j)%is_Shape
1024 <
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
1022 >    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1023 >       call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1024 >            do_pot)
1025      endif
1026  
1027 +    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
1028 +       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1029 +            pot, A, f, t, do_pot)
1030 +    endif
1031 +
1032 +    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1033 +       call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1034 +            pot, A, f, t, do_pot)
1035 +    endif
1036 +    
1037    end subroutine do_pair
1038  
1039    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
# Line 1040 | Line 1051 | contains
1051      real ( kind = dp )                :: r, rc
1052      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1053  
1054 <    logical :: is_EAM_i, is_EAM_j
1044 <
1045 <    integer :: me_i, me_j
1046 <
1047 <
1048 <    r = sqrt(rijsq)
1049 <    if (SIM_uses_molecular_cutoffs) then
1050 <       rc = sqrt(rcijsq)
1051 <    else
1052 <       rc = r
1053 <    endif
1054 >    integer :: me_i, me_j, iMap
1055  
1055 –
1056   #ifdef IS_MPI  
1057      me_i = atid_row(i)
1058      me_j = atid_col(j)  
# Line 1061 | Line 1061 | contains
1061      me_j = atid(j)  
1062   #endif
1063  
1064 <    if (FF_uses_EAM .and. SIM_uses_EAM) then
1065 <
1066 <       if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1067 <            call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1064 >    iMap = InteractionMap(me_i, me_j)%InteractionHash
1065  
1066 +    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1067 +            call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1068      endif
1069 <
1069 >    
1070    end subroutine do_prepair
1071  
1072  

Diff Legend

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