--- trunk/OOPSE/libmdtools/do_Forces.F90 2003/03/31 21:50:59 438 +++ trunk/OOPSE/libmdtools/do_Forces.F90 2003/07/17 19:25:51 631 @@ -4,7 +4,7 @@ !! @author Charles F. Vardeman II !! @author Matthew Meineke -!! @version $Id: do_Forces.F90,v 1.5 2003-03-31 21:50:59 chuckv Exp $, $Date: 2003-03-31 21:50:59 $, $Name: not supported by cvs2svn $, $Revision: 1.5 $ +!! @version $Id: do_Forces.F90,v 1.22 2003-07-17 19:25:51 chuckv Exp $, $Date: 2003-07-17 19:25:51 $, $Name: not supported by cvs2svn $, $Revision: 1.22 $ 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 @@ -167,7 +188,7 @@ contains 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 @@ -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,9 +207,13 @@ contains integer :: listerror, error integer :: localError + 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) @@ -196,10 +221,7 @@ contains nlocal = getNlocal() natoms = nlocal #endif - - call getRcut(rcut,rc2=rcutsq) - call getRlist(rlist,rlistsq) - + call check_initialization(localError) if ( localError .ne. 0 ) then error = -1 @@ -229,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 @@ -245,7 +267,7 @@ contains !! save current configuration, construct neighbor list, !! and calculate forces - call saveNeighborList(q) + call saveNeighborList(nlocal, q) neighborListSize = size(list) nlist = 0 @@ -259,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 @@ -275,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) - 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 +320,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 @@ -312,7 +333,7 @@ contains ! save current configuration, contruct neighbor list, ! and calculate forces - call saveNeighborList(q) + call saveNeighborList(natoms, q) neighborListSize = size(list) @@ -328,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 @@ -344,10 +365,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 @@ -368,7 +388,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 @@ -411,20 +431,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 @@ -472,16 +492,15 @@ contains #ifdef IS_MPI if (do_pot) then - write(*,*) "Fortran is on pot:, pot, pot_local ", pot,pot_local - 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 @@ -496,13 +515,13 @@ contains 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 @@ -518,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) @@ -576,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 @@ -655,7 +755,6 @@ contains 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) @@ -699,7 +798,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 @@ -709,14 +808,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. & @@ -725,7 +828,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