--- trunk/OOPSE/libmdtools/calc_reaction_field.F90	2004/01/30 15:01:09	999
+++ trunk/OOPSE/libmdtools/calc_reaction_field.F90	2004/05/11 21:31:15	1160
@@ -97,14 +97,15 @@ contains
     
   end subroutine createMomentMap  
 
-  subroutine accumulate_rf(atom1, atom2, rij, u_l)
+  subroutine accumulate_rf(atom1, atom2, rij, u_l, taper)
 
     integer, intent(in) :: atom1, atom2
     real (kind = dp), intent(in) :: rij
     real (kind = dp), dimension(3,nLocal) :: u_l    
 
     integer :: me1, me2
-    real (kind = dp) :: taper, mu1, mu2
+    real (kind = dp), intent(in) :: taper
+    real (kind = dp):: mu1, mu2
     real (kind = dp), dimension(3) :: ul1
     real (kind = dp), dimension(3) :: ul2   
 
@@ -124,59 +125,51 @@ contains
        end if
     endif
     
-    if (rij.le.rrf) then
-          
-       if (rij.lt.rt) then
-          taper = 1.0d0
-       else
-  !        write(*,*) 'rf in taper region'
-          taper = (rrf + 2.0d0*rij - 3.0d0*rt)*(rrf-rij)**2/ ((rrf-rt)**3)
-       endif
        
 #ifdef IS_MPI
-       me1 = atid_Row(atom1)
-       ul1(1) = u_l_Row(1,atom1)
-       ul1(2) = u_l_Row(2,atom1)
-       ul1(3) = u_l_Row(3,atom1)
-       
-       me2 = atid_Col(atom2)
-       ul2(1) = u_l_Col(1,atom2)
-       ul2(2) = u_l_Col(2,atom2)
-       ul2(3) = u_l_Col(3,atom2)
+    me1 = atid_Row(atom1)
+    ul1(1) = u_l_Row(1,atom1)
+    ul1(2) = u_l_Row(2,atom1)
+    ul1(3) = u_l_Row(3,atom1)
+    
+    me2 = atid_Col(atom2)
+    ul2(1) = u_l_Col(1,atom2)
+    ul2(2) = u_l_Col(2,atom2)
+    ul2(3) = u_l_Col(3,atom2)
 #else
-       me1 = atid(atom1)
-       ul1(1) = u_l(1,atom1)
-       ul1(2) = u_l(2,atom1)
-       ul1(3) = u_l(3,atom1)
-       
-       me2 = atid(atom2)
-       ul2(1) = u_l(1,atom2)
-       ul2(2) = u_l(2,atom2)
-       ul2(3) = u_l(3,atom2)
-#endif
-       
-       mu1 = MomentMap(me1)%dipole_moment
-       mu2 = MomentMap(me2)%dipole_moment
-       
+    me1 = atid(atom1)
+    ul1(1) = u_l(1,atom1)
+    ul1(2) = u_l(2,atom1)
+    ul1(3) = u_l(3,atom1)
+    
+    me2 = atid(atom2)
+    ul2(1) = u_l(1,atom2)
+    ul2(2) = u_l(2,atom2)
+    ul2(3) = u_l(3,atom2)
+#endif
+    
+    mu1 = MomentMap(me1)%dipole_moment
+    mu2 = MomentMap(me2)%dipole_moment
+    
 #ifdef IS_MPI
-       rf_Row(1,atom1) = rf_Row(1,atom1) + ul2(1)*mu2*taper
-       rf_Row(2,atom1) = rf_Row(2,atom1) + ul2(2)*mu2*taper
-       rf_Row(3,atom1) = rf_Row(3,atom1) + ul2(3)*mu2*taper
-       
-       rf_Col(1,atom2) = rf_Col(1,atom2) + ul1(1)*mu1*taper
-       rf_Col(2,atom2) = rf_Col(2,atom2) + ul1(2)*mu1*taper
-       rf_Col(3,atom2) = rf_Col(3,atom2) + ul1(3)*mu1*taper
+    rf_Row(1,atom1) = rf_Row(1,atom1) + ul2(1)*mu2*taper
+    rf_Row(2,atom1) = rf_Row(2,atom1) + ul2(2)*mu2*taper
+    rf_Row(3,atom1) = rf_Row(3,atom1) + ul2(3)*mu2*taper
+    
+    rf_Col(1,atom2) = rf_Col(1,atom2) + ul1(1)*mu1*taper
+    rf_Col(2,atom2) = rf_Col(2,atom2) + ul1(2)*mu1*taper
+    rf_Col(3,atom2) = rf_Col(3,atom2) + ul1(3)*mu1*taper
 #else
-       rf(1,atom1) = rf(1,atom1) + ul2(1)*mu2*taper
-       rf(2,atom1) = rf(2,atom1) + ul2(2)*mu2*taper
-       rf(3,atom1) = rf(3,atom1) + ul2(3)*mu2*taper
-       
-       rf(1,atom2) = rf(1,atom2) + ul1(1)*mu1*taper
-       rf(2,atom2) = rf(2,atom2) + ul1(2)*mu1*taper
-       rf(3,atom2) = rf(3,atom2) + ul1(3)*mu1*taper     
+    rf(1,atom1) = rf(1,atom1) + ul2(1)*mu2*taper
+    rf(2,atom1) = rf(2,atom1) + ul2(2)*mu2*taper
+    rf(3,atom1) = rf(3,atom1) + ul2(3)*mu2*taper
+    
+    rf(1,atom2) = rf(1,atom2) + ul1(1)*mu1*taper
+    rf(2,atom2) = rf(2,atom2) + ul1(2)*mu1*taper
+    rf(3,atom2) = rf(3,atom2) + ul1(3)*mu1*taper     
 #endif
-       
-    endif
+    
+    
     return  
   end subroutine accumulate_rf
 
@@ -238,11 +231,11 @@ contains
     return
   end subroutine reaction_field_final
   
-  subroutine rf_correct_forces(atom1, atom2, d, rij, u_l, f, do_stress)
+  subroutine rf_correct_forces(atom1, atom2, d, rij, u_l, taper, f, do_stress)
     
     integer, intent(in) :: atom1, atom2
     real(kind=dp), dimension(3), intent(in) :: d
-    real(kind=dp), intent(in) :: rij
+    real(kind=dp), intent(in) :: rij, taper
     real( kind = dp ), dimension(3,nLocal) :: u_l
     real( kind = dp ), dimension(3,nLocal) :: f
     logical, intent(in) :: do_stress