ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/UseTheForce/doForces.F90
(Generate patch)

Comparing trunk/OOPSE-2.0/src/UseTheForce/doForces.F90 (file contents):
Revision 2270 by gezelter, Tue Aug 9 22:33:37 2005 UTC vs.
Revision 2317 by chrisfen, Wed Sep 21 17:20:14 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.28 2005-08-09 22:33:37 gezelter Exp $, $Date: 2005-08-09 22:33:37 $, $Name: not supported by cvs2svn $, $Revision: 1.28 $
48 > !! @version $Id: doForces.F90,v 1.48 2005-09-21 17:20:10 chrisfen Exp $, $Date: 2005-09-21 17:20:10 $, $Name: not supported by cvs2svn $, $Revision: 1.48 $
49  
50  
51   module doForces
# Line 58 | Line 58 | module doForces
58    use lj
59    use sticky
60    use electrostatic_module
61 <  use reaction_field
61 >  use reaction_field_module
62    use gb_pair
63    use shapes
64    use vector_class
# 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 + #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
79  
80 +
81    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
82    INTEGER, PARAMETER:: PAIR_LOOP    = 2
83  
# Line 83 | Line 86 | module doForces
86    logical, save :: haveSaneForceField = .false.
87    logical, save :: haveInteractionHash = .false.
88    logical, save :: haveGtypeCutoffMap = .false.
89 +  logical, save :: haveDefaultCutoffs = .false.
90 +  logical, save :: haveRlist = .false.
91  
92    logical, save :: FF_uses_DirectionalAtoms
93    logical, save :: FF_uses_Dipoles
94    logical, save :: FF_uses_GayBerne
95    logical, save :: FF_uses_EAM
91 –  logical, save :: FF_uses_RF
96  
97    logical, save :: SIM_uses_DirectionalAtoms
98    logical, save :: SIM_uses_EAM
95 –  logical, save :: SIM_uses_RF
99    logical, save :: SIM_requires_postpair_calc
100    logical, save :: SIM_requires_prepair_calc
101    logical, save :: SIM_uses_PBC
102  
103 +  integer, save :: electrostaticSummationMethod
104 +
105    public :: init_FF
106 +  public :: setDefaultCutoffs
107    public :: do_force_loop
108    public :: createInteractionHash
109    public :: createGtypeCutoffMap
110 +  public :: getStickyCut
111 +  public :: getStickyPowerCut
112 +  public :: getGayBerneCut
113 +  public :: getEAMCut
114 +  public :: getShapeCut
115  
116   #ifdef PROFILE
117    public :: getforcetime
# Line 122 | Line 133 | module doForces
133       real(kind=dp) :: rlistsq
134    end type gtypeCutoffs
135    type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
136 +
137 +  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
138 +  real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
139    
140   contains
141  
# Line 146 | Line 160 | contains
160      logical :: j_is_GB
161      logical :: j_is_EAM
162      logical :: j_is_Shape
163 <    
163 >    real(kind=dp) :: myRcut
164 >
165      status = 0  
166  
167      if (.not. associated(atypes)) then
# Line 179 | Line 194 | contains
194         call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
195         call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
196  
182 –       if (i_is_LJ) then
183 –          thisCut = getDefaultLJCutoff(i)
184 –          if (thisCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisCut
185 –       endif
186 –
187 –
188 –
197         do j = i, nAtypes
198  
199            iHash = 0
# Line 238 | Line 246 | contains
246      haveInteractionHash = .true.
247    end subroutine createInteractionHash
248  
249 <  subroutine createGtypeCutoffMap(defaultRcut, defaultSkinThickness, stat)
249 >  subroutine createGtypeCutoffMap(stat)
250  
251 <    real(kind=dp), intent(in), optional :: defaultRCut, defaultSkinThickness
252 <    integer, intent(out) :: stat
251 >    integer, intent(out), optional :: stat
252 >    logical :: i_is_LJ
253 >    logical :: i_is_Elect
254 >    logical :: i_is_Sticky
255 >    logical :: i_is_StickyP
256 >    logical :: i_is_GB
257 >    logical :: i_is_EAM
258 >    logical :: i_is_Shape
259 >    logical :: GtypeFound
260  
261 <    integer :: myStatus, nAtypes
261 >    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
262 >    integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
263 >    integer :: nGroupsInRow
264 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
265 >    real(kind=dp) :: biggestAtypeCutoff
266  
267      stat = 0
268      if (.not. haveInteractionHash) then
# Line 254 | Line 273 | contains
273            return
274         endif
275      endif
276 <
276 > #ifdef IS_MPI
277 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
278 > #endif
279      nAtypes = getSize(atypes)
280 <
280 > ! Set all of the initial cutoffs to zero.
281 >    atypeMaxCutoff = 0.0_dp
282      do i = 1, nAtypes
283 <      
284 <       atypeMaxCutoff(i) =
283 >       if (SimHasAtype(i)) then    
284 >          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
285 >          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
286 >          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
287 >          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
288 >          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
289 >          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
290 >          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
291 >          
292 >
293 >          if (haveDefaultCutoffs) then
294 >             atypeMaxCutoff(i) = defaultRcut
295 >          else
296 >             if (i_is_LJ) then          
297 >                thisRcut = getSigma(i) * 2.5_dp
298 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
299 >             endif
300 >             if (i_is_Elect) then
301 >                thisRcut = defaultRcut
302 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
303 >             endif
304 >             if (i_is_Sticky) then
305 >                thisRcut = getStickyCut(i)
306 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
307 >             endif
308 >             if (i_is_StickyP) then
309 >                thisRcut = getStickyPowerCut(i)
310 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
311 >             endif
312 >             if (i_is_GB) then
313 >                thisRcut = getGayBerneCut(i)
314 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
315 >             endif
316 >             if (i_is_EAM) then
317 >                thisRcut = getEAMCut(i)
318 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
319 >             endif
320 >             if (i_is_Shape) then
321 >                thisRcut = getShapeCut(i)
322 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
323 >             endif
324 >          endif
325 >          
326 >          
327 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
328 >             biggestAtypeCutoff = atypeMaxCutoff(i)
329 >          endif
330 >
331 >       endif
332 >    enddo
333 >  
334 >    nGroupTypes = 0
335 >    
336 >    istart = 1
337 > #ifdef IS_MPI
338 >    iend = nGroupsInRow
339 > #else
340 >    iend = nGroups
341 > #endif
342 >    
343 >    !! allocate the groupToGtype and gtypeMaxCutoff here.
344 >    if(.not.allocated(groupToGtype)) then
345 >       allocate(groupToGtype(iend))
346 >       allocate(groupMaxCutoff(iend))
347 >       allocate(gtypeMaxCutoff(iend))
348 >       groupMaxCutoff = 0.0_dp
349 >       gtypeMaxCutoff = 0.0_dp
350 >    endif
351 >    !! first we do a single loop over the cutoff groups to find the
352 >    !! largest cutoff for any atypes present in this group.  We also
353 >    !! create gtypes at this point.
354 >    
355 >    tol = 1.0d-6
356 >    
357 >    do i = istart, iend      
358 >       n_in_i = groupStartRow(i+1) - groupStartRow(i)
359 >       groupMaxCutoff(i) = 0.0_dp
360 >       do ia = groupStartRow(i), groupStartRow(i+1)-1
361 >          atom1 = groupListRow(ia)
362 > #ifdef IS_MPI
363 >          me_i = atid_row(atom1)
364 > #else
365 >          me_i = atid(atom1)
366 > #endif          
367 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
368 >             groupMaxCutoff(i)=atypeMaxCutoff(me_i)
369 >          endif          
370 >       enddo
371  
372 +       if (nGroupTypes.eq.0) then
373 +          nGroupTypes = nGroupTypes + 1
374 +          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
375 +          groupToGtype(i) = nGroupTypes
376 +       else
377 +          GtypeFound = .false.
378 +          do g = 1, nGroupTypes
379 +             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
380 +                groupToGtype(i) = g
381 +                GtypeFound = .true.
382 +             endif
383 +          enddo
384 +          if (.not.GtypeFound) then            
385 +             nGroupTypes = nGroupTypes + 1
386 +             gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
387 +             groupToGtype(i) = nGroupTypes
388 +          endif
389 +       endif
390 +    enddo    
391 +
392 +    !! allocate the gtypeCutoffMap here.
393 +    allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
394 +    !! then we do a double loop over all the group TYPES to find the cutoff
395 +    !! map between groups of two types
396      
397 +    do i = 1, nGroupTypes
398 +       do j = 1, nGroupTypes
399 +      
400 +          select case(cutoffPolicy)
401 +          case(TRADITIONAL_CUTOFF_POLICY)
402 +             thisRcut = maxval(gtypeMaxCutoff)
403 +          case(MIX_CUTOFF_POLICY)
404 +             thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
405 +          case(MAX_CUTOFF_POLICY)
406 +             thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
407 +          case default
408 +             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
409 +             return
410 +          end select
411 +          gtypeCutoffMap(i,j)%rcut = thisRcut
412 +          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
413 +          skin = defaultRlist - defaultRcut
414 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
415  
416 +          ! sanity check
417  
418 +          if (haveDefaultCutoffs) then
419 +             if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
420 +                call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
421 +             endif
422 +          endif
423 +       enddo
424 +    enddo
425  
426 <     haveGtypeCutoffMap = .true.
426 >    haveGtypeCutoffMap = .true.
427     end subroutine createGtypeCutoffMap
428  
429 +   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
430 +     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
431 +     integer, intent(in) :: cutPolicy
432 +
433 +     defaultRcut = defRcut
434 +     defaultRsw = defRsw
435 +     defaultRlist = defRlist
436 +     cutoffPolicy = cutPolicy
437 +
438 +     haveDefaultCutoffs = .true.
439 +   end subroutine setDefaultCutoffs
440 +
441 +   subroutine setCutoffPolicy(cutPolicy)
442 +
443 +     integer, intent(in) :: cutPolicy
444 +     cutoffPolicy = cutPolicy
445 +     call createGtypeCutoffMap()
446 +   end subroutine setCutoffPolicy
447 +    
448 +    
449    subroutine setSimVariables()
450      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
451      SIM_uses_EAM = SimUsesEAM()
274 –    SIM_uses_RF = SimUsesRF()
452      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
453      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
454      SIM_uses_PBC = SimUsesPBC()
# Line 312 | Line 489 | contains
489         call setSimVariables()
490      endif
491  
492 <    if (.not. haveRlist) then
493 <       write(default_error, *) 'rList has not been set in doForces!'
494 <       error = -1
495 <       return
496 <    endif
492 >  !  if (.not. haveRlist) then
493 >  !     write(default_error, *) 'rList has not been set in doForces!'
494 >  !     error = -1
495 >  !     return
496 >  !  endif
497  
498      if (.not. haveNeighborList) then
499         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 341 | Line 518 | contains
518    end subroutine doReadyCheck
519  
520  
521 <  subroutine init_FF(use_RF_c, thisStat)
345 <
346 <    logical, intent(in) :: use_RF_c
521 >  subroutine init_FF(thisESM, thisStat)
522  
523 +    integer, intent(in) :: thisESM
524      integer, intent(out) :: thisStat  
525      integer :: my_status, nMatches
526      integer, pointer :: MatchList(:) => null()
# Line 353 | Line 529 | contains
529      !! assume things are copacetic, unless they aren't
530      thisStat = 0
531  
532 <    !! Fortran's version of a cast:
357 <    FF_uses_RF = use_RF_c
532 >    electrostaticSummationMethod = thisESM
533  
534      !! init_FF is called *after* all of the atom types have been
535      !! defined in atype_module using the new_atype subroutine.
# Line 385 | Line 560 | contains
560  
561      haveSaneForceField = .true.
562  
563 <    !! check to make sure the FF_uses_RF setting makes sense
563 >    !! check to make sure the reaction field setting makes sense
564  
565      if (FF_uses_Dipoles) then
566 <       if (FF_uses_RF) then
566 >       if (electrostaticSummationMethod == REACTION_FIELD) then
567            dielect = getDielect()
568            call initialize_rf(dielect)
569         endif
570      else
571 <       if (FF_uses_RF) then          
571 >       if (electrostaticSummationMethod == REACTION_FIELD) then
572            write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
573            thisStat = -1
574            haveSaneForceField = .false.
# Line 575 | Line 750 | contains
750   #endif
751         outer: do i = istart, iend
752  
578 – #ifdef IS_MPI
579 –             me_i = atid_row(i)
580 – #else
581 –             me_i = atid(i)
582 – #endif
583 –
753            if (update_nlist) point(i) = nlist + 1
754  
755            n_in_i = groupStartRow(i+1) - groupStartRow(i)
# Line 615 | Line 784 | contains
784               me_j = atid(j)
785               call get_interatomic_vector(q_group(:,i), &
786                    q_group(:,j), d_grp, rgrpsq)
787 < #endif
787 > #endif      
788  
789 <             if (rgrpsq < InteractionHash(me_i,me_j)%rListsq) then
789 >             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
790                  if (update_nlist) then
791                     nlist = nlist + 1
792  
# Line 817 | Line 986 | contains
986  
987      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
988  
989 <       if (FF_uses_RF .and. SIM_uses_RF) then
989 >       if (electrostaticSummationMethod == REACTION_FIELD) then
990  
991   #ifdef IS_MPI
992            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 901 | Line 1070 | contains
1070      real ( kind = dp ), intent(inout) :: rijsq
1071      real ( kind = dp )                :: r
1072      real ( kind = dp ), intent(inout) :: d(3)
904 –    real ( kind = dp ) :: ebalance
1073      integer :: me_i, me_j
1074  
1075      integer :: iHash
# Line 928 | Line 1096 | contains
1096         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1097              pot, eFrame, f, t, do_pot)
1098  
1099 <       if (FF_uses_RF .and. SIM_uses_RF) then
1099 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1100  
1101            ! CHECK ME (RF needs to know about all electrostatic types)
1102            call accumulate_rf(i, j, r, eFrame, sw)
# Line 990 | Line 1158 | contains
1158      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1159  
1160      integer :: me_i, me_j, iHash
1161 +
1162 +    r = sqrt(rijsq)
1163  
1164   #ifdef IS_MPI  
1165      me_i = atid_row(i)
# Line 1208 | Line 1378 | contains
1378  
1379    function FF_RequiresPostpairCalc() result(doesit)
1380      logical :: doesit
1381 <    doesit = FF_uses_RF
1381 >    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1382    end function FF_RequiresPostpairCalc
1383  
1384   #ifdef PROFILE

Diff Legend

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