ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/calc_reaction_field.F90
(Generate patch)

Comparing trunk/OOPSE/libmdtools/calc_reaction_field.F90 (file contents):
Revision 394 by gezelter, Mon Mar 24 21:55:34 2003 UTC vs.
Revision 611 by gezelter, Tue Jul 15 17:10:50 2003 UTC

# Line 3 | Line 3 | module reaction_field
3    use definitions
4    use atype_module
5    use vector_class
6 +  use simulation
7   #ifdef IS_MPI
8    use mpiSimulation
9   #endif
# Line 41 | Line 42 | contains
42  
43      integer, intent(in) :: atom1, atom2
44      real (kind = dp), intent(in) :: rij
45 <    real (kind = dp), dimension(:,:) :: u_l    
45 >    real (kind = dp), dimension(3,getNlocal()) :: u_l    
46  
47      integer :: me1, me2
48      real (kind = dp) :: taper, mu1, mu2
# Line 113 | Line 114 | contains
114      
115      integer, intent(in) :: atom1
116      real(kind=dp), intent(in) :: mu1
117 <    real(kind=dp), dimension(:,:) :: u_l
117 >    real(kind=dp), dimension(3,getNlocal()) :: u_l
118      
119      !! should work for both MPI and non-MPI version since this is not pairwise.
120      rf(1,atom1) = rf(1,atom1) + u_l(1,atom1)*mu1
# Line 129 | Line 130 | contains
130      real (kind=dp), intent(in) :: mu1
131      real (kind=dp), intent(inout) :: rfpot
132      logical, intent(in) :: do_pot
133 <    real (kind = dp), dimension(:,:) :: u_l    
134 <    real (kind = dp), dimension(:,:) :: t
133 >    real (kind = dp), dimension(3,getNlocal()) :: u_l    
134 >    real (kind = dp), dimension(3,getNlocal()) :: t
135  
136      if (.not.rf_initialized) then
137         write(default_error,*) 'Reaction field not initialized!'
# Line 161 | Line 162 | contains
162      integer, intent(in) :: atom1, atom2
163      real(kind=dp), dimension(3), intent(in) :: d
164      real(kind=dp), intent(in) :: rij
165 <    real( kind = dp ), dimension(:,:) :: u_l
166 <    real( kind = dp ), dimension(:,:) :: f
165 >    real( kind = dp ), dimension(3,getNlocal()) :: u_l
166 >    real( kind = dp ), dimension(3,getNlocal()) :: f
167      logical, intent(in) :: do_stress
168      
169      real (kind = dp), dimension(3) :: ul1
# Line 235 | Line 236 | contains
236   #endif
237        
238         if (do_stress) then          
239 <          tau_Temp(1) = tau_Temp(1) + dudx * d(1)
240 <          tau_Temp(2) = tau_Temp(2) + dudx * d(2)
241 <          tau_Temp(3) = tau_Temp(3) + dudx * d(3)
242 <          tau_Temp(4) = tau_Temp(4) + dudy * d(1)
243 <          tau_Temp(5) = tau_Temp(5) + dudy * d(2)
244 <          tau_Temp(6) = tau_Temp(6) + dudy * d(3)
245 <          tau_Temp(7) = tau_Temp(7) + dudz * d(1)
246 <          tau_Temp(8) = tau_Temp(8) + dudz * d(2)
247 <          tau_Temp(9) = tau_Temp(9) + dudz * d(3)
248 <          virial_Temp = virial_Temp + (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
239 >          if (molMembershipList(atom1) .ne. molMembershipList(atom2)) then
240 >
241 >             ! because the d vector is the rj - ri vector, and
242 >             ! because dudx, dudy, and dudz are the
243 >             ! (positive) force on atom i (negative on atom j) we need
244 >             ! a negative sign here:
245 >
246 >             tau_Temp(1) = tau_Temp(1) - d(1) * dudx
247 >             tau_Temp(2) = tau_Temp(2) - d(1) * dudy
248 >             tau_Temp(3) = tau_Temp(3) - d(1) * dudz
249 >             tau_Temp(4) = tau_Temp(4) - d(2) * dudx
250 >             tau_Temp(5) = tau_Temp(5) - d(2) * dudy
251 >             tau_Temp(6) = tau_Temp(6) - d(2) * dudz
252 >             tau_Temp(7) = tau_Temp(7) - d(3) * dudx
253 >             tau_Temp(8) = tau_Temp(8) - d(3) * dudy
254 >             tau_Temp(9) = tau_Temp(9) - d(3) * dudz
255 >             virial_Temp = virial_Temp + &
256 >                  (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
257 >          endif
258         endif      
259      endif
260      

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines