--- trunk/OOPSE/libmdtools/do_Forces.F90 2003/07/16 21:30:56 626 +++ trunk/OOPSE/libmdtools/do_Forces.F90 2003/08/26 20:37:30 726 @@ -4,7 +4,7 @@ !! @author Charles F. Vardeman II !! @author Matthew Meineke -!! @version $Id: do_Forces.F90,v 1.21 2003-07-16 21:30:55 mmeineke Exp $, $Date: 2003-07-16 21:30:55 $, $Name: not supported by cvs2svn $, $Revision: 1.21 $ +!! @version $Id: do_Forces.F90,v 1.30 2003-08-26 20:37:30 tim Exp $, $Date: 2003-08-26 20:37:30 $, $Name: not supported by cvs2svn $, $Revision: 1.30 $ module do_Forces use force_globals @@ -18,6 +18,8 @@ module do_Forces use reaction_field use gb_pair use vector_class + use eam + use status #ifdef IS_MPI use mpiSimulation #endif @@ -43,6 +45,14 @@ contains public :: do_force_loop public :: setRlistDF +#ifdef PROFILE + real(kind = dp) :: forceTime + real(kind = dp) :: forceTimeInitial, forceTimeFinal + real(kind = dp) :: globalForceTime + real(kind = dp) :: maxForceTime + integer, save :: nloops = 0 +#endif + contains subroutine setRlistDF( this_rlist ) @@ -136,11 +146,22 @@ contains if (FF_uses_sticky) then call check_sticky_FF(my_status) + if (my_status /= 0) then + thisStat = -1 + return + end if + endif + + + if (FF_uses_EAM) then + call init_EAM_FF(my_status) if (my_status /= 0) then thisStat = -1 return end if endif + + if (FF_uses_GB) then call check_gb_pair_FF(my_status) @@ -161,10 +182,11 @@ contains return endif endif + havePolicies = .true. if( haveRlist ) do_forces_initialized = .true. - + end subroutine init_FF @@ -192,6 +214,7 @@ contains real( kind = DP ) :: pot_local integer :: nrow integer :: ncol + integer :: nprocs #endif integer :: nlocal integer :: natoms @@ -221,9 +244,10 @@ contains nlocal = getNlocal() natoms = nlocal #endif - + call check_initialization(localError) if ( localError .ne. 0 ) then + call handleError("do_force_loop","Not Initialized") error = -1 return end if @@ -232,6 +256,7 @@ contains do_pot = do_pot_c do_stress = do_stress_c + ! Gather all information needed by all force loops: #ifdef IS_MPI @@ -248,7 +273,13 @@ 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, listSkin, update_nlist) @@ -256,11 +287,160 @@ contains !! do_prepair_loop_if_needed !! if_mpi_scatter_stuff_from_prepair !! if_mpi_gather_stuff_from_prepair_to_main_loop + +!--------------------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 + +#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 + !! 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 @@ -398,7 +578,14 @@ 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 @@ -512,7 +699,39 @@ contains endif #endif - + +#ifdef PROFILE + if (do_pot) then + +#ifdef IS_MPI + + + call printCommTime() + + call mpi_allreduce(forceTime,globalForceTime,1,MPI_DOUBLE_PRECISION, & + mpi_sum,mpi_comm_world,mpi_err) + + call mpi_allreduce(forceTime,maxForceTime,1,MPI_DOUBLE_PRECISION, & + MPI_MAX,mpi_comm_world,mpi_err) + + call mpi_comm_size( MPI_COMM_WORLD, nprocs,mpi_err) + + if (getMyNode() == 0) then + write(*,*) "Total processor time spent in force calculations is: ", globalForceTime + write(*,*) "Total Time spent in force loop per processor is: ", globalforceTime/nprocs + write(*,*) "Maximum force time on any processor is: ", maxForceTime + end if +#else + write(*,*) "Time spent in force loop is: ", forceTime +#endif + + + endif + +#endif + + + end subroutine do_force_loop subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot) @@ -531,6 +750,7 @@ 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 @@ -564,7 +784,6 @@ contains call getElementProperty(atypes, me_j, "is_DP", is_DP_j) if ( is_DP_i .and. is_DP_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 @@ -589,6 +808,7 @@ 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) @@ -598,11 +818,83 @@ contains endif 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 do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress) + 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 @@ -660,6 +952,7 @@ contains error = 0 ! Make sure we are properly initialized. if (.not. do_forces_initialized) then + write(*,*) "Forces not initialized" error = -1 return endif @@ -707,6 +1000,15 @@ 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 @@ -819,4 +1121,6 @@ end module do_Forces doesit = FF_uses_RF end function FF_RequiresPostpairCalc +!! This cleans componets of force arrays belonging only to fortran + end module do_Forces