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 2306 by chrisfen, Fri Sep 16 20:37:05 2005 UTC vs.
Revision 2354 by chuckv, Tue Oct 11 22:00:30 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.45 2005-09-16 20:36:55 chrisfen Exp $, $Date: 2005-09-16 20:36:55 $, $Name: not supported by cvs2svn $, $Revision: 1.45 $
48 > !! @version $Id: doForces.F90,v 1.53 2005-10-11 22:00:30 chuckv Exp $, $Date: 2005-10-11 22:00:30 $, $Name: not supported by cvs2svn $, $Revision: 1.53 $
49  
50  
51   module doForces
# Line 75 | Line 75 | module doForces
75   #include "UseTheForce/fSwitchingFunction.h"
76   #include "UseTheForce/fCutoffPolicy.h"
77   #include "UseTheForce/DarkSide/fInteractionMap.h"
78 + #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
79  
80  
81    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
# Line 85 | Line 86 | module doForces
86    logical, save :: haveSaneForceField = .false.
87    logical, save :: haveInteractionHash = .false.
88    logical, save :: haveGtypeCutoffMap = .false.
89 +  logical, save :: haveDefaultCutoffs = .false.
90    logical, save :: haveRlist = .false.
91  
92    logical, save :: FF_uses_DirectionalAtoms
# Line 122 | Line 124 | module doForces
124    ! Bit hash to determine pair-pair interactions.
125    integer, dimension(:,:), allocatable :: InteractionHash
126    real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
127 <  real(kind=dp), dimension(:), allocatable :: groupMaxCutoff
128 <  integer, dimension(:), allocatable :: groupToGtype
129 <  real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff
127 >  real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
128 >  real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
129 >
130 >  integer, dimension(:), allocatable, target :: groupToGtypeRow
131 >  integer, dimension(:), pointer :: groupToGtypeCol => null()
132 >
133 >  real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
134 >  real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
135    type ::gtypeCutoffs
136       real(kind=dp) :: rcut
137       real(kind=dp) :: rcutsq
# Line 134 | Line 141 | module doForces
141  
142    integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
143    real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
144 <  real(kind=dp),save :: rcuti
144 >  real(kind=dp),save :: listSkin
145    
146   contains
147  
# Line 177 | Line 184 | contains
184      end if
185  
186      if (.not. allocated(InteractionHash)) then
187 +       allocate(InteractionHash(nAtypes,nAtypes))
188 +    else
189 +       deallocate(InteractionHash)
190         allocate(InteractionHash(nAtypes,nAtypes))
191      endif
192  
193      if (.not. allocated(atypeMaxCutoff)) then
194         allocate(atypeMaxCutoff(nAtypes))
195 +    else
196 +       deallocate(atypeMaxCutoff)
197 +       allocate(atypeMaxCutoff(nAtypes))
198      endif
199          
200      do i = 1, nAtypes
# Line 258 | Line 271 | contains
271      logical :: GtypeFound
272  
273      integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
274 <    integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
274 >    integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
275      integer :: nGroupsInRow
276 <    real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
276 >    integer :: nGroupsInCol
277 >    integer :: nGroupTypesRow,nGroupTypesCol
278 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin
279      real(kind=dp) :: biggestAtypeCutoff
280  
281      stat = 0
# Line 274 | Line 289 | contains
289      endif
290   #ifdef IS_MPI
291      nGroupsInRow = getNgroupsInRow(plan_group_row)
292 +    nGroupsInCol = getNgroupsInCol(plan_group_col)
293   #endif
294      nAtypes = getSize(atypes)
295   ! Set all of the initial cutoffs to zero.
# Line 289 | Line 305 | contains
305            call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
306            
307  
308 <          if (i_is_LJ) then
309 <             thisRcut = getSigma(i) * 2.5_dp
310 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
311 <          endif
312 <          if (i_is_Elect) then
313 <             thisRcut = defaultRcut
314 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
315 <          endif
316 <          if (i_is_Sticky) then
317 <             thisRcut = getStickyCut(i)
318 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
319 <          endif
320 <          if (i_is_StickyP) then
321 <             thisRcut = getStickyPowerCut(i)
322 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
323 <          endif
324 <          if (i_is_GB) then
325 <             thisRcut = getGayBerneCut(i)
326 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
327 <          endif
328 <          if (i_is_EAM) then
329 <             thisRcut = getEAMCut(i)
330 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
331 <          endif
332 <          if (i_is_Shape) then
333 <             thisRcut = getShapeCut(i)
334 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
308 >          if (haveDefaultCutoffs) then
309 >             atypeMaxCutoff(i) = defaultRcut
310 >          else
311 >             if (i_is_LJ) then          
312 >                thisRcut = getSigma(i) * 2.5_dp
313 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
314 >             endif
315 >             if (i_is_Elect) then
316 >                thisRcut = defaultRcut
317 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
318 >             endif
319 >             if (i_is_Sticky) then
320 >                thisRcut = getStickyCut(i)
321 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
322 >             endif
323 >             if (i_is_StickyP) then
324 >                thisRcut = getStickyPowerCut(i)
325 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
326 >             endif
327 >             if (i_is_GB) then
328 >                thisRcut = getGayBerneCut(i)
329 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
330 >             endif
331 >             if (i_is_EAM) then
332 >                thisRcut = getEAMCut(i)
333 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
334 >             endif
335 >             if (i_is_Shape) then
336 >                thisRcut = getShapeCut(i)
337 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
338 >             endif
339            endif
340            
341 +          
342            if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
343               biggestAtypeCutoff = atypeMaxCutoff(i)
344            endif
345 +
346         endif
347      enddo
348    
349 <    nGroupTypes = 0
349 >
350      
351      istart = 1
352 +    jstart = 1
353   #ifdef IS_MPI
354      iend = nGroupsInRow
355 +    jend = nGroupsInCol
356   #else
357      iend = nGroups
358 +    jend = nGroups
359   #endif
360      
361      !! allocate the groupToGtype and gtypeMaxCutoff here.
362 <    if(.not.allocated(groupToGtype)) then
363 <       allocate(groupToGtype(iend))
364 <       allocate(groupMaxCutoff(iend))
365 <       allocate(gtypeMaxCutoff(iend))
366 <       groupMaxCutoff = 0.0_dp
367 <       gtypeMaxCutoff = 0.0_dp
362 >    if(.not.allocated(groupToGtypeRow)) then
363 >     !  allocate(groupToGtype(iend))
364 >       allocate(groupToGtypeRow(iend))
365 >    else
366 >       deallocate(groupToGtypeRow)
367 >       allocate(groupToGtypeRow(iend))
368      endif
369 +    if(.not.allocated(groupMaxCutoffRow)) then
370 +       allocate(groupMaxCutoffRow(iend))
371 +    else
372 +       deallocate(groupMaxCutoffRow)
373 +       allocate(groupMaxCutoffRow(iend))
374 +    end if
375 +
376 +    if(.not.allocated(gtypeMaxCutoffRow)) then
377 +       allocate(gtypeMaxCutoffRow(iend))
378 +    else
379 +       deallocate(gtypeMaxCutoffRow)
380 +       allocate(gtypeMaxCutoffRow(iend))
381 +    endif
382 +
383 +
384 + #ifdef IS_MPI
385 +       ! We only allocate new storage if we are in MPI because Ncol /= Nrow
386 +    if(.not.associated(groupToGtypeCol)) then
387 +       allocate(groupToGtypeCol(jend))
388 +    else
389 +       deallocate(groupToGtypeCol)
390 +       allocate(groupToGtypeCol(jend))
391 +    end if
392 +
393 +    if(.not.associated(groupToGtypeCol)) then
394 +       allocate(groupToGtypeCol(jend))
395 +    else
396 +       deallocate(groupToGtypeCol)
397 +       allocate(groupToGtypeCol(jend))
398 +    end if
399 +    if(.not.associated(gtypeMaxCutoffCol)) then
400 +       allocate(gtypeMaxCutoffCol(jend))
401 +    else
402 +       deallocate(gtypeMaxCutoffCol)      
403 +       allocate(gtypeMaxCutoffCol(jend))
404 +    end if
405 +
406 +       groupMaxCutoffCol = 0.0_dp
407 +       gtypeMaxCutoffCol = 0.0_dp
408 +
409 + #endif
410 +       groupMaxCutoffRow = 0.0_dp
411 +       gtypeMaxCutoffRow = 0.0_dp
412 +
413 +
414      !! first we do a single loop over the cutoff groups to find the
415      !! largest cutoff for any atypes present in this group.  We also
416      !! create gtypes at this point.
417      
418      tol = 1.0d-6
419 <    
419 >    nGroupTypesRow = 0
420 >
421      do i = istart, iend      
422         n_in_i = groupStartRow(i+1) - groupStartRow(i)
423 <       groupMaxCutoff(i) = 0.0_dp
423 >       groupMaxCutoffRow(i) = 0.0_dp
424         do ia = groupStartRow(i), groupStartRow(i+1)-1
425            atom1 = groupListRow(ia)
426   #ifdef IS_MPI
# Line 357 | Line 428 | contains
428   #else
429            me_i = atid(atom1)
430   #endif          
431 <          if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
432 <             groupMaxCutoff(i)=atypeMaxCutoff(me_i)
431 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
432 >             groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
433            endif          
434         enddo
435  
436 <       if (nGroupTypes.eq.0) then
437 <          nGroupTypes = nGroupTypes + 1
438 <          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
439 <          groupToGtype(i) = nGroupTypes
436 >       if (nGroupTypesRow.eq.0) then
437 >          nGroupTypesRow = nGroupTypesRow + 1
438 >          gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
439 >          groupToGtypeRow(i) = nGroupTypesRow
440         else
441            GtypeFound = .false.
442 <          do g = 1, nGroupTypes
443 <             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
444 <                groupToGtype(i) = g
442 >          do g = 1, nGroupTypesRow
443 >             if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
444 >                groupToGtypeRow(i) = g
445                  GtypeFound = .true.
446               endif
447            enddo
448            if (.not.GtypeFound) then            
449 <             nGroupTypes = nGroupTypes + 1
450 <             gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
451 <             groupToGtype(i) = nGroupTypes
449 >             nGroupTypesRow = nGroupTypesRow + 1
450 >             gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
451 >             groupToGtypeRow(i) = nGroupTypesRow
452            endif
453         endif
454      enddo    
455  
456 + #ifdef IS_MPI
457 +    do j = jstart, jend      
458 +       n_in_j = groupStartCol(j+1) - groupStartCol(j)
459 +       groupMaxCutoffCol(j) = 0.0_dp
460 +       do ja = groupStartCol(j), groupStartCol(j+1)-1
461 +          atom1 = groupListCol(ja)
462 +
463 +          me_j = atid_col(atom1)
464 +
465 +          if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
466 +             groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
467 +          endif          
468 +       enddo
469 +
470 +       if (nGroupTypesCol.eq.0) then
471 +          nGroupTypesCol = nGroupTypesCol + 1
472 +          gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
473 +          groupToGtypeCol(j) = nGroupTypesCol
474 +       else
475 +          GtypeFound = .false.
476 +          do g = 1, nGroupTypesCol
477 +             if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
478 +                groupToGtypeCol(j) = g
479 +                GtypeFound = .true.
480 +             endif
481 +          enddo
482 +          if (.not.GtypeFound) then            
483 +             nGroupTypesCol = nGroupTypesCol + 1
484 +             gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
485 +             groupToGtypeCol(j) = nGroupTypesCol
486 +          endif
487 +       endif
488 +    enddo    
489 +
490 + #else
491 + ! Set pointers to information we just found
492 +    nGroupTypesCol = nGroupTypesRow
493 +    groupToGtypeCol => groupToGtypeRow
494 +    gtypeMaxCutoffCol => gtypeMaxCutoffRow
495 +    groupMaxCutoffCol => groupMaxCutoffRow
496 + #endif
497 +
498 +
499 +
500 +
501 +
502      !! allocate the gtypeCutoffMap here.
503 <    allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
503 >    allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
504      !! then we do a double loop over all the group TYPES to find the cutoff
505      !! map between groups of two types
506 <    
507 <    do i = 1, nGroupTypes
508 <       do j = 1, nGroupTypes
506 >    tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
507 >
508 >    do i = 1, nGroupTypesRow
509 >       do j = 1, nGroupTypesCol
510        
511            select case(cutoffPolicy)
512            case(TRADITIONAL_CUTOFF_POLICY)
513 <             thisRcut = maxval(gtypeMaxCutoff)
513 >             thisRcut = tradRcut
514            case(MIX_CUTOFF_POLICY)
515 <             thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
515 >             thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
516            case(MAX_CUTOFF_POLICY)
517 <             thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
517 >             thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
518            case default
519               call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
520               return
# Line 404 | Line 522 | contains
522            gtypeCutoffMap(i,j)%rcut = thisRcut
523            gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
524            skin = defaultRlist - defaultRcut
525 +          listSkin = skin ! set neighbor list skin thickness
526            gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
527  
528 +          ! sanity check
529 +
530 +          if (haveDefaultCutoffs) then
531 +             if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
532 +                call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
533 +             endif
534 +          endif
535         enddo
536      enddo
537 +    if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
538 +    if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
539 +    if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
540 + #ifdef IS_MPI
541 +    if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
542 +    if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
543 + #endif
544 +    groupMaxCutoffCol => null()
545 +    gtypeMaxCutoffCol => null()
546      
547      haveGtypeCutoffMap = .true.
548     end subroutine createGtypeCutoffMap
# Line 420 | Line 555 | contains
555       defaultRsw = defRsw
556       defaultRlist = defRlist
557       cutoffPolicy = cutPolicy
558 <     rcuti = 1.0_dp / defaultRcut
558 >
559 >     haveDefaultCutoffs = .true.
560     end subroutine setDefaultCutoffs
561  
562     subroutine setCutoffPolicy(cutPolicy)
# Line 503 | Line 639 | contains
639    end subroutine doReadyCheck
640  
641  
642 <  subroutine init_FF(thisESM, dampingAlpha, thisStat)
642 >  subroutine init_FF(thisESM, thisStat)
643  
644      integer, intent(in) :: thisESM
509 –    real(kind=dp), intent(in) :: dampingAlpha
645      integer, intent(out) :: thisStat  
646      integer :: my_status, nMatches
647      integer, pointer :: MatchList(:) => null()
# Line 549 | Line 684 | contains
684      !! check to make sure the reaction field setting makes sense
685  
686      if (FF_uses_Dipoles) then
687 <       if (electrostaticSummationMethod == 3) then
687 >       if (electrostaticSummationMethod == REACTION_FIELD) then
688            dielect = getDielect()
689            call initialize_rf(dielect)
690         endif
691      else
692 <       if (electrostaticSummationMethod == 3) then
692 >       if (electrostaticSummationMethod == REACTION_FIELD) then
693            write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
694            thisStat = -1
695            haveSaneForceField = .false.
# Line 645 | Line 780 | contains
780      integer :: propPack_i, propPack_j
781      integer :: loopStart, loopEnd, loop
782      integer :: iHash
783 <    real(kind=dp) :: listSkin = 1.0  
783 >  
784  
785      !! initialize local variables  
786  
# Line 770 | Line 905 | contains
905               me_j = atid(j)
906               call get_interatomic_vector(q_group(:,i), &
907                    q_group(:,j), d_grp, rgrpsq)
908 < #endif
908 > #endif      
909  
910 <             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
910 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
911                  if (update_nlist) then
912                     nlist = nlist + 1
913  
# Line 893 | Line 1028 | contains
1028                  endif
1029               end if
1030            enddo
1031 +
1032         enddo outer
1033  
1034         if (update_nlist) then
# Line 972 | Line 1108 | contains
1108  
1109      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1110  
1111 <       if (electrostaticSummationMethod == 3) then
1111 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1112  
1113   #ifdef IS_MPI
1114            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 1082 | Line 1218 | contains
1218         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1219              pot, eFrame, f, t, do_pot)
1220  
1221 <       if (electrostaticSummationMethod == 3) then
1221 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1222  
1223            ! CHECK ME (RF needs to know about all electrostatic types)
1224            call accumulate_rf(i, j, r, eFrame, sw)
# Line 1364 | Line 1500 | contains
1500  
1501    function FF_RequiresPostpairCalc() result(doesit)
1502      logical :: doesit
1503 <    if (electrostaticSummationMethod == 3) doesit = .true.
1503 >    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1504    end function FF_RequiresPostpairCalc
1505  
1506   #ifdef PROFILE

Diff Legend

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