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

Comparing trunk/OOPSE/libmdtools/calc_dipole_dipole.F90 (file contents):
Revision 486 by mmeineke, Thu Apr 10 16:22:00 2003 UTC vs.
Revision 843 by mmeineke, Wed Oct 29 20:41:39 2003 UTC

# Line 10 | Line 10 | module dipole_dipole
10   #endif
11    implicit none
12  
13 +  PRIVATE
14    real(kind=dp), save :: rrf = 0.0
15    real(kind=dp), save :: rt  = 0.0
16 <  real(kind=dp), save :: rrfsq = 0.0
16 <  real(kind=dp), save :: pre = 0.0
16 >   real(kind=dp), save :: pre = 0.0
17    logical, save :: dipole_initialized = .false.
18  
19 +
20 +  public::setCutoffsDipole
21 +  public::do_dipole_pair
22 +
23   contains
24      
25 <  subroutine initialize_dipole(this_rrf, this_rt)
25 >  subroutine setCutoffsDipole(this_rrf, this_rt)
26      real(kind=dp), intent(in) :: this_rrf, this_rt
27      rrf = this_rrf
28      rt = this_rt    
29 <    rrfsq = rrf * rrf
30 <    ! pre converts from mu in units of debye to kcal/mol
29 >
30 >      ! pre converts from mu in units of debye to kcal/mol
31      pre = 14.38362_dp
32  
33      dipole_initialized = .true.
34      
35      return
36 <  end subroutine initialize_dipole
36 >  end subroutine setCutoffsDipole
37  
34
38    subroutine do_dipole_pair(atom1, atom2, d, rij, r2, pot, u_l, f, t, &
39         do_pot, do_stress)
40      
41      logical :: do_pot, do_stress
42  
43 <    integer atom1, atom2, me1, me2
43 >    integer atom1, atom2, me1, me2, id1, id2
44      real(kind=dp) :: rij, mu1, mu2
45      real(kind=dp) :: dfact1, dfact2, dip2, r2, r3, r5
46      real(kind=dp) :: dudx, dudy, dudz, dudu1x, dudu1y, dudu1z
# Line 173 | Line 176 | contains
176         t(2,atom2) = t(2,atom2) - ul2(3)*dudu2x + ul2(1)*dudu2z
177         t(3,atom2) = t(3,atom2) - ul2(1)*dudu2y + ul2(2)*dudu2x
178   #endif
179 <      
179 >
180         if (do_stress) then          
181  
182 <          if (molMembershipList(atom1) .ne. molMembershipList(atom2)) then
182 > #ifdef IS_MPI
183 >          id1 = tagRow(atom1)
184 >          id2 = tagColumn(atom2)
185 > #else
186 >          id1 = atom1
187 >          id2 = atom2
188 > #endif                
189 >          if (molMembershipList(id1) .ne. molMembershipList(id2)) then
190  
191 <             tau_Temp(1) = tau_Temp(1) + dudx * d(1)
192 <             tau_Temp(2) = tau_Temp(2) + dudx * d(2)
193 <             tau_Temp(3) = tau_Temp(3) + dudx * d(3)
194 <             tau_Temp(4) = tau_Temp(4) + dudy * d(1)
195 <             tau_Temp(5) = tau_Temp(5) + dudy * d(2)
196 <             tau_Temp(6) = tau_Temp(6) + dudy * d(3)
197 <             tau_Temp(7) = tau_Temp(7) + dudz * d(1)
198 <             tau_Temp(8) = tau_Temp(8) + dudz * d(2)
199 <             tau_Temp(9) = tau_Temp(9) + dudz * d(3)
191 >             ! because the d vector is the rj - ri vector, and
192 >             ! because dudx, dudy, dudz are the (positive) force on
193 >             ! atom i (negative on atom j) we need a negative sign here:
194 >
195 >             tau_Temp(1) = tau_Temp(1) - d(1) * dudx
196 >             tau_Temp(2) = tau_Temp(2) - d(1) * dudy
197 >             tau_Temp(3) = tau_Temp(3) - d(1) * dudz
198 >             tau_Temp(4) = tau_Temp(4) - d(2) * dudx
199 >             tau_Temp(5) = tau_Temp(5) - d(2) * dudy
200 >             tau_Temp(6) = tau_Temp(6) - d(2) * dudz
201 >             tau_Temp(7) = tau_Temp(7) - d(3) * dudx
202 >             tau_Temp(8) = tau_Temp(8) - d(3) * dudy
203 >             tau_Temp(9) = tau_Temp(9) - d(3) * dudz
204 >
205               virial_Temp = virial_Temp + &
206                    (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
207  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines