--- trunk/OOPSE/libmdtools/do_Forces.F90 2004/01/05 21:00:05 894 +++ trunk/OOPSE/libmdtools/do_Forces.F90 2004/01/05 22:12:11 895 @@ -4,7 +4,7 @@ !! @author Charles F. Vardeman II !! @author Matthew Meineke -!! @version $Id: do_Forces.F90,v 1.41 2004-01-05 21:00:05 chuckv Exp $, $Date: 2004-01-05 21:00:05 $, $Name: not supported by cvs2svn $, $Revision: 1.41 $ +!! @version $Id: do_Forces.F90,v 1.42 2004-01-05 22:12:11 chuckv Exp $, $Date: 2004-01-05 22:12:11 $, $Name: not supported by cvs2svn $, $Revision: 1.42 $ module do_Forces use force_globals @@ -51,6 +51,9 @@ module do_Forces real :: forceTimeInitial, forceTimeFinal integer :: nLoops #endif + + logical, allocatable :: propertyMapI(:,:) + logical, allocatable :: propertyMapJ(:,:) contains @@ -204,6 +207,7 @@ contains real ( kind = dp ), dimension(3,getNlocal()) :: f !! Torsion array provided by C, dimensioned by getNlocal real( kind = dp ), dimension(3,getNlocal()) :: t + !! Stress Tensor real( kind = dp), dimension(9) :: tau real ( kind = dp ) :: pot @@ -255,7 +259,77 @@ contains do_pot = do_pot_c do_stress = do_stress_c - + +#ifdef IS_MPI + if (.not.allocated(propertyMapI)) then + allocate(propertyMapI(5,getNrow()) + endif + + do i = 1, nrow + me_i = atid_row(i) +#else + if (.not.allocated(propertyMapI)) then + allocate(propertyMapI(5,getNlocal()) + 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,getNcol()) + endif + + do j = 1, ncol + me_j = atid_col(j) +#else + if (.not.allocated(propertyMapJ)) then + allocate(propertyMapJ(5,getNlocal()) + 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 +525,7 @@ contains nlist = 0 do i = 1, nrow + point(i) = nlist + 1 inner: do j = 1, ncol @@ -739,35 +814,17 @@ contains me_j = atid(j) #endif - - call getElementProperty(atypes, me_i, "propertyPack", propPack_i) - call getElementProperty(atypes, me_j, "propertyPack", propPack_j) - ! unpack the properties - - if (iand(propPack_i, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) is_LJ_i = .true. - if (iand(propPack_i, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) is_DP_i = .true. - if (iand(propPack_i, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) is_Sticky_i = .true. - if (iand(propPack_i, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) is_GB_i = .true. - if (iand(propPack_i, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) is_EAM_i = .true. - - if (iand(propPack_j, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) is_LJ_j = .true. - if (iand(propPack_j, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) is_DP_j = .true. - if (iand(propPack_j, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) is_Sticky_j = .true. - if (iand(propPack_j, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) is_GB_j = .true. - if (iand(propPack_j, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) is_EAM_j = .true. - - if (FF_uses_LJ .and. SimUsesLJ()) then - 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 - 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 @@ -780,7 +837,7 @@ contains if (FF_uses_Sticky .and. SimUsesSticky()) then - 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 @@ -789,7 +846,7 @@ contains if (FF_uses_GB .and. SimUsesGB()) then - 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 @@ -800,8 +857,9 @@ contains if (FF_uses_EAM .and. SimUsesEAM()) then - if ( is_EAM_i .and. is_EAM_j ) & - call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress) + 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