--- trunk/OOPSE/libmdtools/calc_reaction_field.F90	2003/03/24 21:55:34	394
+++ trunk/OOPSE/libmdtools/calc_reaction_field.F90	2003/07/15 17:10:50	611
@@ -3,6 +3,7 @@ module reaction_field
   use definitions
   use atype_module
   use vector_class
+  use simulation
 #ifdef IS_MPI
   use mpiSimulation
 #endif
@@ -41,7 +42,7 @@ contains
 
     integer, intent(in) :: atom1, atom2
     real (kind = dp), intent(in) :: rij
-    real (kind = dp), dimension(:,:) :: u_l    
+    real (kind = dp), dimension(3,getNlocal()) :: u_l    
 
     integer :: me1, me2
     real (kind = dp) :: taper, mu1, mu2
@@ -113,7 +114,7 @@ contains
     
     integer, intent(in) :: atom1
     real(kind=dp), intent(in) :: mu1
-    real(kind=dp), dimension(:,:) :: u_l
+    real(kind=dp), dimension(3,getNlocal()) :: u_l
     
     !! should work for both MPI and non-MPI version since this is not pairwise.
     rf(1,atom1) = rf(1,atom1) + u_l(1,atom1)*mu1
@@ -129,8 +130,8 @@ contains
     real (kind=dp), intent(in) :: mu1
     real (kind=dp), intent(inout) :: rfpot
     logical, intent(in) :: do_pot
-    real (kind = dp), dimension(:,:) :: u_l    
-    real (kind = dp), dimension(:,:) :: t
+    real (kind = dp), dimension(3,getNlocal()) :: u_l    
+    real (kind = dp), dimension(3,getNlocal()) :: t
 
     if (.not.rf_initialized) then
        write(default_error,*) 'Reaction field not initialized!'
@@ -161,8 +162,8 @@ contains
     integer, intent(in) :: atom1, atom2
     real(kind=dp), dimension(3), intent(in) :: d
     real(kind=dp), intent(in) :: rij
-    real( kind = dp ), dimension(:,:) :: u_l
-    real( kind = dp ), dimension(:,:) :: f
+    real( kind = dp ), dimension(3,getNlocal()) :: u_l
+    real( kind = dp ), dimension(3,getNlocal()) :: f
     logical, intent(in) :: do_stress
     
     real (kind = dp), dimension(3) :: ul1
@@ -235,16 +236,25 @@ contains
 #endif
        
        if (do_stress) then          
-          tau_Temp(1) = tau_Temp(1) + dudx * d(1)
-          tau_Temp(2) = tau_Temp(2) + dudx * d(2)
-          tau_Temp(3) = tau_Temp(3) + dudx * d(3)
-          tau_Temp(4) = tau_Temp(4) + dudy * d(1)
-          tau_Temp(5) = tau_Temp(5) + dudy * d(2)
-          tau_Temp(6) = tau_Temp(6) + dudy * d(3)
-          tau_Temp(7) = tau_Temp(7) + dudz * d(1)
-          tau_Temp(8) = tau_Temp(8) + dudz * d(2)
-          tau_Temp(9) = tau_Temp(9) + dudz * d(3)
-          virial_Temp = virial_Temp + (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
+          if (molMembershipList(atom1) .ne. molMembershipList(atom2)) then
+
+             ! because the d vector is the rj - ri vector, and 
+             ! because dudx, dudy, and dudz are the
+             ! (positive) force on atom i (negative on atom j) we need
+             ! a negative sign here:
+
+             tau_Temp(1) = tau_Temp(1) - d(1) * dudx
+             tau_Temp(2) = tau_Temp(2) - d(1) * dudy
+             tau_Temp(3) = tau_Temp(3) - d(1) * dudz
+             tau_Temp(4) = tau_Temp(4) - d(2) * dudx
+             tau_Temp(5) = tau_Temp(5) - d(2) * dudy
+             tau_Temp(6) = tau_Temp(6) - d(2) * dudz
+             tau_Temp(7) = tau_Temp(7) - d(3) * dudx
+             tau_Temp(8) = tau_Temp(8) - d(3) * dudy
+             tau_Temp(9) = tau_Temp(9) - d(3) * dudz
+             virial_Temp = virial_Temp + &
+                  (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
+          endif
        endif       
     endif