--- trunk/OOPSE/libmdtools/do_Forces.F90 2003/07/02 21:26:55 572 +++ 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.17 2003-07-02 21:26:55 mmeineke Exp $, $Date: 2003-07-02 21:26:55 $, $Name: not supported by cvs2svn $, $Revision: 1.17 $ +!! @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 @@ -150,8 +162,9 @@ contains endif endif - do_forces_initialized = .true. - + havePolicies = .true. + if( haveRlist ) do_forces_initialized = .true. + end subroutine init_FF @@ -185,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 @@ -194,6 +207,9 @@ contains integer :: listerror, error integer :: localError + real(kind=dp) :: listSkin = 1.0 + + !! initialize local variables #ifdef IS_MPI @@ -206,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 @@ -238,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 @@ -268,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 @@ -284,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 @@ -337,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 @@ -353,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 @@ -525,8 +536,6 @@ 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) @@ -589,9 +598,56 @@ 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