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 2284 by gezelter, Wed Sep 7 19:18:54 2005 UTC vs.
Revision 2298 by chuckv, Thu Sep 15 02:48:43 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.37 2005-09-07 19:18:54 gezelter Exp $, $Date: 2005-09-07 19:18:54 $, $Name: not supported by cvs2svn $, $Revision: 1.37 $
48 > !! @version $Id: doForces.F90,v 1.43 2005-09-15 02:48:43 chuckv Exp $, $Date: 2005-09-15 02:48:43 $, $Name: not supported by cvs2svn $, $Revision: 1.43 $
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 74 | Line 74 | module doForces
74   #define __FORTRAN90
75   #include "UseTheForce/fSwitchingFunction.h"
76   #include "UseTheForce/fCutoffPolicy.h"
77 + #include "UseTheForce/fCoulombicCorrection.h"
78   #include "UseTheForce/DarkSide/fInteractionMap.h"
79  
80  
# Line 135 | Line 136 | module doForces
136    type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
137  
138    integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
139 +  integer, save :: coulombicCorrection = NONE
140    real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
141 +  real(kind=dp),save :: rcuti
142    
143   contains
144  
# Line 256 | Line 259 | contains
259      logical :: i_is_GB
260      logical :: i_is_EAM
261      logical :: i_is_Shape
262 +    logical :: GtypeFound
263  
264      integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
265      integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
266 +    integer :: nGroupsInRow
267      real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
268      real(kind=dp) :: biggestAtypeCutoff
269  
# Line 271 | Line 276 | contains
276            return
277         endif
278      endif
279 <
279 > #ifdef IS_MPI
280 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
281 > #endif
282      nAtypes = getSize(atypes)
283 <    
283 > ! Set all of the initial cutoffs to zero.
284 >    atypeMaxCutoff = 0.0_dp
285      do i = 1, nAtypes
286         if (SimHasAtype(i)) then    
287            call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
# Line 284 | Line 292 | contains
292            call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
293            call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
294            
295 <          atypeMaxCutoff(i) = 0.0_dp
295 >
296            if (i_is_LJ) then
297               thisRcut = getSigma(i) * 2.5_dp
298               if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
# Line 334 | Line 342 | contains
342         allocate(groupToGtype(iend))
343         allocate(groupMaxCutoff(iend))
344         allocate(gtypeMaxCutoff(iend))
345 +       groupMaxCutoff = 0.0_dp
346 +       gtypeMaxCutoff = 0.0_dp
347      endif
348      !! first we do a single loop over the cutoff groups to find the
349      !! largest cutoff for any atypes present in this group.  We also
# Line 353 | Line 363 | contains
363   #endif          
364            if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
365               groupMaxCutoff(i)=atypeMaxCutoff(me_i)
366 <          endif
366 >          endif          
367         enddo
368 +
369         if (nGroupTypes.eq.0) then
370            nGroupTypes = nGroupTypes + 1
371            gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
372            groupToGtype(i) = nGroupTypes
373         else
374 +          GtypeFound = .false.
375            do g = 1, nGroupTypes
376 <             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).gt.tol) then
365 <                nGroupTypes = nGroupTypes + 1
366 <                gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
367 <                groupToGtype(i) = nGroupTypes
368 <             else
376 >             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
377                  groupToGtype(i) = g
378 +                GtypeFound = .true.
379               endif
380            enddo
381 +          if (.not.GtypeFound) then            
382 +             nGroupTypes = nGroupTypes + 1
383 +             gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
384 +             groupToGtype(i) = nGroupTypes
385 +          endif
386         endif
387 <    enddo
388 <    
387 >    enddo    
388 >
389      !! allocate the gtypeCutoffMap here.
390      allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
391      !! then we do a double loop over all the group TYPES to find the cutoff
# Line 400 | Line 414 | contains
414      enddo
415      
416      haveGtypeCutoffMap = .true.
417 <    
418 <  end subroutine createGtypeCutoffMap
419 <  
420 <  subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
421 <    real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
408 <    integer, intent(in) :: cutPolicy
409 <    
410 <    defaultRcut = defRcut
411 <    defaultRsw = defRsw
412 <    defaultRlist = defRlist
413 <    cutoffPolicy = cutPolicy
414 <  end subroutine setDefaultCutoffs
415 <  
416 <  subroutine setCutoffPolicy(cutPolicy)
417 >   end subroutine createGtypeCutoffMap
418 >
419 >   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
420 >     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
421 >     integer, intent(in) :: cutPolicy
422  
423 +     defaultRcut = defRcut
424 +     defaultRsw = defRsw
425 +     defaultRlist = defRlist
426 +     cutoffPolicy = cutPolicy
427 +     rcuti = 1.0_dp / defaultRcut
428 +   end subroutine setDefaultCutoffs
429 +
430 +   subroutine setCutoffPolicy(cutPolicy)
431 +
432       integer, intent(in) :: cutPolicy
433       cutoffPolicy = cutPolicy
434       call createGtypeCutoffMap()
421 –
435     end subroutine setCutoffPolicy
436      
437      
438    subroutine setSimVariables()
439      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
440      SIM_uses_EAM = SimUsesEAM()
428 –    SIM_uses_RF = SimUsesRF()
441      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
442      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
443      SIM_uses_PBC = SimUsesPBC()
444 +    SIM_uses_RF = SimUsesRF()
445  
446      haveSIMvariables = .true.
447  
# Line 495 | Line 508 | contains
508    end subroutine doReadyCheck
509  
510  
511 <  subroutine init_FF(use_RF, use_UW, use_DW, thisStat)
511 >  subroutine init_FF(use_RF, correctionMethod, dampingAlpha, thisStat)
512  
513      logical, intent(in) :: use_RF
514 <    logical, intent(in) :: use_UW
515 <    logical, intent(in) :: use_DW
514 >    integer, intent(in) :: correctionMethod
515 >    real(kind=dp), intent(in) :: dampingAlpha
516      integer, intent(out) :: thisStat  
517      integer :: my_status, nMatches
505 –    integer :: corrMethod
518      integer, pointer :: MatchList(:) => null()
519      real(kind=dp) :: rcut, rrf, rt, dielect
520  
# Line 513 | Line 525 | contains
525      FF_uses_RF = use_RF
526  
527      !! set the electrostatic correction method
528 <    if (use_UW) then
528 >    select case(coulombicCorrection)
529 >    case(NONE)
530 >       corrMethod = 0
531 >    case(UNDAMPED_WOLF)
532         corrMethod = 1
533 <    elseif (use_DW) then
533 >    case(WOLF)
534         corrMethod = 2
535 <    else
536 <       corrMethod = 0
537 <    endif
538 <    
535 >    case (REACTION_FIELD)
536 >       corrMethod = 3
537 >    case default
538 >       call handleError("init_FF", "Unknown Coulombic Correction Method")
539 >       return
540 >    end select
541 >        
542      !! init_FF is called *after* all of the atom types have been
543      !! defined in atype_module using the new_atype subroutine.
544      !!
# Line 558 | Line 576 | contains
576            call initialize_rf(dielect)
577         endif
578      else
579 <       if (FF_uses_RF) then          
579 >       if ((corrMethod == 3) .or. FF_uses_RF) then
580            write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
581            thisStat = -1
582            haveSaneForceField = .false.
# Line 976 | Line 994 | contains
994  
995      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
996  
997 <       if (FF_uses_RF .and. SIM_uses_RF) then
997 >       if ((FF_uses_RF .and. SIM_uses_RF) .or. (corrMethod == 3)) then
998  
999   #ifdef IS_MPI
1000            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 1060 | Line 1078 | contains
1078      real ( kind = dp ), intent(inout) :: rijsq
1079      real ( kind = dp )                :: r
1080      real ( kind = dp ), intent(inout) :: d(3)
1063 –    real ( kind = dp ) :: ebalance
1081      integer :: me_i, me_j
1082  
1083      integer :: iHash
# Line 1085 | Line 1102 | contains
1102  
1103      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1104         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1105 <            pot, eFrame, f, t, do_pot, corrMethod)
1105 >            pot, eFrame, f, t, do_pot, corrMethod, rcuti)
1106  
1107 <       if (FF_uses_RF .and. SIM_uses_RF) then
1107 >       if ((FF_uses_RF .and. SIM_uses_RF) .or. (corrMethod == 3)) then
1108  
1109            ! CHECK ME (RF needs to know about all electrostatic types)
1110            call accumulate_rf(i, j, r, eFrame, sw)
# Line 1150 | Line 1167 | contains
1167  
1168      integer :: me_i, me_j, iHash
1169  
1170 +    r = sqrt(rijsq)
1171 +
1172   #ifdef IS_MPI  
1173      me_i = atid_row(i)
1174      me_j = atid_col(j)  
# Line 1368 | Line 1387 | contains
1387    function FF_RequiresPostpairCalc() result(doesit)
1388      logical :: doesit
1389      doesit = FF_uses_RF
1390 +    if (corrMethod == 3) doesit = .true.
1391    end function FF_RequiresPostpairCalc
1392  
1393   #ifdef PROFILE

Diff Legend

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