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 626 by mmeineke, Wed Jul 16 21:30:56 2003 UTC vs.
Revision 648 by chuckv, Wed Jul 23 22:13:59 2003 UTC

# Line 4 | Line 4
4  
5   !! @author Charles F. Vardeman II
6   !! @author Matthew Meineke
7 < !! @version $Id: do_Forces.F90,v 1.21 2003-07-16 21:30:55 mmeineke Exp $, $Date: 2003-07-16 21:30:55 $, $Name: not supported by cvs2svn $, $Revision: 1.21 $
7 > !! @version $Id: do_Forces.F90,v 1.24 2003-07-23 22:13:59 chuckv Exp $, $Date: 2003-07-23 22:13:59 $, $Name: not supported by cvs2svn $, $Revision: 1.24 $
8  
9   module do_Forces
10    use force_globals
# Line 256 | Line 256 | contains
256         !! do_prepair_loop_if_needed
257         !! if_mpi_scatter_stuff_from_prepair
258         !! if_mpi_gather_stuff_from_prepair_to_main_loop
259 +
260 + !--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
261 + #ifdef IS_MPI
262 +    
263 +    if (update_nlist) then
264 +      
265 +       !! save current configuration, construct neighbor list,
266 +       !! and calculate forces
267 +       call saveNeighborList(nlocal, q)
268 +      
269 +       neighborListSize = size(list)
270 +       nlist = 0      
271 +      
272 +       do i = 1, nrow
273 +          point(i) = nlist + 1
274 +          
275 +          prepair_inner: do j = 1, ncol
276 +            
277 +             if (skipThisPair(i,j)) cycle prepair_inner
278 +            
279 +             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
280 +            
281 +             if (rijsq < rlistsq) then            
282 +                
283 +                nlist = nlist + 1
284 +                
285 +                if (nlist > neighborListSize) then
286 +                   call expandNeighborList(nlocal, listerror)
287 +                   if (listerror /= 0) then
288 +                      error = -1
289 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
290 +                      return
291 +                   end if
292 +                   neighborListSize = size(list)
293 +                endif
294 +                
295 +                list(nlist) = j
296 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)                      
297 +             endif
298 +          enddo prepair_inner
299 +       enddo
300 +
301 +       point(nrow + 1) = nlist + 1
302 +      
303 +    else  !! (of update_check)
304 +
305 +       ! use the list to find the neighbors
306 +       do i = 1, nrow
307 +          JBEG = POINT(i)
308 +          JEND = POINT(i+1) - 1
309 +          ! check thiat molecule i has neighbors
310 +          if (jbeg .le. jend) then
311 +            
312 +             do jnab = jbeg, jend
313 +                j = list(jnab)
314 +
315 +                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
316 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
317 +                     u_l, A, f, t, pot_local)
318 +
319 +             enddo
320 +          endif
321 +       enddo
322 +    endif
323 +    
324 + #else
325 +    
326 +    if (update_nlist) then
327 +      
328 +       ! save current configuration, contruct neighbor list,
329 +       ! and calculate forces
330 +       call saveNeighborList(natoms, q)
331 +      
332 +       neighborListSize = size(list)
333 +  
334 +       nlist = 0
335 +      
336 +       do i = 1, natoms-1
337 +          point(i) = nlist + 1
338 +          
339 +          prepair_inner: do j = i+1, natoms
340 +            
341 +             if (skipThisPair(i,j))  cycle prepair_inner
342 +                          
343 +             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
344 +          
345 +
346 +             if (rijsq < rlistsq) then
347 +                
348 +                nlist = nlist + 1
349 +              
350 +                if (nlist > neighborListSize) then
351 +                   call expandNeighborList(natoms, listerror)
352 +                   if (listerror /= 0) then
353 +                      error = -1
354 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
355 +                      return
356 +                   end if
357 +                   neighborListSize = size(list)
358 +                endif
359 +                
360 +                list(nlist) = j
361 +                
362 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
363 +                        u_l, A, f, t, pot)
364 +                
365 +             endif
366 +          enddo prepair_inner
367 +       enddo
368 +      
369 +       point(natoms) = nlist + 1
370 +      
371 +    else !! (update)
372 +      
373 +       ! use the list to find the neighbors
374 +       do i = 1, natoms-1
375 +          JBEG = POINT(i)
376 +          JEND = POINT(i+1) - 1
377 +          ! check thiat molecule i has neighbors
378 +          if (jbeg .le. jend) then
379 +            
380 +             do jnab = jbeg, jend
381 +                j = list(jnab)
382 +
383 +                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
384 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
385 +                     u_l, A, f, t, pot)
386 +
387 +             enddo
388 +          endif
389 +       enddo
390 +    endif    
391 + #endif
392 +    !! Do rest of preforce calculations
393 +   call do_preforce(nlocal,pot)
394      else
395 <       !! See if we need to update neighbor lists
395 >       !! See if we need to update neighbor lists for non pre-pair
396         call checkNeighborList(nlocal, q, listSkin, update_nlist)  
397      endif
398 <    
398 >
399 >
400 >
401 >
402 >
403 > !---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>
404 >
405 >
406 >
407 >
408 >  
409   #ifdef IS_MPI
410      
411      if (update_nlist) then
# Line 531 | Line 676 | contains
676      logical :: is_LJ_i, is_LJ_j
677      logical :: is_DP_i, is_DP_j
678      logical :: is_GB_i, is_GB_j
679 +    logical :: is_EAM_i,is_EAM_j
680      logical :: is_Sticky_i, is_Sticky_j
681      integer :: me_i, me_j
682  
# Line 599 | Line 745 | contains
745      endif
746      
747  
748 +  
749 +   if (FF_uses_EAM .and. SimUsesEAM()) then
750 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
751 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
752 +      
753 +      if ( is_EAM_i .and. is_EAM_j ) &
754 +           call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
755 +   endif
756 +
757  
758 +
759 +
760    end subroutine do_pair
761 +
762 +
763 +
764 +  subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
765 +   real( kind = dp ) :: pot
766 +   real( kind = dp ), dimension(3,getNlocal()) :: u_l
767 +   real (kind=dp), dimension(9,getNlocal()) :: A
768 +   real (kind=dp), dimension(3,getNlocal()) :: f
769 +   real (kind=dp), dimension(3,getNlocal()) :: t
770 +  
771 +   logical, intent(inout) :: do_pot, do_stress
772 +   integer, intent(in) :: i, j
773 +   real ( kind = dp ), intent(inout)    :: rijsq
774 +   real ( kind = dp )                :: r
775 +   real ( kind = dp ), intent(inout) :: d(3)
776 +  
777 +   logical :: is_EAM_i, is_EAM_j
778 +  
779 +   integer :: me_i, me_j
780 +  
781 +   r = sqrt(rijsq)
782 +  
783 + #ifdef IS_MPI
784 +   if (tagRow(i) .eq. tagColumn(j)) then
785 +      write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
786 +   endif
787 +  
788 +   me_i = atid_row(i)
789 +   me_j = atid_col(j)
790 +  
791 + #else
792 +  
793 +   me_i = atid(i)
794 +   me_j = atid(j)
795 +  
796 + #endif
797 +  
798 +   if (FF_uses_EAM .and. SimUsesEAM()) then
799 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
800 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
801 +      
802 +      if ( is_EAM_i .and. is_EAM_j ) &
803 +           call calc_EAM_prepair_rho(i, j, d, r, rijsq )
804 +   endif
805 +  end subroutine do_prepair
806  
807  
808 +
809 +
810 +  subroutine do_preforce(nlocal,pot)
811 +    integer :: nlocal
812 +    real( kind = dp ) :: pot
813 +
814 +   if (FF_uses_EAM .and. SimUsesEAM()) then
815 +      call calc_EAM_preforce_Frho(nlocal,pot)
816 +   endif
817 +
818 +
819 +  end subroutine do_preforce
820 +  
821 +  
822    subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
823      
824      real (kind = dp), dimension(3) :: q_i

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines