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 2270 by gezelter, Tue Aug 9 22:33:37 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.28 2005-08-09 22:33:37 gezelter Exp $, $Date: 2005-08-09 22:33:37 $, $Name: not supported by cvs2svn $, $Revision: 1.28 $
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 73 | Line 73 | module doForces
73  
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 +
81    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
82    INTEGER, PARAMETER:: PAIR_LOOP    = 2
83  
# Line 83 | 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 97 | 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 122 | Line 134 | module doForces
134       real(kind=dp) :: rlistsq
135    end type gtypeCutoffs
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 146 | 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 179 | Line 197 | contains
197         call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
198         call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
199  
182 –       if (i_is_LJ) then
183 –          thisCut = getDefaultLJCutoff(i)
184 –          if (thisCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisCut
185 –       endif
186 –
187 –
188 –
200         do j = i, nAtypes
201  
202            iHash = 0
# Line 238 | Line 249 | contains
249      haveInteractionHash = .true.
250    end subroutine createInteractionHash
251  
252 <  subroutine createGtypeCutoffMap(defaultRcut, defaultSkinThickness, stat)
252 >  subroutine createGtypeCutoffMap(stat)
253  
254 <    real(kind=dp), intent(in), optional :: defaultRCut, defaultSkinThickness
255 <    integer, intent(out) :: stat
254 >    integer, intent(out), optional :: stat
255 >    logical :: i_is_LJ
256 >    logical :: i_is_Elect
257 >    logical :: i_is_Sticky
258 >    logical :: i_is_StickyP
259 >    logical :: i_is_GB
260 >    logical :: i_is_EAM
261 >    logical :: i_is_Shape
262 >    logical :: GtypeFound
263  
264 <    integer :: myStatus, nAtypes
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 254 | 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 <      
287 <       atypeMaxCutoff(i) =
288 <
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)
290 >          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
291 >          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
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
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 >          
325 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
326 >             biggestAtypeCutoff = atypeMaxCutoff(i)
327 >          endif
328 >       endif
329 >    enddo
330 >  
331 >    nGroupTypes = 0
332      
333 +    istart = 1
334 + #ifdef IS_MPI
335 +    iend = nGroupsInRow
336 + #else
337 +    iend = nGroups
338 + #endif
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 +    
352 +    tol = 1.0d-6
353 +    
354 +    do i = istart, iend      
355 +       n_in_i = groupStartRow(i+1) - groupStartRow(i)
356 +       groupMaxCutoff(i) = 0.0_dp
357 +       do ia = groupStartRow(i), groupStartRow(i+1)-1
358 +          atom1 = groupListRow(ia)
359 + #ifdef IS_MPI
360 +          me_i = atid_row(atom1)
361 + #else
362 +          me_i = atid(atom1)
363 + #endif          
364 +          if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
365 +             groupMaxCutoff(i)=atypeMaxCutoff(me_i)
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)).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  
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 +    
394 +    do i = 1, nGroupTypes
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
412  
413 <     haveGtypeCutoffMap = .true.
413 >       enddo
414 >    enddo
415 >    
416 >    haveGtypeCutoffMap = .true.
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()
435 +   end subroutine setCutoffPolicy
436 +    
437 +    
438    subroutine setSimVariables()
439      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
440      SIM_uses_EAM = SimUsesEAM()
274 –    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 312 | 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 341 | Line 508 | contains
508    end subroutine doReadyCheck
509  
510  
511 <  subroutine init_FF(use_RF_c, thisStat)
511 >  subroutine init_FF(use_RF, correctionMethod, dampingAlpha, thisStat)
512  
513 <    logical, intent(in) :: use_RF_c
514 <
513 >    logical, intent(in) :: use_RF
514 >    integer, intent(in) :: correctionMethod
515 >    real(kind=dp), intent(in) :: dampingAlpha
516      integer, intent(out) :: thisStat  
517      integer :: my_status, nMatches
518      integer, pointer :: MatchList(:) => null()
# Line 354 | Line 522 | contains
522      thisStat = 0
523  
524      !! Fortran's version of a cast:
525 <    FF_uses_RF = use_RF_c
525 >    FF_uses_RF = use_RF
526  
527 +    !! set the electrostatic correction method
528 +    select case(coulombicCorrection)
529 +    case(NONE)
530 +       corrMethod = 0
531 +    case(UNDAMPED_WOLF)
532 +       corrMethod = 1
533 +    case(WOLF)
534 +       corrMethod = 2
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 393 | 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 574 | Line 757 | contains
757         iend = nGroups - 1
758   #endif
759         outer: do i = istart, iend
577 –
578 – #ifdef IS_MPI
579 –             me_i = atid_row(i)
580 – #else
581 –             me_i = atid(i)
582 – #endif
760  
761            if (update_nlist) point(i) = nlist + 1
762  
# Line 617 | Line 794 | contains
794                    q_group(:,j), d_grp, rgrpsq)
795   #endif
796  
797 <             if (rgrpsq < InteractionHash(me_i,me_j)%rListsq) then
797 >             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
798                  if (update_nlist) then
799                     nlist = nlist + 1
800  
# Line 817 | 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 901 | Line 1078 | contains
1078      real ( kind = dp ), intent(inout) :: rijsq
1079      real ( kind = dp )                :: r
1080      real ( kind = dp ), intent(inout) :: d(3)
904 –    real ( kind = dp ) :: ebalance
1081      integer :: me_i, me_j
1082  
1083      integer :: iHash
# Line 926 | 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)
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 991 | 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 1209 | 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