--- trunk/OOPSE/libmdtools/calc_eam.F90	2003/07/30 21:17:01	657
+++ trunk/OOPSE/libmdtools/calc_eam.F90	2003/08/07 00:47:33	669
@@ -15,7 +15,7 @@ module eam
   logical, save :: EAM_FF_initialized = .false.
   integer, save :: EAM_Mixing_Policy
   real(kind = dp), save :: EAM_rcut
-  real(kind = dp), save :: EAM_rcut_orig
+  logical, save :: haveRcut = .false.
 
   character(len = statusMsgSize) :: errMesg
   integer :: eam_err
@@ -51,23 +51,23 @@ module eam
 
 
   !! Arrays for derivatives used in force calculation
-  real( kind = dp), dimension(:), allocatable :: frho
-  real( kind = dp), dimension(:), allocatable :: rho
+  real( kind = dp),save, dimension(:), allocatable :: frho
+  real( kind = dp),save, dimension(:), allocatable :: rho
 
-  real( kind = dp), dimension(:), allocatable :: dfrhodrho
-  real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho
+  real( kind = dp),save, dimension(:), allocatable :: dfrhodrho
+  real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho
 
 
 !! Arrays for MPI storage
 #ifdef IS_MPI
-  real( kind = dp), dimension(:), allocatable :: dfrhodrho_col
-  real( kind = dp), dimension(:), allocatable :: dfrhodrho_row
-  real( kind = dp), dimension(:), allocatable :: frho_row
-  real( kind = dp), dimension(:), allocatable :: frho_col
-  real( kind = dp), dimension(:), allocatable :: rho_row
-  real( kind = dp), dimension(:), allocatable :: rho_col
-  real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho_col
-  real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho_row
+  real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_col
+  real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_row
+  real( kind = dp),save, dimension(:), allocatable :: frho_row
+  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 :: d2frhodrhodrho_col
+  real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_row
 #endif
 
   type, private :: EAMTypeList
@@ -117,7 +117,7 @@ contains
 
     status = 0
 
-    write(*,*) "Adding new eamtype: ",eam_ident
+
     !! Assume that atypes has already been set and get the total number of types in atypes
     !! Also assume that every member of atypes is a EAM model.
    
@@ -222,14 +222,14 @@ contains
        enddo
 
 
-       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
-       !! find the smallest rcut for any eam atype
-       do i = 2, EAMList%currentAddition 
-          current_rcut_max = min(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
-       end do
+!       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
+       !! find the smallest rcut for any eam atype
+!       do i = 2, EAMList%currentAddition 
+!          current_rcut_max =max(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
+!       end do
 
-       EAM_rcut = current_rcut_max
-       EAM_rcut_orig = current_rcut_max
+!       EAM_rcut = current_rcut_max
+!       EAM_rcut_orig = current_rcut_max
 !       do i = 1, EAMList%currentAddition
 !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
 !       end do
@@ -350,22 +350,22 @@ contains
 
   end subroutine allocateEAM
 
+!! C sets rcut to be the largest cutoff of any atype 
+!! present in this simulation. Doesn't include all atypes
+!! sim knows about, just those in the simulation.
   subroutine setCutoffEAM(rcut, status)
     real(kind=dp) :: rcut
     integer :: status
+    status = 0
 
-    if (rcut < EAM_rcut) then
-       EAM_rcut = rcut
-    endif
+    EAM_rcut = rcut
 
-
   end subroutine setCutoffEAM
 
 
 
   subroutine clean_EAM()
-
-! clean non-IS_MPI first
+   ! clean non-IS_MPI first
     frho = 0.0_dp
     rho  = 0.0_dp
     dfrhodrho = 0.0_dp
@@ -484,10 +484,15 @@ contains
     integer :: myid_atom2
 
 ! check to see if we need to be cleaned at the start of a force loop
-    if (cleanme) call clean_EAM
-    cleanme = .false.
     
+    if (cleanme) then
+       call clean_EAM
+       cleanme = .false.
+    end if
+      
 
+    
+
 #ifdef IS_MPI
     myid_atom1 = atid_Row(atom1)
     myid_atom2 = atid_Col(atom2)
@@ -532,6 +537,8 @@ contains
 #endif
        endif
 
+     
+
   end subroutine calc_eam_prepair_rho
 
 
@@ -549,7 +556,7 @@ contains
     integer :: n_rho_points
     ! reset clean forces to be true at top of calc rho.
     cleanme = .true.
-
+  
 !! Scatter the electron density from  pre-pair calculation back to local atoms
 #ifdef IS_MPI
     call scatter(rho_row,rho,plan_row,eam_err)
@@ -594,8 +601,8 @@ contains
        pot = pot + u
     enddo
 
+  
 
-
 #ifdef IS_MPI
     !! communicate f(rho) and derivatives back into row and column arrays
     call gather(frho,frho_row,plan_row, eam_err)
@@ -661,12 +668,11 @@ contains
 !Local Variables
     
 
-    
+
     phab = 0.0E0_DP
     dvpdr = 0.0E0_DP
     d2vpdrdr = 0.0E0_DP
-    
-    
+
     if (rij .lt. EAM_rcut) then
 #ifdef IS_MPI
 !!!!! FIX ME
@@ -753,11 +759,11 @@ contains
 #ifdef IS_MPI
        dudr = drhojdr*dfrhodrho_row(atom1)+drhoidr*dfrhodrho_col(atom2) &
             + dvpdr
-    
+
 #else
        dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) &
             + dvpdr
-
+      ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
 #endif
 
        fx = dudr * drdx
@@ -779,8 +785,11 @@ contains
        f_Col(2,atom2) = f_Col(2,atom2) - fy
        f_Col(3,atom2) = f_Col(3,atom2) - fz
 #else
-       if(do_pot) pot = pot + phab
-       
+
+       if(do_pot) then
+          pot = pot + phab
+       end if
+
        f(1,atom1) = f(1,atom1) + fx
        f(2,atom1) = f(2,atom1) + fy
        f(3,atom1) = f(3,atom1) + fz