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 2280 by gezelter, Thu Sep 1 20:17:55 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.34 2005-09-01 20:17:55 gezelter Exp $, $Date: 2005-09-01 20:17:55 $, $Name: not supported by cvs2svn $, $Revision: 1.34 $
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          
286 >       if (SimHasAtype(i)) then    
287            call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
288            call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
289            call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
# 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 +
296            if (i_is_LJ) then
297               thisRcut = getSigma(i) * 2.5_dp
298               if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
# Line 327 | Line 336 | contains
336   #else
337      iend = nGroups
338   #endif
339 <
339 >    
340      !! allocate the groupToGtype and gtypeMaxCutoff here.
341 +    if(.not.allocated(groupToGtype)) then
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
350 +    !! create gtypes at this point.
351      
333 –    !! first we do a single loop over the cutoff groups to find the largest cutoff for any atypes
334 –    !! present in this group.   We also create gtypes at this point.
352      tol = 1.0d-6
353      
354      do i = istart, iend      
# Line 346 | 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
358 <                nGroupTypes = nGroupTypes + 1
359 <                gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
360 <                groupToGtype(i) = nGroupTypes
361 <             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
367 <    
368 <    !! allocate the gtypeCutoffMap here.
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
392      !! map between groups of two types
393      
# Line 374 | Line 395 | contains
395         do j = 1, nGroupTypes
396        
397            select case(cutoffPolicy)
398 <             case(TRADITIONAL_CUTOFF_POLICY)
399 <                thisRcut = maxval(gtypeMaxCutoff)
400 <             case(MIX_CUTOFF_POLICY)
401 <                thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
402 <             case(MAX_CUTOFF_POLICY)
403 <                thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
404 <             case default
405 <                call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
406 <                return
407 <          end select      
408 <         gtypeCutoffMap(i,j)%rcut = thisRcut
409 <         gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
410 <         skin = defaultRlist - defaultRcut
411 <         gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
398 >          case(TRADITIONAL_CUTOFF_POLICY)
399 >             thisRcut = maxval(gtypeMaxCutoff)
400 >          case(MIX_CUTOFF_POLICY)
401 >             thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
402 >          case(MAX_CUTOFF_POLICY)
403 >             thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
404 >          case default
405 >             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
406 >             return
407 >          end select
408 >          gtypeCutoffMap(i,j)%rcut = thisRcut
409 >          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
410 >          skin = defaultRlist - defaultRcut
411 >          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
412 >
413         enddo
414      enddo
415      
# Line 402 | Line 424 | contains
424       defaultRsw = defRsw
425       defaultRlist = defRlist
426       cutoffPolicy = cutPolicy
427 +     rcuti = 1.0_dp / defaultRcut
428     end subroutine setDefaultCutoffs
429  
430     subroutine setCutoffPolicy(cutPolicy)
# Line 409 | Line 432 | contains
432       integer, intent(in) :: cutPolicy
433       cutoffPolicy = cutPolicy
434       call createGtypeCutoffMap()
412 –
435     end subroutine setCutoffPolicy
436      
437      
438    subroutine setSimVariables()
439      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
440      SIM_uses_EAM = SimUsesEAM()
419 –    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 457 | Line 479 | contains
479         call setSimVariables()
480      endif
481  
482 <    if (.not. haveRlist) then
483 <       write(default_error, *) 'rList has not been set in doForces!'
484 <       error = -1
485 <       return
486 <    endif
482 >  !  if (.not. haveRlist) then
483 >  !     write(default_error, *) 'rList has not been set in doForces!'
484 >  !     error = -1
485 >  !     return
486 >  !  endif
487  
488      if (.not. haveNeighborList) then
489         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 486 | 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
496 –    integer :: corrMethod
518      integer, pointer :: MatchList(:) => null()
519      real(kind=dp) :: rcut, rrf, rt, dielect
520  
# Line 504 | 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 549 | 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 967 | 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 1051 | Line 1078 | contains
1078      real ( kind = dp ), intent(inout) :: rijsq
1079      real ( kind = dp )                :: r
1080      real ( kind = dp ), intent(inout) :: d(3)
1054 –    real ( kind = dp ) :: ebalance
1081      integer :: me_i, me_j
1082  
1083      integer :: iHash
# Line 1076 | 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 1141 | 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 1359 | 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