ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/UseTheForce/doForces.F90
(Generate patch)

Comparing trunk/OOPSE-2.0/src/UseTheForce/doForces.F90 (file contents):
Revision 2274 by gezelter, Wed Aug 17 15:26:42 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.30 2005-08-17 15:26:37 gezelter Exp $, $Date: 2005-08-17 15:26:37 $, $Name: not supported by cvs2svn $, $Revision: 1.30 $
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 85 | 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 99 | 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 152 | 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 237 | Line 246 | contains
246      haveInteractionHash = .true.
247    end subroutine createInteractionHash
248  
249 <  subroutine createGtypeCutoffMap()
249 >  subroutine createGtypeCutoffMap(stat)
250  
251 +    integer, intent(out), optional :: stat
252      logical :: i_is_LJ
253      logical :: i_is_Elect
254      logical :: i_is_Sticky
# Line 246 | Line 256 | contains
256      logical :: i_is_GB
257      logical :: i_is_EAM
258      logical :: i_is_Shape
259 +    logical :: GtypeFound
260  
261 <    integer :: myStatus, nAtypes
262 <    real(kind=dp):: thisSigma, bigSigma
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 263 | Line 276 | contains
276      nAtypes = getSize(atypes)
277      
278      do i = 1, nAtypes
279 <       if (SimHasAtype(i)) then          
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)
# Line 272 | Line 285 | contains
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
# Line 284 | Line 298 | contains
298               thisRcut = getStickyCut(i)
299               if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
300            endif
301 <          if (i_is_StickyPower) then
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_GayBerne) then
305 >          if (i_is_GB) then
306               thisRcut = getGayBerneCut(i)
307               if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
308            endif
# Line 306 | Line 320 | contains
320            endif
321         endif
322      enddo
323 <
323 >  
324 >    nGroupTypes = 0
325 >    
326      istart = 1
327   #ifdef IS_MPI
328      iend = nGroupsInRow
329   #else
330      iend = nGroups
331   #endif
332 <    outer: do i = istart, iend
333 <
334 <       n_in_i = groupStartRow(i+1) - groupStartRow(i)
335 <
336 < #ifdef IS_MPI
337 <             jstart = 1
338 <             jend = nGroupsInCol
339 < #else
340 <             jstart = i+1
341 <             jend = nGroups
342 < #endif
343 <
344 <
345 <
346 <            
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 <     haveGtypeCutoffMap = .true.
405 <   end subroutine createGtypeCutoffMap
406 <
407 <   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
408 <     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
409 <     integer, intent(in) :: cutPolicy
410 <
411 <     defaultRcut = defRcut
412 <     defaultRsw = defRsw
413 <     defaultRlist = defRlist
414 <     cutoffPolicy = cutPolicy
415 <   end subroutine setDefaultCutoffs
416 <
417 <   subroutine setCutoffPolicy(cutPolicy)
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 <   end subroutine setDefaultCutoffs
427 >   end subroutine setCutoffPolicy
428      
429      
430    subroutine setSimVariables()
# Line 397 | 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 426 | Line 500 | contains
500    end subroutine doReadyCheck
501  
502  
503 <  subroutine init_FF(use_RF_c, thisStat)
503 >  subroutine init_FF(use_RF, use_UW, use_DW, thisStat)
504  
505 <    logical, intent(in) :: use_RF_c
506 <
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 439 | 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 696 | 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 1005 | 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