--- trunk/OOPSE-4/src/UseTheForce/doForces.F90 2005/09/07 19:18:54 2284 +++ trunk/OOPSE-4/src/UseTheForce/doForces.F90 2005/09/19 23:21:46 2310 @@ -45,7 +45,7 @@ !! @author Charles F. Vardeman II !! @author Matthew Meineke -!! @version $Id: doForces.F90,v 1.37 2005-09-07 19:18:54 gezelter Exp $, $Date: 2005-09-07 19:18:54 $, $Name: not supported by cvs2svn $, $Revision: 1.37 $ +!! @version $Id: doForces.F90,v 1.47 2005-09-19 23:21:46 chrisfen Exp $, $Date: 2005-09-19 23:21:46 $, $Name: not supported by cvs2svn $, $Revision: 1.47 $ module doForces @@ -58,7 +58,7 @@ module doForces use lj use sticky use electrostatic_module - use reaction_field + use reaction_field_module use gb_pair use shapes use vector_class @@ -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 @@ -91,16 +92,14 @@ module doForces 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 @@ -256,9 +255,11 @@ contains logical :: i_is_GB logical :: i_is_EAM logical :: i_is_Shape + logical :: GtypeFound integer :: myStatus, nAtypes, i, j, istart, iend, jstart, jend integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes + integer :: nGroupsInRow real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin real(kind=dp) :: biggestAtypeCutoff @@ -271,9 +272,12 @@ contains return endif endif - +#ifdef IS_MPI + nGroupsInRow = getNgroupsInRow(plan_group_row) +#endif nAtypes = getSize(atypes) - +! Set all of the initial cutoffs to zero. + atypeMaxCutoff = 0.0_dp do i = 1, nAtypes if (SimHasAtype(i)) then call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ) @@ -284,7 +288,7 @@ contains call getElementProperty(atypes, i, "is_EAM", i_is_EAM) call getElementProperty(atypes, i, "is_Shape", i_is_Shape) - atypeMaxCutoff(i) = 0.0_dp + if (i_is_LJ) then thisRcut = getSigma(i) * 2.5_dp if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut @@ -334,6 +338,8 @@ contains allocate(groupToGtype(iend)) allocate(groupMaxCutoff(iend)) allocate(gtypeMaxCutoff(iend)) + groupMaxCutoff = 0.0_dp + gtypeMaxCutoff = 0.0_dp endif !! first we do a single loop over the cutoff groups to find the !! largest cutoff for any atypes present in this group. We also @@ -353,25 +359,29 @@ contains #endif if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then groupMaxCutoff(i)=atypeMaxCutoff(me_i) - endif + endif enddo + if (nGroupTypes.eq.0) then nGroupTypes = nGroupTypes + 1 gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i) groupToGtype(i) = nGroupTypes else + GtypeFound = .false. do g = 1, nGroupTypes - if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).gt.tol) then - nGroupTypes = nGroupTypes + 1 - gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i) - groupToGtype(i) = nGroupTypes - else + if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then groupToGtype(i) = g + GtypeFound = .true. endif enddo + if (.not.GtypeFound) then + nGroupTypes = nGroupTypes + 1 + gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i) + groupToGtype(i) = nGroupTypes + endif endif - enddo - + enddo + !! allocate the gtypeCutoffMap here. allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes)) !! then we do a double loop over all the group TYPES to find the cutoff @@ -400,32 +410,29 @@ contains enddo haveGtypeCutoffMap = .true. - - end subroutine createGtypeCutoffMap - - subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy) - real(kind=dp),intent(in) :: defRcut, defRsw, defRlist - integer, intent(in) :: cutPolicy - - defaultRcut = defRcut - defaultRsw = defRsw - defaultRlist = defRlist - cutoffPolicy = cutPolicy - end subroutine setDefaultCutoffs - - subroutine setCutoffPolicy(cutPolicy) + end subroutine createGtypeCutoffMap + subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy) + real(kind=dp),intent(in) :: defRcut, defRsw, defRlist integer, intent(in) :: cutPolicy + + defaultRcut = defRcut + defaultRsw = defRsw + defaultRlist = defRlist cutoffPolicy = cutPolicy - call createGtypeCutoffMap() + end subroutine setDefaultCutoffs + subroutine setCutoffPolicy(cutPolicy) + + integer, intent(in) :: cutPolicy + cutoffPolicy = cutPolicy + call createGtypeCutoffMap() end subroutine setCutoffPolicy subroutine setSimVariables() SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms() SIM_uses_EAM = SimUsesEAM() - SIM_uses_RF = SimUsesRF() SIM_requires_postpair_calc = SimRequiresPostpairCalc() SIM_requires_prepair_calc = SimRequiresPrepairCalc() SIM_uses_PBC = SimUsesPBC() @@ -495,32 +502,19 @@ contains end subroutine doReadyCheck - subroutine init_FF(use_RF, use_UW, use_DW, thisStat) + subroutine init_FF(thisESM, thisStat) - logical, intent(in) :: use_RF - logical, intent(in) :: use_UW - logical, intent(in) :: use_DW + integer, intent(in) :: thisESM integer, intent(out) :: thisStat integer :: my_status, nMatches - integer :: corrMethod 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 + electrostaticSummationMethod = thisESM - !! set the electrostatic correction method - if (use_UW) then - corrMethod = 1 - elseif (use_DW) then - corrMethod = 2 - else - corrMethod = 0 - endif - !! init_FF is called *after* all of the atom types have been !! defined in atype_module using the new_atype subroutine. !! @@ -550,15 +544,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 (FF_uses_RF) then + if (electrostaticSummationMethod == REACTION_FIELD) then write(default_error,*) 'Using Reaction Field with no dipoles? Huh?' thisStat = -1 haveSaneForceField = .false. @@ -976,7 +970,7 @@ contains if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then - if (FF_uses_RF .and. SIM_uses_RF) then + if (electrostaticSummationMethod == REACTION_FIELD) then #ifdef IS_MPI call scatter(rf_Row,rf,plan_atom_row_3d) @@ -1060,7 +1054,6 @@ contains real ( kind = dp ), intent(inout) :: rijsq real ( kind = dp ) :: r real ( kind = dp ), intent(inout) :: d(3) - real ( kind = dp ) :: ebalance integer :: me_i, me_j integer :: iHash @@ -1085,9 +1078,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) + pot, eFrame, f, t, do_pot) - if (FF_uses_RF .and. SIM_uses_RF) then + if (electrostaticSummationMethod == REACTION_FIELD) then ! CHECK ME (RF needs to know about all electrostatic types) call accumulate_rf(i, j, r, eFrame, sw) @@ -1150,6 +1143,8 @@ contains integer :: me_i, me_j, iHash + r = sqrt(rijsq) + #ifdef IS_MPI me_i = atid_row(i) me_j = atid_col(j) @@ -1367,7 +1362,7 @@ contains function FF_RequiresPostpairCalc() result(doesit) logical :: doesit - doesit = FF_uses_RF + if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true. end function FF_RequiresPostpairCalc #ifdef PROFILE