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 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.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.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
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 260 | Line 274 | contains
274      endif
275  
276      nAtypes = getSize(atypes)
277 <
277 >    
278      do i = 1, nAtypes
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 <          thisCut = getSigma(i) * DEFAULT_SIGMA_MULTIPLIER
289 <          if (thisCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisCut
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 <       if (i_is_Elect) then
323 <          thisCut =
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 <     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 346 | 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 375 | 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 388 | 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 608 | Line 744 | contains
744         iend = nGroups - 1
745   #endif
746         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
747  
748            if (update_nlist) point(i) = nlist + 1
749  
# Line 651 | 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 960 | 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