--- trunk/OOPSE-4/src/UseTheForce/doForces.F90 2005/09/15 22:05:21 2301 +++ trunk/OOPSE-4/src/UseTheForce/doForces.F90 2005/10/11 22:00:30 2354 @@ -45,7 +45,7 @@ !! @author Charles F. Vardeman II !! @author Matthew Meineke -!! @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 $ +!! @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 $ module doForces @@ -75,6 +75,7 @@ module doForces #include "UseTheForce/fSwitchingFunction.h" #include "UseTheForce/fCutoffPolicy.h" #include "UseTheForce/DarkSide/fInteractionMap.h" +#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h" INTEGER, PARAMETER:: PREPAIR_LOOP = 1 @@ -85,22 +86,21 @@ module doForces logical, save :: haveSaneForceField = .false. logical, save :: haveInteractionHash = .false. logical, save :: haveGtypeCutoffMap = .false. + logical, save :: haveDefaultCutoffs = .false. logical, save :: haveRlist = .false. logical, save :: FF_uses_DirectionalAtoms logical, save :: FF_uses_Dipoles logical, save :: FF_uses_GayBerne logical, save :: FF_uses_EAM - logical, save :: FF_uses_RF logical, save :: SIM_uses_DirectionalAtoms logical, save :: SIM_uses_EAM - logical, save :: SIM_uses_RF logical, save :: SIM_requires_postpair_calc logical, save :: SIM_requires_prepair_calc logical, save :: SIM_uses_PBC - integer, save :: corrMethod + integer, save :: electrostaticSummationMethod public :: init_FF public :: setDefaultCutoffs @@ -124,9 +124,14 @@ module doForces ! Bit hash to determine pair-pair interactions. integer, dimension(:,:), allocatable :: InteractionHash real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff - real(kind=dp), dimension(:), allocatable :: groupMaxCutoff - integer, dimension(:), allocatable :: groupToGtype - real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff + real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow + real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol + + integer, dimension(:), allocatable, target :: groupToGtypeRow + integer, dimension(:), pointer :: groupToGtypeCol => null() + + real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow + real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol type ::gtypeCutoffs real(kind=dp) :: rcut real(kind=dp) :: rcutsq @@ -136,7 +141,7 @@ module doForces integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist - real(kind=dp),save :: rcuti + real(kind=dp),save :: listSkin contains @@ -180,10 +185,16 @@ contains if (.not. allocated(InteractionHash)) then allocate(InteractionHash(nAtypes,nAtypes)) + else + deallocate(InteractionHash) + allocate(InteractionHash(nAtypes,nAtypes)) endif if (.not. allocated(atypeMaxCutoff)) then allocate(atypeMaxCutoff(nAtypes)) + else + deallocate(atypeMaxCutoff) + allocate(atypeMaxCutoff(nAtypes)) endif do i = 1, nAtypes @@ -260,9 +271,11 @@ contains logical :: GtypeFound integer :: myStatus, nAtypes, i, j, istart, iend, jstart, jend - integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes + integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j integer :: nGroupsInRow - real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin + integer :: nGroupsInCol + integer :: nGroupTypesRow,nGroupTypesCol + real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin real(kind=dp) :: biggestAtypeCutoff stat = 0 @@ -276,6 +289,7 @@ contains endif #ifdef IS_MPI nGroupsInRow = getNgroupsInRow(plan_group_row) + nGroupsInCol = getNgroupsInCol(plan_group_col) #endif nAtypes = getSize(atypes) ! Set all of the initial cutoffs to zero. @@ -291,67 +305,122 @@ contains call getElementProperty(atypes, i, "is_Shape", i_is_Shape) - if (i_is_LJ) then - thisRcut = getSigma(i) * 2.5_dp - if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut - endif - if (i_is_Elect) then - thisRcut = defaultRcut - if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut - endif - if (i_is_Sticky) then - thisRcut = getStickyCut(i) - if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut - endif - if (i_is_StickyP) then - thisRcut = getStickyPowerCut(i) - if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut - endif - if (i_is_GB) then - thisRcut = getGayBerneCut(i) - if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut - endif - if (i_is_EAM) then - thisRcut = getEAMCut(i) - if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut - endif - if (i_is_Shape) then - thisRcut = getShapeCut(i) - if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut + if (haveDefaultCutoffs) then + atypeMaxCutoff(i) = defaultRcut + else + if (i_is_LJ) then + thisRcut = getSigma(i) * 2.5_dp + if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut + endif + if (i_is_Elect) then + thisRcut = defaultRcut + if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut + endif + if (i_is_Sticky) then + thisRcut = getStickyCut(i) + if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut + endif + if (i_is_StickyP) then + thisRcut = getStickyPowerCut(i) + if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut + endif + if (i_is_GB) then + thisRcut = getGayBerneCut(i) + if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut + endif + if (i_is_EAM) then + thisRcut = getEAMCut(i) + if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut + endif + if (i_is_Shape) then + thisRcut = getShapeCut(i) + if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut + endif endif + if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then biggestAtypeCutoff = atypeMaxCutoff(i) endif + endif enddo - nGroupTypes = 0 + istart = 1 + jstart = 1 #ifdef IS_MPI iend = nGroupsInRow + jend = nGroupsInCol #else iend = nGroups + jend = nGroups #endif !! allocate the groupToGtype and gtypeMaxCutoff here. - if(.not.allocated(groupToGtype)) then - allocate(groupToGtype(iend)) - allocate(groupMaxCutoff(iend)) - allocate(gtypeMaxCutoff(iend)) - groupMaxCutoff = 0.0_dp - gtypeMaxCutoff = 0.0_dp + if(.not.allocated(groupToGtypeRow)) then + ! allocate(groupToGtype(iend)) + allocate(groupToGtypeRow(iend)) + else + deallocate(groupToGtypeRow) + allocate(groupToGtypeRow(iend)) endif + if(.not.allocated(groupMaxCutoffRow)) then + allocate(groupMaxCutoffRow(iend)) + else + deallocate(groupMaxCutoffRow) + allocate(groupMaxCutoffRow(iend)) + end if + + if(.not.allocated(gtypeMaxCutoffRow)) then + allocate(gtypeMaxCutoffRow(iend)) + else + deallocate(gtypeMaxCutoffRow) + allocate(gtypeMaxCutoffRow(iend)) + endif + + +#ifdef IS_MPI + ! We only allocate new storage if we are in MPI because Ncol /= Nrow + if(.not.associated(groupToGtypeCol)) then + allocate(groupToGtypeCol(jend)) + else + deallocate(groupToGtypeCol) + allocate(groupToGtypeCol(jend)) + end if + + if(.not.associated(groupToGtypeCol)) then + allocate(groupToGtypeCol(jend)) + else + deallocate(groupToGtypeCol) + allocate(groupToGtypeCol(jend)) + end if + if(.not.associated(gtypeMaxCutoffCol)) then + allocate(gtypeMaxCutoffCol(jend)) + else + deallocate(gtypeMaxCutoffCol) + allocate(gtypeMaxCutoffCol(jend)) + end if + + groupMaxCutoffCol = 0.0_dp + gtypeMaxCutoffCol = 0.0_dp + +#endif + groupMaxCutoffRow = 0.0_dp + gtypeMaxCutoffRow = 0.0_dp + + !! first we do a single loop over the cutoff groups to find the !! largest cutoff for any atypes present in this group. We also !! create gtypes at this point. tol = 1.0d-6 - + nGroupTypesRow = 0 + do i = istart, iend n_in_i = groupStartRow(i+1) - groupStartRow(i) - groupMaxCutoff(i) = 0.0_dp + groupMaxCutoffRow(i) = 0.0_dp do ia = groupStartRow(i), groupStartRow(i+1)-1 atom1 = groupListRow(ia) #ifdef IS_MPI @@ -359,46 +428,93 @@ contains #else me_i = atid(atom1) #endif - if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then - groupMaxCutoff(i)=atypeMaxCutoff(me_i) + if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then + groupMaxCutoffRow(i)=atypeMaxCutoff(me_i) endif enddo - if (nGroupTypes.eq.0) then - nGroupTypes = nGroupTypes + 1 - gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i) - groupToGtype(i) = nGroupTypes + if (nGroupTypesRow.eq.0) then + nGroupTypesRow = nGroupTypesRow + 1 + gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i) + groupToGtypeRow(i) = nGroupTypesRow else GtypeFound = .false. - do g = 1, nGroupTypes - if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then - groupToGtype(i) = g + do g = 1, nGroupTypesRow + if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then + groupToGtypeRow(i) = g GtypeFound = .true. endif enddo if (.not.GtypeFound) then - nGroupTypes = nGroupTypes + 1 - gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i) - groupToGtype(i) = nGroupTypes + nGroupTypesRow = nGroupTypesRow + 1 + gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i) + groupToGtypeRow(i) = nGroupTypesRow endif endif enddo +#ifdef IS_MPI + do j = jstart, jend + n_in_j = groupStartCol(j+1) - groupStartCol(j) + groupMaxCutoffCol(j) = 0.0_dp + do ja = groupStartCol(j), groupStartCol(j+1)-1 + atom1 = groupListCol(ja) + + me_j = atid_col(atom1) + + if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then + groupMaxCutoffCol(j)=atypeMaxCutoff(me_j) + endif + enddo + + if (nGroupTypesCol.eq.0) then + nGroupTypesCol = nGroupTypesCol + 1 + gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j) + groupToGtypeCol(j) = nGroupTypesCol + else + GtypeFound = .false. + do g = 1, nGroupTypesCol + if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then + groupToGtypeCol(j) = g + GtypeFound = .true. + endif + enddo + if (.not.GtypeFound) then + nGroupTypesCol = nGroupTypesCol + 1 + gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j) + groupToGtypeCol(j) = nGroupTypesCol + endif + endif + enddo + +#else +! Set pointers to information we just found + nGroupTypesCol = nGroupTypesRow + groupToGtypeCol => groupToGtypeRow + gtypeMaxCutoffCol => gtypeMaxCutoffRow + groupMaxCutoffCol => groupMaxCutoffRow +#endif + + + + + !! allocate the gtypeCutoffMap here. - allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes)) + allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol)) !! then we do a double loop over all the group TYPES to find the cutoff !! map between groups of two types - - do i = 1, nGroupTypes - do j = 1, nGroupTypes + tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol)) + + do i = 1, nGroupTypesRow + do j = 1, nGroupTypesCol select case(cutoffPolicy) case(TRADITIONAL_CUTOFF_POLICY) - thisRcut = maxval(gtypeMaxCutoff) + thisRcut = tradRcut case(MIX_CUTOFF_POLICY) - thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j)) + thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j)) case(MAX_CUTOFF_POLICY) - thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j)) + thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j)) case default call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy") return @@ -406,10 +522,27 @@ contains gtypeCutoffMap(i,j)%rcut = thisRcut gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut skin = defaultRlist - defaultRcut + listSkin = skin ! set neighbor list skin thickness gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2 + ! sanity check + + if (haveDefaultCutoffs) then + if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then + call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff") + endif + endif enddo enddo + if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow) + if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow) + if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff) +#ifdef IS_MPI + if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol) + if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol) +#endif + groupMaxCutoffCol => null() + gtypeMaxCutoffCol => null() haveGtypeCutoffMap = .true. end subroutine createGtypeCutoffMap @@ -422,7 +555,8 @@ contains defaultRsw = defRsw defaultRlist = defRlist cutoffPolicy = cutPolicy - rcuti = 1.0_dp / defaultRcut + + haveDefaultCutoffs = .true. end subroutine setDefaultCutoffs subroutine setCutoffPolicy(cutPolicy) @@ -439,7 +573,6 @@ contains SIM_requires_postpair_calc = SimRequiresPostpairCalc() SIM_requires_prepair_calc = SimRequiresPrepairCalc() SIM_uses_PBC = SimUsesPBC() - SIM_uses_RF = SimUsesRF() haveSIMvariables = .true. @@ -506,11 +639,9 @@ contains end subroutine doReadyCheck - subroutine init_FF(use_RF, correctionMethod, dampingAlpha, thisStat) + subroutine init_FF(thisESM, thisStat) - logical, intent(in) :: use_RF - integer, intent(in) :: correctionMethod - real(kind=dp), intent(in) :: dampingAlpha + integer, intent(in) :: thisESM integer, intent(out) :: thisStat integer :: my_status, nMatches integer, pointer :: MatchList(:) => null() @@ -519,10 +650,8 @@ contains !! assume things are copacetic, unless they aren't thisStat = 0 - !! Fortran's version of a cast: - FF_uses_RF = use_RF + electrostaticSummationMethod = thisESM - !! init_FF is called *after* all of the atom types have been !! defined in atype_module using the new_atype subroutine. !! @@ -552,15 +681,15 @@ contains haveSaneForceField = .true. - !! check to make sure the FF_uses_RF setting makes sense + !! check to make sure the reaction field setting makes sense if (FF_uses_Dipoles) then - if (FF_uses_RF) then + if (electrostaticSummationMethod == REACTION_FIELD) then dielect = getDielect() call initialize_rf(dielect) endif else - if ((corrMethod == 3) .or. FF_uses_RF) then + if (electrostaticSummationMethod == REACTION_FIELD) then write(default_error,*) 'Using Reaction Field with no dipoles? Huh?' thisStat = -1 haveSaneForceField = .false. @@ -651,7 +780,7 @@ contains integer :: propPack_i, propPack_j integer :: loopStart, loopEnd, loop integer :: iHash - real(kind=dp) :: listSkin = 1.0 + !! initialize local variables @@ -776,9 +905,9 @@ contains me_j = atid(j) call get_interatomic_vector(q_group(:,i), & q_group(:,j), d_grp, rgrpsq) -#endif +#endif - if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then + if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then if (update_nlist) then nlist = nlist + 1 @@ -899,6 +1028,7 @@ contains endif end if enddo + enddo outer if (update_nlist) then @@ -978,7 +1108,7 @@ contains if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then - if ((FF_uses_RF .and. SIM_uses_RF) .or. (corrMethod == 3)) then + if (electrostaticSummationMethod == REACTION_FIELD) then #ifdef IS_MPI call scatter(rf_Row,rf,plan_atom_row_3d) @@ -1086,9 +1216,9 @@ contains if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, & - pot, eFrame, f, t, do_pot, corrMethod, rcuti) + pot, eFrame, f, t, do_pot) - if ((FF_uses_RF .and. SIM_uses_RF) .or. (corrMethod == 3)) then + if (electrostaticSummationMethod == REACTION_FIELD) then ! CHECK ME (RF needs to know about all electrostatic types) call accumulate_rf(i, j, r, eFrame, sw) @@ -1370,8 +1500,7 @@ contains function FF_RequiresPostpairCalc() result(doesit) logical :: doesit - doesit = FF_uses_RF - if (corrMethod == 3) doesit = .true. + if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true. end function FF_RequiresPostpairCalc #ifdef PROFILE