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 2301 by gezelter, Thu Sep 15 22:05:21 2005 UTC vs.
Revision 2350 by chuckv, Mon Oct 10 21:20:46 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.44 2005-09-15 22:05:17 gezelter Exp $, $Date: 2005-09-15 22:05:17 $, $Name: not supported by cvs2svn $, $Revision: 1.44 $
48 > !! @version $Id: doForces.F90,v 1.51 2005-10-10 21:20:46 chuckv Exp $, $Date: 2005-10-10 21:20:46 $, $Name: not supported by cvs2svn $, $Revision: 1.51 $
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
93    logical, save :: FF_uses_Dipoles
94    logical, save :: FF_uses_GayBerne
95    logical, save :: FF_uses_EAM
94 –  logical, save :: FF_uses_RF
96  
97    logical, save :: SIM_uses_DirectionalAtoms
98    logical, save :: SIM_uses_EAM
98 –  logical, save :: SIM_uses_RF
99    logical, save :: SIM_requires_postpair_calc
100    logical, save :: SIM_requires_prepair_calc
101    logical, save :: SIM_uses_PBC
102  
103 <  integer, save :: corrMethod
103 >  integer, save :: electrostaticSummationMethod
104  
105    public :: init_FF
106    public :: setDefaultCutoffs
# Line 124 | 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 136 | 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 260 | Line 265 | contains
265      logical :: GtypeFound
266  
267      integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
268 <    integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
268 >    integer :: n_in_i, me_i, ia, g, atom1
269      integer :: nGroupsInRow
270 <    real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
270 >    integer :: nGroupsInCol
271 >    integer :: nGroupTypesRow,nGroupTypesCol
272 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin
273      real(kind=dp) :: biggestAtypeCutoff
274  
275      stat = 0
# Line 276 | Line 283 | contains
283      endif
284   #ifdef IS_MPI
285      nGroupsInRow = getNgroupsInRow(plan_group_row)
286 +    nGroupsInCol = getNgroupsInCol(plan_group_col)
287   #endif
288      nAtypes = getSize(atypes)
289   ! Set all of the initial cutoffs to zero.
# Line 291 | Line 299 | contains
299            call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
300            
301  
302 <          if (i_is_LJ) then
303 <             thisRcut = getSigma(i) * 2.5_dp
304 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
305 <          endif
306 <          if (i_is_Elect) then
307 <             thisRcut = defaultRcut
308 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
309 <          endif
310 <          if (i_is_Sticky) then
311 <             thisRcut = getStickyCut(i)
312 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
313 <          endif
314 <          if (i_is_StickyP) then
315 <             thisRcut = getStickyPowerCut(i)
316 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
317 <          endif
318 <          if (i_is_GB) then
319 <             thisRcut = getGayBerneCut(i)
320 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
321 <          endif
322 <          if (i_is_EAM) then
323 <             thisRcut = getEAMCut(i)
324 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
325 <          endif
326 <          if (i_is_Shape) then
327 <             thisRcut = getShapeCut(i)
328 <             if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
302 >          if (haveDefaultCutoffs) then
303 >             atypeMaxCutoff(i) = defaultRcut
304 >          else
305 >             if (i_is_LJ) then          
306 >                thisRcut = getSigma(i) * 2.5_dp
307 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
308 >             endif
309 >             if (i_is_Elect) then
310 >                thisRcut = defaultRcut
311 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
312 >             endif
313 >             if (i_is_Sticky) then
314 >                thisRcut = getStickyCut(i)
315 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
316 >             endif
317 >             if (i_is_StickyP) then
318 >                thisRcut = getStickyPowerCut(i)
319 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
320 >             endif
321 >             if (i_is_GB) then
322 >                thisRcut = getGayBerneCut(i)
323 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
324 >             endif
325 >             if (i_is_EAM) then
326 >                thisRcut = getEAMCut(i)
327 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
328 >             endif
329 >             if (i_is_Shape) then
330 >                thisRcut = getShapeCut(i)
331 >                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
332 >             endif
333            endif
334            
335 +          
336            if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
337               biggestAtypeCutoff = atypeMaxCutoff(i)
338            endif
339 +
340         endif
341      enddo
342    
343 <    nGroupTypes = 0
343 >
344      
345      istart = 1
346 +    jstart = 1
347   #ifdef IS_MPI
348      iend = nGroupsInRow
349 +    jend = nGroupsInCol
350   #else
351      iend = nGroups
352 +    jend = nGroups
353   #endif
354      
355      !! allocate the groupToGtype and gtypeMaxCutoff here.
356 <    if(.not.allocated(groupToGtype)) then
357 <       allocate(groupToGtype(iend))
358 <       allocate(groupMaxCutoff(iend))
359 <       allocate(gtypeMaxCutoff(iend))
360 <       groupMaxCutoff = 0.0_dp
361 <       gtypeMaxCutoff = 0.0_dp
356 >    if(.not.allocated(groupToGtypeRow)) then
357 >     !  allocate(groupToGtype(iend))
358 >       allocate(groupToGtypeRow(iend))
359 >    else
360 >       deallocate(groupToGtypeRow)
361 >       allocate(groupToGtypeRow(iend))
362      endif
363 +    if(.not.allocated(groupMaxCutoffRow)) then
364 +       allocate(groupMaxCutoffRow(iend))
365 +    else
366 +       deallocate(groupMaxCutoffRow)
367 +       allocate(groupMaxCutoffRow(iend))
368 +    end if
369 +
370 +    if(.not.allocated(gtypeMaxCutoffRow)) then
371 +       allocate(gtypeMaxCutoffRow(iend))
372 +    else
373 +       deallocate(gtypeMaxCutoffRow)
374 +       allocate(gtypeMaxCutoffRow(iend))
375 +    endif
376 +
377 +
378 + #ifdef IS_MPI
379 +       ! We only allocate new storage if we are in MPI because Ncol /= Nrow
380 +    if(.not.allocated(groupToGtypeCol)) then
381 +       allocate(groupToGtypeCol(jend))
382 +    else
383 +       deallocate(groupToGtypeCol)
384 +       allocate(groupToGtypeCol(jend))
385 +    end if
386 +
387 +    if(.not.allocated(groupToGtypeCol)) then
388 +       allocate(groupToGtypeCol(jend))
389 +    else
390 +       deallocate(groupToGtypeCol)
391 +       allocate(groupToGtypeCol(jend))
392 +    end if
393 +    if(.not.allocated(gtypeMaxCutoffCol)) then
394 +       allocate(gtypeMaxCutoffCol(jend))
395 +    else
396 +       deallocate(gtypeMaxCutoffCol)      
397 +       allocate(gtypeMaxCutoffCol(jend))
398 +    end if
399 +
400 +       groupMaxCutoffCol = 0.0_dp
401 +       gtypeMaxCutoffCol = 0.0_dp
402 +
403 + #endif
404 +       groupMaxCutoffRow = 0.0_dp
405 +       gtypeMaxCutoffRow = 0.0_dp
406 +
407 +
408      !! first we do a single loop over the cutoff groups to find the
409      !! largest cutoff for any atypes present in this group.  We also
410      !! create gtypes at this point.
411      
412      tol = 1.0d-6
413 <    
413 >    nGroupTypesRow = 0
414 >
415      do i = istart, iend      
416         n_in_i = groupStartRow(i+1) - groupStartRow(i)
417 <       groupMaxCutoff(i) = 0.0_dp
417 >       groupMaxCutoffRow(i) = 0.0_dp
418         do ia = groupStartRow(i), groupStartRow(i+1)-1
419            atom1 = groupListRow(ia)
420   #ifdef IS_MPI
# Line 359 | Line 422 | contains
422   #else
423            me_i = atid(atom1)
424   #endif          
425 <          if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
426 <             groupMaxCutoff(i)=atypeMaxCutoff(me_i)
425 >          if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
426 >             groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
427            endif          
428         enddo
429  
430 <       if (nGroupTypes.eq.0) then
431 <          nGroupTypes = nGroupTypes + 1
432 <          gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
433 <          groupToGtype(i) = nGroupTypes
430 >       if (nGroupTypesRow.eq.0) then
431 >          nGroupTypesRow = nGroupTypesRow + 1
432 >          gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
433 >          groupToGtypeRow(i) = nGroupTypesRow
434         else
435            GtypeFound = .false.
436 <          do g = 1, nGroupTypes
437 <             if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
438 <                groupToGtype(i) = g
436 >          do g = 1, nGroupTypesRow
437 >             if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
438 >                groupToGtypeRow(i) = g
439                  GtypeFound = .true.
440               endif
441            enddo
442            if (.not.GtypeFound) then            
443 <             nGroupTypes = nGroupTypes + 1
444 <             gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
445 <             groupToGtype(i) = nGroupTypes
443 >             nGroupTypesRow = nGroupTypesRow + 1
444 >             gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
445 >             groupToGtypeRow(i) = nGroupTypesRow
446            endif
447         endif
448      enddo    
449  
450 + #ifdef IS_MPI
451 +    do j = jstart, jend      
452 +       n_in_j = groupStartCol(j+1) - groupStartCol(j)
453 +       groupMaxCutoffCol(j) = 0.0_dp
454 +       do ja = groupStartCol(j), groupStartCol(j+1)-1
455 +          atom1 = groupListCol(ja)
456 +
457 +          me_j = atid_col(atom1)
458 +
459 +          if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
460 +             groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
461 +          endif          
462 +       enddo
463 +
464 +       if (nGroupTypesCol.eq.0) then
465 +          nGroupTypesCol = nGroupTypesCol + 1
466 +          gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
467 +          groupToGtypeCol(j) = nGroupTypesCol
468 +       else
469 +          GtypeFound = .false.
470 +          do g = 1, nGroupTypesCol
471 +             if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
472 +                groupToGtypeCol(j) = g
473 +                GtypeFound = .true.
474 +             endif
475 +          enddo
476 +          if (.not.GtypeFound) then            
477 +             nGroupTypesCol = nGroupTypesCol + 1
478 +             gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
479 +             groupToGtypeCol(j) = nGroupTypesCol
480 +          endif
481 +       endif
482 +    enddo    
483 +
484 + #else
485 + ! Set pointers to information we just found
486 +    nGroupTypesCol = nGroupTypesRow
487 +    groupToGtypeCol => groupToGtypeRow
488 +    gtypeMaxCutoffCol => gtypeMaxCutoffRow
489 +    groupMaxCutoffCol => groupMaxCutoffRow
490 + #endif
491 +
492 +
493 +
494 +
495 +
496      !! allocate the gtypeCutoffMap here.
497 <    allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
497 >    allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
498      !! then we do a double loop over all the group TYPES to find the cutoff
499      !! map between groups of two types
500 <    
501 <    do i = 1, nGroupTypes
502 <       do j = 1, nGroupTypes
500 >    tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
501 >
502 >    do i = 1, nGroupTypesRow
503 >       do j = 1, nGroupTypesCol
504        
505            select case(cutoffPolicy)
506            case(TRADITIONAL_CUTOFF_POLICY)
507 <             thisRcut = maxval(gtypeMaxCutoff)
507 >             thisRcut = tradRcut
508            case(MIX_CUTOFF_POLICY)
509 <             thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
509 >             thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
510            case(MAX_CUTOFF_POLICY)
511 <             thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
511 >             thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
512            case default
513               call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
514               return
# Line 406 | Line 516 | contains
516            gtypeCutoffMap(i,j)%rcut = thisRcut
517            gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
518            skin = defaultRlist - defaultRcut
519 +          listSkin = skin ! set neighbor list skin thickness
520            gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
521  
522 +          ! sanity check
523 +
524 +          if (haveDefaultCutoffs) then
525 +             if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
526 +                call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
527 +             endif
528 +          endif
529         enddo
530      enddo
531 +    if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
532 +    if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
533 +    if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
534 + #ifdef IS_MPI
535 +    if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
536 +    if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
537 + #endif
538 +    groupMaxCutoffCol => null()
539 +    gtypeMaxCutoffCol => null()
540      
541      haveGtypeCutoffMap = .true.
542     end subroutine createGtypeCutoffMap
# Line 422 | Line 549 | contains
549       defaultRsw = defRsw
550       defaultRlist = defRlist
551       cutoffPolicy = cutPolicy
552 <     rcuti = 1.0_dp / defaultRcut
552 >
553 >     haveDefaultCutoffs = .true.
554     end subroutine setDefaultCutoffs
555  
556     subroutine setCutoffPolicy(cutPolicy)
# Line 439 | Line 567 | contains
567      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
568      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
569      SIM_uses_PBC = SimUsesPBC()
442 –    SIM_uses_RF = SimUsesRF()
570  
571      haveSIMvariables = .true.
572  
# Line 506 | Line 633 | contains
633    end subroutine doReadyCheck
634  
635  
636 <  subroutine init_FF(use_RF, correctionMethod, dampingAlpha, thisStat)
636 >  subroutine init_FF(thisESM, thisStat)
637  
638 <    logical, intent(in) :: use_RF
512 <    integer, intent(in) :: correctionMethod
513 <    real(kind=dp), intent(in) :: dampingAlpha
638 >    integer, intent(in) :: thisESM
639      integer, intent(out) :: thisStat  
640      integer :: my_status, nMatches
641      integer, pointer :: MatchList(:) => null()
# Line 519 | Line 644 | contains
644      !! assume things are copacetic, unless they aren't
645      thisStat = 0
646  
647 <    !! Fortran's version of a cast:
523 <    FF_uses_RF = use_RF
647 >    electrostaticSummationMethod = thisESM
648  
525 –        
649      !! init_FF is called *after* all of the atom types have been
650      !! defined in atype_module using the new_atype subroutine.
651      !!
# Line 552 | Line 675 | contains
675  
676      haveSaneForceField = .true.
677  
678 <    !! check to make sure the FF_uses_RF setting makes sense
678 >    !! check to make sure the reaction field setting makes sense
679  
680      if (FF_uses_Dipoles) then
681 <       if (FF_uses_RF) then
681 >       if (electrostaticSummationMethod == REACTION_FIELD) then
682            dielect = getDielect()
683            call initialize_rf(dielect)
684         endif
685      else
686 <       if ((corrMethod == 3) .or. FF_uses_RF) then
686 >       if (electrostaticSummationMethod == REACTION_FIELD) then
687            write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
688            thisStat = -1
689            haveSaneForceField = .false.
# Line 651 | Line 774 | contains
774      integer :: propPack_i, propPack_j
775      integer :: loopStart, loopEnd, loop
776      integer :: iHash
777 <    real(kind=dp) :: listSkin = 1.0  
777 >  
778  
779      !! initialize local variables  
780  
# Line 776 | Line 899 | contains
899               me_j = atid(j)
900               call get_interatomic_vector(q_group(:,i), &
901                    q_group(:,j), d_grp, rgrpsq)
902 < #endif
902 > #endif      
903  
904 <             if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
904 >             if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
905                  if (update_nlist) then
906                     nlist = nlist + 1
907  
# Line 899 | Line 1022 | contains
1022                  endif
1023               end if
1024            enddo
1025 +
1026         enddo outer
1027  
1028         if (update_nlist) then
# Line 978 | Line 1102 | contains
1102  
1103      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1104  
1105 <       if ((FF_uses_RF .and. SIM_uses_RF) .or. (corrMethod == 3)) then
1105 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1106  
1107   #ifdef IS_MPI
1108            call scatter(rf_Row,rf,plan_atom_row_3d)
# Line 1086 | Line 1210 | contains
1210  
1211      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1212         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1213 <            pot, eFrame, f, t, do_pot, corrMethod, rcuti)
1213 >            pot, eFrame, f, t, do_pot)
1214  
1215 <       if ((FF_uses_RF .and. SIM_uses_RF) .or. (corrMethod == 3)) then
1215 >       if (electrostaticSummationMethod == REACTION_FIELD) then
1216  
1217            ! CHECK ME (RF needs to know about all electrostatic types)
1218            call accumulate_rf(i, j, r, eFrame, sw)
# Line 1370 | Line 1494 | contains
1494  
1495    function FF_RequiresPostpairCalc() result(doesit)
1496      logical :: doesit
1497 <    doesit = FF_uses_RF
1374 <    if (corrMethod == 3) doesit = .true.
1497 >    if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1498    end function FF_RequiresPostpairCalc
1499  
1500   #ifdef PROFILE

Diff Legend

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