--- trunk/OOPSE-4/src/UseTheForce/doForces.F90 2005/11/01 19:09:30 2402 +++ trunk/OOPSE-4/src/UseTheForce/doForces.F90 2006/06/05 18:24:45 2787 @@ -45,7 +45,7 @@ !! @author Charles F. Vardeman II !! @author Matthew Meineke -!! @version $Id: doForces.F90,v 1.64 2005-11-01 19:09:23 chrisfen Exp $, $Date: 2005-11-01 19:09:23 $, $Name: not supported by cvs2svn $, $Revision: 1.64 $ +!! @version $Id: doForces.F90,v 1.83 2006-06-05 18:24:45 gezelter Exp $, $Date: 2006-06-05 18:24:45 $, $Name: not supported by cvs2svn $, $Revision: 1.83 $ module doForces @@ -62,6 +62,7 @@ module doForces use shapes use vector_class use eam + use suttonchen use status #ifdef IS_MPI use mpiSimulation @@ -71,12 +72,10 @@ module doForces PRIVATE #define __FORTRAN90 -#include "UseTheForce/fSwitchingFunction.h" #include "UseTheForce/fCutoffPolicy.h" #include "UseTheForce/DarkSide/fInteractionMap.h" #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h" - INTEGER, PARAMETER:: PREPAIR_LOOP = 1 INTEGER, PARAMETER:: PAIR_LOOP = 2 @@ -86,31 +85,41 @@ module doForces logical, save :: haveInteractionHash = .false. logical, save :: haveGtypeCutoffMap = .false. logical, save :: haveDefaultCutoffs = .false. - logical, save :: haveRlist = .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_SC + logical, save :: FF_uses_MEAM + logical, save :: SIM_uses_DirectionalAtoms logical, save :: SIM_uses_EAM + 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 :: 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 @@ -138,15 +147,11 @@ 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 :: listSkin - + contains - subroutine createInteractionHash(status) + subroutine createInteractionHash() integer :: nAtypes - integer, intent(out) :: status integer :: i integer :: j integer :: iHash @@ -158,6 +163,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 @@ -165,20 +172,19 @@ 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 @@ -204,6 +210,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 @@ -217,6 +225,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) @@ -238,6 +248,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) @@ -257,9 +271,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 @@ -267,6 +280,7 @@ 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 @@ -274,17 +288,11 @@ contains integer :: nGroupsInRow integer :: nGroupsInCol integer :: nGroupTypesRow,nGroupTypesCol - real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin + 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) @@ -302,7 +310,7 @@ 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 (haveDefaultCutoffs) then atypeMaxCutoff(i) = defaultRcut @@ -335,17 +343,18 @@ contains 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 (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then biggestAtypeCutoff = atypeMaxCutoff(i) endif endif enddo - - istart = 1 jstart = 1 @@ -389,11 +398,11 @@ contains allocate(groupToGtypeCol(jend)) end if - if(.not.associated(groupToGtypeCol)) then - allocate(groupToGtypeCol(jend)) + if(.not.associated(groupMaxCutoffCol)) then + allocate(groupMaxCutoffCol(jend)) else - deallocate(groupToGtypeCol) - allocate(groupToGtypeCol(jend)) + deallocate(groupMaxCutoffCol) + allocate(groupMaxCutoffCol(jend)) end if if(.not.associated(gtypeMaxCutoffCol)) then allocate(gtypeMaxCutoffCol(jend)) @@ -414,9 +423,9 @@ contains !! largest cutoff for any atypes present in this group. We also !! create gtypes at this point. - tol = 1.0d-6 + tol = 1.0e-6_dp nGroupTypesRow = 0 - + nGroupTypesCol = 0 do i = istart, iend n_in_i = groupStartRow(i+1) - groupStartRow(i) groupMaxCutoffRow(i) = 0.0_dp @@ -431,7 +440,6 @@ contains groupMaxCutoffRow(i)=atypeMaxCutoff(me_i) endif enddo - if (nGroupTypesRow.eq.0) then nGroupTypesRow = nGroupTypesRow + 1 gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i) @@ -494,17 +502,13 @@ contains groupMaxCutoffCol => groupMaxCutoffRow #endif - - - - !! allocate the gtypeCutoffMap here. 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 tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol)) - do i = 1, nGroupTypesRow + do i = 1, nGroupTypesRow do j = 1, nGroupTypesCol select case(cutoffPolicy) @@ -519,11 +523,17 @@ contains return end select gtypeCutoffMap(i,j)%rcut = thisRcut + + if (thisRcut.gt.largestRcut) largestRcut = thisRcut + gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut - skin = defaultRlist - defaultRcut - listSkin = skin ! set neighbor list skin thickness - 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 @@ -533,6 +543,7 @@ contains endif enddo enddo + if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow) if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow) if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff) @@ -546,87 +557,129 @@ contains 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 - - haveDefaultCutoffs = .true. - 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 + + localError = 0 + call setLJDefaultCutoff( defaultRcut, defaultDoShift ) + call setElectrostaticCutoffRadius( defaultRcut, defaultRsw ) + call setCutoffEAM( defaultRcut ) + call setCutoffSC( defaultRcut ) + call set_switch(defaultRsw, defaultRcut) + call setHmatDangerousRcutValue(defaultRcut) + + haveDefaultCutoffs = .true. + haveGtypeCutoffMap = .false. - subroutine setCutoffPolicy(cutPolicy) + 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() + 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() + SIM_uses_SC = SimUsesSC() + + haveSIMvariables = .true. + + return + end subroutine setSimVariables + subroutine doReadyCheck(error) integer, intent(out) :: error - integer :: myStatus 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(largestRcut, largestRcut) + call setHmatDangerousRcutValue(largestRcut) + call setCutoffEAM(largestRcut) + call setCutoffSC(largestRcut) + VisitCutoffsAfterComputing = .false. + endif + if (.not. haveSIMvariables) then call setSimVariables() endif - ! if (.not. haveRlist) then - ! write(default_error, *) 'rList has not been set in doForces!' - ! error = -1 - ! return - ! endif - if (.not. haveNeighborList) then write(default_error, *) 'neighbor list has not been initialized in doForces!' error = -1 return end if - + if (.not. haveSaneForceField) then write(default_error, *) 'Force Field is not sane in doForces!' error = -1 return end if - + #ifdef IS_MPI if (.not. isMPISimSet()) then write(default_error,*) "ERROR: mpiSimulation has not been initialized!" @@ -638,9 +691,8 @@ contains end subroutine doReadyCheck - subroutine init_FF(thisESM, thisStat) + subroutine init_FF(thisStat) - integer, intent(in) :: thisESM integer, intent(out) :: thisStat integer :: my_status, nMatches integer, pointer :: MatchList(:) => null() @@ -648,8 +700,6 @@ contains !! assume things are copacetic, unless they aren't thisStat = 0 - electrostaticSummationMethod = thisESM - !! init_FF is called *after* all of the atom types have been !! defined in atype_module using the new_atype subroutine. !! @@ -660,6 +710,7 @@ contains FF_uses_Dipoles = .false. FF_uses_GayBerne = .false. FF_uses_EAM = .false. + FF_uses_SC = .false. call getMatchingElementList(atypes, "is_Directional", .true., & nMatches, MatchList) @@ -676,7 +727,10 @@ contains call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList) if (nMatches .gt. 0) FF_uses_EAM = .true. + call getMatchingElementList(atypes, "is_SC", .true., nMatches, MatchList) + if (nMatches .gt. 0) FF_uses_SC = .true. + haveSaneForceField = .true. if (FF_uses_EAM) then @@ -746,6 +800,7 @@ contains 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 @@ -819,10 +874,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 @@ -902,108 +957,116 @@ contains list(nlist) = j endif + + if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then - 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) -#else - call get_interatomic_vector(q(:,atom1), & - q(:,atom2), d_atm, ratmsq) -#endif - endif - - if (loop .eq. PREPAIR_LOOP) then + rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut + if (loop .eq. PAIR_LOOP) then + vij = 0.0_dp + fij(1) = 0.0_dp + fij(2) = 0.0_dp + fij(3) = 0.0_dp + endif + + call get_switch(rgrpsq, sw, dswdr,rgrp, 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) = d_grp(1) + d_atm(2) = d_grp(2) + d_atm(3) = d_grp(3) + ratmsq = rgrpsq + else +#ifdef IS_MPI + 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) +#endif + 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, d_grp, rgrp) + 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, & - d_grp, rgrp) + 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) = fij(1) + fpair(1) + fij(2) = fij(2) + fpair(2) + fij(3) = fij(3) + fpair(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 @@ -1108,7 +1171,7 @@ contains endif -!!$ if (electrostaticSummationMethod.eq.REACTION_FIELD) then + if (electrostaticSummationMethod.eq.REACTION_FIELD) then ! loop over the excludes to accumulate RF stuff we've ! left out of the normal pair loop @@ -1118,11 +1181,9 @@ contains ! 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) + call get_interatomic_vector(q(:,i), q(:,j), d_atm, ratmsq) + rVal = sqrt(ratmsq) + call get_switch(ratmsq, sw, dswdr, rVal,in_switching_region) #ifdef IS_MPI call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, & vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot) @@ -1132,22 +1193,34 @@ contains #endif endif enddo -!!$ endif + endif enddo endif #ifdef IS_MPI if (do_pot) then +#ifdef SINGLE_PRECISION + call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_real,mpi_sum, & + mpi_comm_world,mpi_err) +#else call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, & mpi_comm_world,mpi_err) +#endif endif if (do_stress) then +#ifdef SINGLE_PRECISION + call mpi_allreduce(tau_Temp, tau, 9,mpi_real,mpi_sum, & + mpi_comm_world,mpi_err) + call mpi_allreduce(virial_Temp, virial,1,mpi_real,mpi_sum, & + mpi_comm_world,mpi_err) +#else 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 endif #else @@ -1162,14 +1235,11 @@ contains end subroutine do_force_loop subroutine do_pair(i, j, rijsq, d, sw, do_pot, & - eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp) -!!$ subroutine do_pair(i, j, rijsq, d, sw, do_pot, & -!!$ eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, felec) + eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut) real( kind = dp ) :: vpair, sw real( kind = dp ), dimension(LR_POT_TYPES) :: pot real( kind = dp ), dimension(3) :: fpair - real( kind = dp ), dimension(3) :: felec real( kind = dp ), dimension(nLocal) :: mfact real( kind = dp ), dimension(9,nLocal) :: eFrame real( kind = dp ), dimension(9,nLocal) :: A @@ -1182,14 +1252,18 @@ contains 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 + real ( kind = dp ) :: a_k, b_k, c_k, d_k, dx integer :: me_i, me_j + integer :: k integer :: iHash r = sqrt(rijsq) - vpair = 0.0d0 - fpair(1:3) = 0.0d0 + + vpair = 0.0_dp + fpair(1:3) = 0.0_dp #ifdef IS_MPI me_i = atid_row(i) @@ -1202,15 +1276,13 @@ 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, & + 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, & + call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, & pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot) -!!$ call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, & -!!$ pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot, felec) endif if ( iand(iHash, STICKY_PAIR).ne.0 ) then @@ -1229,7 +1301,7 @@ contains endif if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then - call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, & + call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, & pot(VDW_POT), A, f, t, do_pot) endif @@ -1247,10 +1319,15 @@ contains call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, & 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 ) :: sw @@ -1262,14 +1339,14 @@ 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) integer :: me_i, me_j, iHash r = sqrt(rijsq) - + #ifdef IS_MPI me_i = atid_row(i) me_j = atid_col(j) @@ -1281,8 +1358,12 @@ 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 @@ -1294,8 +1375,9 @@ contains if (FF_uses_EAM .and. SIM_uses_EAM) then 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 end subroutine do_preforce @@ -1307,46 +1389,59 @@ contains real( kind = dp ) :: d(3), scaled(3) integer i - d(1:3) = q_j(1:3) - q_i(1:3) + d(1) = q_j(1) - q_i(1) + d(2) = q_j(2) - q_i(2) + d(3) = q_j(3) - q_i(3) ! Wrap back into periodic box if necessary if ( SIM_uses_PBC ) then if( .not.boxIsOrthorhombic ) then ! calc the scaled coordinates. + ! scaled = matmul(HmatInv, d) - scaled = matmul(HmatInv, d) - + scaled(1) = HmatInv(1,1)*d(1) + HmatInv(1,2)*d(2) + HmatInv(1,3)*d(3) + scaled(2) = HmatInv(2,1)*d(1) + HmatInv(2,2)*d(2) + HmatInv(2,3)*d(3) + scaled(3) = HmatInv(3,1)*d(1) + HmatInv(3,2)*d(2) + HmatInv(3,3)*d(3) + ! wrap the scaled coordinates - scaled = scaled - anint(scaled) + scaled(1) = scaled(1) - anint(scaled(1), kind=dp) + scaled(2) = scaled(2) - anint(scaled(2), kind=dp) + scaled(3) = scaled(3) - anint(scaled(3), kind=dp) - ! calc the wrapped real coordinates from the wrapped scaled ! coordinates + ! d = matmul(Hmat,scaled) + d(1)= Hmat(1,1)*scaled(1) + Hmat(1,2)*scaled(2) + Hmat(1,3)*scaled(3) + d(2)= Hmat(2,1)*scaled(1) + Hmat(2,2)*scaled(2) + Hmat(2,3)*scaled(3) + d(3)= Hmat(3,1)*scaled(1) + Hmat(3,2)*scaled(2) + Hmat(3,3)*scaled(3) - d = matmul(Hmat,scaled) - else ! calc the scaled coordinates. - do i = 1, 3 - scaled(i) = d(i) * HmatInv(i,i) + scaled(1) = d(1) * HmatInv(1,1) + scaled(2) = d(2) * HmatInv(2,2) + scaled(3) = d(3) * HmatInv(3,3) + + ! wrap the scaled coordinates + + scaled(1) = scaled(1) - anint(scaled(1), kind=dp) + scaled(2) = scaled(2) - anint(scaled(2), kind=dp) + scaled(3) = scaled(3) - anint(scaled(3), kind=dp) - ! wrap the scaled coordinates + ! calc the wrapped real coordinates from the wrapped scaled + ! coordinates - scaled(i) = scaled(i) - anint(scaled(i)) + d(1) = scaled(1)*Hmat(1,1) + d(2) = scaled(2)*Hmat(2,2) + d(3) = scaled(3)*Hmat(3,3) - ! calc the wrapped real coordinates from the wrapped scaled - ! coordinates - - d(i) = scaled(i)*Hmat(i,i) - enddo endif endif - r_sq = dot_product(d,d) + r_sq = d(1)*d(1) + d(2)*d(2) + d(3)*d(3) end subroutine get_interatomic_vector @@ -1384,7 +1479,6 @@ contains call clean_EAM() endif - rf = 0.0_dp tau_Temp = 0.0_dp virial_Temp = 0.0_dp end subroutine zero_work_arrays @@ -1478,7 +1572,8 @@ 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 #ifdef PROFILE