--- trunk/OOPSE/libmdtools/calc_eam.F90	2003/08/27 16:25:11	730
+++ trunk/OOPSE/libmdtools/calc_eam.F90	2004/05/07 21:35:05	1150
@@ -66,6 +66,7 @@ module eam
   real( kind = dp),save, dimension(:), allocatable :: frho_col
   real( kind = dp),save, dimension(:), allocatable :: rho_row
   real( kind = dp),save, dimension(:), allocatable :: rho_col
+  real( kind = dp),save, dimension(:), allocatable :: rho_tmp
   real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_col
   real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_row
 #endif
@@ -78,7 +79,7 @@ module eam
   end type EAMTypeList
 
 
-  type (eamTypeList) :: EAMList
+  type (eamTypeList), save :: EAMList
 
   !! standard eam stuff  
 
@@ -164,14 +165,15 @@ contains
     integer :: alloc_stat
     integer :: number_r, number_rho
 
+
+    status = 0
     if (EAMList%currentAddition == 0) then
        call handleError("init_EAM_FF","No members in EAMList")
        status = -1
        return
     end if
 
-
-
+ 
        do i = 1, EAMList%currentAddition
 
 ! Build array of r values
@@ -221,7 +223,6 @@ contains
                0.0E0_DP, 0.0E0_DP, 'N')
        enddo
 
-
 !       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
        !! find the smallest rcut for any eam atype
 !       do i = 2, EAMList%currentAddition 
@@ -234,12 +235,14 @@ contains
 !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
 !       end do
        !! Allocate arrays for force calculation
-          call allocateEAM(alloc_stat)
-          if (alloc_stat /= 0 ) then
-             status = -1
-             return
-          endif
-
+      
+       call allocateEAM(alloc_stat)
+       if (alloc_stat /= 0 ) then
+          write(*,*) "allocateEAM failed"
+          status = -1
+          return
+       endif
+    
   end subroutine init_EAM_FF
 
 !! routine checks to see if array is allocated, deallocates array if allocated
@@ -247,7 +250,6 @@ contains
   subroutine allocateEAM(status)
     integer, intent(out) :: status
 
-    integer :: nlocal
 #ifdef IS_MPI
     integer :: nrow
     integer :: ncol
@@ -255,8 +257,7 @@ contains
     integer :: alloc_stat
 
 
-    nlocal = getNlocal()
-
+    status = 0
 #ifdef IS_MPI
     nrow = getNrow(plan_row)
     ncol = getNcol(plan_col)
@@ -290,7 +291,15 @@ contains
     end if
     
 #ifdef IS_MPI
+
+    if (allocated(rho_tmp)) deallocate(rho_tmp)
+    allocate(rho_tmp(nlocal),stat=alloc_stat)
+    if (alloc_stat /= 0) then 
+       status = -1
+       return
+    end if
 
+
     if (allocated(frho_row)) deallocate(frho_row)
     allocate(frho_row(nrow),stat=alloc_stat)
     if (alloc_stat /= 0) then 
@@ -374,6 +383,7 @@ contains
     frho_col = 0.0_dp
     rho_row  = 0.0_dp
     rho_col  = 0.0_dp
+    rho_tmp  = 0.0_dp
     dfrhodrho_row = 0.0_dp
     dfrhodrho_col = 0.0_dp
 #endif
@@ -522,7 +532,7 @@ contains
                EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, &
                r, rho_j_at_i,drho,d2rho)
 
-
+    
        
        
 #ifdef  IS_MPI
@@ -533,6 +543,10 @@ contains
        endif
 
 
+
+
+
+
   end subroutine calc_eam_prepair_rho
 
 
@@ -558,14 +572,17 @@ contains
       write(errMsg,*) " Error scattering rho_row into rho"
       call handleError(RoutineName,errMesg)
    endif       
-    call scatter(rho_col,rho,plan_col,eam_err)
+    call scatter(rho_col,rho_tmp,plan_col,eam_err)
     if (eam_err /= 0 ) then
       write(errMsg,*) " Error scattering rho_col into rho"
       call handleError(RoutineName,errMesg)
    endif
+
+      rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
 #endif
 
 
+
 !! Calculate F(rho) and derivative 
     do atom = 1, nlocal
        me = atid(atom)
@@ -634,12 +651,13 @@ contains
 
 
   !! Does EAM pairwise Force calculation.  
-  subroutine do_eam_pair(atom1,atom2,d,rij,r2,pot,f,do_pot,do_stress)
+  subroutine do_eam_pair(atom1, atom2, d, rij, r2, sw, vpair, pot, f, &
+       do_pot, do_stress)
     !Arguments    
     integer, intent(in) ::  atom1, atom2
     real( kind = dp ), intent(in) :: rij, r2
-    real( kind = dp ) :: pot
-    real( kind = dp ), dimension(3,getNlocal()) :: f
+    real( kind = dp ) :: pot, sw, vpair
+    real( kind = dp ), dimension(3,nLocal) :: f
     real( kind = dp ), intent(in), dimension(3) :: d
     logical, intent(in) :: do_pot, do_stress
     
@@ -771,6 +789,7 @@ contains
           pot_Row(atom1) = pot_Row(atom1) + phab*0.5
           pot_Col(atom2) = pot_Col(atom2) + phab*0.5
        end if
+       vpair = vpair + phab
 
        f_Row(1,atom1) = f_Row(1,atom1) + fx
        f_Row(2,atom1) = f_Row(2,atom1) + fy
@@ -784,6 +803,7 @@ contains
        if(do_pot) then
           pot = pot + phab
        end if
+       vpair = vpair + phab
 
        f(1,atom1) = f(1,atom1) + fx
        f(2,atom1) = f(2,atom1) + fy