ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/do_Forces.F90
(Generate patch)

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents):
Revision 891 by mmeineke, Fri Dec 19 20:36:35 2003 UTC vs.
Revision 898 by chuckv, Mon Jan 5 22:49:14 2004 UTC

# Line 4 | Line 4
4  
5   !! @author Charles F. Vardeman II
6   !! @author Matthew Meineke
7 < !! @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 $
7 > !! @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 $
8  
9   module do_Forces
10    use force_globals
# Line 52 | Line 52 | contains
52    integer :: nLoops
53   #endif
54  
55 +  logical, allocatable :: propertyMapI(:,:)
56 +  logical, allocatable :: propertyMapJ(:,:)
57 +
58   contains
59  
60    subroutine setRlistDF( this_rlist )
# Line 175 | Line 178 | contains
178      endif
179      if (.not. do_forces_initialized) then
180         !! Create neighbor lists
181 <       call expandNeighborList(getNlocal(), my_status)
181 >       call expandNeighborList(nLocal, my_status)
182         if (my_Status /= 0) then
183            write(default_error,*) "SimSetup: ExpandNeighborList returned error."
184            thisStat = -1
# Line 195 | Line 198 | contains
198    subroutine do_force_loop(q, A, u_l, f, t, tau, pot, do_pot_c, do_stress_c, &
199         error)
200      !! Position array provided by C, dimensioned by getNlocal
201 <    real ( kind = dp ), dimension(3,getNlocal()) :: q
201 >    real ( kind = dp ), dimension(3,nLocal) :: q
202      !! Rotation Matrix for each long range particle in simulation.
203 <    real( kind = dp), dimension(9,getNlocal()) :: A    
203 >    real( kind = dp), dimension(9,nLocal) :: A    
204      !! Unit vectors for dipoles (lab frame)
205 <    real( kind = dp ), dimension(3,getNlocal()) :: u_l
205 >    real( kind = dp ), dimension(3,nLocal) :: u_l
206      !! Force array provided by C, dimensioned by getNlocal
207 <    real ( kind = dp ), dimension(3,getNlocal()) :: f
207 >    real ( kind = dp ), dimension(3,nLocal) :: f
208      !! Torsion array provided by C, dimensioned by getNlocal
209 <    real( kind = dp ), dimension(3,getNlocal()) :: t    
209 >    real( kind = dp ), dimension(3,nLocal) :: t    
210 >
211      !! Stress Tensor
212      real( kind = dp), dimension(9) :: tau  
213      real ( kind = dp ) :: pot
# Line 216 | Line 220 | contains
220      integer :: ncol
221      integer :: nprocs
222   #endif
219 –    integer :: nlocal
223      integer :: natoms    
224      logical :: update_nlist  
225      integer :: i, j, jbeg, jend, jnab
# Line 224 | Line 227 | contains
227      real( kind = DP ) ::  rijsq
228      real(kind=dp),dimension(3) :: d
229      real(kind=dp) :: rfpot, mu_i, virial
230 <    integer :: me_i
230 >    integer :: me_i, me_j
231      logical :: is_dp_i
232      integer :: neighborListSize
233      integer :: listerror, error
234      integer :: localError
235 +    integer :: propPack_i, propPack_j
236  
237      real(kind=dp) :: listSkin = 1.0  
238  
# Line 236 | Line 240 | contains
240  
241   #ifdef IS_MPI
242      pot_local = 0.0_dp
239 –    nlocal = getNlocal()
243      nrow   = getNrow(plan_row)
244      ncol   = getNcol(plan_col)
245   #else
243 –    nlocal = getNlocal()
246      natoms = nlocal
247   #endif
248  
# Line 255 | Line 257 | contains
257      do_pot = do_pot_c
258      do_stress = do_stress_c
259  
260 <
260 >
261 > #ifdef IS_MPI
262 >    if (.not.allocated(propertyMapI)) then
263 >       allocate(propertyMapI(5,nrow))
264 >    endif
265 >
266 >    do i = 1, nrow
267 >       me_i = atid_row(i)
268 > #else
269 >    if (.not.allocated(propertyMapI)) then
270 >       allocate(propertyMapI(5,nlocal))
271 >    endif
272 >
273 >    do i = 1, natoms
274 >       me_i = atid(i)
275 > #endif
276 >      
277 >       propertyMapI(1:5,i) = .false.
278 >
279 >       call getElementProperty(atypes, me_i, "propertyPack", propPack_i)
280 >    
281 >       ! unpack the properties
282 >      
283 >       if (iand(propPack_i, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) &
284 >            propertyMapI(1, i) = .true.
285 >       if (iand(propPack_i, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) &
286 >            propertyMapI(2, i) = .true.
287 >       if (iand(propPack_i, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) &
288 >            propertyMapI(3, i) = .true.
289 >       if (iand(propPack_i, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) &
290 >            propertyMapI(4, i) = .true.
291 >       if (iand(propPack_i, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) &
292 >            propertyMapI(5, i) = .true.
293 >
294 >    end do
295 >
296 > #ifdef IS_MPI
297 >    if (.not.allocated(propertyMapJ)) then
298 >       allocate(propertyMapJ(5,ncol))
299 >    endif
300 >
301 >    do j = 1, ncol
302 >       me_j = atid_col(j)
303 > #else
304 >    if (.not.allocated(propertyMapJ)) then
305 >       allocate(propertyMapJ(5,nlocal))
306 >    endif
307 >
308 >    do j = 1, natoms
309 >       me_j = atid(j)
310 > #endif
311 >      
312 >       propertyMapJ(1:5,j) = .false.
313 >
314 >       call getElementProperty(atypes, me_j, "propertyPack", propPack_j)
315 >    
316 >       ! unpack the properties
317 >      
318 >       if (iand(propPack_j, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) &
319 >            propertyMapJ(1, j) = .true.
320 >       if (iand(propPack_j, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) &
321 >            propertyMapJ(2, j) = .true.
322 >       if (iand(propPack_j, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) &
323 >            propertyMapJ(3, j) = .true.
324 >       if (iand(propPack_j, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) &
325 >            propertyMapJ(4, j) = .true.
326 >       if (iand(propPack_j, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) &
327 >            propertyMapJ(5, j) = .true.
328 >
329 >    end do
330 >
331      ! Gather all information needed by all force loops:
332      
333   #ifdef IS_MPI    
# Line 451 | Line 523 | contains
523         nlist = 0      
524        
525         do i = 1, nrow
526 +
527            point(i) = nlist + 1
528            
529            inner: do j = 1, ncol
# Line 645 | Line 718 | contains
718            end do
719   #endif
720            
721 <          do i = 1, getNlocal()
721 >          do i = 1, nLocal
722  
723               rfpot = 0.0_DP
724   #ifdef IS_MPI
# Line 705 | Line 778 | contains
778    subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
779  
780      real( kind = dp ) :: pot
781 <    real( kind = dp ), dimension(3,getNlocal()) :: u_l
782 <    real (kind=dp), dimension(9,getNlocal()) :: A
783 <    real (kind=dp), dimension(3,getNlocal()) :: f
784 <    real (kind=dp), dimension(3,getNlocal()) :: t
781 >    real( kind = dp ), dimension(3,nLocal) :: u_l
782 >    real (kind=dp), dimension(9,nLocal) :: A
783 >    real (kind=dp), dimension(3,nLocal) :: f
784 >    real (kind=dp), dimension(3,nLocal) :: t
785  
786      logical, intent(inout) :: do_pot, do_stress
787      integer, intent(in) :: i, j
# Line 721 | Line 794 | contains
794      logical :: is_EAM_i,is_EAM_j
795      logical :: is_Sticky_i, is_Sticky_j
796      integer :: me_i, me_j
797 <
797 >    integer :: propPack_i
798 >    integer :: propPack_j
799      r = sqrt(rijsq)
800  
801   #ifdef IS_MPI
# Line 738 | Line 812 | contains
812      me_j = atid(j)
813  
814   #endif
815 <
815 >    
816      if (FF_uses_LJ .and. SimUsesLJ()) then
743 –       call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i)
744 –       call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j)
817  
818 <       if ( is_LJ_i .and. is_LJ_j ) &
818 >       if ( propertyMapI(1, me_i) .and. propertyMapJ(1, me_j) ) &
819              call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
820 +
821      endif
822  
823      if (FF_uses_dipoles .and. SimUsesDipoles()) then
751 –       call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
752 –       call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
824        
825 <       if ( is_DP_i .and. is_DP_j ) then
825 >       if ( propertyMapI(2, me_i) .and. propertyMapJ(2, me_j)) then
826            call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
827                 do_pot, do_stress)
828            if (FF_uses_RF .and. SimUsesRF()) then
# Line 764 | Line 835 | contains
835  
836      if (FF_uses_Sticky .and. SimUsesSticky()) then
837  
838 <       call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i)
768 <       call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j)
769 <
770 <       if ( is_Sticky_i .and. is_Sticky_j ) then
838 >       if ( propertyMapI(3, me_i) .and. propertyMapJ(3, me_j)) then
839            call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
840                 do_pot, do_stress)
841         endif
# Line 775 | Line 843 | contains
843  
844  
845      if (FF_uses_GB .and. SimUsesGB()) then
778 –
779 –
780 –       call getElementProperty(atypes, me_i, "is_GB", is_GB_i)
781 –       call getElementProperty(atypes, me_j, "is_GB", is_GB_j)
846        
847 <       if ( is_GB_i .and. is_GB_j ) then
847 >       if ( propertyMapI(4, me_i) .and. propertyMapJ(4, me_j)) then
848            call do_gb_pair(i, j, d, r, rijsq, u_l, pot, f, t, &
849                 do_pot, do_stress)          
850         endif
851 +
852      endif
853      
854  
855    
856     if (FF_uses_EAM .and. SimUsesEAM()) then
792 –      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
793 –      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
857        
858 <      if ( is_EAM_i .and. is_EAM_j ) &
859 <           call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
860 <   endif
798 <
858 >      if ( propertyMapI(5, me_i) .and. propertyMapJ(5, me_j)) then
859 >         call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
860 >      endif
861  
862 <
862 >   endif
863  
864    end subroutine do_pair
865  
# Line 805 | Line 867 | contains
867  
868    subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
869     real( kind = dp ) :: pot
870 <   real( kind = dp ), dimension(3,getNlocal()) :: u_l
871 <   real (kind=dp), dimension(9,getNlocal()) :: A
872 <   real (kind=dp), dimension(3,getNlocal()) :: f
873 <   real (kind=dp), dimension(3,getNlocal()) :: t
870 >   real( kind = dp ), dimension(3,nLocal) :: u_l
871 >   real (kind=dp), dimension(9,nLocal) :: A
872 >   real (kind=dp), dimension(3,nLocal) :: f
873 >   real (kind=dp), dimension(3,nLocal) :: t
874    
875     logical, intent(inout) :: do_pot, do_stress
876     integer, intent(in) :: i, j

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines