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 2309 by chrisfen, Sun Sep 18 20:45:38 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.46 2005-09-18 20:45:38 chrisfen Exp $, $Date: 2005-09-18 20:45:38 $, $Name: not supported by cvs2svn $, $Revision: 1.46 $
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
92    logical, save :: FF_uses_GayBerne
93    logical, save :: FF_uses_EAM
91 –  logical, save :: FF_uses_RF
94  
95    logical, save :: SIM_uses_DirectionalAtoms
96    logical, save :: SIM_uses_EAM
95 –  logical, save :: SIM_uses_RF
97    logical, save :: SIM_requires_postpair_calc
98    logical, save :: SIM_requires_prepair_calc
99    logical, save :: SIM_uses_PBC
100  
101 +  integer, save :: electrostaticSummationMethod
102 +
103    public :: init_FF
104 +  public :: setDefaultCutoffs
105    public :: do_force_loop
106    public :: createInteractionHash
107    public :: createGtypeCutoffMap
108 +  public :: getStickyCut
109 +  public :: getStickyPowerCut
110 +  public :: getGayBerneCut
111 +  public :: getEAMCut
112 +  public :: getShapeCut
113  
114   #ifdef PROFILE
115    public :: getforcetime
# Line 122 | Line 131 | module doForces
131       real(kind=dp) :: rlistsq
132    end type gtypeCutoffs
133    type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
134 +
135 +  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
136 +  real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
137    
138   contains
139  
# Line 146 | Line 158 | contains
158      logical :: j_is_GB
159      logical :: j_is_EAM
160      logical :: j_is_Shape
161 <    
161 >    real(kind=dp) :: myRcut
162 >
163      status = 0  
164  
165      if (.not. associated(atypes)) then
# Line 179 | Line 192 | contains
192         call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
193         call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
194  
182 –       if (i_is_LJ) then
183 –          thisCut = getDefaultLJCutoff(i)
184 –          if (thisCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisCut
185 –       endif
186 –
187 –
188 –
195         do j = i, nAtypes
196  
197            iHash = 0
# Line 238 | Line 244 | contains
244      haveInteractionHash = .true.
245    end subroutine createInteractionHash
246  
247 <  subroutine createGtypeCutoffMap(defaultRcut, defaultSkinThickness, stat)
242 <
243 <    real(kind=dp), intent(in), optional :: defaultRCut, defaultSkinThickness
244 <    integer, intent(out) :: stat
247 >  subroutine createGtypeCutoffMap(stat)
248  
249 <    integer :: myStatus, nAtypes
249 >    integer, intent(out), optional :: stat
250 >    logical :: i_is_LJ
251 >    logical :: i_is_Elect
252 >    logical :: i_is_Sticky
253 >    logical :: i_is_StickyP
254 >    logical :: i_is_GB
255 >    logical :: i_is_EAM
256 >    logical :: i_is_Shape
257 >    logical :: GtypeFound
258  
259 +    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
260 +    integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
261 +    integer :: nGroupsInRow
262 +    real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
263 +    real(kind=dp) :: biggestAtypeCutoff
264 +
265      stat = 0
266      if (.not. haveInteractionHash) then
267         call createInteractionHash(myStatus)      
# Line 254 | Line 271 | contains
271            return
272         endif
273      endif
274 <
274 > #ifdef IS_MPI
275 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
276 > #endif
277      nAtypes = getSize(atypes)
278 <
278 > ! Set all of the initial cutoffs to zero.
279 >    atypeMaxCutoff = 0.0_dp
280      do i = 1, nAtypes
281 <      
282 <       atypeMaxCutoff(i) =
283 <
281 >       if (SimHasAtype(i)) then    
282 >          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
283 >          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
284 >          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
285 >          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
286 >          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
287 >          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
288 >          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
289 >          
290 >
291 >          if (i_is_LJ) then
292 >             thisRcut = getSigma(i) * 2.5_dp
293 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
294 >          endif
295 >          if (i_is_Elect) then
296 >             thisRcut = defaultRcut
297 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
298 >          endif
299 >          if (i_is_Sticky) then
300 >             thisRcut = getStickyCut(i)
301 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
302 >          endif
303 >          if (i_is_StickyP) then
304 >             thisRcut = getStickyPowerCut(i)
305 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
306 >          endif
307 >          if (i_is_GB) then
308 >             thisRcut = getGayBerneCut(i)
309 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
310 >          endif
311 >          if (i_is_EAM) then
312 >             thisRcut = getEAMCut(i)
313 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
314 >          endif
315 >          if (i_is_Shape) then
316 >             thisRcut = getShapeCut(i)
317 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
318 >          endif
319 >          
320 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
321 >             biggestAtypeCutoff = atypeMaxCutoff(i)
322 >          endif
323 >       endif
324 >    enddo
325 >  
326 >    nGroupTypes = 0
327 >    
328 >    istart = 1
329 > #ifdef IS_MPI
330 >    iend = nGroupsInRow
331 > #else
332 >    iend = nGroups
333 > #endif
334 >    
335 >    !! allocate the groupToGtype and gtypeMaxCutoff here.
336 >    if(.not.allocated(groupToGtype)) then
337 >       allocate(groupToGtype(iend))
338 >       allocate(groupMaxCutoff(iend))
339 >       allocate(gtypeMaxCutoff(iend))
340 >       groupMaxCutoff = 0.0_dp
341 >       gtypeMaxCutoff = 0.0_dp
342 >    endif
343 >    !! first we do a single loop over the cutoff groups to find the
344 >    !! largest cutoff for any atypes present in this group.  We also
345 >    !! create gtypes at this point.
346 >    
347 >    tol = 1.0d-6
348      
349 +    do i = istart, iend      
350 +       n_in_i = groupStartRow(i+1) - groupStartRow(i)
351 +       groupMaxCutoff(i) = 0.0_dp
352 +       do ia = groupStartRow(i), groupStartRow(i+1)-1
353 +          atom1 = groupListRow(ia)
354 + #ifdef IS_MPI
355 +          me_i = atid_row(atom1)
356 + #else
357 +          me_i = atid(atom1)
358 + #endif          
359 +          if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
360 +             groupMaxCutoff(i)=atypeMaxCutoff(me_i)
361 +          endif          
362 +       enddo
363  
364 +       if (nGroupTypes.eq.0) then
365 +          nGroupTypes = nGroupTypes + 1
366 +          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
367 +          groupToGtype(i) = nGroupTypes
368 +       else
369 +          GtypeFound = .false.
370 +          do g = 1, nGroupTypes
371 +             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
372 +                groupToGtype(i) = g
373 +                GtypeFound = .true.
374 +             endif
375 +          enddo
376 +          if (.not.GtypeFound) then            
377 +             nGroupTypes = nGroupTypes + 1
378 +             gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
379 +             groupToGtype(i) = nGroupTypes
380 +          endif
381 +       endif
382 +    enddo    
383  
384 +    !! allocate the gtypeCutoffMap here.
385 +    allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
386 +    !! then we do a double loop over all the group TYPES to find the cutoff
387 +    !! map between groups of two types
388 +    
389 +    do i = 1, nGroupTypes
390 +       do j = 1, nGroupTypes
391 +      
392 +          select case(cutoffPolicy)
393 +          case(TRADITIONAL_CUTOFF_POLICY)
394 +             thisRcut = maxval(gtypeMaxCutoff)
395 +          case(MIX_CUTOFF_POLICY)
396 +             thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
397 +          case(MAX_CUTOFF_POLICY)
398 +             thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
399 +          case default
400 +             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
401 +             return
402 +          end select
403 +          gtypeCutoffMap(i,j)%rcut = thisRcut
404 +          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
405 +          skin = defaultRlist - defaultRcut
406 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
407  
408 <     haveGtypeCutoffMap = .true.
408 >       enddo
409 >    enddo
410 >    
411 >    haveGtypeCutoffMap = .true.
412     end subroutine createGtypeCutoffMap
413 +
414 +   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
415 +     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
416 +     integer, intent(in) :: cutPolicy
417  
418 +     defaultRcut = defRcut
419 +     defaultRsw = defRsw
420 +     defaultRlist = defRlist
421 +     cutoffPolicy = cutPolicy
422 +   end subroutine setDefaultCutoffs
423 +
424 +   subroutine setCutoffPolicy(cutPolicy)
425 +
426 +     integer, intent(in) :: cutPolicy
427 +     cutoffPolicy = cutPolicy
428 +     call createGtypeCutoffMap()
429 +   end subroutine setCutoffPolicy
430 +    
431 +    
432    subroutine setSimVariables()
433      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
434      SIM_uses_EAM = SimUsesEAM()
274 –    SIM_uses_RF = SimUsesRF()
435      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
436      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
437      SIM_uses_PBC = SimUsesPBC()
# Line 312 | Line 472 | contains
472         call setSimVariables()
473      endif
474  
475 <    if (.not. haveRlist) then
476 <       write(default_error, *) 'rList has not been set in doForces!'
477 <       error = -1
478 <       return
479 <    endif
475 >  !  if (.not. haveRlist) then
476 >  !     write(default_error, *) 'rList has not been set in doForces!'
477 >  !     error = -1
478 >  !     return
479 >  !  endif
480  
481      if (.not. haveNeighborList) then
482         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 341 | Line 501 | contains
501    end subroutine doReadyCheck
502  
503  
504 <  subroutine init_FF(use_RF_c, thisStat)
345 <
346 <    logical, intent(in) :: use_RF_c
504 >  subroutine init_FF(thisESM, thisStat)
505  
506 +    integer, intent(in) :: thisESM
507 +    real(kind=dp), intent(in) :: dampingAlpha
508      integer, intent(out) :: thisStat  
509      integer :: my_status, nMatches
510      integer, pointer :: MatchList(:) => null()
# Line 353 | Line 513 | contains
513      !! assume things are copacetic, unless they aren't
514      thisStat = 0
515  
516 <    !! Fortran's version of a cast:
357 <    FF_uses_RF = use_RF_c
516 >    electrostaticSummationMethod = thisESM
517  
518      !! init_FF is called *after* all of the atom types have been
519      !! defined in atype_module using the new_atype subroutine.
# Line 385 | Line 544 | contains
544  
545      haveSaneForceField = .true.
546  
547 <    !! check to make sure the FF_uses_RF setting makes sense
547 >    !! check to make sure the reaction field setting makes sense
548  
549      if (FF_uses_Dipoles) then
550 <       if (FF_uses_RF) then
550 >       if (electrostaticSummationMethod == 3) then
551            dielect = getDielect()
552            call initialize_rf(dielect)
553         endif
554      else
555 <       if (FF_uses_RF) then          
555 >       if (electrostaticSummationMethod == 3) then
556            write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
557            thisStat = -1
558            haveSaneForceField = .false.
# Line 575 | Line 734 | contains
734   #endif
735         outer: do i = istart, iend
736  
578 – #ifdef IS_MPI
579 –             me_i = atid_row(i)
580 – #else
581 –             me_i = atid(i)
582 – #endif
583 –
737            if (update_nlist) point(i) = nlist + 1
738  
739            n_in_i = groupStartRow(i+1) - groupStartRow(i)
# Line 617 | Line 770 | contains
770                    q_group(:,j), d_grp, rgrpsq)
771   #endif
772  
773 <             if (rgrpsq < InteractionHash(me_i,me_j)%rListsq) then
773 >             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
774                  if (update_nlist) then
775                     nlist = nlist + 1
776  
# Line 817 | Line 970 | contains
970  
971      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
972  
973 <       if (FF_uses_RF .and. SIM_uses_RF) then
973 >       if (electrostaticSummationMethod == 3) then
974  
975   #ifdef IS_MPI
976            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 901 | Line 1054 | contains
1054      real ( kind = dp ), intent(inout) :: rijsq
1055      real ( kind = dp )                :: r
1056      real ( kind = dp ), intent(inout) :: d(3)
904 –    real ( kind = dp ) :: ebalance
1057      integer :: me_i, me_j
1058  
1059      integer :: iHash
# Line 928 | Line 1080 | contains
1080         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1081              pot, eFrame, f, t, do_pot)
1082  
1083 <       if (FF_uses_RF .and. SIM_uses_RF) then
1083 >       if (electrostaticSummationMethod == 3) then
1084  
1085            ! CHECK ME (RF needs to know about all electrostatic types)
1086            call accumulate_rf(i, j, r, eFrame, sw)
# Line 991 | Line 1143 | contains
1143  
1144      integer :: me_i, me_j, iHash
1145  
1146 +    r = sqrt(rijsq)
1147 +
1148   #ifdef IS_MPI  
1149      me_i = atid_row(i)
1150      me_j = atid_col(j)  
# Line 1208 | Line 1362 | contains
1362  
1363    function FF_RequiresPostpairCalc() result(doesit)
1364      logical :: doesit
1365 <    doesit = FF_uses_RF
1365 >    if (electrostaticSummationMethod == 3) doesit = .true.
1366    end function FF_RequiresPostpairCalc
1367  
1368   #ifdef PROFILE

Diff Legend

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