--- trunk/OOPSE/libmdtools/do_Forces.F90 2003/12/19 20:36:35 891 +++ trunk/OOPSE/libmdtools/do_Forces.F90 2004/01/05 22:49:14 898 @@ -4,7 +4,7 @@ !! @author Charles F. Vardeman II !! @author Matthew Meineke -!! @version $Id: do_Forces.F90,v 1.40 2003-12-19 20:36:35 mmeineke Exp $, $Date: 2003-12-19 20:36:35 $, $Name: not supported by cvs2svn $, $Revision: 1.40 $ +!! @version $Id: do_Forces.F90,v 1.44 2004-01-05 22:49:14 chuckv Exp $, $Date: 2004-01-05 22:49:14 $, $Name: not supported by cvs2svn $, $Revision: 1.44 $ module do_Forces use force_globals @@ -52,6 +52,9 @@ contains integer :: nLoops #endif + logical, allocatable :: propertyMapI(:,:) + logical, allocatable :: propertyMapJ(:,:) + contains subroutine setRlistDF( this_rlist ) @@ -175,7 +178,7 @@ contains endif if (.not. do_forces_initialized) then !! Create neighbor lists - call expandNeighborList(getNlocal(), my_status) + call expandNeighborList(nLocal, my_status) if (my_Status /= 0) then write(default_error,*) "SimSetup: ExpandNeighborList returned error." thisStat = -1 @@ -195,15 +198,16 @@ contains subroutine do_force_loop(q, A, u_l, f, t, tau, pot, do_pot_c, do_stress_c, & error) !! Position array provided by C, dimensioned by getNlocal - real ( kind = dp ), dimension(3,getNlocal()) :: q + real ( kind = dp ), dimension(3,nLocal) :: q !! Rotation Matrix for each long range particle in simulation. - real( kind = dp), dimension(9,getNlocal()) :: A + real( kind = dp), dimension(9,nLocal) :: A !! Unit vectors for dipoles (lab frame) - real( kind = dp ), dimension(3,getNlocal()) :: u_l + real( kind = dp ), dimension(3,nLocal) :: u_l !! Force array provided by C, dimensioned by getNlocal - real ( kind = dp ), dimension(3,getNlocal()) :: f + real ( kind = dp ), dimension(3,nLocal) :: f !! Torsion array provided by C, dimensioned by getNlocal - real( kind = dp ), dimension(3,getNlocal()) :: t + real( kind = dp ), dimension(3,nLocal) :: t + !! Stress Tensor real( kind = dp), dimension(9) :: tau real ( kind = dp ) :: pot @@ -216,7 +220,6 @@ contains integer :: ncol integer :: nprocs #endif - integer :: nlocal integer :: natoms logical :: update_nlist integer :: i, j, jbeg, jend, jnab @@ -224,11 +227,12 @@ contains 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 @@ -236,11 +240,9 @@ contains #ifdef IS_MPI pot_local = 0.0_dp - nlocal = getNlocal() nrow = getNrow(plan_row) ncol = getNcol(plan_col) #else - nlocal = getNlocal() natoms = nlocal #endif @@ -255,7 +257,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 @@ -451,6 +523,7 @@ contains nlist = 0 do i = 1, nrow + point(i) = nlist + 1 inner: do j = 1, ncol @@ -645,7 +718,7 @@ contains end do #endif - do i = 1, getNlocal() + do i = 1, nLocal rfpot = 0.0_DP #ifdef IS_MPI @@ -705,10 +778,10 @@ contains 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(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 + real( kind = dp ), dimension(3,nLocal) :: u_l + real (kind=dp), dimension(9,nLocal) :: A + real (kind=dp), dimension(3,nLocal) :: f + real (kind=dp), dimension(3,nLocal) :: t logical, intent(inout) :: do_pot, do_stress integer, intent(in) :: i, j @@ -721,7 +794,8 @@ contains 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 @@ -738,20 +812,17 @@ 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 @@ -764,10 +835,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 @@ -775,29 +843,23 @@ 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 - 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 - + 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 @@ -805,10 +867,10 @@ contains 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 + real( kind = dp ), dimension(3,nLocal) :: u_l + real (kind=dp), dimension(9,nLocal) :: A + real (kind=dp), dimension(3,nLocal) :: f + real (kind=dp), dimension(3,nLocal) :: t logical, intent(inout) :: do_pot, do_stress integer, intent(in) :: i, j