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 2273 by gezelter, Thu Aug 11 21:04:03 2005 UTC vs.
Revision 2290 by gezelter, Wed Sep 14 20:28:05 2005 UTC

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

Diff Legend

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