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 2260 by chuckv, Mon Jun 27 22:21:37 2005 UTC vs.
Revision 2268 by gezelter, Fri Jul 29 19:38:27 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
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 $
48 > !! @version $Id: doForces.F90,v 1.26 2005-07-29 19:38:27 gezelter Exp $, $Date: 2005-07-29 19:38:27 $, $Name: not supported by cvs2svn $, $Revision: 1.26 $
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 122 | Line 122 | module doForces
122    public :: init_FF
123    public :: do_force_loop
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 130 | Line 135 | module doForces
135    integer :: nLoops
136   #endif
137  
133 –  type :: Properties
134 –     logical :: is_Directional   = .false.
135 –     logical :: is_LennardJones  = .false.
136 –     logical :: is_Electrostatic = .false.
137 –     logical :: is_Charge        = .false.
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.
145 –     logical :: is_FLARB         = .false.
146 –  end type Properties
147 –
148 –  type(Properties), dimension(:),allocatable :: PropertyMap
149 –
150 –
151 –  
138    type, public :: Interaction
139       integer :: InteractionHash
140 <     real(kind=dp) :: rCut
140 >     real(kind=dp) :: rCut = 0.0_dp
141 >     real(kind=dp) :: rCutSq = 0.0_dp    
142 >     real(kind=dp) :: rListSq = 0.0_dp
143    end type Interaction
144    
145 <  type(Interaction), public, dimension(:,:), allocatable :: InteractionMap
145 >  type(Interaction), dimension(:,:),allocatable :: InteractionMap
146    
159 –  !public :: addInteraction
160 –  !public :: setInteractionHash
161 –  !public :: getInteractionHash
162 –  public :: createInteractionMap
147  
148 +  
149   contains
150  
151  
152    subroutine createInteractionMap(status)
153      integer :: nAtypes
154 <    integer :: status
154 >    integer, intent(out) :: status
155      integer :: i
156      integer :: j
157      integer :: ihash
158      real(kind=dp) :: myRcut
159 < ! Test Types
159 >    !! Test Types
160      logical :: i_is_LJ
161      logical :: i_is_Elect
162      logical :: i_is_Sticky
# Line 187 | Line 172 | contains
172      logical :: j_is_EAM
173      logical :: j_is_Shape
174      
175 <    
175 >    status = 0  
176 >
177      if (.not. associated(atypes)) then
178         call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
179         status = -1
# Line 228 | Line 214 | contains
214            call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
215  
216            if (i_is_LJ .and. j_is_LJ) then
217 <             iHash = ior(iHash, LJ_PAIR)
218 <            
217 >             iHash = ior(iHash, LJ_PAIR)            
218 >          endif
219 >          
220 >          if (i_is_Elect .and. j_is_Elect) then
221 >             iHash = ior(iHash, ELECTROSTATIC_PAIR)
222 >          endif
223 >          
224 >          if (i_is_Sticky .and. j_is_Sticky) then
225 >             iHash = ior(iHash, STICKY_PAIR)
226 >          endif
227  
228 <
228 >          if (i_is_StickyP .and. j_is_StickyP) then
229 >             iHash = ior(iHash, STICKYPOWER_PAIR)
230            endif
231  
232 <
233 <
234 <          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)
232 >          if (i_is_EAM .and. j_is_EAM) then
233 >             iHash = ior(iHash, EAM_PAIR)
234 >          endif
235  
236            if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
237            if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
# Line 257 | Line 248 | contains
248         end do
249  
250      end do
251 +
252 +    haveInteractionMap = .true.
253    end subroutine createInteractionMap
254  
255 +  ! Query each potential and return the cutoff for that potential. We
256 +  ! build the neighbor list based on the largest cutoff value for that
257 +  ! atype. Each potential can decide whether to calculate the force for
258 +  ! that atype based upon it's own cutoff.
259 +  
260 +  subroutine createRcuts(defaultRcut, defaultSkinThickness, stat)
261  
262 +    real(kind=dp), intent(in), optional :: defaultRCut, defaultSkinThickness
263 +    integer :: iMap
264 +    integer :: map_i,map_j
265 +    real(kind=dp) :: thisRCut = 0.0_dp
266 +    real(kind=dp) :: actualCutoff = 0.0_dp
267 +    integer, intent(out) :: stat
268 +    integer :: nAtypes
269 +    integer :: myStatus
270  
271 +    stat = 0
272 +    if (.not. haveInteractionMap) then
273 +
274 +       call createInteractionMap(myStatus)
275 +
276 +       if (myStatus .ne. 0) then
277 +          write(default_error, *) 'createInteractionMap failed in doForces!'
278 +          stat = -1
279 +          return
280 +       endif
281 +    endif
282 +
283 +    nAtypes = getSize(atypes)
284 +    !! If we pass a default rcut, set all atypes to that cutoff distance
285 +    if(present(defaultRList)) then
286 +       InteractionMap(:,:)%rCut = defaultRCut
287 +       InteractionMap(:,:)%rCutSq = defaultRCut*defaultRCut
288 +       InteractionMap(:,:)%rListSq = (defaultRCut+defaultSkinThickness)**2
289 +       haveRlist = .true.
290 +       return
291 +    end if
292 +
293 +    do map_i = 1,nAtypes
294 +       do map_j = map_i,nAtypes
295 +          iMap = InteractionMap(map_i, map_j)%InteractionHash
296 +          
297 +          if ( iand(iMap, LJ_PAIR).ne.0 ) then
298 +             ! thisRCut = getLJCutOff(map_i,map_j)
299 +             if (thisRcut > actualCutoff) actualCutoff = thisRcut
300 +          endif
301 +          
302 +          if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
303 +             ! thisRCut = getElectrostaticCutOff(map_i,map_j)
304 +             if (thisRcut > actualCutoff) actualCutoff = thisRcut
305 +          endif
306 +          
307 +          if ( iand(iMap, STICKY_PAIR).ne.0 ) then
308 +             ! thisRCut = getStickyCutOff(map_i,map_j)
309 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
310 +           endif
311 +          
312 +           if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
313 +              ! thisRCut = getStickyPowerCutOff(map_i,map_j)
314 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
315 +           endif
316 +          
317 +           if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
318 +              ! thisRCut = getGayberneCutOff(map_i,map_j)
319 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
320 +           endif
321 +          
322 +           if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
323 + !              thisRCut = getGaybrneLJCutOff(map_i,map_j)
324 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
325 +           endif
326 +          
327 +           if ( iand(iMap, EAM_PAIR).ne.0 ) then      
328 + !              thisRCut = getEAMCutOff(map_i,map_j)
329 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
330 +           endif
331 +          
332 +           if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
333 + !              thisRCut = getShapeCutOff(map_i,map_j)
334 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
335 +           endif
336 +          
337 +           if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
338 + !              thisRCut = getShapeLJCutOff(map_i,map_j)
339 +              if (thisRcut > actualCutoff) actualCutoff = thisRcut
340 +           endif
341 +           InteractionMap(map_i, map_j)%rCut = actualCutoff
342 +           InteractionMap(map_i, map_j)%rCutSq = actualCutoff * actualCutoff
343 +           InteractionMap(map_i, map_j)%rListSq = (actualCutoff + skinThickness)**2
344 +
345 +           InteractionMap(map_j, map_i)%rCut = InteractionMap(map_i, map_j)%rCut
346 +           InteractionMap(map_j, map_i)%rCutSq = InteractionMap(map_i, map_j)%rCutSq
347 +           InteractionMap(map_j, map_i)%rListSq = InteractionMap(map_i, map_j)%rListSq
348 +        end do
349 +     end do
350 +     ! now the groups
351 +
352 +
353 +
354 +     haveRlist = .true.
355 +  end subroutine createRcuts
356 +
357 +
358   !!! THIS GOES AWAY FOR SIZE DEPENDENT CUTOFF
359   !!$  subroutine setRlistDF( this_rlist )
360   !!$
# Line 272 | Line 366 | contains
366   !!$    haveRlist = .true.
367   !!$
368   !!$  end subroutine setRlistDF
275 –
276 –  subroutine createPropertyMap(status)
277 –    integer :: nAtypes
278 –    integer :: status
279 –    integer :: i
280 –    logical :: thisProperty
281 –    real (kind=DP) :: thisDPproperty
282 –
283 –    status = 0
284 –
285 –    nAtypes = getSize(atypes)
286 –
287 –    if (nAtypes == 0) then
288 –       status = -1
289 –       return
290 –    end if
291 –
292 –    if (.not. allocated(PropertyMap)) then
293 –       allocate(PropertyMap(nAtypes))
294 –    endif
295 –
296 –    do i = 1, nAtypes
297 –       call getElementProperty(atypes, i, "is_Directional", thisProperty)
298 –       PropertyMap(i)%is_Directional = thisProperty
299 –
300 –       call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
301 –       PropertyMap(i)%is_LennardJones = thisProperty
302 –
303 –       call getElementProperty(atypes, i, "is_Electrostatic", thisProperty)
304 –       PropertyMap(i)%is_Electrostatic = thisProperty
305 –
306 –       call getElementProperty(atypes, i, "is_Charge", thisProperty)
307 –       PropertyMap(i)%is_Charge = thisProperty
308 –
309 –       call getElementProperty(atypes, i, "is_Dipole", thisProperty)
310 –       PropertyMap(i)%is_Dipole = thisProperty
311 –
312 –       call getElementProperty(atypes, i, "is_Quadrupole", thisProperty)
313 –       PropertyMap(i)%is_Quadrupole = thisProperty
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
369  
321 –       call getElementProperty(atypes, i, "is_GayBerne", thisProperty)
322 –       PropertyMap(i)%is_GayBerne = thisProperty
370  
324 –       call getElementProperty(atypes, i, "is_EAM", thisProperty)
325 –       PropertyMap(i)%is_EAM = thisProperty
326 –
327 –       call getElementProperty(atypes, i, "is_Shape", thisProperty)
328 –       PropertyMap(i)%is_Shape = thisProperty
329 –
330 –       call getElementProperty(atypes, i, "is_FLARB", thisProperty)
331 –       PropertyMap(i)%is_FLARB = thisProperty
332 –    end do
333 –
334 –    havePropertyMap = .true.
335 –
336 –  end subroutine createPropertyMap
337 –
371    subroutine setSimVariables()
372      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
373      SIM_uses_LennardJones = SimUsesLennardJones()
# Line 364 | Line 397 | contains
397  
398      error = 0
399  
400 <    if (.not. havePropertyMap) then
401 <
402 <       myStatus = 0
403 <
404 <       call createPropertyMap(myStatus)
372 <
400 >    if (.not. haveInteractionMap) then
401 >      
402 >       myStatus = 0      
403 >       call createInteractionMap(myStatus)
404 >      
405         if (myStatus .ne. 0) then
406 <          write(default_error, *) 'createPropertyMap failed in doForces!'
406 >          write(default_error, *) 'createInteractionMap failed in doForces!'
407            error = -1
408            return
409         endif
# Line 629 | Line 661 | contains
661      integer :: localError
662      integer :: propPack_i, propPack_j
663      integer :: loopStart, loopEnd, loop
664 <
664 >    integer :: iMap
665      real(kind=dp) :: listSkin = 1.0  
666  
667      !! initialize local variables  
# Line 721 | Line 753 | contains
753   #endif
754         outer: do i = istart, iend
755  
756 + #ifdef IS_MPI
757 +             me_i = atid_row(i)
758 + #else
759 +             me_i = atid(i)
760 + #endif
761 +
762            if (update_nlist) point(i) = nlist + 1
763  
764            n_in_i = groupStartRow(i+1) - groupStartRow(i)
# Line 748 | Line 786 | contains
786               endif
787  
788   #ifdef IS_MPI
789 +             me_j = atid_col(j)
790               call get_interatomic_vector(q_group_Row(:,i), &
791                    q_group_Col(:,j), d_grp, rgrpsq)
792   #else
793 +             me_j = atid(j)
794               call get_interatomic_vector(q_group(:,i), &
795                    q_group(:,j), d_grp, rgrpsq)
796   #endif
797  
798 <             if (rgrpsq < rlistsq) then
798 >             if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then
799                  if (update_nlist) then
800                     nlist = nlist + 1
801  
# Line 973 | Line 1013 | contains
1013   #else
1014               me_i = atid(i)
1015   #endif
1016 <
1017 <             if (PropertyMap(me_i)%is_Dipole) then
1016 >             iMap = InteractionMap(me_i, me_j)%InteractionHash
1017 >            
1018 >             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
1019  
1020                  mu_i = getDipoleMoment(me_i)
1021  
# Line 1065 | Line 1106 | contains
1106         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1107              pot, eFrame, f, t, do_pot)
1108  
1109 <       if (FF_uses_dipoles .and. SIM_uses_dipoles) then                
1110 <          if ( PropertyMap(me_i)%is_Dipole .and. &
1111 <               PropertyMap(me_j)%is_Dipole) then
1112 <             if (FF_uses_RF .and. SIM_uses_RF) then
1113 <                call accumulate_rf(i, j, r, eFrame, sw)
1073 <                call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1074 <             endif
1075 <          endif
1109 >       if (FF_uses_RF .and. SIM_uses_RF) then
1110 >
1111 >          ! CHECK ME (RF needs to know about all electrostatic types)
1112 >          call accumulate_rf(i, j, r, eFrame, sw)
1113 >          call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
1114         endif
1115 +
1116      endif
1117  
1118      if ( iand(iMap, STICKY_PAIR).ne.0 ) then
# Line 1092 | Line 1131 | contains
1131      endif
1132      
1133      if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
1134 <       call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1135 <            pot, A, f, t, do_pot)
1134 > !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1135 > !           pot, A, f, t, do_pot)
1136      endif
1137  
1138      if ( iand(iMap, EAM_PAIR).ne.0 ) then      

Diff Legend

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