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 2286 by gezelter, Wed Sep 7 22:08:39 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.39 2005-09-07 22:08:39 gezelter Exp $, $Date: 2005-09-07 22:08:39 $, $Name: not supported by cvs2svn $, $Revision: 1.39 $
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 >    logical :: GtypeFound
260  
261 <    integer :: myStatus, nAtypes
261 >    integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
262 >    integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
263 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
264 >    real(kind=dp) :: biggestAtypeCutoff
265  
266      stat = 0
267      if (.not. haveInteractionHash) then
# Line 256 | Line 274 | contains
274      endif
275  
276      nAtypes = getSize(atypes)
277 <
277 >    
278      do i = 1, nAtypes
279 <      
280 <       atypeMaxCutoff(i) =
279 >       if (SimHasAtype(i)) then    
280 >          call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
281 >          call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
282 >          call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
283 >          call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
284 >          call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
285 >          call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
286 >          call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
287 >          
288 >          atypeMaxCutoff(i) = 0.0_dp
289 >          if (i_is_LJ) then
290 >             thisRcut = getSigma(i) * 2.5_dp
291 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
292 >          endif
293 >          if (i_is_Elect) then
294 >             thisRcut = defaultRcut
295 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
296 >          endif
297 >          if (i_is_Sticky) then
298 >             thisRcut = getStickyCut(i)
299 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
300 >          endif
301 >          if (i_is_StickyP) then
302 >             thisRcut = getStickyPowerCut(i)
303 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
304 >          endif
305 >          if (i_is_GB) then
306 >             thisRcut = getGayBerneCut(i)
307 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
308 >          endif
309 >          if (i_is_EAM) then
310 >             thisRcut = getEAMCut(i)
311 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
312 >          endif
313 >          if (i_is_Shape) then
314 >             thisRcut = getShapeCut(i)
315 >             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
316 >          endif
317 >          
318 >          if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
319 >             biggestAtypeCutoff = atypeMaxCutoff(i)
320 >          endif
321 >       endif
322 >    enddo
323 >  
324 >    nGroupTypes = 0
325 >    
326 >    istart = 1
327 > #ifdef IS_MPI
328 >    iend = nGroupsInRow
329 > #else
330 >    iend = nGroups
331 > #endif
332 >    
333 >    !! allocate the groupToGtype and gtypeMaxCutoff here.
334 >    if(.not.allocated(groupToGtype)) then
335 >       allocate(groupToGtype(iend))
336 >       allocate(groupMaxCutoff(iend))
337 >       allocate(gtypeMaxCutoff(iend))
338 >    endif
339 >    !! first we do a single loop over the cutoff groups to find the
340 >    !! largest cutoff for any atypes present in this group.  We also
341 >    !! create gtypes at this point.
342 >    
343 >    tol = 1.0d-6
344 >    
345 >    do i = istart, iend      
346 >       n_in_i = groupStartRow(i+1) - groupStartRow(i)
347 >       groupMaxCutoff(i) = 0.0_dp
348 >       do ia = groupStartRow(i), groupStartRow(i+1)-1
349 >          atom1 = groupListRow(ia)
350 > #ifdef IS_MPI
351 >          me_i = atid_row(atom1)
352 > #else
353 >          me_i = atid(atom1)
354 > #endif          
355 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
356 >             groupMaxCutoff(i)=atypeMaxCutoff(me_i)
357 >          endif          
358 >       enddo
359  
360 +       if (nGroupTypes.eq.0) then
361 +          nGroupTypes = nGroupTypes + 1
362 +          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
363 +          groupToGtype(i) = nGroupTypes
364 +       else
365 +          GtypeFound = .false.
366 +          do g = 1, nGroupTypes
367 +             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
368 +                groupToGtype(i) = g
369 +                GtypeFound = .true.
370 +             endif
371 +          enddo
372 +          if (.not.GtypeFound) then            
373 +             nGroupTypes = nGroupTypes + 1
374 +             gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
375 +             groupToGtype(i) = nGroupTypes
376 +          endif
377 +       endif
378 +    enddo    
379 +
380 +    !! allocate the gtypeCutoffMap here.
381 +    allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
382 +    !! then we do a double loop over all the group TYPES to find the cutoff
383 +    !! map between groups of two types
384      
385 +    do i = 1, nGroupTypes
386 +       do j = 1, nGroupTypes
387 +      
388 +          select case(cutoffPolicy)
389 +          case(TRADITIONAL_CUTOFF_POLICY)
390 +             thisRcut = maxval(gtypeMaxCutoff)
391 +          case(MIX_CUTOFF_POLICY)
392 +             thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
393 +          case(MAX_CUTOFF_POLICY)
394 +             thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
395 +          case default
396 +             call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
397 +             return
398 +          end select
399 +          gtypeCutoffMap(i,j)%rcut = thisRcut
400 +          gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
401 +          skin = defaultRlist - defaultRcut
402 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
403  
404 +       enddo
405 +    enddo
406 +    
407 +    haveGtypeCutoffMap = .true.
408 +    
409 +  end subroutine createGtypeCutoffMap
410 +  
411 +  subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
412 +    real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
413 +    integer, intent(in) :: cutPolicy
414 +    
415 +    defaultRcut = defRcut
416 +    defaultRsw = defRsw
417 +    defaultRlist = defRlist
418 +    cutoffPolicy = cutPolicy
419 +  end subroutine setDefaultCutoffs
420 +  
421 +  subroutine setCutoffPolicy(cutPolicy)
422  
423 +     integer, intent(in) :: cutPolicy
424 +     cutoffPolicy = cutPolicy
425 +     call createGtypeCutoffMap()
426  
427 <     haveGtypeCutoffMap = .true.
428 <   end subroutine createGtypeCutoffMap
429 <
427 >   end subroutine setCutoffPolicy
428 >    
429 >    
430    subroutine setSimVariables()
431      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
432      SIM_uses_EAM = SimUsesEAM()
# Line 312 | Line 471 | contains
471         call setSimVariables()
472      endif
473  
474 <    if (.not. haveRlist) then
475 <       write(default_error, *) 'rList has not been set in doForces!'
476 <       error = -1
477 <       return
478 <    endif
474 >  !  if (.not. haveRlist) then
475 >  !     write(default_error, *) 'rList has not been set in doForces!'
476 >  !     error = -1
477 >  !     return
478 >  !  endif
479  
480      if (.not. haveNeighborList) then
481         write(default_error, *) 'neighbor list has not been initialized in doForces!'
# Line 341 | Line 500 | contains
500    end subroutine doReadyCheck
501  
502  
503 <  subroutine init_FF(use_RF_c, thisStat)
345 <
346 <    logical, intent(in) :: use_RF_c
503 >  subroutine init_FF(use_RF, use_UW, use_DW, thisStat)
504  
505 +    logical, intent(in) :: use_RF
506 +    logical, intent(in) :: use_UW
507 +    logical, intent(in) :: use_DW
508      integer, intent(out) :: thisStat  
509      integer :: my_status, nMatches
510 +    integer :: corrMethod
511      integer, pointer :: MatchList(:) => null()
512      real(kind=dp) :: rcut, rrf, rt, dielect
513  
# Line 354 | Line 515 | contains
515      thisStat = 0
516  
517      !! Fortran's version of a cast:
518 <    FF_uses_RF = use_RF_c
518 >    FF_uses_RF = use_RF
519  
520 +    !! set the electrostatic correction method
521 +    if (use_UW) then
522 +       corrMethod = 1
523 +    elseif (use_DW) then
524 +       corrMethod = 2
525 +    else
526 +       corrMethod = 0
527 +    endif
528 +    
529      !! init_FF is called *after* all of the atom types have been
530      !! defined in atype_module using the new_atype subroutine.
531      !!
# Line 574 | Line 744 | contains
744         iend = nGroups - 1
745   #endif
746         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
747  
748            if (update_nlist) point(i) = nlist + 1
749  
# Line 617 | Line 781 | contains
781                    q_group(:,j), d_grp, rgrpsq)
782   #endif
783  
784 <             if (rgrpsq < InteractionHash(me_i,me_j)%rListsq) then
784 >             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
785                  if (update_nlist) then
786                     nlist = nlist + 1
787  
# Line 926 | Line 1090 | contains
1090  
1091      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1092         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1093 <            pot, eFrame, f, t, do_pot)
1093 >            pot, eFrame, f, t, do_pot, corrMethod)
1094  
1095         if (FF_uses_RF .and. SIM_uses_RF) then
1096  

Diff Legend

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