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 2261 by gezelter, Tue Jun 28 13:58:45 2005 UTC vs.
Revision 2273 by gezelter, Thu Aug 11 21:04:03 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
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 $
48 > !! @version $Id: doForces.F90,v 1.29 2005-08-11 21:04:03 gezelter Exp $, $Date: 2005-08-11 21:04:03 $, $Name: not supported by cvs2svn $, $Revision: 1.29 $
49  
50  
51   module doForces
# Line 73 | Line 73 | module doForces
73  
74   #define __FORTRAN90
75   #include "UseTheForce/fSwitchingFunction.h"
76 + #include "UseTheForce/fCutoffPolicy.h"
77   #include "UseTheForce/DarkSide/fInteractionMap.h"
78  
79 +
80    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
81    INTEGER, PARAMETER:: PAIR_LOOP    = 2
82  
81 –  logical, save :: haveRlist = .false.
83    logical, save :: haveNeighborList = .false.
84    logical, save :: haveSIMvariables = .false.
85    logical, save :: haveSaneForceField = .false.
86 +  logical, save :: haveInteractionHash = .false.
87 +  logical, save :: haveGtypeCutoffMap = .false.
88  
89    logical, save :: FF_uses_DirectionalAtoms
87 –  logical, save :: FF_uses_LennardJones
88 –  logical, save :: FF_uses_Electrostatics
89 –  logical, save :: FF_uses_Charges
90    logical, save :: FF_uses_Dipoles
91 –  logical, save :: FF_uses_Quadrupoles
92 –  logical, save :: FF_uses_Sticky
93 –  logical, save :: FF_uses_StickyPower
91    logical, save :: FF_uses_GayBerne
92    logical, save :: FF_uses_EAM
96 –  logical, save :: FF_uses_Shapes
97 –  logical, save :: FF_uses_FLARB
93    logical, save :: FF_uses_RF
94  
95    logical, save :: SIM_uses_DirectionalAtoms
101 –  logical, save :: SIM_uses_LennardJones
102 –  logical, save :: SIM_uses_Electrostatics
103 –  logical, save :: SIM_uses_Charges
104 –  logical, save :: SIM_uses_Dipoles
105 –  logical, save :: SIM_uses_Quadrupoles
106 –  logical, save :: SIM_uses_Sticky
107 –  logical, save :: SIM_uses_StickyPower
108 –  logical, save :: SIM_uses_GayBerne
96    logical, save :: SIM_uses_EAM
110 –  logical, save :: SIM_uses_Shapes
111 –  logical, save :: SIM_uses_FLARB
97    logical, save :: SIM_uses_RF
98    logical, save :: SIM_requires_postpair_calc
99    logical, save :: SIM_requires_prepair_calc
100    logical, save :: SIM_uses_PBC
116 –  logical, save :: SIM_uses_molecular_cutoffs
101  
118 –  !!!GO AWAY---------
119 –  !!!!!real(kind=dp), save :: rlist, rlistsq
120 –
102    public :: init_FF
103 +  public :: setDefaultCutoffs
104    public :: do_force_loop
105 < !  public :: setRlistDF
105 >  public :: createInteractionHash
106 >  public :: createGtypeCutoffMap
107  
108   #ifdef PROFILE
109    public :: getforcetime
# Line 128 | Line 111 | module doForces
111    real :: forceTimeInitial, forceTimeFinal
112    integer :: nLoops
113   #endif
114 +  
115 +  !! Variables for cutoff mapping and interaction mapping
116 +  ! Bit hash to determine pair-pair interactions.
117 +  integer, dimension(:,:), allocatable :: InteractionHash
118 +  real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
119 +  real(kind=dp), dimension(:), allocatable :: groupMaxCutoff
120 +  integer, dimension(:), allocatable :: groupToGtype
121 +  real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff
122 +  type ::gtypeCutoffs
123 +     real(kind=dp) :: rcut
124 +     real(kind=dp) :: rcutsq
125 +     real(kind=dp) :: rlistsq
126 +  end type gtypeCutoffs
127 +  type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
128  
129 <  type, public :: Interaction
130 <     integer :: InteractionHash
134 <     real(kind=dp) :: rCut
135 <  end type Interaction
129 >  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
130 >  real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
131    
137 –  type(Interaction), dimension(:,:),allocatable :: InteractionMap
138 –  
139 –  !public :: addInteraction
140 –  !public :: setInteractionHash
141 –  !public :: getInteractionHash
142 –  public :: createInteractionMap
143 –  
132   contains
133  
134 <
147 <  subroutine createInteractionMap(status)
134 >  subroutine createInteractionHash(status)
135      integer :: nAtypes
136 <    integer :: status
136 >    integer, intent(out) :: status
137      integer :: i
138      integer :: j
139 <    integer :: ihash
140 <    real(kind=dp) :: myRcut
154 < ! Test Types
139 >    integer :: iHash
140 >    !! Test Types
141      logical :: i_is_LJ
142      logical :: i_is_Elect
143      logical :: i_is_Sticky
# Line 167 | Line 153 | contains
153      logical :: j_is_EAM
154      logical :: j_is_Shape
155      
156 <    
156 >    status = 0  
157 >
158      if (.not. associated(atypes)) then
159 <       call handleError("atype", "atypes was not present before call of createDefaultInteractionMap!")
159 >       call handleError("atype", "atypes was not present before call of createInteractionHash!")
160         status = -1
161         return
162      endif
# Line 181 | Line 168 | contains
168         return
169      end if
170  
171 <    if (.not. allocated(InteractionMap)) then
172 <       allocate(InteractionMap(nAtypes,nAtypes))
171 >    if (.not. allocated(InteractionHash)) then
172 >       allocate(InteractionHash(nAtypes,nAtypes))
173      endif
174 +
175 +    if (.not. allocated(atypeMaxCutoff)) then
176 +       allocate(atypeMaxCutoff(nAtypes))
177 +    endif
178          
179      do i = 1, nAtypes
180         call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
# Line 236 | Line 227 | contains
227            if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
228  
229  
230 <          InteractionMap(i,j)%InteractionHash = iHash
231 <          InteractionMap(j,i)%InteractionHash = iHash
230 >          InteractionHash(i,j) = iHash
231 >          InteractionHash(j,i) = iHash
232  
233         end do
234  
235      end do
236 <  end subroutine createInteractionMap
236 >
237 >    haveInteractionHash = .true.
238 >  end subroutine createInteractionHash
239  
240 +  subroutine createGtypeCutoffMap()
241  
242 +    logical :: i_is_LJ
243 +    logical :: i_is_Elect
244 +    logical :: i_is_Sticky
245 +    logical :: i_is_StickyP
246 +    logical :: i_is_GB
247 +    logical :: i_is_EAM
248 +    logical :: i_is_Shape
249  
250 < !!! 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
250 >    integer :: myStatus, nAtypes
251  
252 +    stat = 0
253 +    if (.not. haveInteractionHash) then
254 +       call createInteractionHash(myStatus)      
255 +       if (myStatus .ne. 0) then
256 +          write(default_error, *) 'createInteractionHash failed in doForces!'
257 +          stat = -1
258 +          return
259 +       endif
260 +    endif
261  
262 +    nAtypes = getSize(atypes)
263 +
264 +    do i = 1, nAtypes
265 +       call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
266 +       call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
267 +       call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
268 +       call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
269 +       call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
270 +       call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
271 +       call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
272 +      
273 +       if (i_is_LJ) then
274 +          thisCut = getSigma(i) * DEFAULT_SIGMA_MULTIPLIER
275 +          if (thisCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisCut
276 +       endif
277 +       if (i_is_Elect) then
278 +          thisCut =
279 +    
280 +
281 +
282 +
283 +     haveGtypeCutoffMap = .true.
284 +   end subroutine createGtypeCutoffMap
285 +
286 +   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
287 +     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
288 +     integer, intent(in) :: cutPolicy
289 +
290 +     defaultRcut = defRcut
291 +     defaultRsw = defRsw
292 +     defaultRlist = defRlist
293 +     cutoffPolicy = cutPolicy
294 +   end subroutine setDefaultCutoffs
295 +
296 +   subroutine setCutoffPolicy(cutPolicy)
297 +
298 +     integer, intent(in) :: cutPolicy
299 +     cutoffPolicy = cutPolicy
300 +     call createGtypeCutoffMap()
301 +
302 +   end subroutine setDefaultCutoffs
303 +    
304 +    
305    subroutine setSimVariables()
306      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
264 –    SIM_uses_LennardJones = SimUsesLennardJones()
265 –    SIM_uses_Electrostatics = SimUsesElectrostatics()
266 –    SIM_uses_Charges = SimUsesCharges()
267 –    SIM_uses_Dipoles = SimUsesDipoles()
268 –    SIM_uses_Sticky = SimUsesSticky()
269 –    SIM_uses_StickyPower = SimUsesStickyPower()
270 –    SIM_uses_GayBerne = SimUsesGayBerne()
307      SIM_uses_EAM = SimUsesEAM()
272 –    SIM_uses_Shapes = SimUsesShapes()
273 –    SIM_uses_FLARB = SimUsesFLARB()
308      SIM_uses_RF = SimUsesRF()
309      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
310      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
# Line 288 | Line 322 | contains
322  
323      error = 0
324  
325 <    if (.not. haveInteractionMap) then
325 >    if (.not. haveInteractionHash) then      
326 >       myStatus = 0      
327 >       call createInteractionHash(myStatus)      
328 >       if (myStatus .ne. 0) then
329 >          write(default_error, *) 'createInteractionHash failed in doForces!'
330 >          error = -1
331 >          return
332 >       endif
333 >    endif
334  
335 <       myStatus = 0
336 <
337 <       call createInteractionMap(myStatus)
296 <
335 >    if (.not. haveGtypeCutoffMap) then        
336 >       myStatus = 0      
337 >       call createGtypeCutoffMap(myStatus)      
338         if (myStatus .ne. 0) then
339 <          write(default_error, *) 'createInteractionMap failed in doForces!'
339 >          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
340            error = -1
341            return
342         endif
# Line 356 | Line 397 | contains
397      !! interactions are used by the force field.    
398  
399      FF_uses_DirectionalAtoms = .false.
359 –    FF_uses_LennardJones = .false.
360 –    FF_uses_Electrostatics = .false.
361 –    FF_uses_Charges = .false.    
400      FF_uses_Dipoles = .false.
363 –    FF_uses_Sticky = .false.
364 –    FF_uses_StickyPower = .false.
401      FF_uses_GayBerne = .false.
402      FF_uses_EAM = .false.
367 –    FF_uses_Shapes = .false.
368 –    FF_uses_FLARB = .false.
403  
404      call getMatchingElementList(atypes, "is_Directional", .true., &
405           nMatches, MatchList)
406      if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
407  
374 –    call getMatchingElementList(atypes, "is_LennardJones", .true., &
375 –         nMatches, MatchList)
376 –    if (nMatches .gt. 0) FF_uses_LennardJones = .true.
377 –
378 –    call getMatchingElementList(atypes, "is_Electrostatic", .true., &
379 –         nMatches, MatchList)
380 –    if (nMatches .gt. 0) then
381 –       FF_uses_Electrostatics = .true.
382 –    endif
383 –
384 –    call getMatchingElementList(atypes, "is_Charge", .true., &
385 –         nMatches, MatchList)
386 –    if (nMatches .gt. 0) then
387 –       FF_uses_Charges = .true.  
388 –       FF_uses_Electrostatics = .true.
389 –    endif
390 –
408      call getMatchingElementList(atypes, "is_Dipole", .true., &
409           nMatches, MatchList)
410 <    if (nMatches .gt. 0) then
394 <       FF_uses_Dipoles = .true.
395 <       FF_uses_Electrostatics = .true.
396 <       FF_uses_DirectionalAtoms = .true.
397 <    endif
398 <
399 <    call getMatchingElementList(atypes, "is_Quadrupole", .true., &
400 <         nMatches, MatchList)
401 <    if (nMatches .gt. 0) then
402 <       FF_uses_Quadrupoles = .true.
403 <       FF_uses_Electrostatics = .true.
404 <       FF_uses_DirectionalAtoms = .true.
405 <    endif
406 <
407 <    call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
408 <         MatchList)
409 <    if (nMatches .gt. 0) then
410 <       FF_uses_Sticky = .true.
411 <       FF_uses_DirectionalAtoms = .true.
412 <    endif
413 <
414 <    call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
415 <         MatchList)
416 <    if (nMatches .gt. 0) then
417 <       FF_uses_StickyPower = .true.
418 <       FF_uses_DirectionalAtoms = .true.
419 <    endif
410 >    if (nMatches .gt. 0) FF_uses_Dipoles = .true.
411      
412      call getMatchingElementList(atypes, "is_GayBerne", .true., &
413           nMatches, MatchList)
414 <    if (nMatches .gt. 0) then
424 <       FF_uses_GayBerne = .true.
425 <       FF_uses_DirectionalAtoms = .true.
426 <    endif
414 >    if (nMatches .gt. 0) FF_uses_GayBerne = .true.
415  
416      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
417      if (nMatches .gt. 0) FF_uses_EAM = .true.
418  
431 –    call getMatchingElementList(atypes, "is_Shape", .true., &
432 –         nMatches, MatchList)
433 –    if (nMatches .gt. 0) then
434 –       FF_uses_Shapes = .true.
435 –       FF_uses_DirectionalAtoms = .true.
436 –    endif
419  
438 –    call getMatchingElementList(atypes, "is_FLARB", .true., &
439 –         nMatches, MatchList)
440 –    if (nMatches .gt. 0) FF_uses_FLARB = .true.
441 –
442 –    !! Assume sanity (for the sake of argument)
420      haveSaneForceField = .true.
421  
422      !! check to make sure the FF_uses_RF setting makes sense
423  
424 <    if (FF_uses_dipoles) then
424 >    if (FF_uses_Dipoles) then
425         if (FF_uses_RF) then
426            dielect = getDielect()
427            call initialize_rf(dielect)
# Line 457 | Line 434 | contains
434            return
435         endif
436      endif
460 –
461 –    !sticky module does not contain check_sticky_FF anymore
462 –    !if (FF_uses_sticky) then
463 –    !   call check_sticky_FF(my_status)
464 –    !   if (my_status /= 0) then
465 –    !      thisStat = -1
466 –    !      haveSaneForceField = .false.
467 –    !      return
468 –    !   end if
469 –    !endif
437  
438      if (FF_uses_EAM) then
439         call init_EAM_FF(my_status)
# Line 487 | Line 454 | contains
454         endif
455      endif
456  
490 –    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
491 –    endif
492 –
457      if (.not. haveNeighborList) then
458         !! Create neighbor lists
459         call expandNeighborList(nLocal, my_status)
# Line 553 | Line 517 | contains
517      integer :: localError
518      integer :: propPack_i, propPack_j
519      integer :: loopStart, loopEnd, loop
520 <
520 >    integer :: iHash
521      real(kind=dp) :: listSkin = 1.0  
522  
523      !! initialize local variables  
# Line 645 | Line 609 | contains
609   #endif
610         outer: do i = istart, iend
611  
612 + #ifdef IS_MPI
613 +             me_i = atid_row(i)
614 + #else
615 +             me_i = atid(i)
616 + #endif
617 +
618            if (update_nlist) point(i) = nlist + 1
619  
620            n_in_i = groupStartRow(i+1) - groupStartRow(i)
# Line 672 | Line 642 | contains
642               endif
643  
644   #ifdef IS_MPI
645 +             me_j = atid_col(j)
646               call get_interatomic_vector(q_group_Row(:,i), &
647                    q_group_Col(:,j), d_grp, rgrpsq)
648   #else
649 +             me_j = atid(j)
650               call get_interatomic_vector(q_group(:,i), &
651                    q_group(:,j), d_grp, rgrpsq)
652   #endif
653  
654 <             if (rgrpsq < rlistsq) then
654 >             if (rgrpsq < InteractionHash(me_i,me_j)%rListsq) then
655                  if (update_nlist) then
656                     nlist = nlist + 1
657  
# Line 897 | Line 869 | contains
869   #else
870               me_i = atid(i)
871   #endif
872 <             iMap = InteractionMap(me_i, me_j)%InteractionHash
872 >             iHash = InteractionHash(me_i,me_j)
873              
874 <             if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
874 >             if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
875  
876                  mu_i = getDipoleMoment(me_i)
877  
# Line 966 | Line 938 | contains
938      real ( kind = dp ) :: ebalance
939      integer :: me_i, me_j
940  
941 <    integer :: iMap
941 >    integer :: iHash
942  
943      r = sqrt(rijsq)
944      vpair = 0.0d0
# Line 980 | Line 952 | contains
952      me_j = atid(j)
953   #endif
954  
955 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
955 >    iHash = InteractionHash(me_i, me_j)
956  
957 <    if ( iand(iMap, LJ_PAIR).ne.0 ) then
957 >    if ( iand(iHash, LJ_PAIR).ne.0 ) then
958         call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
959      endif
960  
961 <    if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
961 >    if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
962         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
963              pot, eFrame, f, t, do_pot)
964  
# Line 999 | Line 971 | contains
971  
972      endif
973  
974 <    if ( iand(iMap, STICKY_PAIR).ne.0 ) then
974 >    if ( iand(iHash, STICKY_PAIR).ne.0 ) then
975         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
976              pot, A, f, t, do_pot)
977      endif
978  
979 <    if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
979 >    if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
980         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
981              pot, A, f, t, do_pot)
982      endif
983  
984 <    if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
984 >    if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
985         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
986              pot, A, f, t, do_pot)
987      endif
988      
989 <    if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
990 <       call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
991 <            pot, A, f, t, do_pot)
989 >    if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
990 > !      call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
991 > !           pot, A, f, t, do_pot)
992      endif
993  
994 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
994 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
995         call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
996              do_pot)
997      endif
998  
999 <    if ( iand(iMap, SHAPE_PAIR).ne.0 ) then      
999 >    if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1000         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1001              pot, A, f, t, do_pot)
1002      endif
1003  
1004 <    if ( iand(iMap, SHAPE_LJ).ne.0 ) then      
1004 >    if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1005         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1006              pot, A, f, t, do_pot)
1007      endif
# Line 1051 | Line 1023 | contains
1023      real ( kind = dp )                :: r, rc
1024      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1025  
1026 <    integer :: me_i, me_j, iMap
1026 >    integer :: me_i, me_j, iHash
1027  
1028   #ifdef IS_MPI  
1029      me_i = atid_row(i)
# Line 1061 | Line 1033 | contains
1033      me_j = atid(j)  
1034   #endif
1035  
1036 <    iMap = InteractionMap(me_i, me_j)%InteractionHash
1036 >    iHash = InteractionHash(me_i, me_j)
1037  
1038 <    if ( iand(iMap, EAM_PAIR).ne.0 ) then      
1038 >    if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1039              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1040      endif
1041      
# Line 1260 | Line 1232 | contains
1232  
1233    function FF_UsesDirectionalAtoms() result(doesit)
1234      logical :: doesit
1235 <    doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1264 <         FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1265 <         FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1235 >    doesit = FF_uses_DirectionalAtoms
1236    end function FF_UsesDirectionalAtoms
1237  
1238    function FF_RequiresPrepairCalc() result(doesit)

Diff Legend

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