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 2301 by gezelter, Thu Sep 15 22:05:21 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.44 2005-09-15 22:05:17 gezelter Exp $, $Date: 2005-09-15 22:05:17 $, $Name: not supported by cvs2svn $, $Revision: 1.44 $
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 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 :: rcuti
139    
140   contains
141  
# Line 291 | Line 290 | contains
290            call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
291            
292  
293 <          if (i_is_LJ) then
294 <             thisRcut = getSigma(i) * 2.5_dp
295 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
296 <          endif
297 <          if (i_is_Elect) then
298 <             thisRcut = defaultRcut
299 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
300 <          endif
301 <          if (i_is_Sticky) then
302 <             thisRcut = getStickyCut(i)
303 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
304 <          endif
305 <          if (i_is_StickyP) then
306 <             thisRcut = getStickyPowerCut(i)
307 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
308 <          endif
309 <          if (i_is_GB) then
310 <             thisRcut = getGayBerneCut(i)
311 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
312 <          endif
313 <          if (i_is_EAM) then
314 <             thisRcut = getEAMCut(i)
315 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
316 <          endif
317 <          if (i_is_Shape) then
318 <             thisRcut = getShapeCut(i)
319 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
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    
# Line 408 | Line 413 | contains
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 <    
425 >
426      haveGtypeCutoffMap = .true.
427     end subroutine createGtypeCutoffMap
428  
# Line 422 | Line 434 | contains
434       defaultRsw = defRsw
435       defaultRlist = defRlist
436       cutoffPolicy = cutPolicy
437 <     rcuti = 1.0_dp / defaultRcut
437 >
438 >     haveDefaultCutoffs = .true.
439     end subroutine setDefaultCutoffs
440  
441     subroutine setCutoffPolicy(cutPolicy)
# Line 439 | Line 452 | contains
452      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
453      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
454      SIM_uses_PBC = SimUsesPBC()
442 –    SIM_uses_RF = SimUsesRF()
455  
456      haveSIMvariables = .true.
457  
# Line 506 | Line 518 | contains
518    end subroutine doReadyCheck
519  
520  
521 <  subroutine init_FF(use_RF, correctionMethod, dampingAlpha, thisStat)
521 >  subroutine init_FF(thisESM, thisStat)
522  
523 <    logical, intent(in) :: use_RF
512 <    integer, intent(in) :: correctionMethod
513 <    real(kind=dp), intent(in) :: dampingAlpha
523 >    integer, intent(in) :: thisESM
524      integer, intent(out) :: thisStat  
525      integer :: my_status, nMatches
526      integer, pointer :: MatchList(:) => null()
# Line 519 | Line 529 | contains
529      !! assume things are copacetic, unless they aren't
530      thisStat = 0
531  
532 <    !! Fortran's version of a cast:
523 <    FF_uses_RF = use_RF
532 >    electrostaticSummationMethod = thisESM
533  
525 –        
534      !! init_FF is called *after* all of the atom types have been
535      !! defined in atype_module using the new_atype subroutine.
536      !!
# Line 552 | 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 ((corrMethod == 3) .or. 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 776 | 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 < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
790                  if (update_nlist) then
# Line 978 | Line 986 | contains
986  
987      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
988  
989 <       if ((FF_uses_RF .and. SIM_uses_RF) .or. (corrMethod == 3)) then
989 >       if (electrostaticSummationMethod == REACTION_FIELD) then
990  
991   #ifdef IS_MPI
992            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 1086 | Line 1094 | contains
1094  
1095      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1096         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1097 <            pot, eFrame, f, t, do_pot, corrMethod, rcuti)
1097 >            pot, eFrame, f, t, do_pot)
1098  
1099 <       if ((FF_uses_RF .and. SIM_uses_RF) .or. (corrMethod == 3)) 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 1370 | Line 1378 | contains
1378  
1379    function FF_RequiresPostpairCalc() result(doesit)
1380      logical :: doesit
1381 <    doesit = FF_uses_RF
1374 <    if (corrMethod == 3) doesit = .true.
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