--- trunk/OOPSE-4/src/UseTheForce/doForces.F90 2005/09/18 20:45:38 2309 +++ trunk/OOPSE-4/src/UseTheForce/doForces.F90 2005/10/18 15:01:42 2381 @@ -45,7 +45,7 @@ !! @author Charles F. Vardeman II !! @author Matthew Meineke -!! @version $Id: doForces.F90,v 1.46 2005-09-18 20:45:38 chrisfen Exp $, $Date: 2005-09-18 20:45:38 $, $Name: not supported by cvs2svn $, $Revision: 1.46 $ +!! @version $Id: doForces.F90,v 1.60 2005-10-18 15:01:42 chrisfen Exp $, $Date: 2005-10-18 15:01:42 $, $Name: not supported by cvs2svn $, $Revision: 1.60 $ module doForces @@ -58,8 +58,7 @@ 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 @@ -75,6 +74,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,6 +85,7 @@ 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 @@ -122,9 +123,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,6 +140,7 @@ module doForces integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist + real(kind=dp),save :: listSkin contains @@ -177,10 +184,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 @@ -257,9 +270,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 @@ -273,6 +288,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. @@ -288,67 +304,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 + 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 (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 (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 @@ -356,46 +427,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 @@ -403,10 +521,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 @@ -419,6 +554,8 @@ contains defaultRsw = defRsw defaultRlist = defRlist cutoffPolicy = cutPolicy + + haveDefaultCutoffs = .true. end subroutine setDefaultCutoffs subroutine setCutoffPolicy(cutPolicy) @@ -504,7 +641,6 @@ contains subroutine init_FF(thisESM, thisStat) integer, intent(in) :: thisESM - real(kind=dp), intent(in) :: dampingAlpha integer, intent(out) :: thisStat integer :: my_status, nMatches integer, pointer :: MatchList(:) => null() @@ -543,22 +679,6 @@ contains haveSaneForceField = .true. - - !! check to make sure the reaction field setting makes sense - - if (FF_uses_Dipoles) then - if (electrostaticSummationMethod == 3) then - dielect = getDielect() - call initialize_rf(dielect) - endif - else - if (electrostaticSummationMethod == 3) 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) @@ -570,15 +690,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) @@ -612,13 +723,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 @@ -643,7 +754,7 @@ contains integer :: propPack_i, propPack_j integer :: loopStart, loopEnd, loop integer :: iHash - real(kind=dp) :: listSkin = 1.0 + !! initialize local variables @@ -768,9 +879,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 @@ -891,6 +1002,7 @@ contains endif end if enddo + enddo outer if (update_nlist) then @@ -950,19 +1062,23 @@ 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 @@ -970,7 +1086,7 @@ contains if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then - if (electrostaticSummationMethod == 3) then + if (electrostaticSummationMethod == REACTION_FIELD) then #ifdef IS_MPI call scatter(rf_Row,rf,plan_atom_row_3d) @@ -1001,9 +1117,9 @@ contains !! potential and torques: call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot) #ifdef IS_MPI - pot_local = pot_local + rfpot + pot_local(ELECTROSTATIC_POT) = pot_local(ELECTROSTATIC_POT) + rfpot #else - pot = pot + rfpot + pot(ELECTROSTATIC_POT) = pot(ELECTROSTATIC_POT) + rfpot #endif endif @@ -1015,7 +1131,8 @@ contains #ifdef IS_MPI if (do_pot) then - pot = pot + pot_local + pot(1:LR_POT_TYPES) = pot(1:LR_POT_TYPES) & + + pot_local(1:LR_POT_TYPES) !! we assume the c code will do the allreduce to get the total potential !! we could do it right here if we needed to... endif @@ -1041,7 +1158,8 @@ contains subroutine do_pair(i, j, rijsq, d, sw, do_pot, & eFrame, A, f, t, pot, vpair, fpair) - 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 @@ -1073,14 +1191,15 @@ contains 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, 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) + pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot) - if (electrostaticSummationMethod == 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) @@ -1091,37 +1210,37 @@ contains 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, 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 end subroutine do_pair @@ -1129,7 +1248,8 @@ contains subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, & 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 @@ -1164,10 +1284,10 @@ contains 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 @@ -1362,7 +1482,7 @@ contains function FF_RequiresPostpairCalc() result(doesit) logical :: doesit - if (electrostaticSummationMethod == 3) doesit = .true. + if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true. end function FF_RequiresPostpairCalc #ifdef PROFILE