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 2280 by gezelter, Thu Sep 1 20:17:55 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.34 2005-09-01 20:17:55 gezelter Exp $, $Date: 2005-09-01 20:17:55 $, $Name: not supported by cvs2svn $, $Revision: 1.34 $
49  
50  
51   module doForces
# 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    
140   contains
141  
# Line 146 | Line 160 | contains
160      logical :: j_is_GB
161      logical :: j_is_EAM
162      logical :: j_is_Shape
163 <    
163 >    real(kind=dp) :: myRcut
164 >
165      status = 0  
166  
167      if (.not. associated(atypes)) then
# Line 179 | Line 194 | contains
194         call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
195         call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
196  
182 –       if (i_is_LJ) then
183 –          thisCut = getDefaultLJCutoff(i)
184 –          if (thisCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisCut
185 –       endif
186 –
187 –
188 –
197         do j = i, nAtypes
198  
199            iHash = 0
# Line 238 | Line 246 | contains
246      haveInteractionHash = .true.
247    end subroutine createInteractionHash
248  
249 <  subroutine createGtypeCutoffMap(defaultRcut, defaultSkinThickness, stat)
249 >  subroutine createGtypeCutoffMap(stat)
250  
251 <    real(kind=dp), intent(in), optional :: defaultRCut, defaultSkinThickness
252 <    integer, intent(out) :: stat
251 >    integer, intent(out), optional :: stat
252 >    logical :: i_is_LJ
253 >    logical :: i_is_Elect
254 >    logical :: i_is_Sticky
255 >    logical :: i_is_StickyP
256 >    logical :: i_is_GB
257 >    logical :: i_is_EAM
258 >    logical :: i_is_Shape
259  
260 <    integer :: myStatus, nAtypes
260 >    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
261 >    integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
262 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
263 >    real(kind=dp) :: biggestAtypeCutoff
264  
265      stat = 0
266      if (.not. haveInteractionHash) then
# Line 256 | Line 273 | contains
273      endif
274  
275      nAtypes = getSize(atypes)
276 <
276 >    
277      do i = 1, nAtypes
278 <      
279 <       atypeMaxCutoff(i) =
280 <
278 >       if (SimHasAtype(i)) then          
279 >          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
280 >          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
281 >          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
282 >          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
283 >          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
284 >          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
285 >          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
286 >          
287 >          if (i_is_LJ) then
288 >             thisRcut = getSigma(i) * 2.5_dp
289 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
290 >          endif
291 >          if (i_is_Elect) then
292 >             thisRcut = defaultRcut
293 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
294 >          endif
295 >          if (i_is_Sticky) then
296 >             thisRcut = getStickyCut(i)
297 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
298 >          endif
299 >          if (i_is_StickyP) then
300 >             thisRcut = getStickyPowerCut(i)
301 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
302 >          endif
303 >          if (i_is_GB) then
304 >             thisRcut = getGayBerneCut(i)
305 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
306 >          endif
307 >          if (i_is_EAM) then
308 >             thisRcut = getEAMCut(i)
309 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
310 >          endif
311 >          if (i_is_Shape) then
312 >             thisRcut = getShapeCut(i)
313 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
314 >          endif
315 >          
316 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
317 >             biggestAtypeCutoff = atypeMaxCutoff(i)
318 >          endif
319 >       endif
320 >    enddo
321 >  
322 >    nGroupTypes = 0
323      
324 +    istart = 1
325 + #ifdef IS_MPI
326 +    iend = nGroupsInRow
327 + #else
328 +    iend = nGroups
329 + #endif
330  
331 <
331 >    !! allocate the groupToGtype and gtypeMaxCutoff here.
332 >    
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.
335 >    tol = 1.0d-6
336 >    
337 >    do i = istart, iend      
338 >       n_in_i = groupStartRow(i+1) - groupStartRow(i)
339 >       groupMaxCutoff(i) = 0.0_dp
340 >       do ia = groupStartRow(i), groupStartRow(i+1)-1
341 >          atom1 = groupListRow(ia)
342 > #ifdef IS_MPI
343 >          me_i = atid_row(atom1)
344 > #else
345 >          me_i = atid(atom1)
346 > #endif          
347 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
348 >             groupMaxCutoff(i)=atypeMaxCutoff(me_i)
349 >          endif
350 >       enddo
351 >       if (nGroupTypes.eq.0) then
352 >          nGroupTypes = nGroupTypes + 1
353 >          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
354 >          groupToGtype(i) = nGroupTypes
355 >       else
356 >          do g = 1, nGroupTypes
357 >             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).gt.tol) then
358 >                nGroupTypes = nGroupTypes + 1
359 >                gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
360 >                groupToGtype(i) = nGroupTypes
361 >             else
362 >                groupToGtype(i) = g
363 >             endif
364 >          enddo
365 >       endif
366 >    enddo
367 >    
368 >    !! allocate the gtypeCutoffMap here.
369  
370 <     haveGtypeCutoffMap = .true.
370 >    !! then we do a double loop over all the group TYPES to find the cutoff
371 >    !! map between groups of two types
372 >    
373 >    do i = 1, nGroupTypes
374 >       do j = 1, nGroupTypes
375 >      
376 >          select case(cutoffPolicy)
377 >             case(TRADITIONAL_CUTOFF_POLICY)
378 >                thisRcut = maxval(gtypeMaxCutoff)
379 >             case(MIX_CUTOFF_POLICY)
380 >                thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
381 >             case(MAX_CUTOFF_POLICY)
382 >                thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
383 >             case default
384 >                call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
385 >                return
386 >          end select      
387 >         gtypeCutoffMap(i,j)%rcut = thisRcut
388 >         gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
389 >         skin = defaultRlist - defaultRcut
390 >         gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
391 >       enddo
392 >    enddo
393 >    
394 >    haveGtypeCutoffMap = .true.
395     end subroutine createGtypeCutoffMap
396  
397 +   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
398 +     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
399 +     integer, intent(in) :: cutPolicy
400 +
401 +     defaultRcut = defRcut
402 +     defaultRsw = defRsw
403 +     defaultRlist = defRlist
404 +     cutoffPolicy = cutPolicy
405 +   end subroutine setDefaultCutoffs
406 +
407 +   subroutine setCutoffPolicy(cutPolicy)
408 +
409 +     integer, intent(in) :: cutPolicy
410 +     cutoffPolicy = cutPolicy
411 +     call createGtypeCutoffMap()
412 +
413 +   end subroutine setCutoffPolicy
414 +    
415 +    
416    subroutine setSimVariables()
417      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
418      SIM_uses_EAM = SimUsesEAM()
# Line 341 | Line 486 | contains
486    end subroutine doReadyCheck
487  
488  
489 <  subroutine init_FF(use_RF_c, thisStat)
489 >  subroutine init_FF(use_RF, use_UW, use_DW, thisStat)
490  
491 <    logical, intent(in) :: use_RF_c
492 <
491 >    logical, intent(in) :: use_RF
492 >    logical, intent(in) :: use_UW
493 >    logical, intent(in) :: use_DW
494      integer, intent(out) :: thisStat  
495      integer :: my_status, nMatches
496 +    integer :: corrMethod
497      integer, pointer :: MatchList(:) => null()
498      real(kind=dp) :: rcut, rrf, rt, dielect
499  
# Line 354 | Line 501 | contains
501      thisStat = 0
502  
503      !! Fortran's version of a cast:
504 <    FF_uses_RF = use_RF_c
504 >    FF_uses_RF = use_RF
505  
506 +    !! set the electrostatic correction method
507 +    if (use_UW) then
508 +       corrMethod = 1
509 +    elseif (use_DW) then
510 +       corrMethod = 2
511 +    else
512 +       corrMethod = 0
513 +    endif
514 +    
515      !! init_FF is called *after* all of the atom types have been
516      !! defined in atype_module using the new_atype subroutine.
517      !!
# Line 575 | Line 731 | contains
731   #endif
732         outer: do i = istart, iend
733  
578 – #ifdef IS_MPI
579 –             me_i = atid_row(i)
580 – #else
581 –             me_i = atid(i)
582 – #endif
583 –
734            if (update_nlist) point(i) = nlist + 1
735  
736            n_in_i = groupStartRow(i+1) - groupStartRow(i)
# Line 617 | Line 767 | contains
767                    q_group(:,j), d_grp, rgrpsq)
768   #endif
769  
770 <             if (rgrpsq < InteractionHash(me_i,me_j)%rListsq) then
770 >             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
771                  if (update_nlist) then
772                     nlist = nlist + 1
773  
# Line 926 | Line 1076 | contains
1076  
1077      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1078         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1079 <            pot, eFrame, f, t, do_pot)
1079 >            pot, eFrame, f, t, do_pot, corrMethod)
1080  
1081         if (FF_uses_RF .and. SIM_uses_RF) then
1082  

Diff Legend

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