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 2301 by gezelter, Thu Sep 15 22:05:21 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.44 2005-09-15 22:05:17 gezelter Exp $, $Date: 2005-09-15 22:05:17 $, $Name: not supported by cvs2svn $, $Revision: 1.44 $
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/DarkSide/fInteractionMap.h"
78  
79 +
80    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
81    INTEGER, PARAMETER:: PAIR_LOOP    = 2
82  
# Line 83 | Line 85 | module doForces
85    logical, save :: haveSaneForceField = .false.
86    logical, save :: haveInteractionHash = .false.
87    logical, save :: haveGtypeCutoffMap = .false.
88 +  logical, save :: haveRlist = .false.
89  
90    logical, save :: FF_uses_DirectionalAtoms
91    logical, save :: FF_uses_Dipoles
# Line 97 | Line 100 | module doForces
100    logical, save :: SIM_requires_prepair_calc
101    logical, save :: SIM_uses_PBC
102  
103 +  integer, save :: corrMethod
104 +
105    public :: init_FF
106 +  public :: setDefaultCutoffs
107    public :: do_force_loop
108    public :: createInteractionHash
109    public :: createGtypeCutoffMap
110 +  public :: getStickyCut
111 +  public :: getStickyPowerCut
112 +  public :: getGayBerneCut
113 +  public :: getEAMCut
114 +  public :: getShapeCut
115  
116   #ifdef PROFILE
117    public :: getforcetime
# Line 122 | Line 133 | module doForces
133       real(kind=dp) :: rlistsq
134    end type gtypeCutoffs
135    type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
136 +
137 +  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
138 +  real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
139 +  real(kind=dp),save :: rcuti
140    
141   contains
142  
# Line 146 | Line 161 | contains
161      logical :: j_is_GB
162      logical :: j_is_EAM
163      logical :: j_is_Shape
164 <    
164 >    real(kind=dp) :: myRcut
165 >
166      status = 0  
167  
168      if (.not. associated(atypes)) then
# Line 179 | Line 195 | contains
195         call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
196         call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
197  
182 –       if (i_is_LJ) then
183 –          thisCut = getDefaultLJCutoff(i)
184 –          if (thisCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisCut
185 –       endif
186 –
187 –
188 –
198         do j = i, nAtypes
199  
200            iHash = 0
# Line 238 | Line 247 | contains
247      haveInteractionHash = .true.
248    end subroutine createInteractionHash
249  
250 <  subroutine createGtypeCutoffMap(defaultRcut, defaultSkinThickness, stat)
250 >  subroutine createGtypeCutoffMap(stat)
251  
252 <    real(kind=dp), intent(in), optional :: defaultRCut, defaultSkinThickness
253 <    integer, intent(out) :: stat
252 >    integer, intent(out), optional :: stat
253 >    logical :: i_is_LJ
254 >    logical :: i_is_Elect
255 >    logical :: i_is_Sticky
256 >    logical :: i_is_StickyP
257 >    logical :: i_is_GB
258 >    logical :: i_is_EAM
259 >    logical :: i_is_Shape
260 >    logical :: GtypeFound
261  
262 <    integer :: myStatus, nAtypes
262 >    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
263 >    integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
264 >    integer :: nGroupsInRow
265 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
266 >    real(kind=dp) :: biggestAtypeCutoff
267  
268      stat = 0
269      if (.not. haveInteractionHash) then
# Line 254 | Line 274 | contains
274            return
275         endif
276      endif
277 <
277 > #ifdef IS_MPI
278 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
279 > #endif
280      nAtypes = getSize(atypes)
281 <
281 > ! Set all of the initial cutoffs to zero.
282 >    atypeMaxCutoff = 0.0_dp
283      do i = 1, nAtypes
284 <      
285 <       atypeMaxCutoff(i) =
286 <
284 >       if (SimHasAtype(i)) then    
285 >          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
286 >          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
287 >          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
288 >          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
289 >          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
290 >          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
291 >          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
292 >          
293 >
294 >          if (i_is_LJ) then
295 >             thisRcut = getSigma(i) * 2.5_dp
296 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
297 >          endif
298 >          if (i_is_Elect) then
299 >             thisRcut = defaultRcut
300 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
301 >          endif
302 >          if (i_is_Sticky) then
303 >             thisRcut = getStickyCut(i)
304 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
305 >          endif
306 >          if (i_is_StickyP) then
307 >             thisRcut = getStickyPowerCut(i)
308 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
309 >          endif
310 >          if (i_is_GB) then
311 >             thisRcut = getGayBerneCut(i)
312 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
313 >          endif
314 >          if (i_is_EAM) then
315 >             thisRcut = getEAMCut(i)
316 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
317 >          endif
318 >          if (i_is_Shape) then
319 >             thisRcut = getShapeCut(i)
320 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
321 >          endif
322 >          
323 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
324 >             biggestAtypeCutoff = atypeMaxCutoff(i)
325 >          endif
326 >       endif
327 >    enddo
328 >  
329 >    nGroupTypes = 0
330      
331 +    istart = 1
332 + #ifdef IS_MPI
333 +    iend = nGroupsInRow
334 + #else
335 +    iend = nGroups
336 + #endif
337 +    
338 +    !! allocate the groupToGtype and gtypeMaxCutoff here.
339 +    if(.not.allocated(groupToGtype)) then
340 +       allocate(groupToGtype(iend))
341 +       allocate(groupMaxCutoff(iend))
342 +       allocate(gtypeMaxCutoff(iend))
343 +       groupMaxCutoff = 0.0_dp
344 +       gtypeMaxCutoff = 0.0_dp
345 +    endif
346 +    !! first we do a single loop over the cutoff groups to find the
347 +    !! largest cutoff for any atypes present in this group.  We also
348 +    !! create gtypes at this point.
349 +    
350 +    tol = 1.0d-6
351 +    
352 +    do i = istart, iend      
353 +       n_in_i = groupStartRow(i+1) - groupStartRow(i)
354 +       groupMaxCutoff(i) = 0.0_dp
355 +       do ia = groupStartRow(i), groupStartRow(i+1)-1
356 +          atom1 = groupListRow(ia)
357 + #ifdef IS_MPI
358 +          me_i = atid_row(atom1)
359 + #else
360 +          me_i = atid(atom1)
361 + #endif          
362 +          if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
363 +             groupMaxCutoff(i)=atypeMaxCutoff(me_i)
364 +          endif          
365 +       enddo
366  
367 +       if (nGroupTypes.eq.0) then
368 +          nGroupTypes = nGroupTypes + 1
369 +          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
370 +          groupToGtype(i) = nGroupTypes
371 +       else
372 +          GtypeFound = .false.
373 +          do g = 1, nGroupTypes
374 +             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
375 +                groupToGtype(i) = g
376 +                GtypeFound = .true.
377 +             endif
378 +          enddo
379 +          if (.not.GtypeFound) then            
380 +             nGroupTypes = nGroupTypes + 1
381 +             gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
382 +             groupToGtype(i) = nGroupTypes
383 +          endif
384 +       endif
385 +    enddo    
386  
387 +    !! allocate the gtypeCutoffMap here.
388 +    allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
389 +    !! then we do a double loop over all the group TYPES to find the cutoff
390 +    !! map between groups of two types
391 +    
392 +    do i = 1, nGroupTypes
393 +       do j = 1, nGroupTypes
394 +      
395 +          select case(cutoffPolicy)
396 +          case(TRADITIONAL_CUTOFF_POLICY)
397 +             thisRcut = maxval(gtypeMaxCutoff)
398 +          case(MIX_CUTOFF_POLICY)
399 +             thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
400 +          case(MAX_CUTOFF_POLICY)
401 +             thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
402 +          case default
403 +             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
404 +             return
405 +          end select
406 +          gtypeCutoffMap(i,j)%rcut = thisRcut
407 +          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
408 +          skin = defaultRlist - defaultRcut
409 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
410  
411 <     haveGtypeCutoffMap = .true.
411 >       enddo
412 >    enddo
413 >    
414 >    haveGtypeCutoffMap = .true.
415     end subroutine createGtypeCutoffMap
416 +
417 +   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
418 +     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
419 +     integer, intent(in) :: cutPolicy
420 +
421 +     defaultRcut = defRcut
422 +     defaultRsw = defRsw
423 +     defaultRlist = defRlist
424 +     cutoffPolicy = cutPolicy
425 +     rcuti = 1.0_dp / defaultRcut
426 +   end subroutine setDefaultCutoffs
427  
428 +   subroutine setCutoffPolicy(cutPolicy)
429 +
430 +     integer, intent(in) :: cutPolicy
431 +     cutoffPolicy = cutPolicy
432 +     call createGtypeCutoffMap()
433 +   end subroutine setCutoffPolicy
434 +    
435 +    
436    subroutine setSimVariables()
437      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
438      SIM_uses_EAM = SimUsesEAM()
274 –    SIM_uses_RF = SimUsesRF()
439      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
440      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
441      SIM_uses_PBC = SimUsesPBC()
442 +    SIM_uses_RF = SimUsesRF()
443  
444      haveSIMvariables = .true.
445  
# Line 312 | Line 477 | contains
477         call setSimVariables()
478      endif
479  
480 <    if (.not. haveRlist) then
481 <       write(default_error, *) 'rList has not been set in doForces!'
482 <       error = -1
483 <       return
484 <    endif
480 >  !  if (.not. haveRlist) then
481 >  !     write(default_error, *) 'rList has not been set in doForces!'
482 >  !     error = -1
483 >  !     return
484 >  !  endif
485  
486      if (.not. haveNeighborList) then
487         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 341 | Line 506 | contains
506    end subroutine doReadyCheck
507  
508  
509 <  subroutine init_FF(use_RF_c, thisStat)
345 <
346 <    logical, intent(in) :: use_RF_c
509 >  subroutine init_FF(use_RF, correctionMethod, dampingAlpha, thisStat)
510  
511 +    logical, intent(in) :: use_RF
512 +    integer, intent(in) :: correctionMethod
513 +    real(kind=dp), intent(in) :: dampingAlpha
514      integer, intent(out) :: thisStat  
515      integer :: my_status, nMatches
516      integer, pointer :: MatchList(:) => null()
# Line 354 | Line 520 | contains
520      thisStat = 0
521  
522      !! Fortran's version of a cast:
523 <    FF_uses_RF = use_RF_c
523 >    FF_uses_RF = use_RF
524  
525 +        
526      !! init_FF is called *after* all of the atom types have been
527      !! defined in atype_module using the new_atype subroutine.
528      !!
# Line 393 | Line 560 | contains
560            call initialize_rf(dielect)
561         endif
562      else
563 <       if (FF_uses_RF) then          
563 >       if ((corrMethod == 3) .or. FF_uses_RF) then
564            write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
565            thisStat = -1
566            haveSaneForceField = .false.
# Line 574 | Line 741 | contains
741         iend = nGroups - 1
742   #endif
743         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
744  
745            if (update_nlist) point(i) = nlist + 1
746  
# Line 617 | Line 778 | contains
778                    q_group(:,j), d_grp, rgrpsq)
779   #endif
780  
781 <             if (rgrpsq < InteractionHash(me_i,me_j)%rListsq) then
781 >             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
782                  if (update_nlist) then
783                     nlist = nlist + 1
784  
# Line 817 | Line 978 | contains
978  
979      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
980  
981 <       if (FF_uses_RF .and. SIM_uses_RF) then
981 >       if ((FF_uses_RF .and. SIM_uses_RF) .or. (corrMethod == 3)) then
982  
983   #ifdef IS_MPI
984            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 901 | Line 1062 | contains
1062      real ( kind = dp ), intent(inout) :: rijsq
1063      real ( kind = dp )                :: r
1064      real ( kind = dp ), intent(inout) :: d(3)
904 –    real ( kind = dp ) :: ebalance
1065      integer :: me_i, me_j
1066  
1067      integer :: iHash
# Line 926 | Line 1086 | contains
1086  
1087      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1088         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1089 <            pot, eFrame, f, t, do_pot)
1089 >            pot, eFrame, f, t, do_pot, corrMethod, rcuti)
1090  
1091 <       if (FF_uses_RF .and. SIM_uses_RF) then
1091 >       if ((FF_uses_RF .and. SIM_uses_RF) .or. (corrMethod == 3)) then
1092  
1093            ! CHECK ME (RF needs to know about all electrostatic types)
1094            call accumulate_rf(i, j, r, eFrame, sw)
# Line 991 | Line 1151 | contains
1151  
1152      integer :: me_i, me_j, iHash
1153  
1154 +    r = sqrt(rijsq)
1155 +
1156   #ifdef IS_MPI  
1157      me_i = atid_row(i)
1158      me_j = atid_col(j)  
# Line 1209 | Line 1371 | contains
1371    function FF_RequiresPostpairCalc() result(doesit)
1372      logical :: doesit
1373      doesit = FF_uses_RF
1374 +    if (corrMethod == 3) doesit = .true.
1375    end function FF_RequiresPostpairCalc
1376  
1377   #ifdef PROFILE

Diff Legend

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