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 2259 by gezelter, Mon Jun 27 21:01:36 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.20 2005-06-27 21:01:30 gezelter Exp $, $Date: 2005-06-27 21:01:30 $, $Name: not supported by cvs2svn $, $Revision: 1.20 $
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 81 | Line 81 | module doForces
81    logical, save :: haveRlist = .false.
82    logical, save :: haveNeighborList = .false.
83    logical, save :: haveSIMvariables = .false.
84 –  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 116 | 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 129 | 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.
140 <     logical :: is_StickyPower   = .false.
141 <     logical :: is_GayBerne      = .false.
142 <     logical :: is_EAM           = .false.
143 <     logical :: is_Shape         = .false.
144 <     logical :: is_FLARB         = .false.
145 <  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
148 <
147 >  
148   contains
149  
150 <  subroutine setRlistDF( this_rlist )
151 <
153 <    real(kind=dp) :: this_rlist
154 <
155 <    rlist = this_rlist
156 <    rlistsq = rlist * rlist
157 <
158 <    haveRlist = .true.
159 <
160 <  end subroutine setRlistDF
161 <
162 <  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)
187 <       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
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_Dipole", thisProperty)
216 <       PropertyMap(i)%is_Dipole = 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_Quadrupole", thisProperty)
228 <       PropertyMap(i)%is_Quadrupole = 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_Sticky", thisProperty)
232 <       PropertyMap(i)%is_Sticky = thisProperty
233 <      
204 <       call getElementProperty(atypes, i, "is_StickyPower", thisProperty)
205 <       PropertyMap(i)%is_StickyPower = 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_GayBerne", thisProperty)
236 <       PropertyMap(i)%is_GayBerne = 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_EAM", thisProperty)
240 <       PropertyMap(i)%is_EAM = 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  
213 –       call getElementProperty(atypes, i, "is_Shape", thisProperty)
214 –       PropertyMap(i)%is_Shape = thisProperty
243  
244 <       call getElementProperty(atypes, i, "is_FLARB", thisProperty)
245 <       PropertyMap(i)%is_FLARB = thisProperty
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 250 | 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 515 | 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 738 | contains
738         iend = nGroups - 1
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  
# Line 634 | 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 859 | 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 951 | Line 1092 | contains
1092         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1093              pot, eFrame, f, t, do_pot)
1094  
1095 <       if (FF_uses_dipoles .and. SIM_uses_dipoles) then                
1096 <          if ( PropertyMap(me_i)%is_Dipole .and. &
1097 <               PropertyMap(me_j)%is_Dipole) then
1098 <             if (FF_uses_RF .and. SIM_uses_RF) then
1099 <                call accumulate_rf(i, j, r, eFrame, sw)
959 <                call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
960 <             endif
961 <          endif
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 ( iand(iMap, STICKY_PAIR).ne.0 ) then
# Line 978 | Line 1117 | contains
1117      endif
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)
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 ( iand(iMap, EAM_PAIR).ne.0 ) then      

Diff Legend

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