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 2288 by chuckv, Wed Sep 7 22:44:48 2005 UTC vs.
Revision 2342 by chrisfen, Tue Oct 4 19:32:58 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.40 2005-09-07 22:44:48 chuckv Exp $, $Date: 2005-09-07 22:44:48 $, $Name: not supported by cvs2svn $, $Revision: 1.40 $
48 > !! @version $Id: doForces.F90,v 1.50 2005-10-04 19:32:58 chrisfen Exp $, $Date: 2005-10-04 19:32:58 $, $Name: not supported by cvs2svn $, $Revision: 1.50 $
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
94 –  logical, save :: FF_uses_RF
96  
97    logical, save :: SIM_uses_DirectionalAtoms
98    logical, save :: SIM_uses_EAM
98 –  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 :: corrMethod
103 >  integer, save :: electrostaticSummationMethod
104  
105    public :: init_FF
106    public :: setDefaultCutoffs
# Line 136 | 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 277 | Line 278 | contains
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         if (SimHasAtype(i)) then    
285            call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
# Line 288 | Line 290 | contains
290            call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
291            call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
292            
293 <          atypeMaxCutoff(i) = 0.0_dp
294 <          if (i_is_LJ) then
295 <             thisRcut = getSigma(i) * 2.5_dp
296 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
297 <          endif
298 <          if (i_is_Elect) then
299 <             thisRcut = defaultRcut
300 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
301 <          endif
302 <          if (i_is_Sticky) then
303 <             thisRcut = getStickyCut(i)
304 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
305 <          endif
306 <          if (i_is_StickyP) then
307 <             thisRcut = getStickyPowerCut(i)
308 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
309 <          endif
310 <          if (i_is_GB) then
311 <             thisRcut = getGayBerneCut(i)
312 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
313 <          endif
314 <          if (i_is_EAM) then
315 <             thisRcut = getEAMCut(i)
316 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
317 <          endif
318 <          if (i_is_Shape) then
319 <             thisRcut = getShapeCut(i)
320 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
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      enddo
334    
# Line 338 | Line 346 | contains
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
# Line 402 | Line 412 | contains
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 +          ! 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 <    
427 >
428      haveGtypeCutoffMap = .true.
429 <    
412 <  end subroutine createGtypeCutoffMap
413 <  
414 <  subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
415 <    real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
416 <    integer, intent(in) :: cutPolicy
417 <    
418 <    defaultRcut = defRcut
419 <    defaultRsw = defRsw
420 <    defaultRlist = defRlist
421 <    cutoffPolicy = cutPolicy
422 <  end subroutine setDefaultCutoffs
423 <  
424 <  subroutine setCutoffPolicy(cutPolicy)
429 >   end subroutine createGtypeCutoffMap
430  
431 +   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
432 +     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
433       integer, intent(in) :: cutPolicy
434 +
435 +     defaultRcut = defRcut
436 +     defaultRsw = defRsw
437 +     defaultRlist = defRlist
438       cutoffPolicy = cutPolicy
428 –     call createGtypeCutoffMap()
439  
440 +     haveDefaultCutoffs = .true.
441 +   end subroutine setDefaultCutoffs
442 +
443 +   subroutine setCutoffPolicy(cutPolicy)
444 +
445 +     integer, intent(in) :: cutPolicy
446 +     cutoffPolicy = cutPolicy
447 +     call createGtypeCutoffMap()
448     end subroutine setCutoffPolicy
449      
450      
451    subroutine setSimVariables()
452      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
453      SIM_uses_EAM = SimUsesEAM()
436 –    SIM_uses_RF = SimUsesRF()
454      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
455      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
456      SIM_uses_PBC = SimUsesPBC()
# Line 503 | Line 520 | contains
520    end subroutine doReadyCheck
521  
522  
523 <  subroutine init_FF(use_RF, use_UW, use_DW, thisStat)
523 >  subroutine init_FF(thisESM, thisStat)
524  
525 <    logical, intent(in) :: use_RF
509 <    logical, intent(in) :: use_UW
510 <    logical, intent(in) :: use_DW
525 >    integer, intent(in) :: thisESM
526      integer, intent(out) :: thisStat  
527      integer :: my_status, nMatches
513 –    integer :: corrMethod
528      integer, pointer :: MatchList(:) => null()
529      real(kind=dp) :: rcut, rrf, rt, dielect
530  
531      !! assume things are copacetic, unless they aren't
532      thisStat = 0
533  
534 <    !! Fortran's version of a cast:
521 <    FF_uses_RF = use_RF
534 >    electrostaticSummationMethod = thisESM
535  
523 –    !! set the electrostatic correction method
524 –    if (use_UW) then
525 –       corrMethod = 1
526 –    elseif (use_DW) then
527 –       corrMethod = 2
528 –    else
529 –       corrMethod = 0
530 –    endif
531 –    
536      !! init_FF is called *after* all of the atom types have been
537      !! defined in atype_module using the new_atype subroutine.
538      !!
# Line 558 | 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 657 | 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 782 | 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 < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
792                  if (update_nlist) then
# Line 905 | Line 909 | contains
909                  endif
910               end if
911            enddo
912 +
913         enddo outer
914  
915         if (update_nlist) then
# Line 984 | Line 989 | contains
989  
990      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
991  
992 <       if (FF_uses_RF .and. SIM_uses_RF) then
992 >       if (electrostaticSummationMethod == REACTION_FIELD) then
993  
994   #ifdef IS_MPI
995            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 1068 | Line 1073 | contains
1073      real ( kind = dp ), intent(inout) :: rijsq
1074      real ( kind = dp )                :: r
1075      real ( kind = dp ), intent(inout) :: d(3)
1071 –    real ( kind = dp ) :: ebalance
1076      integer :: me_i, me_j
1077  
1078      integer :: iHash
# Line 1093 | Line 1097 | contains
1097  
1098      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1099         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1100 <            pot, eFrame, f, t, do_pot, corrMethod)
1100 >            pot, eFrame, f, t, do_pot)
1101  
1102 <       if (FF_uses_RF .and. SIM_uses_RF) then
1102 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1103  
1104            ! CHECK ME (RF needs to know about all electrostatic types)
1105            call accumulate_rf(i, j, r, eFrame, sw)
# Line 1157 | Line 1161 | contains
1161      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1162  
1163      integer :: me_i, me_j, iHash
1164 +
1165 +    r = sqrt(rijsq)
1166  
1167   #ifdef IS_MPI  
1168      me_i = atid_row(i)
# Line 1375 | Line 1381 | contains
1381  
1382    function FF_RequiresPostpairCalc() result(doesit)
1383      logical :: doesit
1384 <    doesit = FF_uses_RF
1384 >    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1385    end function FF_RequiresPostpairCalc
1386  
1387   #ifdef PROFILE

Diff Legend

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