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 2274 by gezelter, Wed Aug 17 15:26:42 2005 UTC vs.
Revision 2295 by chrisfen, Thu Sep 15 00:13:15 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.30 2005-08-17 15:26:37 gezelter Exp $, $Date: 2005-08-17 15:26:37 $, $Name: not supported by cvs2svn $, $Revision: 1.30 $
48 > !! @version $Id: doForces.F90,v 1.42 2005-09-15 00:13:14 chrisfen Exp $, $Date: 2005-09-15 00:13:14 $, $Name: not supported by cvs2svn $, $Revision: 1.42 $
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 85 | Line 86 | module doForces
86    logical, save :: haveSaneForceField = .false.
87    logical, save :: haveInteractionHash = .false.
88    logical, save :: haveGtypeCutoffMap = .false.
89 +  logical, save :: haveRlist = .false.
90  
91    logical, save :: FF_uses_DirectionalAtoms
92    logical, save :: FF_uses_Dipoles
# Line 99 | Line 101 | module doForces
101    logical, save :: SIM_requires_prepair_calc
102    logical, save :: SIM_uses_PBC
103  
104 +  integer, save :: corrMethod
105 +
106    public :: init_FF
107    public :: setDefaultCutoffs
108    public :: do_force_loop
109    public :: createInteractionHash
110    public :: createGtypeCutoffMap
111 +  public :: getStickyCut
112 +  public :: getStickyPowerCut
113 +  public :: getGayBerneCut
114 +  public :: getEAMCut
115 +  public :: getShapeCut
116  
117   #ifdef PROFILE
118    public :: getforcetime
# Line 127 | 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 152 | Line 163 | contains
163      logical :: j_is_GB
164      logical :: j_is_EAM
165      logical :: j_is_Shape
166 <    
166 >    real(kind=dp) :: myRcut
167 >
168      status = 0  
169  
170      if (.not. associated(atypes)) then
# Line 237 | Line 249 | contains
249      haveInteractionHash = .true.
250    end subroutine createInteractionHash
251  
252 <  subroutine createGtypeCutoffMap()
252 >  subroutine createGtypeCutoffMap(stat)
253  
254 +    integer, intent(out), optional :: stat
255      logical :: i_is_LJ
256      logical :: i_is_Elect
257      logical :: i_is_Sticky
# Line 246 | 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
265 <    real(kind=dp):: thisSigma, bigSigma
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  
270      stat = 0
271      if (.not. haveInteractionHash) then
# Line 259 | 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      
284      do i = 1, nAtypes
285 <       if (SimHasAtype(i)) then          
285 >       if (SimHasAtype(i)) then    
286            call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
287            call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
288            call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
# Line 272 | Line 291 | contains
291            call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
292            call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
293            
294 +          atypeMaxCutoff(i) = 0.0_dp
295            if (i_is_LJ) then
296               thisRcut = getSigma(i) * 2.5_dp
297               if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
# Line 284 | Line 304 | contains
304               thisRcut = getStickyCut(i)
305               if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
306            endif
307 <          if (i_is_StickyPower) then
307 >          if (i_is_StickyP) then
308               thisRcut = getStickyPowerCut(i)
309               if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
310            endif
311 <          if (i_is_GayBerne) then
311 >          if (i_is_GB) then
312               thisRcut = getGayBerneCut(i)
313               if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
314            endif
# Line 306 | Line 326 | contains
326            endif
327         endif
328      enddo
329 <
329 >  
330 >    nGroupTypes = 0
331 >    
332      istart = 1
333   #ifdef IS_MPI
334      iend = nGroupsInRow
335   #else
336      iend = nGroups
337   #endif
338 <    outer: do i = istart, iend
339 <
338 >    
339 >    !! allocate the groupToGtype and gtypeMaxCutoff here.
340 >    if(.not.allocated(groupToGtype)) then
341 >       allocate(groupToGtype(iend))
342 >       allocate(groupMaxCutoff(iend))
343 >       allocate(gtypeMaxCutoff(iend))
344 >    endif
345 >    !! first we do a single loop over the cutoff groups to find the
346 >    !! largest cutoff for any atypes present in this group.  We also
347 >    !! create gtypes at this point.
348 >    
349 >    tol = 1.0d-6
350 >    
351 >    do i = istart, iend      
352         n_in_i = groupStartRow(i+1) - groupStartRow(i)
353 <
353 >       groupMaxCutoff(i) = 0.0_dp
354 >       do ia = groupStartRow(i), groupStartRow(i+1)-1
355 >          atom1 = groupListRow(ia)
356   #ifdef IS_MPI
357 <             jstart = 1
322 <             jend = nGroupsInCol
357 >          me_i = atid_row(atom1)
358   #else
359 <             jstart = i+1
360 <             jend = nGroups
361 < #endif
359 >          me_i = atid(atom1)
360 > #endif          
361 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
362 >             groupMaxCutoff(i)=atypeMaxCutoff(me_i)
363 >          endif          
364 >       enddo
365  
366 <
367 <
368 <            
366 >       if (nGroupTypes.eq.0) then
367 >          nGroupTypes = nGroupTypes + 1
368 >          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
369 >          groupToGtype(i) = nGroupTypes
370 >       else
371 >          GtypeFound = .false.
372 >          do g = 1, nGroupTypes
373 >             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
374 >                groupToGtype(i) = g
375 >                GtypeFound = .true.
376 >             endif
377 >          enddo
378 >          if (.not.GtypeFound) then            
379 >             nGroupTypes = nGroupTypes + 1
380 >             gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
381 >             groupToGtype(i) = nGroupTypes
382 >          endif
383 >       endif
384 >    enddo    
385  
386 +    !! allocate the gtypeCutoffMap here.
387 +    allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
388 +    !! then we do a double loop over all the group TYPES to find the cutoff
389 +    !! map between groups of two types
390 +    
391 +    do i = 1, nGroupTypes
392 +       do j = 1, nGroupTypes
393 +      
394 +          select case(cutoffPolicy)
395 +          case(TRADITIONAL_CUTOFF_POLICY)
396 +             thisRcut = maxval(gtypeMaxCutoff)
397 +          case(MIX_CUTOFF_POLICY)
398 +             thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
399 +          case(MAX_CUTOFF_POLICY)
400 +             thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
401 +          case default
402 +             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
403 +             return
404 +          end select
405 +          gtypeCutoffMap(i,j)%rcut = thisRcut
406 +          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
407 +          skin = defaultRlist - defaultRcut
408 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
409  
410 <
411 <     haveGtypeCutoffMap = .true.
410 >       enddo
411 >    enddo
412 >    
413 >    haveGtypeCutoffMap = .true.
414     end subroutine createGtypeCutoffMap
415 <
415 >
416     subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
417       real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
418       integer, intent(in) :: cutPolicy
# Line 342 | Line 421 | contains
421       defaultRsw = defRsw
422       defaultRlist = defRlist
423       cutoffPolicy = cutPolicy
424 +     rcuti = 1.0_dp / defaultRcut
425     end subroutine setDefaultCutoffs
426  
427     subroutine setCutoffPolicy(cutPolicy)
# Line 349 | Line 429 | contains
429       integer, intent(in) :: cutPolicy
430       cutoffPolicy = cutPolicy
431       call createGtypeCutoffMap()
432 <
353 <   end subroutine setDefaultCutoffs
432 >   end subroutine setCutoffPolicy
433      
434      
435    subroutine setSimVariables()
436      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
437      SIM_uses_EAM = SimUsesEAM()
359 –    SIM_uses_RF = SimUsesRF()
438      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
439      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
440      SIM_uses_PBC = SimUsesPBC()
441 +    SIM_uses_RF = SimUsesRF()
442  
443      haveSIMvariables = .true.
444  
# Line 397 | Line 476 | contains
476         call setSimVariables()
477      endif
478  
479 <    if (.not. haveRlist) then
480 <       write(default_error, *) 'rList has not been set in doForces!'
481 <       error = -1
482 <       return
483 <    endif
479 >  !  if (.not. haveRlist) then
480 >  !     write(default_error, *) 'rList has not been set in doForces!'
481 >  !     error = -1
482 >  !     return
483 >  !  endif
484  
485      if (.not. haveNeighborList) then
486         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 426 | Line 505 | contains
505    end subroutine doReadyCheck
506  
507  
508 <  subroutine init_FF(use_RF_c, thisStat)
508 >  subroutine init_FF(use_RF, correctionMethod, dampingAlpha, thisStat)
509  
510 <    logical, intent(in) :: use_RF_c
511 <
510 >    logical, intent(in) :: use_RF
511 >    integer, intent(in) :: correctionMethod
512 >    real(kind=dp), intent(in) :: dampingAlpha
513      integer, intent(out) :: thisStat  
514      integer :: my_status, nMatches
515      integer, pointer :: MatchList(:) => null()
# Line 439 | Line 519 | contains
519      thisStat = 0
520  
521      !! Fortran's version of a cast:
522 <    FF_uses_RF = use_RF_c
522 >    FF_uses_RF = use_RF
523  
524 +    !! set the electrostatic correction method
525 +    select case(coulombicCorrection)
526 +    case(NONE)
527 +       corrMethod = 0
528 +    case(UNDAMPED_WOLF)
529 +       corrMethod = 1
530 +    case(WOLF)
531 +       corrMethod = 2
532 +    case (REACTION_FIELD)
533 +       corrMethod = 3
534 +    case default
535 +       call handleError("init_FF", "Unknown Coulombic Correction Method")
536 +       return
537 +    end select
538 +        
539      !! init_FF is called *after* all of the atom types have been
540      !! defined in atype_module using the new_atype subroutine.
541      !!
# Line 478 | Line 573 | contains
573            call initialize_rf(dielect)
574         endif
575      else
576 <       if (FF_uses_RF) then          
576 >       if ((corrMethod == 3) .or. FF_uses_RF) then
577            write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
578            thisStat = -1
579            haveSaneForceField = .false.
# Line 696 | Line 791 | contains
791                    q_group(:,j), d_grp, rgrpsq)
792   #endif
793  
794 <             if (rgrpsq < InteractionHash(me_i,me_j)%rListsq) then
794 >             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
795                  if (update_nlist) then
796                     nlist = nlist + 1
797  
# Line 896 | Line 991 | contains
991  
992      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
993  
994 <       if (FF_uses_RF .and. SIM_uses_RF) then
994 >       if ((FF_uses_RF .and. SIM_uses_RF) .or. (corrMethod == 3)) then
995  
996   #ifdef IS_MPI
997            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 980 | Line 1075 | contains
1075      real ( kind = dp ), intent(inout) :: rijsq
1076      real ( kind = dp )                :: r
1077      real ( kind = dp ), intent(inout) :: d(3)
983 –    real ( kind = dp ) :: ebalance
1078      integer :: me_i, me_j
1079  
1080      integer :: iHash
# Line 1005 | Line 1099 | contains
1099  
1100      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1101         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1102 <            pot, eFrame, f, t, do_pot)
1102 >            pot, eFrame, f, t, do_pot, corrMethod, rcuti)
1103  
1104 <       if (FF_uses_RF .and. SIM_uses_RF) then
1104 >       if ((FF_uses_RF .and. SIM_uses_RF) .or. (corrMethod == 3)) then
1105  
1106            ! CHECK ME (RF needs to know about all electrostatic types)
1107            call accumulate_rf(i, j, r, eFrame, sw)
# Line 1070 | Line 1164 | contains
1164  
1165      integer :: me_i, me_j, iHash
1166  
1167 +    r = sqrt(rijsq)
1168 +
1169   #ifdef IS_MPI  
1170      me_i = atid_row(i)
1171      me_j = atid_col(j)  
# Line 1288 | Line 1384 | contains
1384    function FF_RequiresPostpairCalc() result(doesit)
1385      logical :: doesit
1386      doesit = FF_uses_RF
1387 +    if (corrMethod == 3) doesit = .true.
1388    end function FF_RequiresPostpairCalc
1389  
1390   #ifdef PROFILE

Diff Legend

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