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 2273 by gezelter, Thu Aug 11 21:04:03 2005 UTC vs.
Revision 2323 by chuckv, Fri Sep 23 20:31:02 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
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 $
48 > !! @version $Id: doForces.F90,v 1.49 2005-09-23 20:31:02 chuckv Exp $, $Date: 2005-09-23 20:31:02 $, $Name: not supported by cvs2svn $, $Revision: 1.49 $
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 75 | Line 75 | module doForces
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
# Line 85 | 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
93  logical, save :: FF_uses_RF
96  
97    logical, save :: SIM_uses_DirectionalAtoms
98    logical, save :: SIM_uses_EAM
97  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 128 | Line 136 | module doForces
136  
137    integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
138    real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
139 +  real(kind=dp),save :: listSkin
140    
141   contains
142  
# Line 152 | Line 161 | contains
161      logical :: j_is_GB
162      logical :: j_is_EAM
163      logical :: j_is_Shape
164 <    
164 >    real(kind=dp) :: myRcut
165 >
166      status = 0  
167  
168      if (.not. associated(atypes)) then
# Line 237 | Line 247 | contains
247      haveInteractionHash = .true.
248    end subroutine createInteractionHash
249  
250 <  subroutine createGtypeCutoffMap()
250 >  subroutine createGtypeCutoffMap(stat)
251  
252 +    integer, intent(out), optional :: stat
253      logical :: i_is_LJ
254      logical :: i_is_Elect
255      logical :: i_is_Sticky
# Line 246 | Line 257 | contains
257      logical :: i_is_GB
258      logical :: i_is_EAM
259      logical :: i_is_Shape
260 +    logical :: GtypeFound
261  
262 <    integer :: myStatus, nAtypes
262 >    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
263 >    integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
264 >    integer :: nGroupsInRow
265 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
266 >    real(kind=dp) :: biggestAtypeCutoff
267  
268      stat = 0
269      if (.not. haveInteractionHash) then
# Line 258 | Line 274 | contains
274            return
275         endif
276      endif
277 <
277 > #ifdef IS_MPI
278 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
279 > #endif
280      nAtypes = getSize(atypes)
281 <
281 > ! Set all of the initial cutoffs to zero.
282 >    atypeMaxCutoff = 0.0_dp
283      do i = 1, nAtypes
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 <       if (i_is_LJ) then
293 <          thisCut = getSigma(i) * DEFAULT_SIGMA_MULTIPLIER
294 <          if (thisCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisCut
284 >       if (SimHasAtype(i)) then    
285 >          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
286 >          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
287 >          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
288 >          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
289 >          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
290 >          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
291 >          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
292 >          
293 >
294 >          if (haveDefaultCutoffs) then
295 >             atypeMaxCutoff(i) = defaultRcut
296 >          else
297 >             if (i_is_LJ) then          
298 >                thisRcut = getSigma(i) * 2.5_dp
299 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
300 >             endif
301 >             if (i_is_Elect) then
302 >                thisRcut = defaultRcut
303 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
304 >             endif
305 >             if (i_is_Sticky) then
306 >                thisRcut = getStickyCut(i)
307 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
308 >             endif
309 >             if (i_is_StickyP) then
310 >                thisRcut = getStickyPowerCut(i)
311 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
312 >             endif
313 >             if (i_is_GB) then
314 >                thisRcut = getGayBerneCut(i)
315 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
316 >             endif
317 >             if (i_is_EAM) then
318 >                thisRcut = getEAMCut(i)
319 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
320 >             endif
321 >             if (i_is_Shape) then
322 >                thisRcut = getShapeCut(i)
323 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
324 >             endif
325 >          endif
326 >          
327 >          
328 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
329 >             biggestAtypeCutoff = atypeMaxCutoff(i)
330 >          endif
331 >
332         endif
333 <       if (i_is_Elect) then
334 <          thisCut =
335 <    
333 >    enddo
334 >  
335 >    nGroupTypes = 0
336 >    
337 >    istart = 1
338 > #ifdef IS_MPI
339 >    iend = nGroupsInRow
340 > #else
341 >    iend = nGroups
342 > #endif
343 >    
344 >    !! allocate the groupToGtype and gtypeMaxCutoff here.
345 >    if(.not.allocated(groupToGtype)) then
346 >       allocate(groupToGtype(iend))
347 >       allocate(groupMaxCutoff(iend))
348 >       allocate(gtypeMaxCutoff(iend))
349 >       groupMaxCutoff = 0.0_dp
350 >       gtypeMaxCutoff = 0.0_dp
351 >    endif
352 >    !! first we do a single loop over the cutoff groups to find the
353 >    !! largest cutoff for any atypes present in this group.  We also
354 >    !! create gtypes at this point.
355 >    
356 >    tol = 1.0d-6
357 >    
358 >    do i = istart, iend      
359 >       n_in_i = groupStartRow(i+1) - groupStartRow(i)
360 >       groupMaxCutoff(i) = 0.0_dp
361 >       do ia = groupStartRow(i), groupStartRow(i+1)-1
362 >          atom1 = groupListRow(ia)
363 > #ifdef IS_MPI
364 >          me_i = atid_row(atom1)
365 > #else
366 >          me_i = atid(atom1)
367 > #endif          
368 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
369 >             groupMaxCutoff(i)=atypeMaxCutoff(me_i)
370 >          endif          
371 >       enddo
372  
373 +       if (nGroupTypes.eq.0) then
374 +          nGroupTypes = nGroupTypes + 1
375 +          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
376 +          groupToGtype(i) = nGroupTypes
377 +       else
378 +          GtypeFound = .false.
379 +          do g = 1, nGroupTypes
380 +             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
381 +                groupToGtype(i) = g
382 +                GtypeFound = .true.
383 +             endif
384 +          enddo
385 +          if (.not.GtypeFound) then            
386 +             nGroupTypes = nGroupTypes + 1
387 +             gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
388 +             groupToGtype(i) = nGroupTypes
389 +          endif
390 +       endif
391 +    enddo    
392  
393 +    !! allocate the gtypeCutoffMap here.
394 +    allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
395 +    !! then we do a double loop over all the group TYPES to find the cutoff
396 +    !! map between groups of two types
397 +    
398 +    do i = 1, nGroupTypes
399 +       do j = 1, nGroupTypes
400 +      
401 +          select case(cutoffPolicy)
402 +          case(TRADITIONAL_CUTOFF_POLICY)
403 +             thisRcut = maxval(gtypeMaxCutoff)
404 +          case(MIX_CUTOFF_POLICY)
405 +             thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
406 +          case(MAX_CUTOFF_POLICY)
407 +             thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
408 +          case default
409 +             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
410 +             return
411 +          end select
412 +          gtypeCutoffMap(i,j)%rcut = thisRcut
413 +          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
414 +          skin = defaultRlist - defaultRcut
415 +          listSkin = skin ! set neighbor list skin thickness
416 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
417  
418 <     haveGtypeCutoffMap = .true.
418 >          ! sanity check
419 >
420 >          if (haveDefaultCutoffs) then
421 >             if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
422 >                call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
423 >             endif
424 >          endif
425 >       enddo
426 >    enddo
427 >
428 >    haveGtypeCutoffMap = .true.
429     end subroutine createGtypeCutoffMap
430 <
430 >
431     subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
432       real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
433       integer, intent(in) :: cutPolicy
# Line 291 | Line 436 | contains
436       defaultRsw = defRsw
437       defaultRlist = defRlist
438       cutoffPolicy = cutPolicy
439 +
440 +     haveDefaultCutoffs = .true.
441     end subroutine setDefaultCutoffs
442  
443     subroutine setCutoffPolicy(cutPolicy)
# Line 298 | Line 445 | contains
445       integer, intent(in) :: cutPolicy
446       cutoffPolicy = cutPolicy
447       call createGtypeCutoffMap()
448 <
302 <   end subroutine setDefaultCutoffs
448 >   end subroutine setCutoffPolicy
449      
450      
451    subroutine setSimVariables()
452      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
453      SIM_uses_EAM = SimUsesEAM()
308    SIM_uses_RF = SimUsesRF()
454      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
455      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
456      SIM_uses_PBC = SimUsesPBC()
# Line 346 | Line 491 | contains
491         call setSimVariables()
492      endif
493  
494 <    if (.not. haveRlist) then
495 <       write(default_error, *) 'rList has not been set in doForces!'
496 <       error = -1
497 <       return
498 <    endif
494 >  !  if (.not. haveRlist) then
495 >  !     write(default_error, *) 'rList has not been set in doForces!'
496 >  !     error = -1
497 >  !     return
498 >  !  endif
499  
500      if (.not. haveNeighborList) then
501         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 375 | Line 520 | contains
520    end subroutine doReadyCheck
521  
522  
523 <  subroutine init_FF(use_RF_c, thisStat)
523 >  subroutine init_FF(thisESM, thisStat)
524  
525 <    logical, intent(in) :: use_RF_c
381 <
525 >    integer, intent(in) :: thisESM
526      integer, intent(out) :: thisStat  
527      integer :: my_status, nMatches
528      integer, pointer :: MatchList(:) => null()
# Line 387 | Line 531 | contains
531      !! assume things are copacetic, unless they aren't
532      thisStat = 0
533  
534 <    !! Fortran's version of a cast:
391 <    FF_uses_RF = use_RF_c
534 >    electrostaticSummationMethod = thisESM
535  
536      !! init_FF is called *after* all of the atom types have been
537      !! defined in atype_module using the new_atype subroutine.
# Line 419 | Line 562 | contains
562  
563      haveSaneForceField = .true.
564  
565 <    !! check to make sure the FF_uses_RF setting makes sense
565 >    !! check to make sure the reaction field setting makes sense
566  
567      if (FF_uses_Dipoles) then
568 <       if (FF_uses_RF) then
568 >       if (electrostaticSummationMethod == REACTION_FIELD) then
569            dielect = getDielect()
570            call initialize_rf(dielect)
571         endif
572      else
573 <       if (FF_uses_RF) then          
573 >       if (electrostaticSummationMethod == REACTION_FIELD) then
574            write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
575            thisStat = -1
576            haveSaneForceField = .false.
# Line 518 | Line 661 | contains
661      integer :: propPack_i, propPack_j
662      integer :: loopStart, loopEnd, loop
663      integer :: iHash
664 <    real(kind=dp) :: listSkin = 1.0  
664 >  
665  
666      !! initialize local variables  
667  
# Line 609 | Line 752 | contains
752   #endif
753         outer: do i = istart, iend
754  
612 #ifdef IS_MPI
613             me_i = atid_row(i)
614 #else
615             me_i = atid(i)
616 #endif
617
755            if (update_nlist) point(i) = nlist + 1
756  
757            n_in_i = groupStartRow(i+1) - groupStartRow(i)
# Line 649 | Line 786 | contains
786               me_j = atid(j)
787               call get_interatomic_vector(q_group(:,i), &
788                    q_group(:,j), d_grp, rgrpsq)
789 < #endif
789 > #endif      
790  
791 <             if (rgrpsq < InteractionHash(me_i,me_j)%rListsq) then
791 >             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
792                  if (update_nlist) then
793                     nlist = nlist + 1
794  
# Line 851 | Line 988 | contains
988  
989      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
990  
991 <       if (FF_uses_RF .and. SIM_uses_RF) then
991 >       if (electrostaticSummationMethod == REACTION_FIELD) then
992  
993   #ifdef IS_MPI
994            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 935 | Line 1072 | contains
1072      real ( kind = dp ), intent(inout) :: rijsq
1073      real ( kind = dp )                :: r
1074      real ( kind = dp ), intent(inout) :: d(3)
938    real ( kind = dp ) :: ebalance
1075      integer :: me_i, me_j
1076  
1077      integer :: iHash
# Line 962 | Line 1098 | contains
1098         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1099              pot, eFrame, f, t, do_pot)
1100  
1101 <       if (FF_uses_RF .and. SIM_uses_RF) then
1101 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1102  
1103            ! CHECK ME (RF needs to know about all electrostatic types)
1104            call accumulate_rf(i, j, r, eFrame, sw)
# Line 1024 | Line 1160 | contains
1160      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1161  
1162      integer :: me_i, me_j, iHash
1163 +
1164 +    r = sqrt(rijsq)
1165  
1166   #ifdef IS_MPI  
1167      me_i = atid_row(i)
# Line 1242 | Line 1380 | contains
1380  
1381    function FF_RequiresPostpairCalc() result(doesit)
1382      logical :: doesit
1383 <    doesit = FF_uses_RF
1383 >    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1384    end function FF_RequiresPostpairCalc
1385  
1386   #ifdef PROFILE

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines