--- trunk/OOPSE-4/src/UseTheForce/doForces.F90 2005/06/27 21:01:36 2259 +++ trunk/OOPSE-4/src/UseTheForce/doForces.F90 2005/08/09 19:40:56 2269 @@ -45,7 +45,7 @@ !! @author Charles F. Vardeman II !! @author Matthew Meineke -!! @version $Id: doForces.F90,v 1.20 2005-06-27 21:01:30 gezelter Exp $, $Date: 2005-06-27 21:01:30 $, $Name: not supported by cvs2svn $, $Revision: 1.20 $ +!! @version $Id: doForces.F90,v 1.27 2005-08-09 19:40:56 chuckv Exp $, $Date: 2005-08-09 19:40:56 $, $Name: not supported by cvs2svn $, $Revision: 1.27 $ module doForces @@ -81,8 +81,8 @@ module doForces logical, save :: haveRlist = .false. logical, save :: haveNeighborList = .false. logical, save :: haveSIMvariables = .false. - logical, save :: havePropertyMap = .false. logical, save :: haveSaneForceField = .false. + logical, save :: haveInteractionMap = .false. logical, save :: FF_uses_DirectionalAtoms logical, save :: FF_uses_LennardJones @@ -116,11 +116,15 @@ module doForces logical, save :: SIM_uses_PBC logical, save :: SIM_uses_molecular_cutoffs - real(kind=dp), save :: rlist, rlistsq public :: init_FF public :: do_force_loop - public :: setRlistDF +! public :: setRlistDF + !public :: addInteraction + !public :: setInteractionHash + !public :: getInteractionHash + public :: createInteractionMap + public :: createGroupCutoffs #ifdef PROFILE public :: getforcetime @@ -128,99 +132,239 @@ module doForces real :: forceTimeInitial, forceTimeFinal integer :: nLoops #endif + +!! Variables for cutoff mapping and interaction mapping +! Bit hash to determine pair-pair interactions. + integer, dimension(:,:),allocatable :: InteractionHash +!! Cuttoffs in OOPSE are handled on a Group-Group pair basis. +! Largest cutoff for atypes for all potentials + real(kind=dp), dimension(:), allocatable :: atypeMaxCuttoff +! Largest cutoff for groups + real(kind=dp), dimension(:), allocatable :: groupMaxCutoff +! Group to Gtype transformation Map + integer,dimension(:), allocatable :: groupToGtype +! Group Type Max Cutoff + real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff +! GroupType definition + type ::gtype + real(kind=dp) :: rcut ! Group Cutoff + real(kind=dp) :: rcutsq ! Group Cutoff Squared + real(kind=dp) :: rlistsq ! List cutoff Squared + end type gtype - type :: Properties - logical :: is_Directional = .false. - logical :: is_LennardJones = .false. - logical :: is_Electrostatic = .false. - logical :: is_Charge = .false. - logical :: is_Dipole = .false. - logical :: is_Quadrupole = .false. - logical :: is_Sticky = .false. - logical :: is_StickyPower = .false. - logical :: is_GayBerne = .false. - logical :: is_EAM = .false. - logical :: is_Shape = .false. - logical :: is_FLARB = .false. - end type Properties - - type(Properties), dimension(:),allocatable :: PropertyMap - + type(gtype), dimension(:,:), allocatable :: gtypeCutoffMap + contains - subroutine setRlistDF( this_rlist ) - - real(kind=dp) :: this_rlist - - rlist = this_rlist - rlistsq = rlist * rlist - - haveRlist = .true. - - end subroutine setRlistDF - - subroutine createPropertyMap(status) + + subroutine createInteractionMap(status) integer :: nAtypes - integer :: status + integer, intent(out) :: status integer :: i - logical :: thisProperty - real (kind=DP) :: thisDPproperty + integer :: j + integer :: ihash + real(kind=dp) :: myRcut + !! Test Types + logical :: i_is_LJ + logical :: i_is_Elect + logical :: i_is_Sticky + logical :: i_is_StickyP + logical :: i_is_GB + logical :: i_is_EAM + logical :: i_is_Shape + logical :: j_is_LJ + logical :: j_is_Elect + logical :: j_is_Sticky + logical :: j_is_StickyP + logical :: j_is_GB + logical :: j_is_EAM + logical :: j_is_Shape + + status = 0 - status = 0 - + if (.not. associated(atypes)) then + call handleError("atype", "atypes was not present before call of createDefaultInteractionHash!") + status = -1 + return + endif + nAtypes = getSize(atypes) - + if (nAtypes == 0) then status = -1 return end if - if (.not. allocated(PropertyMap)) then - allocate(PropertyMap(nAtypes)) + if (.not. allocated(InteractionHash)) then + allocate(InteractionHash(nAtypes,nAtypes)) endif - + do i = 1, nAtypes - call getElementProperty(atypes, i, "is_Directional", thisProperty) - PropertyMap(i)%is_Directional = thisProperty + call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ) + call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect) + call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky) + call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP) + 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_LennardJones", thisProperty) - PropertyMap(i)%is_LennardJones = thisProperty + do j = i, nAtypes - call getElementProperty(atypes, i, "is_Electrostatic", thisProperty) - PropertyMap(i)%is_Electrostatic = thisProperty + iHash = 0 + myRcut = 0.0_dp - call getElementProperty(atypes, i, "is_Charge", thisProperty) - PropertyMap(i)%is_Charge = thisProperty + call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ) + call getElementProperty(atypes, j, "is_Electrostatic", j_is_Elect) + call getElementProperty(atypes, j, "is_Sticky", j_is_Sticky) + call getElementProperty(atypes, j, "is_StickyPower", j_is_StickyP) + 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, i, "is_Dipole", thisProperty) - PropertyMap(i)%is_Dipole = thisProperty + if (i_is_LJ .and. j_is_LJ) then + iHash = ior(iHash, LJ_PAIR) + endif + + if (i_is_Elect .and. j_is_Elect) then + iHash = ior(iHash, ELECTROSTATIC_PAIR) + endif + + if (i_is_Sticky .and. j_is_Sticky) then + iHash = ior(iHash, STICKY_PAIR) + endif - call getElementProperty(atypes, i, "is_Quadrupole", thisProperty) - PropertyMap(i)%is_Quadrupole = thisProperty + if (i_is_StickyP .and. j_is_StickyP) then + iHash = ior(iHash, STICKYPOWER_PAIR) + endif - call getElementProperty(atypes, i, "is_Sticky", thisProperty) - PropertyMap(i)%is_Sticky = thisProperty - - call getElementProperty(atypes, i, "is_StickyPower", thisProperty) - PropertyMap(i)%is_StickyPower = thisProperty + if (i_is_EAM .and. j_is_EAM) then + iHash = ior(iHash, EAM_PAIR) + endif - call getElementProperty(atypes, i, "is_GayBerne", thisProperty) - PropertyMap(i)%is_GayBerne = thisProperty + 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) - call getElementProperty(atypes, i, "is_EAM", thisProperty) - PropertyMap(i)%is_EAM = thisProperty + if (i_is_Shape .and. j_is_Shape) iHash = ior(iHash, SHAPE_PAIR) + if (i_is_Shape .and. j_is_LJ) iHash = ior(iHash, SHAPE_LJ) + if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ) - call getElementProperty(atypes, i, "is_Shape", thisProperty) - PropertyMap(i)%is_Shape = thisProperty + + InteractionHash(i,j) = iHash + InteractionHash(j,i) = iHash - call getElementProperty(atypes, i, "is_FLARB", thisProperty) - PropertyMap(i)%is_FLARB = thisProperty + end do + end do - havePropertyMap = .true. + haveInteractionMap = .true. + end subroutine createInteractionMap - end subroutine createPropertyMap + subroutine createGroupCutoffs(skinThickness,defaultrList,stat) + real(kind=dp), intent(in), optional :: defaultRList + real(kind-dp), intent(in), :: skinThickenss + ! Query each potential and return the cutoff for that potential. We + ! build the neighbor list based on the largest cutoff value for that + ! atype. Each potential can decide whether to calculate the force for + ! that atype based upon it's own cutoff. + + + real(kind=dp), intent(in), optional :: defaultRCut, defaultSkinThickness + + integer :: iMap + integer :: map_i,map_j + real(kind=dp) :: thisRCut = 0.0_dp + real(kind=dp) :: actualCutoff = 0.0_dp + integer, intent(out) :: stat + integer :: nAtypes + integer :: myStatus + + stat = 0 + if (.not. haveInteractionMap) then + + call createInteractionMap(myStatus) + + if (myStatus .ne. 0) then + write(default_error, *) 'createInteractionMap failed in doForces!' + stat = -1 + return + endif + endif + + nAtypes = getSize(atypes) + !! If we pass a default rcut, set all atypes to that cutoff distance + if(present(defaultRList)) then + InteractionMap(:,:)%rCut = defaultRCut + InteractionMap(:,:)%rCutSq = defaultRCut*defaultRCut + InteractionMap(:,:)%rListSq = (defaultRCut+defaultSkinThickness)**2 + haveRlist = .true. + return + end if + + do map_i = 1,nAtypes + do map_j = map_i,nAtypes + iMap = InteractionMap(map_i, map_j)%InteractionHash + + if ( iand(iMap, LJ_PAIR).ne.0 ) then + ! thisRCut = getLJCutOff(map_i,map_j) + if (thisRcut > actualCutoff) actualCutoff = thisRcut + endif + + if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then + ! thisRCut = getElectrostaticCutOff(map_i,map_j) + if (thisRcut > actualCutoff) actualCutoff = thisRcut + endif + + if ( iand(iMap, STICKY_PAIR).ne.0 ) then + ! thisRCut = getStickyCutOff(map_i,map_j) + if (thisRcut > actualCutoff) actualCutoff = thisRcut + endif + + if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then + ! thisRCut = getStickyPowerCutOff(map_i,map_j) + if (thisRcut > actualCutoff) actualCutoff = thisRcut + endif + + if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then + ! thisRCut = getGayberneCutOff(map_i,map_j) + if (thisRcut > actualCutoff) actualCutoff = thisRcut + endif + + if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then +! thisRCut = getGaybrneLJCutOff(map_i,map_j) + if (thisRcut > actualCutoff) actualCutoff = thisRcut + endif + + if ( iand(iMap, EAM_PAIR).ne.0 ) then +! thisRCut = getEAMCutOff(map_i,map_j) + if (thisRcut > actualCutoff) actualCutoff = thisRcut + endif + + if ( iand(iMap, SHAPE_PAIR).ne.0 ) then +! thisRCut = getShapeCutOff(map_i,map_j) + if (thisRcut > actualCutoff) actualCutoff = thisRcut + endif + + if ( iand(iMap, SHAPE_LJ).ne.0 ) then +! thisRCut = getShapeLJCutOff(map_i,map_j) + if (thisRcut > actualCutoff) actualCutoff = thisRcut + endif + InteractionMap(map_i, map_j)%rCut = actualCutoff + InteractionMap(map_i, map_j)%rCutSq = actualCutoff * actualCutoff + InteractionMap(map_i, map_j)%rListSq = (actualCutoff + skinThickness)**2 + InteractionMap(map_j, map_i)%rCut = InteractionMap(map_i, map_j)%rCut + InteractionMap(map_j, map_i)%rCutSq = InteractionMap(map_i, map_j)%rCutSq + InteractionMap(map_j, map_i)%rListSq = InteractionMap(map_i, map_j)%rListSq + end do + end do + ! now the groups + + + + haveRlist = .true. + end subroutine createGroupCutoffs + subroutine setSimVariables() SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms() SIM_uses_LennardJones = SimUsesLennardJones() @@ -250,14 +394,13 @@ contains error = 0 - if (.not. havePropertyMap) then - - myStatus = 0 - - call createPropertyMap(myStatus) - + if (.not. haveInteractionMap) then + + myStatus = 0 + call createInteractionMap(myStatus) + if (myStatus .ne. 0) then - write(default_error, *) 'createPropertyMap failed in doForces!' + write(default_error, *) 'createInteractionMap failed in doForces!' error = -1 return endif @@ -515,7 +658,7 @@ contains integer :: localError integer :: propPack_i, propPack_j integer :: loopStart, loopEnd, loop - + integer :: iMap real(kind=dp) :: listSkin = 1.0 !! initialize local variables @@ -606,6 +749,12 @@ contains iend = nGroups - 1 #endif outer: do i = istart, iend + +#ifdef IS_MPI + me_i = atid_row(i) +#else + me_i = atid(i) +#endif if (update_nlist) point(i) = nlist + 1 @@ -634,14 +783,16 @@ contains endif #ifdef IS_MPI + me_j = atid_col(j) call get_interatomic_vector(q_group_Row(:,i), & q_group_Col(:,j), d_grp, rgrpsq) #else + me_j = atid(j) call get_interatomic_vector(q_group(:,i), & q_group(:,j), d_grp, rgrpsq) #endif - if (rgrpsq < rlistsq) then + if (rgrpsq < InteractionMap(me_i,me_j)%rListsq) then if (update_nlist) then nlist = nlist + 1 @@ -859,8 +1010,9 @@ contains #else me_i = atid(i) #endif - - if (PropertyMap(me_i)%is_Dipole) then + iMap = InteractionHash(me_i,me_j) + + if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then mu_i = getDipoleMoment(me_i) @@ -951,15 +1103,13 @@ contains call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, & pot, eFrame, f, t, do_pot) - if (FF_uses_dipoles .and. SIM_uses_dipoles) then - if ( PropertyMap(me_i)%is_Dipole .and. & - PropertyMap(me_j)%is_Dipole) then - if (FF_uses_RF .and. SIM_uses_RF) then - call accumulate_rf(i, j, r, eFrame, sw) - call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair) - endif - endif + if (FF_uses_RF .and. SIM_uses_RF) 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 + endif if ( iand(iMap, STICKY_PAIR).ne.0 ) then @@ -978,8 +1128,8 @@ contains endif if ( iand(iMap, 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_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, & + ! pot, A, f, t, do_pot) endif if ( iand(iMap, EAM_PAIR).ne.0 ) then