--- trunk/OOPSE-4/src/UseTheForce/doForces.F90 2005/12/15 21:43:16 2512 +++ 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.71 2005-12-15 21:43:16 gezelter Exp $, $Date: 2005-12-15 21:43:16 $, $Name: not supported by cvs2svn $, $Revision: 1.71 $ +!! @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 @@ -72,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 @@ -149,6 +147,7 @@ contains end type gtypeCutoffs type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap + contains subroutine createInteractionHash() @@ -281,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 @@ -310,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 @@ -341,6 +341,10 @@ contains 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 @@ -394,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)) @@ -419,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 @@ -575,19 +579,18 @@ contains defaultDoShift = .true. endif - + 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) - + call setCutoffEAM( defaultRcut ) + call setCutoffSC( defaultRcut ) + call set_switch(defaultRsw, defaultRcut) + call setHmatDangerousRcutValue(defaultRcut) + haveDefaultCutoffs = .true. haveGtypeCutoffMap = .false. + end subroutine setCutoffs subroutine cWasLame() @@ -632,6 +635,7 @@ contains SIM_requires_postpair_calc = SimRequiresPostpairCalc() SIM_requires_prepair_calc = SimRequiresPrepairCalc() SIM_uses_PBC = SimUsesPBC() + SIM_uses_SC = SimUsesSC() haveSIMvariables = .true. @@ -640,7 +644,6 @@ contains subroutine doReadyCheck(error) integer, intent(out) :: error - integer :: myStatus error = 0 @@ -653,34 +656,30 @@ contains call createGtypeCutoffMap() endif - if (VisitCutoffsAfterComputing) then - call set_switch(GROUP_SWITCH, largestRcut, largestRcut) + 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!" @@ -711,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) @@ -727,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 @@ -954,19 +957,18 @@ contains list(nlist) = j endif - - - + if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut if (loop .eq. PAIR_LOOP) then - vij = 0.0d0 - fij(1:3) = 0.0d0 + 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, & - group_switch, in_switching_region) + call get_switch(rgrpsq, sw, dswdr,rgrp, in_switching_region) n_in_j = groupStartCol(j+1) - groupStartCol(j) @@ -981,7 +983,9 @@ contains 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) + d_atm(1) = d_grp(1) + d_atm(2) = d_grp(2) + d_atm(3) = d_grp(3) ratmsq = rgrpsq else #ifdef IS_MPI @@ -1014,7 +1018,9 @@ contains d_grp, rgrp, rCut) #endif vij = vij + vpair - fij(1:3) = fij(1:3) + fpair(1:3) + fij(1) = fij(1) + fpair(1) + fij(2) = fij(2) + fpair(2) + fij(3) = fij(3) + fpair(3) endif enddo inner enddo @@ -1175,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) @@ -1196,15 +1200,27 @@ contains #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 @@ -1238,13 +1254,16 @@ contains 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) @@ -1282,7 +1301,7 @@ contains endif if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then - call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, & + call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, & pot(VDW_POT), A, f, t, do_pot) endif @@ -1305,8 +1324,6 @@ contains call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, & pot(METALLIC_POT), f, do_pot) endif - - end subroutine do_pair @@ -1329,7 +1346,7 @@ contains integer :: me_i, me_j, iHash r = sqrt(rijsq) - + #ifdef IS_MPI me_i = atid_row(i) me_j = atid_col(j) @@ -1361,8 +1378,6 @@ contains if (FF_uses_SC .and. SIM_uses_SC) then call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT)) endif - - end subroutine do_preforce @@ -1374,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