--- trunk/OOPSE-4/src/UseTheForce/doForces.F90 2005/09/15 22:05:21 2301 +++ trunk/OOPSE-4/src/UseTheForce/doForces.F90 2005/12/30 21:25:56 2532 @@ -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.73 2005-12-30 21:25:56 tim Exp $, $Date: 2005-12-30 21:25:56 $, $Name: not supported by cvs2svn $, $Revision: 1.73 $ module doForces @@ -58,11 +58,11 @@ module doForces use lj use sticky use electrostatic_module - use reaction_field_module - use gb_pair + use gayberne use shapes use vector_class use eam + use suttonchen use status #ifdef IS_MPI use mpiSimulation @@ -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,33 +86,42 @@ module doForces logical, save :: haveSaneForceField = .false. logical, save :: haveInteractionHash = .false. logical, save :: haveGtypeCutoffMap = .false. - logical, save :: haveRlist = .false. + logical, save :: haveDefaultCutoffs = .false. + logical, save :: haveSkinThickness = .false. + logical, save :: haveElectrostaticSummationMethod = .false. + logical, save :: haveCutoffPolicy = .false. + logical, save :: VisitCutoffsAfterComputing = .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 :: FF_uses_SC + logical, save :: FF_uses_MEAM + logical, save :: SIM_uses_DirectionalAtoms logical, save :: SIM_uses_EAM - logical, save :: SIM_uses_RF + logical, save :: SIM_uses_SC + logical, save :: SIM_uses_MEAM logical, save :: SIM_requires_postpair_calc logical, save :: SIM_requires_prepair_calc logical, save :: SIM_uses_PBC - integer, save :: corrMethod + integer, save :: electrostaticSummationMethod + integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY + real(kind=dp), save :: defaultRcut, defaultRsw, largestRcut + real(kind=dp), save :: skinThickness + logical, save :: defaultDoShift + public :: init_FF - public :: setDefaultCutoffs + public :: setCutoffs + public :: cWasLame + public :: setElectrostaticMethod + public :: setCutoffPolicy + public :: setSkinThickness public :: do_force_loop - public :: createInteractionHash - public :: createGtypeCutoffMap - public :: getStickyCut - public :: getStickyPowerCut - public :: getGayBerneCut - public :: getEAMCut - public :: getShapeCut #ifdef PROFILE public :: getforcetime @@ -124,9 +134,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 @@ -134,15 +149,10 @@ module doForces end type gtypeCutoffs type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap - integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY - real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist - real(kind=dp),save :: rcuti - contains - subroutine createInteractionHash(status) + subroutine createInteractionHash() integer :: nAtypes - integer, intent(out) :: status integer :: i integer :: j integer :: iHash @@ -154,6 +164,8 @@ contains logical :: i_is_GB logical :: i_is_EAM logical :: i_is_Shape + logical :: i_is_SC + logical :: i_is_MEAM logical :: j_is_LJ logical :: j_is_Elect logical :: j_is_Sticky @@ -161,29 +173,34 @@ contains logical :: j_is_GB logical :: j_is_EAM logical :: j_is_Shape + logical :: j_is_SC + logical :: j_is_MEAM real(kind=dp) :: myRcut - status = 0 - if (.not. associated(atypes)) then - call handleError("atype", "atypes was not present before call of createInteractionHash!") - status = -1 + call handleError("doForces", "atypes was not present before call of createInteractionHash!") return endif nAtypes = getSize(atypes) if (nAtypes == 0) then - status = -1 + call handleError("doForces", "nAtypes was zero during call of createInteractionHash!") return end if 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 @@ -194,6 +211,8 @@ contains call getElementProperty(atypes, i, "is_GayBerne", i_is_GB) call getElementProperty(atypes, i, "is_EAM", i_is_EAM) call getElementProperty(atypes, i, "is_Shape", i_is_Shape) + call getElementProperty(atypes, i, "is_SC", i_is_SC) + call getElementProperty(atypes, i, "is_MEAM", i_is_MEAM) do j = i, nAtypes @@ -207,6 +226,8 @@ contains call getElementProperty(atypes, j, "is_GayBerne", j_is_GB) call getElementProperty(atypes, j, "is_EAM", j_is_EAM) call getElementProperty(atypes, j, "is_Shape", j_is_Shape) + call getElementProperty(atypes, j, "is_SC", j_is_SC) + call getElementProperty(atypes, j, "is_MEAM", j_is_MEAM) if (i_is_LJ .and. j_is_LJ) then iHash = ior(iHash, LJ_PAIR) @@ -228,6 +249,10 @@ contains iHash = ior(iHash, EAM_PAIR) endif + if (i_is_SC .and. j_is_SC) then + iHash = ior(iHash, SC_PAIR) + endif + if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR) if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ) if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ) @@ -247,9 +272,8 @@ contains haveInteractionHash = .true. end subroutine createInteractionHash - subroutine createGtypeCutoffMap(stat) + subroutine createGtypeCutoffMap() - integer, intent(out), optional :: stat logical :: i_is_LJ logical :: i_is_Elect logical :: i_is_Sticky @@ -257,25 +281,23 @@ contains logical :: i_is_GB logical :: i_is_EAM logical :: i_is_Shape + logical :: i_is_SC 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 real(kind=dp) :: biggestAtypeCutoff - stat = 0 if (.not. haveInteractionHash) then - call createInteractionHash(myStatus) - if (myStatus .ne. 0) then - write(default_error, *) 'createInteractionHash failed in doForces!' - stat = -1 - return - endif + call createInteractionHash() 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. @@ -289,69 +311,125 @@ contains call getElementProperty(atypes, i, "is_GayBerne", i_is_GB) call getElementProperty(atypes, i, "is_EAM", i_is_EAM) call getElementProperty(atypes, i, "is_Shape", i_is_Shape) - + call getElementProperty(atypes, i, "is_SC", i_is_SC) - 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 + 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 + if (i_is_SC) then + thisRcut = getSCCut(i) + if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut + endif endif - if (i_is_Shape) then - thisRcut = getShapeCut(i) - if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut - 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(groupMaxCutoffCol)) then + allocate(groupMaxCutoffCol(jend)) + else + deallocate(groupMaxCutoffCol) + allocate(groupMaxCutoffCol(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 + nGroupTypesCol = 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,93 +437,212 @@ 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 (nGroupTypesRow.eq.0) then + nGroupTypesRow = nGroupTypesRow + 1 + gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i) + groupToGtypeRow(i) = nGroupTypesRow + else + GtypeFound = .false. + do g = 1, nGroupTypesRow + if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then + groupToGtypeRow(i) = g + GtypeFound = .true. + endif + enddo + if (.not.GtypeFound) then + 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 (nGroupTypes.eq.0) then - nGroupTypes = nGroupTypes + 1 - gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i) - groupToGtype(i) = nGroupTypes + if (nGroupTypesCol.eq.0) then + nGroupTypesCol = nGroupTypesCol + 1 + gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j) + groupToGtypeCol(j) = nGroupTypesCol else GtypeFound = .false. - do g = 1, nGroupTypes - if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then - groupToGtype(i) = g + do g = 1, nGroupTypesCol + if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then + groupToGtypeCol(j) = g GtypeFound = .true. endif enddo if (.not.GtypeFound) then - nGroupTypes = nGroupTypes + 1 - gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i) - groupToGtype(i) = nGroupTypes + 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 end select gtypeCutoffMap(i,j)%rcut = thisRcut + + if (thisRcut.gt.largestRcut) largestRcut = thisRcut + gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut - skin = defaultRlist - defaultRcut - gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2 + if (.not.haveSkinThickness) then + skinThickness = 1.0_dp + endif + + gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**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 - subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy) - real(kind=dp),intent(in) :: defRcut, defRsw, defRlist - integer, intent(in) :: cutPolicy + subroutine setCutoffs(defRcut, defRsw) + real(kind=dp),intent(in) :: defRcut, defRsw + character(len = statusMsgSize) :: errMsg + integer :: localError + defaultRcut = defRcut defaultRsw = defRsw - defaultRlist = defRlist - cutoffPolicy = cutPolicy - rcuti = 1.0_dp / defaultRcut - end subroutine setDefaultCutoffs + + defaultDoShift = .false. + if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then + + write(errMsg, *) & + 'cutoffRadius and switchingRadius are set to the same', newline & + // tab, 'value. OOPSE will use shifted ', newline & + // tab, 'potentials instead of switching functions.' + + call handleInfo("setCutoffs", errMsg) + + defaultDoShift = .true. + + endif - subroutine setCutoffPolicy(cutPolicy) + localError = 0 + call setLJDefaultCutoff( defaultRcut, defaultDoShift ) + call setElectrostaticCutoffRadius( defaultRcut, defaultRsw ) + call setCutoffEAM( defaultRcut, localError) + if (localError /= 0) then + write(errMsg, *) 'An error has occured in setting the EAM cutoff' + call handleError("setCutoffs", errMsg) + end if + call set_switch(GROUP_SWITCH, defaultRsw, defaultRcut) + + haveDefaultCutoffs = .true. + haveGtypeCutoffMap = .false. + end subroutine setCutoffs + subroutine cWasLame() + + VisitCutoffsAfterComputing = .true. + return + + end subroutine cWasLame + + subroutine setCutoffPolicy(cutPolicy) + integer, intent(in) :: cutPolicy + cutoffPolicy = cutPolicy - call createGtypeCutoffMap() - end subroutine setCutoffPolicy - + haveCutoffPolicy = .true. + haveGtypeCutoffMap = .false. - subroutine setSimVariables() - SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms() - SIM_uses_EAM = SimUsesEAM() - SIM_requires_postpair_calc = SimRequiresPostpairCalc() - SIM_requires_prepair_calc = SimRequiresPrepairCalc() - SIM_uses_PBC = SimUsesPBC() - SIM_uses_RF = SimUsesRF() + end subroutine setCutoffPolicy + + subroutine setElectrostaticMethod( thisESM ) - haveSIMvariables = .true. + integer, intent(in) :: thisESM - return - end subroutine setSimVariables + electrostaticSummationMethod = thisESM + haveElectrostaticSummationMethod = .true. + + end subroutine setElectrostaticMethod + subroutine setSkinThickness( thisSkin ) + + real(kind=dp), intent(in) :: thisSkin + + skinThickness = thisSkin + haveSkinThickness = .true. + haveGtypeCutoffMap = .false. + + end subroutine setSkinThickness + + subroutine setSimVariables() + SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms() + SIM_uses_EAM = SimUsesEAM() + SIM_requires_postpair_calc = SimRequiresPostpairCalc() + SIM_requires_prepair_calc = SimRequiresPrepairCalc() + SIM_uses_PBC = SimUsesPBC() + + haveSIMvariables = .true. + + return + end subroutine setSimVariables + subroutine doReadyCheck(error) integer, intent(out) :: error @@ -454,25 +651,19 @@ contains error = 0 if (.not. haveInteractionHash) then - myStatus = 0 - call createInteractionHash(myStatus) - if (myStatus .ne. 0) then - write(default_error, *) 'createInteractionHash failed in doForces!' - error = -1 - return - endif + call createInteractionHash() endif if (.not. haveGtypeCutoffMap) then - myStatus = 0 - call createGtypeCutoffMap(myStatus) - if (myStatus .ne. 0) then - write(default_error, *) 'createGtypeCutoffMap failed in doForces!' - error = -1 - return - endif + call createGtypeCutoffMap() endif + + if (VisitCutoffsAfterComputing) then + call set_switch(GROUP_SWITCH, largestRcut, largestRcut) + endif + + if (.not. haveSIMvariables) then call setSimVariables() endif @@ -506,23 +697,15 @@ contains end subroutine doReadyCheck - subroutine init_FF(use_RF, correctionMethod, dampingAlpha, thisStat) + subroutine init_FF(thisStat) - logical, intent(in) :: use_RF - integer, intent(in) :: correctionMethod - real(kind=dp), intent(in) :: dampingAlpha integer, intent(out) :: thisStat integer :: my_status, nMatches integer, pointer :: MatchList(:) => null() - real(kind=dp) :: rcut, rrf, rt, dielect !! assume things are copacetic, unless they aren't thisStat = 0 - !! Fortran's version of a cast: - FF_uses_RF = use_RF - - !! init_FF is called *after* all of the atom types have been !! defined in atype_module using the new_atype subroutine. !! @@ -552,22 +735,6 @@ contains haveSaneForceField = .true. - !! check to make sure the FF_uses_RF setting makes sense - - if (FF_uses_Dipoles) then - if (FF_uses_RF) then - dielect = getDielect() - call initialize_rf(dielect) - endif - else - if ((corrMethod == 3) .or. FF_uses_RF) then - write(default_error,*) 'Using Reaction Field with no dipoles? Huh?' - thisStat = -1 - haveSaneForceField = .false. - return - endif - endif - if (FF_uses_EAM) then call init_EAM_FF(my_status) if (my_status /= 0) then @@ -578,15 +745,6 @@ contains end if endif - if (FF_uses_GayBerne) then - call check_gb_pair_FF(my_status) - if (my_status .ne. 0) then - thisStat = -1 - haveSaneForceField = .false. - return - endif - endif - if (.not. haveNeighborList) then !! Create neighbor lists call expandNeighborList(nLocal, my_status) @@ -620,13 +778,13 @@ contains !! Stress Tensor real( kind = dp), dimension(9) :: tau - real ( kind = dp ) :: pot + real ( kind = dp ),dimension(LR_POT_TYPES) :: pot logical ( kind = 2) :: do_pot_c, do_stress_c logical :: do_pot logical :: do_stress logical :: in_switching_region #ifdef IS_MPI - real( kind = DP ) :: pot_local + real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local integer :: nAtomsInRow integer :: nAtomsInCol integer :: nprocs @@ -641,8 +799,10 @@ contains integer :: nlist real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij real( kind = DP ) :: sw, dswdr, swderiv, mf + real( kind = DP ) :: rVal real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij real(kind=dp) :: rfpot, mu_i, virial + real(kind=dp):: rCut integer :: me_i, me_j, n_in_i, n_in_j logical :: is_dp_i integer :: neighborListSize @@ -651,7 +811,8 @@ contains integer :: propPack_i, propPack_j integer :: loopStart, loopEnd, loop integer :: iHash - real(kind=dp) :: listSkin = 1.0 + integer :: i1 + !! initialize local variables @@ -715,10 +876,10 @@ contains ! (but only on the first time through): if (loop .eq. loopStart) then #ifdef IS_MPI - call checkNeighborList(nGroupsInRow, q_group_row, listSkin, & + call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, & update_nlist) #else - call checkNeighborList(nGroups, q_group, listSkin, & + call checkNeighborList(nGroups, q_group, skinThickness, & update_nlist) #endif endif @@ -776,9 +937,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 @@ -798,107 +959,113 @@ contains list(nlist) = j endif + - if (loop .eq. PAIR_LOOP) then - vij = 0.0d0 - fij(1:3) = 0.0d0 - endif + + if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then - call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, & - in_switching_region) - - n_in_j = groupStartCol(j+1) - groupStartCol(j) - - do ia = groupStartRow(i), groupStartRow(i+1)-1 - - atom1 = groupListRow(ia) - - inner: do jb = groupStartCol(j), groupStartCol(j+1)-1 - - atom2 = groupListCol(jb) - - if (skipThisPair(atom1, atom2)) cycle inner - - if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then - d_atm(1:3) = d_grp(1:3) - ratmsq = rgrpsq - else + rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut + if (loop .eq. PAIR_LOOP) then + vij = 0.0d0 + fij(1:3) = 0.0d0 + endif + + call get_switch(rgrpsq, sw, dswdr, rgrp, & + group_switch, in_switching_region) + + n_in_j = groupStartCol(j+1) - groupStartCol(j) + + do ia = groupStartRow(i), groupStartRow(i+1)-1 + + atom1 = groupListRow(ia) + + inner: do jb = groupStartCol(j), groupStartCol(j+1)-1 + + atom2 = groupListCol(jb) + + if (skipThisPair(atom1, atom2)) cycle inner + + if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then + d_atm(1:3) = d_grp(1:3) + ratmsq = rgrpsq + else #ifdef IS_MPI - call get_interatomic_vector(q_Row(:,atom1), & - q_Col(:,atom2), d_atm, ratmsq) + call get_interatomic_vector(q_Row(:,atom1), & + q_Col(:,atom2), d_atm, ratmsq) #else - call get_interatomic_vector(q(:,atom1), & - q(:,atom2), d_atm, ratmsq) + call get_interatomic_vector(q(:,atom1), & + q(:,atom2), d_atm, ratmsq) #endif - endif - - if (loop .eq. PREPAIR_LOOP) then + endif + + if (loop .eq. PREPAIR_LOOP) then #ifdef IS_MPI - call do_prepair(atom1, atom2, ratmsq, d_atm, sw, & - rgrpsq, d_grp, do_pot, do_stress, & - eFrame, A, f, t, pot_local) + call do_prepair(atom1, atom2, ratmsq, d_atm, sw, & + rgrpsq, d_grp, rCut, do_pot, do_stress, & + eFrame, A, f, t, pot_local) #else - call do_prepair(atom1, atom2, ratmsq, d_atm, sw, & - rgrpsq, d_grp, do_pot, do_stress, & - eFrame, A, f, t, pot) + call do_prepair(atom1, atom2, ratmsq, d_atm, sw, & + rgrpsq, d_grp, rCut, do_pot, do_stress, & + eFrame, A, f, t, pot) #endif - else + else #ifdef IS_MPI - call do_pair(atom1, atom2, ratmsq, d_atm, sw, & - do_pot, & - eFrame, A, f, t, pot_local, vpair, fpair) + call do_pair(atom1, atom2, ratmsq, d_atm, sw, & + do_pot, eFrame, A, f, t, pot_local, vpair, & + fpair, d_grp, rgrp, rCut) #else - call do_pair(atom1, atom2, ratmsq, d_atm, sw, & - do_pot, & - eFrame, A, f, t, pot, vpair, fpair) + call do_pair(atom1, atom2, ratmsq, d_atm, sw, & + do_pot, eFrame, A, f, t, pot, vpair, fpair, & + d_grp, rgrp, rCut) #endif + vij = vij + vpair + fij(1:3) = fij(1:3) + fpair(1:3) + endif + enddo inner + enddo - vij = vij + vpair - fij(1:3) = fij(1:3) + fpair(1:3) - endif - enddo inner - enddo - - if (loop .eq. PAIR_LOOP) then - if (in_switching_region) then - swderiv = vij*dswdr/rgrp - fij(1) = fij(1) + swderiv*d_grp(1) - fij(2) = fij(2) + swderiv*d_grp(2) - fij(3) = fij(3) + swderiv*d_grp(3) - - do ia=groupStartRow(i), groupStartRow(i+1)-1 - atom1=groupListRow(ia) - mf = mfactRow(atom1) + if (loop .eq. PAIR_LOOP) then + if (in_switching_region) then + swderiv = vij*dswdr/rgrp + fij(1) = fij(1) + swderiv*d_grp(1) + fij(2) = fij(2) + swderiv*d_grp(2) + fij(3) = fij(3) + swderiv*d_grp(3) + + do ia=groupStartRow(i), groupStartRow(i+1)-1 + atom1=groupListRow(ia) + mf = mfactRow(atom1) #ifdef IS_MPI - f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf - f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf - f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf + f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf + f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf + f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf #else - f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf - f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf - f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf + f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf + f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf + f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf #endif - enddo - - do jb=groupStartCol(j), groupStartCol(j+1)-1 - atom2=groupListCol(jb) - mf = mfactCol(atom2) + enddo + + do jb=groupStartCol(j), groupStartCol(j+1)-1 + atom2=groupListCol(jb) + mf = mfactCol(atom2) #ifdef IS_MPI - f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf - f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf - f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf + f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf + f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf + f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf #else - f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf - f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf - f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf + f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf + f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf + f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf #endif - enddo - endif + enddo + endif - if (do_stress) call add_stress_tensor(d_grp, fij) + if (do_stress) call add_stress_tensor(d_grp, fij) + endif endif - end if + endif enddo + enddo outer if (update_nlist) then @@ -958,98 +1125,109 @@ contains if (do_pot) then ! scatter/gather pot_row into the members of my column - call scatter(pot_Row, pot_Temp, plan_atom_row) - + do i = 1,LR_POT_TYPES + call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row) + end do ! scatter/gather pot_local into all other procs ! add resultant to get total pot do i = 1, nlocal - pot_local = pot_local + pot_Temp(i) + pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) & + + pot_Temp(1:LR_POT_TYPES,i) enddo pot_Temp = 0.0_DP - - call scatter(pot_Col, pot_Temp, plan_atom_col) + do i = 1,LR_POT_TYPES + call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col) + end do do i = 1, nlocal - pot_local = pot_local + pot_Temp(i) + pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)& + + pot_Temp(1:LR_POT_TYPES,i) enddo endif #endif - if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then + if (SIM_requires_postpair_calc) then + do i = 1, nlocal + + ! we loop only over the local atoms, so we don't need row and column + ! lookups for the types + + me_i = atid(i) + + ! is the atom electrostatic? See if it would have an + ! electrostatic interaction with itself + iHash = InteractionHash(me_i,me_i) - if ((FF_uses_RF .and. SIM_uses_RF) .or. (corrMethod == 3)) then - + if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then #ifdef IS_MPI - call scatter(rf_Row,rf,plan_atom_row_3d) - call scatter(rf_Col,rf_Temp,plan_atom_col_3d) - do i = 1,nlocal - rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i) - end do -#endif - - do i = 1, nLocal - - rfpot = 0.0_DP -#ifdef IS_MPI - me_i = atid_row(i) + call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), & + t, do_pot) #else - me_i = atid(i) + call self_self(i, eFrame, pot(ELECTROSTATIC_POT), & + t, do_pot) #endif - iHash = InteractionHash(me_i,me_j) + endif + + + if (electrostaticSummationMethod.eq.REACTION_FIELD) then - if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then - - mu_i = getDipoleMoment(me_i) - - !! The reaction field needs to include a self contribution - !! to the field: - call accumulate_self_rf(i, mu_i, eFrame) - !! Get the reaction field contribution to the - !! potential and torques: - call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot) + ! loop over the excludes to accumulate RF stuff we've + ! left out of the normal pair loop + + do i1 = 1, nSkipsForAtom(i) + j = skipsForAtom(i, i1) + + ! prevent overcounting of the skips + if (i.lt.j) then + call get_interatomic_vector(q(:,i), & + q(:,j), d_atm, ratmsq) + rVal = dsqrt(ratmsq) + call get_switch(ratmsq, sw, dswdr, rVal, group_switch, & + in_switching_region) #ifdef IS_MPI - pot_local = pot_local + rfpot + call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, & + vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot) #else - pot = pot + rfpot - + call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, & + vpair, pot(ELECTROSTATIC_POT), f, t, do_pot) #endif - endif - enddo - endif + endif + enddo + endif + enddo endif - - + #ifdef IS_MPI - + if (do_pot) then - pot = pot + pot_local - !! we assume the c code will do the allreduce to get the total potential - !! we could do it right here if we needed to... + call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, & + mpi_comm_world,mpi_err) endif - + if (do_stress) then call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, & mpi_comm_world,mpi_err) call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, & mpi_comm_world,mpi_err) endif - + #else - + if (do_stress) then tau = tau_Temp virial = virial_Temp endif - + #endif - + end subroutine do_force_loop subroutine do_pair(i, j, rijsq, d, sw, do_pot, & - eFrame, A, f, t, pot, vpair, fpair) + eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut) - real( kind = dp ) :: pot, vpair, sw + real( kind = dp ) :: vpair, sw + real( kind = dp ), dimension(LR_POT_TYPES) :: pot real( kind = dp ), dimension(3) :: fpair real( kind = dp ), dimension(nLocal) :: mfact real( kind = dp ), dimension(9,nLocal) :: eFrame @@ -1060,8 +1238,11 @@ contains logical, intent(inout) :: do_pot integer, intent(in) :: i, j real ( kind = dp ), intent(inout) :: rijsq - real ( kind = dp ) :: r + real ( kind = dp ), intent(inout) :: r_grp real ( kind = dp ), intent(inout) :: d(3) + real ( kind = dp ), intent(inout) :: d_grp(3) + real ( kind = dp ), intent(inout) :: rCut + real ( kind = dp ) :: r integer :: me_i, me_j integer :: iHash @@ -1079,65 +1260,66 @@ contains #endif iHash = InteractionHash(me_i, me_j) - + if ( iand(iHash, LJ_PAIR).ne.0 ) then - call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot) + call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, & + pot(VDW_POT), f, do_pot) endif - + 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) - - if ((FF_uses_RF .and. SIM_uses_RF) .or. (corrMethod == 3)) then - - ! CHECK ME (RF needs to know about all electrostatic types) - call accumulate_rf(i, j, r, eFrame, sw) - call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair) - endif - + call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, & + pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot) endif - + if ( iand(iHash, STICKY_PAIR).ne.0 ) then call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, & - pot, A, f, t, do_pot) + pot(HB_POT), A, f, t, do_pot) endif - + if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, & - pot, A, f, t, do_pot) + pot(HB_POT), A, f, t, do_pot) endif - + if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, & - pot, A, f, t, do_pot) + pot(VDW_POT), A, f, t, do_pot) endif if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then - ! call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, & - ! pot, A, f, t, do_pot) + call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, & + pot(VDW_POT), A, f, t, do_pot) endif - + if ( iand(iHash, EAM_PAIR).ne.0 ) then - call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, & - do_pot) + call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, & + pot(METALLIC_POT), f, do_pot) endif - + if ( iand(iHash, SHAPE_PAIR).ne.0 ) then call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, & - pot, A, f, t, do_pot) + pot(VDW_POT), A, f, t, do_pot) endif - + if ( iand(iHash, SHAPE_LJ).ne.0 ) then call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, & - pot, A, f, t, do_pot) + pot(VDW_POT), A, f, t, do_pot) endif + + if ( iand(iHash, SC_PAIR).ne.0 ) then + call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, & + pot(METALLIC_POT), f, do_pot) + endif + + end subroutine do_pair - subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, & + subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, & do_pot, do_stress, eFrame, A, f, t, pot) - real( kind = dp ) :: pot, sw + real( kind = dp ) :: sw + real( kind = dp ), dimension(LR_POT_TYPES) :: pot real( kind = dp ), dimension(9,nLocal) :: eFrame real (kind=dp), dimension(9,nLocal) :: A real (kind=dp), dimension(3,nLocal) :: f @@ -1145,7 +1327,7 @@ contains logical, intent(inout) :: do_pot, do_stress integer, intent(in) :: i, j - real ( kind = dp ), intent(inout) :: rijsq, rcijsq + real ( kind = dp ), intent(inout) :: rijsq, rcijsq, rCut real ( kind = dp ) :: r, rc real ( kind = dp ), intent(inout) :: d(3), dc(3) @@ -1164,18 +1346,25 @@ contains iHash = InteractionHash(me_i, me_j) if ( iand(iHash, EAM_PAIR).ne.0 ) then - call calc_EAM_prepair_rho(i, j, d, r, rijsq ) + call calc_EAM_prepair_rho(i, j, d, r, rijsq) endif + + if ( iand(iHash, SC_PAIR).ne.0 ) then + call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut ) + endif end subroutine do_prepair subroutine do_preforce(nlocal,pot) integer :: nlocal - real( kind = dp ) :: pot + real( kind = dp ),dimension(LR_POT_TYPES) :: pot if (FF_uses_EAM .and. SIM_uses_EAM) then - call calc_EAM_preforce_Frho(nlocal,pot) + call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT)) + endif + if (FF_uses_SC .and. SIM_uses_SC) then + call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT)) endif @@ -1261,17 +1450,12 @@ contains pot_Col = 0.0_dp pot_Temp = 0.0_dp - rf_Row = 0.0_dp - rf_Col = 0.0_dp - rf_Temp = 0.0_dp - #endif if (FF_uses_EAM .and. SIM_uses_EAM) then call clean_EAM() endif - rf = 0.0_dp tau_Temp = 0.0_dp virial_Temp = 0.0_dp end subroutine zero_work_arrays @@ -1365,15 +1549,10 @@ contains function FF_RequiresPrepairCalc() result(doesit) logical :: doesit - doesit = FF_uses_EAM + doesit = FF_uses_EAM .or. FF_uses_SC & + .or. FF_uses_MEAM end function FF_RequiresPrepairCalc - function FF_RequiresPostpairCalc() result(doesit) - logical :: doesit - doesit = FF_uses_RF - if (corrMethod == 3) doesit = .true. - end function FF_RequiresPostpairCalc - #ifdef PROFILE function getforcetime() result(totalforcetime) real(kind=dp) :: totalforcetime