--- trunk/OOPSE/libmdtools/do_Forces.F90 2003/04/07 20:50:46 470 +++ trunk/OOPSE/libmdtools/do_Forces.F90 2003/07/17 19:38:23 634 @@ -4,7 +4,7 @@ !! @author Charles F. Vardeman II !! @author Matthew Meineke -!! @version $Id: do_Forces.F90,v 1.11 2003-04-07 20:50:46 chuckv Exp $, $Date: 2003-04-07 20:50:46 $, $Name: not supported by cvs2svn $, $Revision: 1.11 $ +!! @version $Id: do_Forces.F90,v 1.23 2003-07-17 19:38:23 mmeineke Exp $, $Date: 2003-07-17 19:38:23 $, $Name: not supported by cvs2svn $, $Revision: 1.23 $ module do_Forces use force_globals @@ -17,6 +17,7 @@ module do_Forces use dipole_dipole use reaction_field use gb_pair + use vector_class #ifdef IS_MPI use mpiSimulation #endif @@ -27,7 +28,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 +37,26 @@ 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 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 +104,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 +114,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 @@ -140,10 +152,19 @@ 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 - - do_forces_initialized = .true. - + havePolicies = .true. + if( haveRlist ) do_forces_initialized = .true. + end subroutine init_FF @@ -177,7 +198,7 @@ contains 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 @@ -186,6 +207,9 @@ contains integer :: listerror, error integer :: localError + real(kind=dp) :: listSkin = 1.0 + + !! initialize local variables #ifdef IS_MPI @@ -198,9 +222,6 @@ contains natoms = nlocal #endif - call getRcut(rcut,rc2=rcutsq) - call getRlist(rlist,rlistsq) - call check_initialization(localError) if ( localError .ne. 0 ) then error = -1 @@ -210,7 +231,6 @@ contains do_pot = do_pot_c do_stress = do_stress_c - ! Gather all information needed by all force loops: @@ -231,14 +251,14 @@ contains 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) + call checkNeighborList(nlocal, q, listSkin, update_nlist) endif #ifdef IS_MPI @@ -261,7 +281,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 @@ -277,10 +297,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_local) - endif + call do_pair(i, j, rijsq, d, do_pot, do_stress, & + u_l, A, f, t, pot_local) + endif enddo inner enddo @@ -330,7 +349,7 @@ contains call get_interatomic_vector(q(:,i), q(:,j), d, rijsq) - if (rijsq < rlistsq) then + if (rijsq < rlistsq) then nlist = nlist + 1 @@ -346,10 +365,9 @@ contains list(nlist) = j - if (rijsq < rcutsq) then - call do_pair(i, j, rijsq, d, do_pot, do_stress, & + call do_pair(i, j, rijsq, d, do_pot, do_stress, & u_l, A, f, t, pot) - endif + endif enddo inner enddo @@ -480,7 +498,7 @@ contains endif if (do_stress) then - call mpi_allreduce(tau_Temp, tau,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_Temp, virial,1,mpi_double_precision,mpi_sum, & mpi_comm_world,mpi_err) @@ -519,6 +537,9 @@ contains 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) @@ -577,29 +598,107 @@ contains endif endif + + end subroutine do_pair + + subroutine do_preforce(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(i, j, d, r, rijsq ) + 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 - 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) + if( .not.boxIsOrthorhombic ) then + ! calc the scaled coordinates. + + scaled = matmul(HmatInv, d) + + ! wrap the scaled coordinates + + 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