--- trunk/OOPSE/libmdtools/do_Forces.F90	2003/07/17 19:38:23	634
+++ trunk/OOPSE/libmdtools/do_Forces.F90	2003/07/23 22:13:59	648
@@ -4,7 +4,7 @@
 
 !! @author Charles F. Vardeman II
 !! @author Matthew Meineke
-!! @version $Id: do_Forces.F90,v 1.23 2003-07-17 19:38:23 mmeineke Exp $, $Date: 2003-07-17 19:38:23 $, $Name: not supported by cvs2svn $, $Revision: 1.23 $
+!! @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 $
 
 module do_Forces
   use force_globals
@@ -256,11 +256,156 @@ contains
        !! do_prepair_loop_if_needed
        !! if_mpi_scatter_stuff_from_prepair
        !! if_mpi_gather_stuff_from_prepair_to_main_loop
+
+!--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
+#ifdef IS_MPI
+    
+    if (update_nlist) then
+       
+       !! save current configuration, construct neighbor list, 
+       !! and calculate forces
+       call saveNeighborList(nlocal, q)
+       
+       neighborListSize = size(list)
+       nlist = 0       
+       
+       do i = 1, nrow
+          point(i) = nlist + 1
+          
+          prepair_inner: do j = 1, ncol
+             
+             if (skipThisPair(i,j)) cycle prepair_inner
+             
+             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
+             
+             if (rijsq < rlistsq) then             
+                
+                nlist = nlist + 1
+                
+                if (nlist > neighborListSize) then
+                   call expandNeighborList(nlocal, listerror)
+                   if (listerror /= 0) then
+                      error = -1
+                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
+                      return
+                   end if
+                   neighborListSize = size(list)
+                endif
+                
+                list(nlist) = j
+                call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)                      
+             endif
+          enddo prepair_inner
+       enddo
+
+       point(nrow + 1) = nlist + 1
+       
+    else  !! (of update_check)
+
+       ! use the list to find the neighbors
+       do i = 1, nrow
+          JBEG = POINT(i)
+          JEND = POINT(i+1) - 1
+          ! check thiat molecule i has neighbors
+          if (jbeg .le. jend) then
+             
+             do jnab = jbeg, jend
+                j = list(jnab)
+
+                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
+                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
+                     u_l, A, f, t, pot_local)
+
+             enddo
+          endif
+       enddo
+    endif
+    
+#else
+    
+    if (update_nlist) then
+       
+       ! save current configuration, contruct neighbor list, 
+       ! and calculate forces
+       call saveNeighborList(natoms, q)
+       
+       neighborListSize = size(list)
+   
+       nlist = 0
+       
+       do i = 1, natoms-1
+          point(i) = nlist + 1
+          
+          prepair_inner: do j = i+1, natoms
+             
+             if (skipThisPair(i,j))  cycle prepair_inner
+                          
+             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
+           
+
+             if (rijsq < rlistsq) then
+                
+                nlist = nlist + 1
+              
+                if (nlist > neighborListSize) then
+                   call expandNeighborList(natoms, listerror)
+                   if (listerror /= 0) then
+                      error = -1
+                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
+                      return
+                   end if
+                   neighborListSize = size(list)
+                endif
+                
+                list(nlist) = j
+                
+                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
+                        u_l, A, f, t, pot)
+                
+             endif
+          enddo prepair_inner
+       enddo
+       
+       point(natoms) = nlist + 1
+       
+    else !! (update)
+       
+       ! use the list to find the neighbors
+       do i = 1, natoms-1
+          JBEG = POINT(i)
+          JEND = POINT(i+1) - 1
+          ! check thiat molecule i has neighbors
+          if (jbeg .le. jend) then
+             
+             do jnab = jbeg, jend
+                j = list(jnab)
+
+                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
+                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
+                     u_l, A, f, t, pot)
+
+             enddo
+          endif
+       enddo
+    endif    
+#endif 
+    !! Do rest of preforce calculations
+   call do_preforce(nlocal,pot)
     else
-       !! See if we need to update neighbor lists
+       !! See if we need to update neighbor lists for non pre-pair
        call checkNeighborList(nlocal, q, listSkin, update_nlist)   
     endif
-    
+ 
+
+
+
+
+!---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>
+
+
+
+
+   
 #ifdef IS_MPI
     
     if (update_nlist) then
@@ -531,6 +676,7 @@ contains
     logical :: is_LJ_i, is_LJ_j
     logical :: is_DP_i, is_DP_j
     logical :: is_GB_i, is_GB_j
+    logical :: is_EAM_i,is_EAM_j
     logical :: is_Sticky_i, is_Sticky_j
     integer :: me_i, me_j
 
@@ -598,13 +744,24 @@ contains
        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
+ 
+
 
 
   end subroutine do_pair
 
 
 
-  subroutine do_preforce(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
+  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
@@ -642,9 +799,23 @@ contains
       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 )
+      if ( is_EAM_i .and. is_EAM_j ) &
+           call calc_EAM_prepair_rho(i, j, d, r, rijsq )
+   endif
+  end subroutine do_prepair
+
+
+
+
+  subroutine do_preforce(nlocal,pot)
+    integer :: nlocal
+    real( kind = dp ) :: pot
+
+   if (FF_uses_EAM .and. SimUsesEAM()) then
+      call calc_EAM_preforce_Frho(nlocal,pot)
    endif
+
+
   end subroutine do_preforce