ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/calc_reaction_field.F90
Revision: 901
Committed: Tue Jan 6 19:49:18 2004 UTC (20 years, 6 months ago) by chuckv
File size: 10217 byte(s)
Log Message:
performance fixes in the dipole dipole and reaction field code

File Contents

# User Rev Content
1 mmeineke 377 module reaction_field
2     use force_globals
3     use definitions
4     use atype_module
5     use vector_class
6 chuckv 460 use simulation
7 chuckv 901 use status
8 mmeineke 377 #ifdef IS_MPI
9     use mpiSimulation
10     #endif
11     implicit none
12    
13 mmeineke 626 PRIVATE
14    
15     real(kind=dp), save :: rrf = 1.0_dp
16 mmeineke 377 real(kind=dp), save :: rt
17 mmeineke 626 real(kind=dp), save :: dielect = 1.0_dp
18     real(kind=dp), save :: rrfsq = 1.0_dp
19 mmeineke 377 real(kind=dp), save :: pre
20 chuckv 901 logical, save :: haveCutoffs = .false.
21     logical, save :: haveMomentMap = .false.
22     logical, save :: haveDielectric = .false.
23 mmeineke 377
24 chuckv 901 type :: MomentList
25     real(kind=DP) :: dipole_moment = 0.0_DP
26     end type MomentList
27    
28     type(MomentList), dimension(:),allocatable :: MomentMap
29    
30 mmeineke 626 PUBLIC::initialize_rf
31     PUBLIC::setCutoffsRF
32     PUBLIC::accumulate_rf
33     PUBLIC::accumulate_self_rf
34     PUBLIC::reaction_field_final
35     PUBLIC::rf_correct_forces
36    
37 mmeineke 377 contains
38    
39 mmeineke 626 subroutine initialize_rf(this_dielect)
40     real(kind=dp), intent(in) :: this_dielect
41    
42     dielect = this_dielect
43 mmeineke 377
44 mmeineke 626 pre = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf)
45    
46 chuckv 901 haveDielectric = .true.
47 mmeineke 626
48     return
49     end subroutine initialize_rf
50    
51     subroutine setCutoffsRF( this_rrf, this_rt )
52    
53     real(kind=dp), intent(in) :: this_rrf, this_rt
54    
55 mmeineke 377 rrf = this_rrf
56 mmeineke 626 rt = this_rt
57    
58 mmeineke 377 rrfsq = rrf * rrf
59     pre = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf)
60    
61 chuckv 901 haveCutoffs = .true.
62 gezelter 394
63 mmeineke 626 end subroutine setCutoffsRF
64    
65 chuckv 901 subroutine createMomentMap(status)
66     integer :: nAtypes
67     integer :: status
68     integer :: i
69     real (kind=DP) :: thisDP
70     logical :: thisProperty
71    
72     status = 0
73    
74     nAtypes = getSize(atypes)
75    
76     if (nAtypes == 0) then
77     status = -1
78     return
79     end if
80    
81     if (.not. allocated(MomentMap)) then
82     allocate(MomentMap(nAtypes))
83     endif
84    
85     do i = 1, nAtypes
86    
87     call getElementProperty(atypes, i, "is_DP", thisProperty)
88    
89     if (thisProperty) then
90     call getElementProperty(atypes, i, "dipole_moment", thisDP)
91     MomentMap(i)%dipole_moment = thisDP
92     endif
93    
94     end do
95    
96     haveMomentMap = .true.
97    
98     end subroutine createMomentMap
99    
100 mmeineke 377 subroutine accumulate_rf(atom1, atom2, rij, u_l)
101    
102     integer, intent(in) :: atom1, atom2
103     real (kind = dp), intent(in) :: rij
104 chuckv 898 real (kind = dp), dimension(3,nLocal) :: u_l
105 mmeineke 377
106     integer :: me1, me2
107     real (kind = dp) :: taper, mu1, mu2
108     real (kind = dp), dimension(3) :: ul1
109     real (kind = dp), dimension(3) :: ul2
110    
111 chuckv 901 integer :: localError
112    
113     if ((.not.haveDielectric).or.(.not.haveCutoffs)) then
114 mmeineke 377 write(default_error,*) 'Reaction field not initialized!'
115     return
116     endif
117 chuckv 901
118     if (.not.haveMomentMap) then
119     localError = 0
120     call createMomentMap(localError)
121     if ( localError .ne. 0 ) then
122     call handleError("reaction-field", "MomentMap creation failed!")
123     return
124     end if
125     endif
126 mmeineke 377
127     if (rij.le.rrf) then
128    
129     if (rij.lt.rt) then
130     taper = 1.0d0
131     else
132 gezelter 621 write(*,*) 'rf in taper region'
133 mmeineke 377 taper = (rrf + 2.0d0*rij - 3.0d0*rt)*(rrf-rij)**2/ ((rrf-rt)**3)
134     endif
135    
136     #ifdef IS_MPI
137     me1 = atid_Row(atom1)
138     ul1(1) = u_l_Row(1,atom1)
139     ul1(2) = u_l_Row(2,atom1)
140     ul1(3) = u_l_Row(3,atom1)
141    
142     me2 = atid_Col(atom2)
143     ul2(1) = u_l_Col(1,atom2)
144     ul2(2) = u_l_Col(2,atom2)
145     ul2(3) = u_l_Col(3,atom2)
146     #else
147     me1 = atid(atom1)
148     ul1(1) = u_l(1,atom1)
149     ul1(2) = u_l(2,atom1)
150     ul1(3) = u_l(3,atom1)
151    
152     me2 = atid(atom2)
153     ul2(1) = u_l(1,atom2)
154     ul2(2) = u_l(2,atom2)
155     ul2(3) = u_l(3,atom2)
156     #endif
157    
158 chuckv 901 mu1 = MomentMap(me1)%dipole_moment
159     mu2 = MomentMap(me2)%dipole_moment
160 mmeineke 377
161     #ifdef IS_MPI
162     rf_Row(1,atom1) = rf_Row(1,atom1) + ul2(1)*mu2*taper
163     rf_Row(2,atom1) = rf_Row(2,atom1) + ul2(2)*mu2*taper
164     rf_Row(3,atom1) = rf_Row(3,atom1) + ul2(3)*mu2*taper
165    
166 gezelter 394 rf_Col(1,atom2) = rf_Col(1,atom2) + ul1(1)*mu1*taper
167     rf_Col(2,atom2) = rf_Col(2,atom2) + ul1(2)*mu1*taper
168     rf_Col(3,atom2) = rf_Col(3,atom2) + ul1(3)*mu1*taper
169 mmeineke 377 #else
170     rf(1,atom1) = rf(1,atom1) + ul2(1)*mu2*taper
171     rf(2,atom1) = rf(2,atom1) + ul2(2)*mu2*taper
172     rf(3,atom1) = rf(3,atom1) + ul2(3)*mu2*taper
173    
174 gezelter 394 rf(1,atom2) = rf(1,atom2) + ul1(1)*mu1*taper
175     rf(2,atom2) = rf(2,atom2) + ul1(2)*mu1*taper
176     rf(3,atom2) = rf(3,atom2) + ul1(3)*mu1*taper
177 mmeineke 377 #endif
178    
179     endif
180     return
181     end subroutine accumulate_rf
182    
183     subroutine accumulate_self_rf(atom1, mu1, u_l)
184    
185     integer, intent(in) :: atom1
186     real(kind=dp), intent(in) :: mu1
187 chuckv 898 real(kind=dp), dimension(3,nLocal) :: u_l
188 mmeineke 377
189     !! should work for both MPI and non-MPI version since this is not pairwise.
190     rf(1,atom1) = rf(1,atom1) + u_l(1,atom1)*mu1
191     rf(2,atom1) = rf(2,atom1) + u_l(2,atom1)*mu1
192     rf(3,atom1) = rf(3,atom1) + u_l(3,atom1)*mu1
193    
194     return
195     end subroutine accumulate_self_rf
196    
197     subroutine reaction_field_final(a1, mu1, u_l, rfpot, t, do_pot)
198    
199     integer, intent(in) :: a1
200     real (kind=dp), intent(in) :: mu1
201     real (kind=dp), intent(inout) :: rfpot
202     logical, intent(in) :: do_pot
203 chuckv 898 real (kind = dp), dimension(3,nLocal) :: u_l
204     real (kind = dp), dimension(3,nLocal) :: t
205 mmeineke 377
206 chuckv 901 integer :: localError
207    
208     if ((.not.haveDielectric).or.(.not.haveCutoffs)) then
209 mmeineke 377 write(default_error,*) 'Reaction field not initialized!'
210     return
211     endif
212    
213 chuckv 901 if (.not.haveMomentMap) then
214     localError = 0
215     call createMomentMap(localError)
216     if ( localError .ne. 0 ) then
217     call handleError("reaction-field", "MomentMap creation failed!")
218     return
219     end if
220     endif
221    
222 mmeineke 377 ! compute torques on dipoles:
223     ! pre converts from mu in units of debye to kcal/mol
224    
225 gezelter 394 ! The torque contribution is dipole cross reaction_field
226    
227 mmeineke 377 t(1,a1) = t(1,a1) + pre*mu1*(u_l(2,a1)*rf(3,a1) - u_l(3,a1)*rf(2,a1))
228     t(2,a1) = t(2,a1) + pre*mu1*(u_l(3,a1)*rf(1,a1) - u_l(1,a1)*rf(3,a1))
229     t(3,a1) = t(3,a1) + pre*mu1*(u_l(1,a1)*rf(2,a1) - u_l(2,a1)*rf(1,a1))
230    
231     ! the potential contribution is -1/2 dipole dot reaction_field
232    
233     if (do_pot) then
234     rfpot = rfpot - 0.5d0 * pre * mu1 * &
235     (rf(1,a1)*u_l(1,a1) + rf(2,a1)*u_l(2,a1) + rf(3,a1)*u_l(3,a1))
236     endif
237    
238     return
239     end subroutine reaction_field_final
240    
241     subroutine rf_correct_forces(atom1, atom2, d, rij, u_l, f, do_stress)
242    
243     integer, intent(in) :: atom1, atom2
244     real(kind=dp), dimension(3), intent(in) :: d
245     real(kind=dp), intent(in) :: rij
246 chuckv 898 real( kind = dp ), dimension(3,nLocal) :: u_l
247     real( kind = dp ), dimension(3,nLocal) :: f
248 mmeineke 377 logical, intent(in) :: do_stress
249    
250     real (kind = dp), dimension(3) :: ul1
251     real (kind = dp), dimension(3) :: ul2
252     real (kind = dp) :: dtdr
253     real (kind = dp) :: dudx, dudy, dudz, u1dotu2
254 gezelter 730 integer :: me1, me2, id1, id2
255 mmeineke 377 real (kind = dp) :: mu1, mu2
256    
257 chuckv 901 integer :: localError
258    
259     if ((.not.haveDielectric).or.(.not.haveCutoffs)) then
260 mmeineke 377 write(default_error,*) 'Reaction field not initialized!'
261     return
262     endif
263    
264 chuckv 901 if (.not.haveMomentMap) then
265     localError = 0
266     call createMomentMap(localError)
267     if ( localError .ne. 0 ) then
268     call handleError("reaction-field", "MomentMap creation failed!")
269     return
270     end if
271     endif
272    
273 mmeineke 377 if (rij.le.rrf) then
274    
275     if (rij.lt.rt) then
276     dtdr = 0.0d0
277     else
278 gezelter 621 write(*,*) 'rf correct in taper region'
279 mmeineke 377 dtdr = 6.0d0*(rij*rij - rij*rt - rij*rrf +rrf*rt)/((rrf-rt)**3)
280     endif
281    
282     #ifdef IS_MPI
283     me1 = atid_Row(atom1)
284     ul1(1) = u_l_Row(1,atom1)
285     ul1(2) = u_l_Row(2,atom1)
286     ul1(3) = u_l_Row(3,atom1)
287    
288     me2 = atid_Col(atom2)
289     ul2(1) = u_l_Col(1,atom2)
290     ul2(2) = u_l_Col(2,atom2)
291     ul2(3) = u_l_Col(3,atom2)
292     #else
293     me1 = atid(atom1)
294     ul1(1) = u_l(1,atom1)
295     ul1(2) = u_l(2,atom1)
296     ul1(3) = u_l(3,atom1)
297    
298     me2 = atid(atom2)
299     ul2(1) = u_l(1,atom2)
300     ul2(2) = u_l(2,atom2)
301     ul2(3) = u_l(3,atom2)
302     #endif
303    
304 chuckv 901 mu1 = MomentMap(me1)%dipole_moment
305     mu2 = MomentMap(me2)%dipole_moment
306 mmeineke 377
307     u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3)
308    
309     dudx = - pre*mu1*mu2*u1dotu2*dtdr*d(1)/rij
310     dudy = - pre*mu1*mu2*u1dotu2*dtdr*d(2)/rij
311     dudz = - pre*mu1*mu2*u1dotu2*dtdr*d(3)/rij
312    
313     #ifdef IS_MPI
314     f_Row(1,atom1) = f_Row(1,atom1) + dudx
315     f_Row(2,atom1) = f_Row(2,atom1) + dudy
316     f_Row(3,atom1) = f_Row(3,atom1) + dudz
317    
318     f_Col(1,atom2) = f_Col(1,atom2) - dudx
319     f_Col(2,atom2) = f_Col(2,atom2) - dudy
320     f_Col(3,atom2) = f_Col(3,atom2) - dudz
321     #else
322     f(1,atom1) = f(1,atom1) + dudx
323     f(2,atom1) = f(2,atom1) + dudy
324     f(3,atom1) = f(3,atom1) + dudz
325    
326     f(1,atom2) = f(1,atom2) - dudx
327     f(2,atom2) = f(2,atom2) - dudy
328     f(3,atom2) = f(3,atom2) - dudz
329     #endif
330    
331     if (do_stress) then
332 gezelter 611
333 gezelter 730 #ifdef IS_MPI
334     id1 = tagRow(atom1)
335     id2 = tagColumn(atom2)
336     #else
337     id1 = atom1
338     id2 = atom2
339     #endif
340    
341     if (molMembershipList(id1) .ne. molMembershipList(id2)) then
342    
343 gezelter 611 ! because the d vector is the rj - ri vector, and
344     ! because dudx, dudy, and dudz are the
345     ! (positive) force on atom i (negative on atom j) we need
346     ! a negative sign here:
347    
348     tau_Temp(1) = tau_Temp(1) - d(1) * dudx
349     tau_Temp(2) = tau_Temp(2) - d(1) * dudy
350     tau_Temp(3) = tau_Temp(3) - d(1) * dudz
351     tau_Temp(4) = tau_Temp(4) - d(2) * dudx
352     tau_Temp(5) = tau_Temp(5) - d(2) * dudy
353     tau_Temp(6) = tau_Temp(6) - d(2) * dudz
354     tau_Temp(7) = tau_Temp(7) - d(3) * dudx
355     tau_Temp(8) = tau_Temp(8) - d(3) * dudy
356     tau_Temp(9) = tau_Temp(9) - d(3) * dudz
357 gezelter 483 virial_Temp = virial_Temp + &
358     (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
359     endif
360 mmeineke 377 endif
361     endif
362    
363     return
364     end subroutine rf_correct_forces
365     end module reaction_field