ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/reactionField.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/reactionField.F90 (file contents):
Revision 2203 by chrisfen, Mon Mar 21 20:51:10 2005 UTC vs.
Revision 2204 by gezelter, Fri Apr 15 22:04:00 2005 UTC

# Line 53 | Line 53 | module reaction_field
53    implicit none
54  
55    PRIVATE
56 <  
56 >
57    real(kind=dp), save :: rrf = 1.0_dp
58    real(kind=dp), save :: rt
59    real(kind=dp), save :: dielect = 1.0_dp
# Line 71 | Line 71 | contains
71    PUBLIC::rf_correct_forces
72  
73   contains
74 <  
74 >
75    subroutine initialize_rf(this_dielect)
76      real(kind=dp), intent(in) :: this_dielect
77 <    
77 >
78      dielect = this_dielect
79  
80      pre = pre22 * 2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf)
81 <    
81 >
82      haveDielectric = .true.
83  
84      return
# Line 93 | Line 93 | contains
93  
94      rrfsq = rrf * rrf
95      pre = pre22 * 2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf)
96 <    
96 >
97      haveCutoffs = .true.
98  
99    end subroutine setCutoffsRF
# Line 115 | Line 115 | contains
115      if ((.not.haveDielectric).or.(.not.haveCutoffs)) then
116         write(default_error,*) 'Reaction field not initialized!'
117         return
118 <    endif    
119 <      
118 >    endif
119 >
120   #ifdef IS_MPI
121      me1 = atid_Row(atom1)
122      ul1(1) = eFrame_Row(3,atom1)
123      ul1(2) = eFrame_Row(6,atom1)
124      ul1(3) = eFrame_Row(9,atom1)
125 <    
125 >
126      me2 = atid_Col(atom2)
127      ul2(1) = eFrame_Col(3,atom2)
128      ul2(2) = eFrame_Col(6,atom2)
# Line 132 | Line 132 | contains
132      ul1(1) = eFrame(3,atom1)
133      ul1(2) = eFrame(6,atom1)
134      ul1(3) = eFrame(9,atom1)
135 <    
135 >
136      me2 = atid(atom2)
137      ul2(1) = eFrame(3,atom2)
138      ul2(2) = eFrame(6,atom2)
139      ul2(3) = eFrame(9,atom2)
140   #endif
141 <    
141 >
142      mu1 = getDipoleMoment(me1)
143      mu2 = getDipoleMoment(me2)
144 <    
144 >
145   #ifdef IS_MPI
146      rf_Row(1,atom1) = rf_Row(1,atom1) + ul2(1)*mu2*taper
147      rf_Row(2,atom1) = rf_Row(2,atom1) + ul2(2)*mu2*taper
148      rf_Row(3,atom1) = rf_Row(3,atom1) + ul2(3)*mu2*taper
149 <    
149 >
150      rf_Col(1,atom2) = rf_Col(1,atom2) + ul1(1)*mu1*taper
151      rf_Col(2,atom2) = rf_Col(2,atom2) + ul1(2)*mu1*taper
152      rf_Col(3,atom2) = rf_Col(3,atom2) + ul1(3)*mu1*taper
# Line 154 | Line 154 | contains
154      rf(1,atom1) = rf(1,atom1) + ul2(1)*mu2*taper
155      rf(2,atom1) = rf(2,atom1) + ul2(2)*mu2*taper
156      rf(3,atom1) = rf(3,atom1) + ul2(3)*mu2*taper
157 <    
157 >
158      rf(1,atom2) = rf(1,atom2) + ul1(1)*mu1*taper
159      rf(2,atom2) = rf(2,atom2) + ul1(2)*mu1*taper
160      rf(3,atom2) = rf(3,atom2) + ul1(3)*mu1*taper    
161   #endif
162 <    
163 <    
162 >
163 >
164      return  
165    end subroutine accumulate_rf
166  
167    subroutine accumulate_self_rf(atom1, mu1, eFrame)
168 <    
168 >
169      integer, intent(in) :: atom1
170      real(kind=dp), intent(in) :: mu1
171      real(kind=dp), dimension(9,nLocal) :: eFrame
172 <    
172 >
173      !! should work for both MPI and non-MPI version since this is not pairwise.
174      rf(1,atom1) = rf(1,atom1) + eFrame(3,atom1)*mu1
175      rf(2,atom1) = rf(2,atom1) + eFrame(6,atom1)*mu1
176      rf(3,atom1) = rf(3,atom1) + eFrame(9,atom1)*mu1
177 <        
177 >
178      return
179    end subroutine accumulate_self_rf
180 <  
180 >
181    subroutine reaction_field_final(a1, mu1, eFrame, rfpot, t, do_pot)
182 <            
182 >
183      integer, intent(in) :: a1
184      real (kind=dp), intent(in) :: mu1
185      real (kind=dp), intent(inout) :: rfpot
# Line 196 | Line 196 | contains
196  
197      ! compute torques on dipoles:
198      ! pre converts from mu in units of debye to kcal/mol
199 <    
199 >
200      ! The torque contribution is dipole cross reaction_field  
201  
202      t(1,a1) = t(1,a1) + pre*mu1*(eFrame(6,a1)*rf(3,a1) - eFrame(9,a1)*rf(2,a1))
203      t(2,a1) = t(2,a1) + pre*mu1*(eFrame(9,a1)*rf(1,a1) - eFrame(3,a1)*rf(3,a1))
204      t(3,a1) = t(3,a1) + pre*mu1*(eFrame(3,a1)*rf(2,a1) - eFrame(6,a1)*rf(1,a1))
205 <    
205 >
206      ! the potential contribution is -1/2 dipole dot reaction_field
207 <    
207 >
208      if (do_pot) then
209         rfpot = rfpot - 0.5d0 * pre * mu1 * &
210              (rf(1,a1)*eFrame(3,a1) + rf(2,a1)*eFrame(6,a1) + rf(3,a1)*eFrame(9,a1))
211      endif
212 <
212 >
213      return
214    end subroutine reaction_field_final
215 <  
215 >
216    subroutine rf_correct_forces(atom1, atom2, d, rij, eFrame, taper, f, fpair)
217 <    
217 >
218      integer, intent(in) :: atom1, atom2
219      real(kind=dp), dimension(3), intent(in) :: d
220      real(kind=dp), intent(in) :: rij, taper
221      real( kind = dp ), dimension(9,nLocal) :: eFrame
222      real( kind = dp ), dimension(3,nLocal) :: f
223      real( kind = dp ), dimension(3), intent(inout) :: fpair
224 <    
224 >
225      real (kind = dp), dimension(3) :: ul1
226      real (kind = dp), dimension(3) :: ul2
227      real (kind = dp) :: dtdr
228      real (kind = dp) :: dudx, dudy, dudz, u1dotu2
229      integer :: me1, me2, id1, id2
230      real (kind = dp) :: mu1, mu2
231 <    
231 >
232      integer :: localError
233  
234      if ((.not.haveDielectric).or.(.not.haveCutoffs)) then
# Line 237 | Line 237 | contains
237      endif
238  
239      if (rij.le.rrf) then
240 <      
240 >
241         if (rij.lt.rt) then
242            dtdr = 0.0d0
243         else
244 < !         write(*,*) 'rf correct in taper region'
244 >          !         write(*,*) 'rf correct in taper region'
245            dtdr = 6.0d0*(rij*rij - rij*rt - rij*rrf +rrf*rt)/((rrf-rt)**3)
246         endif
247 <      
247 >
248   #ifdef IS_MPI
249         me1 = atid_Row(atom1)
250         ul1(1) = eFrame_Row(3,atom1)
251         ul1(2) = eFrame_Row(6,atom1)
252         ul1(3) = eFrame_Row(9,atom1)
253 <      
253 >
254         me2 = atid_Col(atom2)
255         ul2(1) = eFrame_Col(3,atom2)
256         ul2(2) = eFrame_Col(6,atom2)
# Line 260 | Line 260 | contains
260         ul1(1) = eFrame(3,atom1)
261         ul1(2) = eFrame(6,atom1)
262         ul1(3) = eFrame(9,atom1)
263 <      
263 >
264         me2 = atid(atom2)
265         ul2(1) = eFrame(3,atom2)
266         ul2(2) = eFrame(6,atom2)
267         ul2(3) = eFrame(9,atom2)
268   #endif
269 <      
269 >
270         mu1 = getDipoleMoment(me1)
271         mu2 = getDipoleMoment(me2)
272 <      
272 >
273         u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3)
274 <      
274 >
275         dudx = - pre*mu1*mu2*u1dotu2*dtdr*d(1)/rij
276         dudy = - pre*mu1*mu2*u1dotu2*dtdr*d(2)/rij
277         dudz = - pre*mu1*mu2*u1dotu2*dtdr*d(3)/rij
278 <      
278 >
279   #ifdef IS_MPI
280         f_Row(1,atom1) = f_Row(1,atom1) + dudx
281         f_Row(2,atom1) = f_Row(2,atom1) + dudy
282         f_Row(3,atom1) = f_Row(3,atom1) + dudz
283 <      
283 >
284         f_Col(1,atom2) = f_Col(1,atom2) - dudx
285         f_Col(2,atom2) = f_Col(2,atom2) - dudy
286         f_Col(3,atom2) = f_Col(3,atom2) - dudz
# Line 288 | Line 288 | contains
288         f(1,atom1) = f(1,atom1) + dudx
289         f(2,atom1) = f(2,atom1) + dudy
290         f(3,atom1) = f(3,atom1) + dudz
291 <      
291 >
292         f(1,atom2) = f(1,atom2) - dudx
293         f(2,atom2) = f(2,atom2) - dudy
294         f(3,atom2) = f(3,atom2) - dudz
# Line 301 | Line 301 | contains
301         id1 = atom1
302         id2 = atom2
303   #endif
304 <      
304 >
305         if (molMembershipList(id1) .ne. molMembershipList(id2)) then
306 <          
306 >
307            fpair(1) = fpair(1) + dudx
308            fpair(2) = fpair(2) + dudy
309            fpair(3) = fpair(3) + dudz
310 <          
310 >
311         endif
312 <      
312 >
313      end if
314      return
315    end subroutine rf_correct_forces

Diff Legend

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