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 460 by chuckv, Fri Apr 4 22:22:30 2003 UTC vs.
Revision 621 by gezelter, Wed Jul 16 02:11:02 2003 UTC

# Line 10 | Line 10 | module dipole_dipole
10   #endif
11    implicit none
12  
13 <  real(kind=dp), save :: rrf
14 <  real(kind=dp), save :: rt  
15 <  real(kind=dp), save :: rrfsq
16 <  real(kind=dp), save :: pre
13 >  real(kind=dp), save :: rrf = 0.0
14 >  real(kind=dp), save :: rt  = 0.0
15 >  real(kind=dp), save :: rrfsq = 0.0
16 >  real(kind=dp), save :: pre = 0.0
17    logical, save :: dipole_initialized = .false.
18  
19   contains
# Line 31 | Line 31 | contains
31      return
32    end subroutine initialize_dipole
33  
34 –
34    subroutine do_dipole_pair(atom1, atom2, d, rij, r2, pot, u_l, f, t, &
35         do_pot, do_stress)
36      
# Line 39 | Line 38 | contains
38  
39      integer atom1, atom2, me1, me2
40      real(kind=dp) :: rij, mu1, mu2
41 <    real(kind=dp) :: dfact1, dfact2, dip2, r2, r3, r5, pre
41 >    real(kind=dp) :: dfact1, dfact2, dip2, r2, r3, r5
42      real(kind=dp) :: dudx, dudy, dudz, dudu1x, dudu1y, dudu1z
43      real(kind=dp) :: dudu2x, dudu2y, dudu2z, rdotu1, rdotu2, u1dotu2
44      real(kind=dp) :: taper, dtdr, vterm
# Line 80 | Line 79 | contains
79      ul2(3) = u_l(3,atom2)
80   #endif
81  
83 –    if( atom1 .eq. 2 )then
84 –       write (*,*) 'ul =', ul1(1), ul1(2), ul1(3)
85 –    endif
82  
87 –    if( atom2 .eq. 2 )then
88 –       write (*,*) 'ul =', ul2(1), ul2(2), ul2(3)
89 –    endif
90 –
91 –
83      call getElementProperty(atypes, me1, "dipole_moment", mu1)
84      call getElementProperty(atypes, me2, "dipole_moment", mu2)
85  
# Line 101 | Line 92 | contains
92            taper = (rrf + 2.0d0*rij - 3.0d0*rt)*(rrf-rij)**2/ ((rrf-rt)**3)
93            dtdr = 6.0d0*(rij*rij - rij*rt - rij*rrf +rrf*rt)/((rrf-rt)**3)
94         endif
95 <      
95 >
96         r3 = r2*rij
97         r5 = r3*r2
98        
99         rdotu1 = d(1)*ul1(1) + d(2)*ul1(2) + d(3)*ul1(3)
100         rdotu2 = d(1)*ul2(1) + d(2)*ul2(2) + d(3)*ul2(3)
101         u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3)
102 <      
102 >
103         dip2 = pre * mu1 * mu2
113 –      
104         dfact1 = 3.0d0*dip2 / r2
105         dfact2 = 3.0d0*dip2 / r5
106 <      
106 >
107         vterm = dip2*((u1dotu2/r3) - 3.0d0*(rdotu1*rdotu2/r5))
108 <      
108 >
109         if (do_pot) then
110   #ifdef IS_MPI
111            pot_row(atom1) = pot_row(atom1) + 0.5d0*vterm*taper
# Line 139 | Line 129 | contains
129              (5.0d0*(rdotu1*rdotu2)/r5)) -  &
130              dfact2*(ul1(3)*rdotu2 + ul2(3)*rdotu1))*taper +  &
131              vterm*dtdr*d(3)/rij
132 <      
132 >
133         dudu1x = (dip2*((ul2(1)/r3) - (3.0d0*d(1)*rdotu2/r5)))*taper
134         dudu1y = (dip2*((ul2(2)/r3) - (3.0d0*d(2)*rdotu2/r5)))*taper
135         dudu1z = (dip2*((ul2(3)/r3) - (3.0d0*d(3)*rdotu2/r5)))*taper
136 <      
136 >
137         dudu2x = (dip2*((ul1(1)/r3) - (3.0d0*d(1)*rdotu1/r5)))*taper
138         dudu2y = (dip2*((ul1(2)/r3) - (3.0d0*d(2)*rdotu1/r5)))*taper
139         dudu2z = (dip2*((ul1(3)/r3) - (3.0d0*d(3)*rdotu1/r5)))*taper
140        
141 <      
141 >
142   #ifdef IS_MPI
143 <       f_Row(1,atom1) = f_Row(1,atom1) - dudx
144 <       f_Row(2,atom1) = f_Row(2,atom1) - dudy
145 <       f_Row(3,atom1) = f_Row(3,atom1) - dudz
143 >       f_Row(1,atom1) = f_Row(1,atom1) + dudx
144 >       f_Row(2,atom1) = f_Row(2,atom1) + dudy
145 >       f_Row(3,atom1) = f_Row(3,atom1) + dudz
146  
147 <       f_Col(1,atom2) = f_Col(1,atom2) + dudx
148 <       f_Col(2,atom2) = f_Col(2,atom2) + dudy
149 <       f_Col(3,atom2) = f_Col(3,atom2) + dudz
147 >       f_Col(1,atom2) = f_Col(1,atom2) - dudx
148 >       f_Col(2,atom2) = f_Col(2,atom2) - dudy
149 >       f_Col(3,atom2) = f_Col(3,atom2) - dudz
150        
151         t_Row(1,atom1) = t_Row(1,atom1) - ul1(2)*dudu1z + ul1(3)*dudu1y
152         t_Row(2,atom1) = t_Row(2,atom1) - ul1(3)*dudu1x + ul1(1)*dudu1z
# Line 166 | Line 156 | contains
156         t_Col(2,atom2) = t_Col(2,atom2) - ul2(3)*dudu2x + ul2(1)*dudu2z
157         t_Col(3,atom2) = t_Col(3,atom2) - ul2(1)*dudu2y + ul2(2)*dudu2x
158   #else
159 <       f(1,atom1) = f(1,atom1) - dudx
160 <       f(2,atom1) = f(2,atom1) - dudy
161 <       f(3,atom1) = f(3,atom1) - dudz
159 >       f(1,atom1) = f(1,atom1) + dudx
160 >       f(2,atom1) = f(2,atom1) + dudy
161 >       f(3,atom1) = f(3,atom1) + dudz
162        
163 <       f(1,atom2) = f(1,atom2) + dudx
164 <       f(2,atom2) = f(2,atom2) + dudy
165 <       f(3,atom2) = f(3,atom2) + dudz
163 >       f(1,atom2) = f(1,atom2) - dudx
164 >       f(2,atom2) = f(2,atom2) - dudy
165 >       f(3,atom2) = f(3,atom2) - dudz
166        
167         t(1,atom1) = t(1,atom1) - ul1(2)*dudu1z + ul1(3)*dudu1y
168         t(2,atom1) = t(2,atom1) - ul1(3)*dudu1x + ul1(1)*dudu1z
# Line 184 | Line 174 | contains
174   #endif
175        
176         if (do_stress) then          
177 <          tau_Temp(1) = tau_Temp(1) + dudx * d(1)
178 <          tau_Temp(2) = tau_Temp(2) + dudx * d(2)
179 <          tau_Temp(3) = tau_Temp(3) + dudx * d(3)
180 <          tau_Temp(4) = tau_Temp(4) + dudy * d(1)
181 <          tau_Temp(5) = tau_Temp(5) + dudy * d(2)
182 <          tau_Temp(6) = tau_Temp(6) + dudy * d(3)
183 <          tau_Temp(7) = tau_Temp(7) + dudz * d(1)
184 <          tau_Temp(8) = tau_Temp(8) + dudz * d(2)
185 <          tau_Temp(9) = tau_Temp(9) + dudz * d(3)
186 <          virial_Temp = virial_Temp + (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
177 >
178 >          if (molMembershipList(atom1) .ne. molMembershipList(atom2)) then
179 >
180 >             ! because the d vector is the rj - ri vector, and
181 >             ! because dudx, dudy, dudz are the (positive) force on
182 >             ! atom i (negative on atom j) we need a negative sign here:
183 >
184 >             tau_Temp(1) = tau_Temp(1) - d(1) * dudx
185 >             tau_Temp(2) = tau_Temp(2) - d(1) * dudy
186 >             tau_Temp(3) = tau_Temp(3) - d(1) * dudz
187 >             tau_Temp(4) = tau_Temp(4) - d(2) * dudx
188 >             tau_Temp(5) = tau_Temp(5) - d(2) * dudy
189 >             tau_Temp(6) = tau_Temp(6) - d(2) * dudz
190 >             tau_Temp(7) = tau_Temp(7) - d(3) * dudx
191 >             tau_Temp(8) = tau_Temp(8) - d(3) * dudy
192 >             tau_Temp(9) = tau_Temp(9) - d(3) * dudz
193 >
194 >             virial_Temp = virial_Temp + &
195 >                  (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
196 >
197 >          endif          
198         endif
199        
200      endif

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines