--- trunk/OOPSE/libmdtools/do_Forces.F90 2003/03/21 17:42:12 378 +++ trunk/OOPSE/libmdtools/do_Forces.F90 2004/01/05 22:18:52 897 @@ -4,7 +4,7 @@ !! @author Charles F. Vardeman II !! @author Matthew Meineke -!! @version $Id: do_Forces.F90,v 1.1.1.1 2003-03-21 17:42:12 mmeineke Exp $, $Date: 2003-03-21 17:42:12 $, $Name: not supported by cvs2svn $, $Revision: 1.1.1.1 $ +!! @version $Id: do_Forces.F90,v 1.43 2004-01-05 22:18:52 chuckv Exp $, $Date: 2004-01-05 22:18:52 $, $Name: not supported by cvs2svn $, $Revision: 1.43 $ module do_Forces use force_globals @@ -17,6 +17,9 @@ module do_Forces use dipole_dipole use reaction_field use gb_pair + use vector_class + use eam + use status #ifdef IS_MPI use mpiSimulation #endif @@ -27,7 +30,8 @@ module do_Forces #define __FORTRAN90 #include "fForceField.h" - logical, save :: do_forces_initialized = .false. + logical, save :: do_forces_initialized = .false., haveRlist = .false. + logical, save :: havePolicies = .false. logical, save :: FF_uses_LJ logical, save :: FF_uses_sticky logical, save :: FF_uses_dipoles @@ -35,11 +39,36 @@ module do_Forces logical, save :: FF_uses_GB logical, save :: FF_uses_EAM + real(kind=dp), save :: rlist, rlistsq + public :: init_FF public :: do_force_loop + public :: setRlistDF +#ifdef PROFILE + public :: getforcetime + real, save :: forceTime = 0 + real :: forceTimeInitial, forceTimeFinal + integer :: nLoops +#endif + + logical, allocatable :: propertyMapI(:,:) + logical, allocatable :: propertyMapJ(:,:) + contains + subroutine setRlistDF( this_rlist ) + + real(kind=dp) :: this_rlist + + rlist = this_rlist + rlistsq = rlist * rlist + + haveRlist = .true. + if( havePolicies ) do_forces_initialized = .true. + + end subroutine setRlistDF + subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat) integer, intent(in) :: LJMIXPOLICY @@ -87,12 +116,9 @@ contains !! check to make sure the FF_uses_RF setting makes sense if (FF_uses_dipoles) then - rrf = getRrf() - rt = getRt() - call initialize_dipole(rrf, rt) if (FF_uses_RF) then dielect = getDielect() - call initialize_rf(rrf, rt, dielect) + call initialize_rf(dielect) endif else if (FF_uses_RF) then @@ -100,17 +126,15 @@ contains thisStat = -1 return endif - endif + endif if (FF_uses_LJ) then - call getRcut(rcut) - select case (LJMIXPOLICY) case (LB_MIXING_RULE) - call init_lj_FF(LB_MIXING_RULE, rcut, my_status) + call init_lj_FF(LB_MIXING_RULE, my_status) case (EXPLICIT_MIXING_RULE) - call init_lj_FF(EXPLICIT_MIXING_RULE, rcut, my_status) + call init_lj_FF(EXPLICIT_MIXING_RULE, my_status) case default write(default_error,*) 'unknown LJ Mixing Policy!' thisStat = -1 @@ -129,6 +153,18 @@ contains return end if endif + + + if (FF_uses_EAM) then + call init_EAM_FF(my_status) + if (my_status /= 0) then + write(*,*) "init_EAM_FF returned a bad status" + thisStat = -1 + return + end if + endif + + if (FF_uses_GB) then call check_gb_pair_FF(my_status) @@ -140,10 +176,20 @@ contains if (FF_uses_GB .and. FF_uses_LJ) then endif + if (.not. do_forces_initialized) then + !! Create neighbor lists + call expandNeighborList(getNlocal(), my_status) + if (my_Status /= 0) then + write(default_error,*) "SimSetup: ExpandNeighborList returned error." + thisStat = -1 + return + endif + endif + + havePolicies = .true. + if( haveRlist ) do_forces_initialized = .true. - do_forces_initialized = .true. - end subroutine init_FF @@ -161,34 +207,40 @@ contains real ( kind = dp ), dimension(3,getNlocal()) :: f !! Torsion array provided by C, dimensioned by getNlocal real( kind = dp ), dimension(3,getNlocal()) :: t + !! Stress Tensor real( kind = dp), dimension(9) :: tau real ( kind = dp ) :: pot logical ( kind = 2) :: do_pot_c, do_stress_c logical :: do_pot logical :: do_stress -#ifdef IS_MPI +#ifdef IS_MPI real( kind = DP ) :: pot_local integer :: nrow integer :: ncol + integer :: nprocs #endif integer :: nlocal integer :: natoms logical :: update_nlist integer :: i, j, jbeg, jend, jnab integer :: nlist - real( kind = DP ) :: rijsq, rlistsq, rcutsq, rlist, rcut + real( kind = DP ) :: rijsq real(kind=dp),dimension(3) :: d real(kind=dp) :: rfpot, mu_i, virial - integer :: me_i + integer :: me_i, me_j logical :: is_dp_i integer :: neighborListSize integer :: listerror, error integer :: localError + integer :: propPack_i, propPack_j + real(kind=dp) :: listSkin = 1.0 + !! initialize local variables #ifdef IS_MPI + pot_local = 0.0_dp nlocal = getNlocal() nrow = getNrow(plan_row) ncol = getNcol(plan_col) @@ -197,11 +249,9 @@ contains natoms = nlocal #endif - call getRcut(rcut,rc2=rcutsq) - call getRlist(rlist,rlistsq) - call check_initialization(localError) if ( localError .ne. 0 ) then + call handleError("do_force_loop","Not Initialized") error = -1 return end if @@ -210,6 +260,77 @@ contains do_pot = do_pot_c do_stress = do_stress_c + +#ifdef IS_MPI + if (.not.allocated(propertyMapI)) then + allocate(propertyMapI(5,nrow)) + endif + + do i = 1, nrow + me_i = atid_row(i) +#else + if (.not.allocated(propertyMapI)) then + allocate(propertyMapI(5,nlocal)) + endif + + do i = 1, natoms + me_i = atid(i) +#endif + + propertyMapI(1:5,i) = .false. + + call getElementProperty(atypes, me_i, "propertyPack", propPack_i) + + ! unpack the properties + + if (iand(propPack_i, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) & + propertyMapI(1, i) = .true. + if (iand(propPack_i, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) & + propertyMapI(2, i) = .true. + if (iand(propPack_i, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) & + propertyMapI(3, i) = .true. + if (iand(propPack_i, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) & + propertyMapI(4, i) = .true. + if (iand(propPack_i, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) & + propertyMapI(5, i) = .true. + + end do + +#ifdef IS_MPI + if (.not.allocated(propertyMapJ)) then + allocate(propertyMapJ(5,ncol)) + endif + + do j = 1, ncol + me_j = atid_col(j) +#else + if (.not.allocated(propertyMapJ)) then + allocate(propertyMapJ(5,nlocal)) + endif + + do j = 1, natoms + me_j = atid(j) +#endif + + propertyMapJ(1:5,j) = .false. + + call getElementProperty(atypes, me_j, "propertyPack", propPack_j) + + ! unpack the properties + + if (iand(propPack_j, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) & + propertyMapJ(1, j) = .true. + if (iand(propPack_j, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) & + propertyMapJ(2, j) = .true. + if (iand(propPack_j, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) & + propertyMapJ(3, j) = .true. + if (iand(propPack_j, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) & + propertyMapJ(4, j) = .true. + if (iand(propPack_j, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) & + propertyMapJ(5, j) = .true. + + end do + ! Gather all information needed by all force loops: #ifdef IS_MPI @@ -226,31 +347,186 @@ contains endif #endif - + +!! Begin force loop timing: +#ifdef PROFILE + call cpu_time(forceTimeInitial) + nloops = nloops + 1 +#endif + if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then !! See if we need to update neighbor lists - call checkNeighborList(nlocal, q, rcut, rlist, update_nlist) + call checkNeighborList(nlocal, q, listSkin, update_nlist) !! if_mpi_gather_stuff_for_prepair !! do_prepair_loop_if_needed !! if_mpi_scatter_stuff_from_prepair !! if_mpi_gather_stuff_from_prepair_to_main_loop - else - !! See if we need to update neighbor lists - call checkNeighborList(nlocal, q, rcut, rlist, update_nlist) + +!--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>> +#ifdef IS_MPI + + if (update_nlist) then + + !! save current configuration, construct neighbor list, + !! and calculate forces + call saveNeighborList(nlocal, q) + + neighborListSize = size(list) + nlist = 0 + + do i = 1, nrow + point(i) = nlist + 1 + + prepair_inner: do j = 1, ncol + + if (skipThisPair(i,j)) cycle prepair_inner + + call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq) + + if (rijsq < rlistsq) then + + nlist = nlist + 1 + + if (nlist > neighborListSize) then + call expandNeighborList(nlocal, listerror) + if (listerror /= 0) then + error = -1 + write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded." + return + end if + neighborListSize = size(list) + endif + + list(nlist) = j + call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local) + endif + enddo prepair_inner + enddo + + point(nrow + 1) = nlist + 1 + + else !! (of update_check) + + ! use the list to find the neighbors + do i = 1, nrow + JBEG = POINT(i) + JEND = POINT(i+1) - 1 + ! check thiat molecule i has neighbors + if (jbeg .le. jend) then + + do jnab = jbeg, jend + j = list(jnab) + + call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq) + call do_prepair(i, j, rijsq, d, do_pot, do_stress, & + u_l, A, f, t, pot_local) + + enddo + endif + enddo endif -#ifdef IS_MPI +#else if (update_nlist) then + ! save current configuration, contruct neighbor list, + ! and calculate forces + call saveNeighborList(natoms, q) + + neighborListSize = size(list) + + nlist = 0 + + do i = 1, natoms-1 + point(i) = nlist + 1 + + prepair_inner: do j = i+1, natoms + + if (skipThisPair(i,j)) cycle prepair_inner + + call get_interatomic_vector(q(:,i), q(:,j), d, rijsq) + + + if (rijsq < rlistsq) then + + + nlist = nlist + 1 + + if (nlist > neighborListSize) then + call expandNeighborList(natoms, listerror) + if (listerror /= 0) then + error = -1 + write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded." + return + end if + neighborListSize = size(list) + endif + + list(nlist) = j + + call do_prepair(i, j, rijsq, d, do_pot, do_stress, & + u_l, A, f, t, pot) + + endif + enddo prepair_inner + enddo + + point(natoms) = nlist + 1 + + else !! (update) + + ! use the list to find the neighbors + do i = 1, natoms-1 + JBEG = POINT(i) + JEND = POINT(i+1) - 1 + ! check thiat molecule i has neighbors + if (jbeg .le. jend) then + + do jnab = jbeg, jend + j = list(jnab) + + call get_interatomic_vector(q(:,i), q(:,j), d, rijsq) + call do_prepair(i, j, rijsq, d, do_pot, do_stress, & + u_l, A, f, t, pot) + + enddo + endif + enddo + endif +#endif + !! Do rest of preforce calculations + !! do necessary preforce calculations + call do_preforce(nlocal,pot) + ! we have already updated the neighbor list set it to false... + update_nlist = .false. + else + !! See if we need to update neighbor lists for non pre-pair + call checkNeighborList(nlocal, q, listSkin, update_nlist) + endif + + + + + +!---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>> + + + + + +#ifdef IS_MPI + + if (update_nlist) then !! save current configuration, construct neighbor list, !! and calculate forces - call saveNeighborList(q) + call saveNeighborList(nlocal, q) neighborListSize = size(list) nlist = 0 do i = 1, nrow + point(i) = nlist + 1 inner: do j = 1, ncol @@ -259,7 +535,7 @@ contains call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq) - if (rijsq < rlistsq) then + if (rijsq < rlistsq) then nlist = nlist + 1 @@ -275,10 +551,9 @@ contains list(nlist) = j - if (rijsq < rcutsq) then - call do_pair(i, j, rijsq, d, do_pot, do_stress, & - u_l, A, f, t,pot) - endif + call do_pair(i, j, rijsq, d, do_pot, do_stress, & + u_l, A, f, t, pot_local) + endif enddo inner enddo @@ -299,7 +574,7 @@ contains call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq) call do_pair(i, j, rijsq, d, do_pot, do_stress, & - u_l, A, f, t,pot) + u_l, A, f, t, pot_local) enddo endif @@ -309,10 +584,10 @@ contains #else if (update_nlist) then - + ! save current configuration, contruct neighbor list, ! and calculate forces - call saveNeighborList(q) + call saveNeighborList(natoms, q) neighborListSize = size(list) @@ -323,11 +598,12 @@ contains inner: do j = i+1, natoms - if (skipThisPair(i,j)) cycle inner - + if (skipThisPair(i,j)) cycle inner + call get_interatomic_vector(q(:,i), q(:,j), d, rijsq) - if (rijsq < rlistsq) then + + if (rijsq < rlistsq) then nlist = nlist + 1 @@ -343,10 +619,9 @@ contains list(nlist) = j - if (rijsq < rcutsq) then - call do_pair(i, j, rijsq, d, do_pot, do_stress, & - u_l, A, f, t,pot) - endif + call do_pair(i, j, rijsq, d, do_pot, do_stress, & + u_l, A, f, t, pot) + endif enddo inner enddo @@ -367,7 +642,7 @@ contains call get_interatomic_vector(q(:,i), q(:,j), d, rijsq) call do_pair(i, j, rijsq, d, do_pot, do_stress, & - u_l, A, f, t,pot) + u_l, A, f, t, pot) enddo endif @@ -377,18 +652,36 @@ contains #endif ! phew, done with main loop. - + +!! Do timing +#ifdef PROFILE + call cpu_time(forceTimeFinal) + forceTime = forceTime + forceTimeFinal - forceTimeInitial +#endif + + #ifdef IS_MPI !!distribute forces - - call scatter(f_Row,f,plan_row3d) + + f_temp = 0.0_dp + call scatter(f_Row,f_temp,plan_row3d) + do i = 1,nlocal + f(1:3,i) = f(1:3,i) + f_temp(1:3,i) + end do + + f_temp = 0.0_dp call scatter(f_Col,f_temp,plan_col3d) do i = 1,nlocal f(1:3,i) = f(1:3,i) + f_temp(1:3,i) end do if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then - call scatter(t_Row,t,plan_row3d) + t_temp = 0.0_dp + call scatter(t_Row,t_temp,plan_row3d) + do i = 1,nlocal + t(1:3,i) = t(1:3,i) + t_temp(1:3,i) + end do + t_temp = 0.0_dp call scatter(t_Col,t_temp,plan_col3d) do i = 1,nlocal @@ -399,20 +692,20 @@ contains if (do_pot) then ! scatter/gather pot_row into the members of my column call scatter(pot_Row, pot_Temp, plan_row) - + ! 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) enddo + + pot_Temp = 0.0_DP - pot_Temp = 0.0_DP - call scatter(pot_Col, pot_Temp, plan_col) do i = 1, nlocal pot_local = pot_local + pot_Temp(i) enddo - + endif #endif @@ -460,15 +753,15 @@ contains #ifdef IS_MPI if (do_pot) then - pot = pot_local + pot = pot + pot_local !! 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 if (do_stress) then - call mpi_allreduce(tau, tau_Temp,9,mpi_double_precision,mpi_sum, & + call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, & mpi_comm_world,mpi_err) - call mpi_allreduce(virial, virial_Temp,1,mpi_double_precision,mpi_sum, & + call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, & mpi_comm_world,mpi_err) endif @@ -478,18 +771,20 @@ contains tau = tau_Temp virial = virial_Temp endif - + #endif + + end subroutine do_force_loop - subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t,pot) + subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot) real( kind = dp ) :: pot - real( kind = dp ), dimension(:,:) :: u_l - real (kind=dp), dimension(:,:) :: A - real (kind=dp), dimension(:,:) :: f - real (kind=dp), dimension(:,:) :: t + real( kind = dp ), dimension(3,getNlocal()) :: u_l + real (kind=dp), dimension(9,getNlocal()) :: A + real (kind=dp), dimension(3,getNlocal()) :: f + real (kind=dp), dimension(3,getNlocal()) :: t logical, intent(inout) :: do_pot, do_stress integer, intent(in) :: i, j @@ -499,12 +794,17 @@ contains logical :: is_LJ_i, is_LJ_j logical :: is_DP_i, is_DP_j logical :: is_GB_i, is_GB_j + logical :: is_EAM_i,is_EAM_j logical :: is_Sticky_i, is_Sticky_j integer :: me_i, me_j - + integer :: propPack_i + integer :: propPack_j r = sqrt(rijsq) #ifdef IS_MPI + if (tagRow(i) .eq. tagColumn(j)) then + write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j) + endif me_i = atid_row(i) me_j = atid_col(j) @@ -515,29 +815,22 @@ contains me_j = atid(j) #endif - + if (FF_uses_LJ .and. SimUsesLJ()) then - call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i) - call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j) - if ( is_LJ_i .and. is_LJ_j ) & + if ( propertyMapI(1, me_i) .and. propertyMapJ(1, me_j) ) & call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress) + endif if (FF_uses_dipoles .and. SimUsesDipoles()) then - call getElementProperty(atypes, me_i, "is_DP", is_DP_i) - call getElementProperty(atypes, me_j, "is_DP", is_DP_j) - if ( is_DP_i .and. is_DP_j ) then - + if ( propertyMapI(2, me_i) .and. propertyMapJ(2, me_j)) then call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, & do_pot, do_stress) - if (FF_uses_RF .and. SimUsesRF()) then - call accumulate_rf(i, j, r, u_l) call rf_correct_forces(i, j, d, r, u_l, f, do_stress) - endif endif @@ -545,10 +838,7 @@ contains if (FF_uses_Sticky .and. SimUsesSticky()) then - call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i) - call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j) - - if ( is_Sticky_i .and. is_Sticky_j ) then + if ( propertyMapI(3, me_i) .and. propertyMapJ(3, me_j)) then call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, & do_pot, do_stress) endif @@ -556,45 +846,146 @@ contains if (FF_uses_GB .and. SimUsesGB()) then - - call getElementProperty(atypes, me_i, "is_GB", is_GB_i) - call getElementProperty(atypes, me_j, "is_GB", is_GB_j) - if ( is_GB_i .and. is_GB_j ) then + if ( propertyMapI(4, me_i) .and. propertyMapJ(4, me_j)) then call do_gb_pair(i, j, d, r, rijsq, u_l, pot, f, t, & do_pot, do_stress) endif + endif + + + if (FF_uses_EAM .and. SimUsesEAM()) then + + if ( propertyMapI(5, me_i) .and. propertyMapJ(5, me_j)) then + call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress) + endif + + endif + end subroutine do_pair + + + + subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot) + real( kind = dp ) :: pot + real( kind = dp ), dimension(3,getNlocal()) :: u_l + real (kind=dp), dimension(9,getNlocal()) :: A + real (kind=dp), dimension(3,getNlocal()) :: f + real (kind=dp), dimension(3,getNlocal()) :: t + + logical, intent(inout) :: do_pot, do_stress + integer, intent(in) :: i, j + real ( kind = dp ), intent(inout) :: rijsq + real ( kind = dp ) :: r + real ( kind = dp ), intent(inout) :: d(3) + + logical :: is_EAM_i, is_EAM_j + + integer :: me_i, me_j + + r = sqrt(rijsq) + +#ifdef IS_MPI + if (tagRow(i) .eq. tagColumn(j)) then + write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j) + endif + + me_i = atid_row(i) + me_j = atid_col(j) + +#else + + me_i = atid(i) + me_j = atid(j) + +#endif + + if (FF_uses_EAM .and. SimUsesEAM()) then + call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i) + call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j) + + if ( is_EAM_i .and. is_EAM_j ) & + call calc_EAM_prepair_rho(i, j, d, r, rijsq ) + endif + end subroutine do_prepair + + + + + subroutine do_preforce(nlocal,pot) + integer :: nlocal + real( kind = dp ) :: pot + + if (FF_uses_EAM .and. SimUsesEAM()) then + call calc_EAM_preforce_Frho(nlocal,pot) + endif + + + end subroutine do_preforce + + subroutine get_interatomic_vector(q_i, q_j, d, r_sq) real (kind = dp), dimension(3) :: q_i real (kind = dp), dimension(3) :: q_j real ( kind = dp ), intent(out) :: r_sq - real( kind = dp ) :: d(3) - real( kind = dp ) :: d_old(3) - d(1:3) = q_i(1:3) - q_j(1:3) - d_old = d + real( kind = dp ) :: d(3), scaled(3) + integer i + + d(1:3) = q_j(1:3) - q_i(1:3) + ! Wrap back into periodic box if necessary if ( SimUsesPBC() ) then + + if( .not.boxIsOrthorhombic ) then + ! calc the scaled coordinates. + + scaled = matmul(HmatInv, d) + + ! wrap the scaled coordinates - d(1:3) = d(1:3) - box(1:3) * sign(1.0_dp,d(1:3)) * & - int(abs(d(1:3)/box(1:3)) + 0.5_dp) + scaled = scaled - anint(scaled) + + ! calc the wrapped real coordinates from the wrapped scaled + ! coordinates + + d = matmul(Hmat,scaled) + + else + ! calc the scaled coordinates. + + do i = 1, 3 + scaled(i) = d(i) * HmatInv(i,i) + + ! wrap the scaled coordinates + + scaled(i) = scaled(i) - anint(scaled(i)) + + ! 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) - + end subroutine get_interatomic_vector - + subroutine check_initialization(error) integer, intent(out) :: error error = 0 ! Make sure we are properly initialized. if (.not. do_forces_initialized) then + write(*,*) "Forces not initialized" error = -1 return endif @@ -642,18 +1033,26 @@ contains #endif + + if (FF_uses_EAM .and. SimUsesEAM()) then + call clean_EAM() + endif + + + + + rf = 0.0_dp tau_Temp = 0.0_dp virial_Temp = 0.0_dp - end subroutine zero_work_arrays function skipThisPair(atom1, atom2) result(skip_it) - integer, intent(in) :: atom1 integer, intent(in), optional :: atom2 logical :: skip_it integer :: unique_id_1, unique_id_2 + integer :: me_i,me_j integer :: i skip_it = .false. @@ -671,7 +1070,7 @@ contains !! in the normal loop, the atom numbers are unique unique_id_1 = atom1 #endif - + !! We were called with only one atom, so just check the global exclude !! list for this atom if (.not. present(atom2)) then @@ -689,7 +1088,7 @@ contains #else unique_id_2 = atom2 #endif - + #ifdef IS_MPI !! this situation should only arise in MPI simulations if (unique_id_1 == unique_id_2) then @@ -699,14 +1098,18 @@ contains !! this prevents us from doing the pair on multiple processors if (unique_id_1 < unique_id_2) then - if (mod(unique_id_1 + unique_id_2,2) == 0) skip_it = .true. - return + if (mod(unique_id_1 + unique_id_2,2) == 0) then + skip_it = .true. + return + endif else - if (mod(unique_id_1 + unique_id_2,2) == 1) skip_it = .true. - return + if (mod(unique_id_1 + unique_id_2,2) == 1) then + skip_it = .true. + return + endif endif #endif - + !! the rest of these situations can happen in all simulations: do i = 1, nExcludes_global if ((excludesGlobal(i) == unique_id_1) .or. & @@ -715,7 +1118,7 @@ contains return endif enddo - + do i = 1, nExcludes_local if (excludesLocal(1,i) == unique_id_1) then if (excludesLocal(2,i) == unique_id_2) then @@ -751,4 +1154,13 @@ end module do_Forces doesit = FF_uses_RF end function FF_RequiresPostpairCalc +#ifdef PROFILE + function getforcetime() result(totalforcetime) + real(kind=dp) :: totalforcetime + totalforcetime = forcetime + end function getforcetime +#endif + +!! This cleans componets of force arrays belonging only to fortran + end module do_Forces