ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/calc_sticky_pair.F90
Revision: 324
Committed: Wed Mar 12 16:38:17 2003 UTC (21 years, 5 months ago) by gezelter
File size: 12655 byte(s)
Log Message: 
MPI bug fix

File Contents

# User Rev Content
1 gezelter 319 !! This Module Calculates forces due to SSD potential and VDW interactions
2     !! [Chandra and Ichiye, J. Chem. Phys. 111, 2701 (1999)].
3    
4     !! This module contains the Public procedures:
5    
6    
7     !! Corresponds to the force field defined in ssd_FF.cpp
8     !! @author Charles F. Vardeman II
9     !! @author Matthew Meineke
10 gezelter 322 !! @author Christopher Fennel
11     !! @author J. Daniel Gezelter
12 gezelter 324 !! @version $Id: calc_sticky_pair.F90,v 1.4 2003-03-12 16:38:17 gezelter Exp $, $Date: 2003-03-12 16:38:17 $, $Name: not supported by cvs2svn $, $Revision: 1.4 $
13 gezelter 319
14 gezelter 322 module sticky_pair
15 gezelter 319
16     use simulation
17     use definitions
18 gezelter 322 use forceGlobals
19 gezelter 319 #ifdef IS_MPI
20     use mpiSimulation
21     #endif
22 gezelter 322
23 gezelter 319 implicit none
24 gezelter 322
25 gezelter 319 PRIVATE
26    
27 gezelter 322 logical, save :: sticky_initialized = .false.
28 gezelter 323 real( kind = dp ), save :: SSD_w0
29 gezelter 322 real( kind = dp ), save :: SSD_v0
30     real( kind = dp ), save :: SSD_rl
31     real( kind = dp ), save :: SSD_ru
32     real( kind = dp ), save :: SSD_rup
33 gezelter 319
34 gezelter 323 public :: check_sticky_FF
35     public :: set_sticky_params
36 gezelter 322 public :: do_sticky_pair
37 gezelter 319
38     contains
39    
40 gezelter 323 subroutine check_sticky_FF(status)
41 gezelter 324 integer :: status
42     status = -1
43 gezelter 323 if (sticky_initialized) status = 0
44     return
45     end subroutine check_sticky_FF
46    
47     subroutine set_sticky_params(sticky_w0, sticky_v0)
48 gezelter 322 real( kind = dp ), intent(in) :: sticky_w0, sticky_v0
49    
50     ! we could pass all 5 parameters if we felt like it...
51 gezelter 323
52 gezelter 322 SSD_w0 = sticky_w0
53     SSD_v0 = sticky_v0
54     SSD_rl = 2.75_DP
55     SSD_ru = 3.35_DP
56     SSD_rup = 4.0_DP
57 gezelter 323
58 gezelter 322 sticky_initialized = .true.
59     return
60 gezelter 323 end subroutine set_sticky_params
61 gezelter 319
62 gezelter 323 subroutine do_sticky_pair(atom1, atom2, d, rij, r2, A, pot, f, t)
63 gezelter 322
64     !! This routine does only the sticky portion of the SSD potential
65     !! [Chandra and Ichiye, J. Chem. Phys. 111, 2701 (1999)].
66     !! The Lennard-Jones and dipolar interaction must be handled separately.
67    
68     !! We assume that the rotation matrices have already been calculated
69     !! and placed in the A array.
70    
71     !! i and j are pointers to the two SSD atoms
72 gezelter 319
73 gezelter 322 integer, intent(in) :: atom1, atom2
74 gezelter 323 real (kind=dp), intent(inout) :: rij, r2
75 gezelter 322 real (kind=dp), dimension(3), intent(in) :: d
76 gezelter 323 real (kind=dp) :: pot
77 gezelter 322 real (kind=dp), dimension(9,getNlocal()) :: A
78     real (kind=dp), dimension(3,getNlocal()) :: f
79     real (kind=dp), dimension(3,getNlocal()) :: t
80 gezelter 319
81 gezelter 322 real (kind=dp) :: xi, yi, zi, xj, yj, zj, xi2, yi2, zi2, xj2, yj2, zj2
82 gezelter 323 real (kind=dp) :: r3, r5, r6, s, sp, dsdr, dspdr
83 gezelter 322 real (kind=dp) :: wi, wj, w, wip, wjp, wp
84     real (kind=dp) :: dwidx, dwidy, dwidz, dwjdx, dwjdy, dwjdz
85     real (kind=dp) :: dwipdx, dwipdy, dwipdz, dwjpdx, dwjpdy, dwjpdz
86     real (kind=dp) :: dwidux, dwiduy, dwiduz, dwjdux, dwjduy, dwjduz
87     real (kind=dp) :: dwipdux, dwipduy, dwipduz, dwjpdux, dwjpduy, dwjpduz
88     real (kind=dp) :: zif, zis, zjf, zjs, uglyi, uglyj
89     real (kind=dp) :: drdx, drdy, drdz
90     real (kind=dp) :: txi, tyi, tzi, txj, tyj, tzj
91     real (kind=dp) :: fxii, fyii, fzii, fxjj, fyjj, fzjj
92     real (kind=dp) :: fxij, fyij, fzij, fxji, fyji, fzji
93     real (kind=dp) :: fxradial, fyradial, fzradial
94    
95     if (.not.sticky_initialized) then
96     write(*,*) 'Sticky forces not initialized!'
97 gezelter 319 return
98     endif
99    
100 gezelter 322 r3 = r2*rij
101     r5 = r3*r2
102    
103     drdx = d(1) / rij
104     drdy = d(2) / rij
105     drdz = d(3) / rij
106    
107 gezelter 319 #ifdef IS_MPI
108 gezelter 322 ! rotate the inter-particle separation into the two different
109     ! body-fixed coordinate systems:
110 gezelter 319
111 gezelter 322 xi = A_row(1,atom1)*d(1) + A_row(2,atom1)*d(2) + A_row(3,atom1)*d(3)
112     yi = A_row(4,atom1)*d(1) + A_row(5,atom1)*d(2) + A_row(6,atom1)*d(3)
113     zi = A_row(7,atom1)*d(1) + A_row(8,atom1)*d(2) + A_row(9,atom1)*d(3)
114 gezelter 319
115 gezelter 322 ! negative sign because this is the vector from j to i:
116    
117     xj = -(A_Col(1,atom2)*d(1) + A_Col(2,atom2)*d(2) + A_Col(3,atom2)*d(3))
118     yj = -(A_Col(4,atom2)*d(1) + A_Col(5,atom2)*d(2) + A_Col(6,atom2)*d(3))
119     zj = -(A_Col(7,atom2)*d(1) + A_Col(8,atom2)*d(2) + A_Col(9,atom2)*d(3))
120 gezelter 319 #else
121 gezelter 322 ! rotate the inter-particle separation into the two different
122     ! body-fixed coordinate systems:
123 gezelter 319
124 gezelter 322 xi = a(1,atom1)*d(1) + a(2,atom1)*d(2) + a(3,atom1)*d(3)
125     yi = a(4,atom1)*d(1) + a(5,atom1)*d(2) + a(6,atom1)*d(3)
126     zi = a(7,atom1)*d(1) + a(8,atom1)*d(2) + a(9,atom1)*d(3)
127    
128     ! negative sign because this is the vector from j to i:
129    
130     xj = -(a(1,atom2)*d(1) + a(2,atom2)*d(2) + a(3,atom2)*d(3))
131     yj = -(a(4,atom2)*d(1) + a(5,atom2)*d(2) + a(6,atom2)*d(3))
132     zj = -(a(7,atom2)*d(1) + a(8,atom2)*d(2) + a(9,atom2)*d(3))
133 gezelter 319 #endif
134    
135 gezelter 322 xi2 = xi*xi
136     yi2 = yi*yi
137     zi2 = zi*zi
138 gezelter 319
139 gezelter 322 xj2 = xj*xj
140     yj2 = yj*yj
141     zj2 = zj*zj
142 gezelter 319
143 gezelter 322 call calc_sw_fnc(rij, s, sp, dsdr, dspdr)
144    
145     wi = 2.0d0*(xi2-yi2)*zi / r3
146     wj = 2.0d0*(xj2-yj2)*zj / r3
147     w = wi+wj
148    
149     zif = zi/rij - 0.6d0
150     zis = zi/rij + 0.8d0
151    
152     zjf = zj/rij - 0.6d0
153     zjs = zj/rij + 0.8d0
154    
155     wip = zif*zif*zis*zis - SSD_w0
156     wjp = zjf*zjf*zjs*zjs - SSD_w0
157     wp = wip + wjp
158 gezelter 319
159 gezelter 322 #ifdef IS_MPI
160     pot_row(atom1) = pot_row(atom1) + 0.25d0*SSD_v0*(s*w + sp*wp)
161     pot_col(atom2) = pot_col(atom2) + 0.25d0*SSD_v0*(s*w + sp*wp)
162     #else
163     pot = pot + 0.5d0*SSD_v0*(s*w + sp*wp)
164     #endif
165 gezelter 319
166 gezelter 322 dwidx = 4.0d0*xi*zi/r3 - 6.0d0*xi*zi*(xi2-yi2)/r5
167     dwidy = - 4.0d0*yi*zi/r3 - 6.0d0*yi*zi*(xi2-yi2)/r5
168     dwidz = 2.0d0*(xi2-yi2)/r3 - 6.0d0*zi2*(xi2-yi2)/r5
169    
170     dwjdx = 4.0d0*xj*zj/r3 - 6.0d0*xj*zj*(xj2-yj2)/r5
171     dwjdy = - 4.0d0*yj*zj/r3 - 6.0d0*yj*zj*(xj2-yj2)/r5
172     dwjdz = 2.0d0*(xj2-yj2)/r3 - 6.0d0*zj2*(xj2-yj2)/r5
173    
174     uglyi = zif*zif*zis + zif*zis*zis
175     uglyj = zjf*zjf*zjs + zjf*zjs*zjs
176    
177     dwipdx = -2.0d0*xi*zi*uglyi/r3
178     dwipdy = -2.0d0*yi*zi*uglyi/r3
179     dwipdz = 2.0d0*(1.0d0/rij - zi2/r3)*uglyi
180    
181     dwjpdx = -2.0d0*xj*zj*uglyj/r3
182     dwjpdy = -2.0d0*yj*zj*uglyj/r3
183     dwjpdz = 2.0d0*(1.0d0/rij - zj2/r3)*uglyj
184    
185     dwidux = 4.0d0*(yi*zi2 + 0.5d0*yi*(xi2-yi2))/r3
186     dwiduy = 4.0d0*(xi*zi2 - 0.5d0*xi*(xi2-yi2))/r3
187     dwiduz = - 8.0d0*xi*yi*zi/r3
188    
189     dwjdux = 4.0d0*(yj*zj2 + 0.5d0*yj*(xj2-yj2))/r3
190     dwjduy = 4.0d0*(xj*zj2 - 0.5d0*xj*(xj2-yj2))/r3
191     dwjduz = - 8.0d0*xj*yj*zj/r3
192    
193     dwipdux = 2.0d0*yi*uglyi/rij
194     dwipduy = -2.0d0*xi*uglyi/rij
195     dwipduz = 0.0d0
196    
197     dwjpdux = 2.0d0*yj*uglyj/rij
198     dwjpduy = -2.0d0*xj*uglyj/rij
199     dwjpduz = 0.0d0
200    
201     ! do the torques first since they are easy:
202     ! remember that these are still in the body fixed axes
203    
204     txi = 0.5d0*SSD_v0*(s*dwidux + sp*dwipdux)
205     tyi = 0.5d0*SSD_v0*(s*dwiduy + sp*dwipduy)
206     tzi = 0.5d0*SSD_v0*(s*dwiduz + sp*dwipduz)
207    
208     txj = 0.5d0*SSD_v0*(s*dwjdux + sp*dwjpdux)
209     tyj = 0.5d0*SSD_v0*(s*dwjduy + sp*dwjpduy)
210     tzj = 0.5d0*SSD_v0*(s*dwjduz + sp*dwjpduz)
211    
212     ! go back to lab frame using transpose of rotation matrix:
213 gezelter 319
214     #ifdef IS_MPI
215 gezelter 322 t_Row(1,atom1) = t_Row(1,atom1) + a_Row(1,atom1)*txi + &
216     a_Row(4,atom1)*tyi + a_Row(7,atom1)*tzi
217     t_Row(2,atom1) = t_Row(2,atom1) + a_Row(2,atom1)*txi + &
218     a_Row(5,atom1)*tyi + a_Row(8,atom1)*tzi
219     t_Row(3,atom1) = t_Row(3,atom1) + a_Row(3,atom1)*txi + &
220     a_Row(6,atom1)*tyi + a_Row(9,atom1)*tzi
221 gezelter 319
222 gezelter 322 t_Col(1,atom2) = t_Col(1,atom2) + a_Col(1,atom2)*txj + &
223     a_Col(4,atom2)*tyj + a_Col(7,atom2)*tzj
224     t_Col(2,atom2) = t_Col(2,atom2) + a_Col(2,atom2)*txj + &
225     a_Col(5,atom2)*tyj + a_Col(8,atom2)*tzj
226     t_Col(3,atom2) = t_Col(3,atom2) + a_Col(3,atom2)*txj + &
227     a_Col(6,atom2)*tyj + a_Col(9,atom2)*tzj
228 gezelter 319 #else
229 gezelter 322 t(1,atom1) = t(1,atom1) + a(1,atom1)*txi + a(4,atom1)*tyi + a(7,atom1)*tzi
230     t(2,atom1) = t(2,atom1) + a(2,atom1)*txi + a(5,atom1)*tyi + a(8,atom1)*tzi
231     t(3,atom1) = t(3,atom1) + a(3,atom1)*txi + a(6,atom1)*tyi + a(9,atom1)*tzi
232    
233     t(1,atom2) = t(1,atom2) + a(1,atom2)*txj + a(4,atom2)*tyj + a(7,atom2)*tzj
234     t(2,atom2) = t(2,atom2) + a(2,atom2)*txj + a(5,atom2)*tyj + a(8,atom2)*tzj
235     t(3,atom2) = t(3,atom2) + a(3,atom2)*txj + a(6,atom2)*tyj + a(9,atom2)*tzj
236     #endif
237     ! Now, on to the forces:
238    
239     ! first rotate the i terms back into the lab frame:
240 gezelter 319
241 gezelter 322 #ifdef IS_MPI
242     fxii = a_Row(1,atom1)*(s*dwidx+sp*dwipdx) + &
243     a_Row(4,atom1)*(s*dwidy+sp*dwipdy) + &
244     a_Row(7,atom1)*(s*dwidz+sp*dwipdz)
245     fyii = a_Row(2,atom1)*(s*dwidx+sp*dwipdx) + &
246     a_Row(5,atom1)*(s*dwidy+sp*dwipdy) + &
247     a_Row(8,atom1)*(s*dwidz+sp*dwipdz)
248     fzii = a_Row(3,atom1)*(s*dwidx+sp*dwipdx) + &
249     a_Row(6,atom1)*(s*dwidy+sp*dwipdy) + &
250     a_Row(9,atom1)*(s*dwidz+sp*dwipdz)
251 gezelter 319
252 gezelter 322 fxjj = a_Col(1,atom2)*(s*dwjdx+sp*dwjpdx) + &
253     a_Col(4,atom2)*(s*dwjdy+sp*dwjpdy) + &
254     a_Col(7,atom2)*(s*dwjdz+sp*dwjpdz)
255     fyjj = a_Col(2,atom2)*(s*dwjdx+sp*dwjpdx) + &
256     a_Col(5,atom2)*(s*dwjdy+sp*dwjpdy) + &
257     a_Col(8,atom2)*(s*dwjdz+sp*dwjpdz)
258     fzjj = a_Col(3,atom2)*(s*dwjdx+sp*dwjpdx)+ &
259     a_Col(6,atom2)*(s*dwjdy+sp*dwjpdy) + &
260     a_Col(9,atom2)*(s*dwjdz+sp*dwjpdz)
261 gezelter 319 #else
262 gezelter 322 fxii = a(1,atom1)*(s*dwidx+sp*dwipdx) + &
263     a(4,atom1)*(s*dwidy+sp*dwipdy) + &
264     a(7,atom1)*(s*dwidz+sp*dwipdz)
265     fyii = a(2,atom1)*(s*dwidx+sp*dwipdx) + &
266     a(5,atom1)*(s*dwidy+sp*dwipdy) + &
267     a(8,atom1)*(s*dwidz+sp*dwipdz)
268     fzii = a(3,atom1)*(s*dwidx+sp*dwipdx) + &
269     a(6,atom1)*(s*dwidy+sp*dwipdy) + &
270     a(9,atom1)*(s*dwidz+sp*dwipdz)
271 gezelter 319
272 gezelter 322 fxjj = a(1,atom2)*(s*dwjdx+sp*dwjpdx) + &
273     a(4,atom2)*(s*dwjdy+sp*dwjpdy) + &
274     a(7,atom2)*(s*dwjdz+sp*dwjpdz)
275     fyjj = a(2,atom2)*(s*dwjdx+sp*dwjpdx) + &
276     a(5,atom2)*(s*dwjdy+sp*dwjpdy) + &
277     a(8,atom2)*(s*dwjdz+sp*dwjpdz)
278     fzjj = a(3,atom2)*(s*dwjdx+sp*dwjpdx)+ &
279     a(6,atom2)*(s*dwjdy+sp*dwjpdy) + &
280     a(9,atom2)*(s*dwjdz+sp*dwjpdz)
281 gezelter 319 #endif
282 gezelter 322
283     fxij = -fxii
284     fyij = -fyii
285     fzij = -fzii
286    
287     fxji = -fxjj
288     fyji = -fyjj
289     fzji = -fzjj
290    
291     ! now assemble these with the radial-only terms:
292 gezelter 319
293 gezelter 322 fxradial = 0.5d0*SSD_v0*(dsdr*drdx*w + dspdr*drdx*wp + fxii + fxji)
294     fyradial = 0.5d0*SSD_v0*(dsdr*drdy*w + dspdr*drdy*wp + fyii + fyji)
295     fzradial = 0.5d0*SSD_v0*(dsdr*drdz*w + dspdr*drdz*wp + fzii + fzji)
296 gezelter 319
297     #ifdef IS_MPI
298 gezelter 322 f_Row(1,atom1) = f_Row(1,atom1) + fxradial
299     f_Row(2,atom1) = f_Row(2,atom1) + fyradial
300     f_Row(3,atom1) = f_Row(3,atom1) + fzradial
301    
302     f_Col(1,atom2) = f_Col(1,atom2) + 0.5d0*SSD_v0*(-dsdr*drdx*w - &
303     dspdr*drdx*wp + fxjj + fxij)
304     f_Col(2,atom2) = f_Col(2,atom2) + 0.5d0*SSD_v0*(-dsdr*drdy*w - &
305     dspdr*drdy*wp + fyjj + fyij)
306     f_Col(3,atom2) = f_Col(3,atom2) + 0.5d0*SSD_v0*(-dsdr*drdz*w - &
307     dspdr*drdz*wp + fzjj + fzij)
308 gezelter 319 #else
309 gezelter 322 f(1,atom1) = f(1,atom1) + fxradial
310     f(2,atom1) = f(2,atom1) + fyradial
311     f(3,atom1) = f(3,atom1) + fzradial
312    
313     f(1,atom2) = f(1,atom2) + 0.5d0*SSD_v0*(-dsdr*drdx*w - dspdr*drdx*wp + &
314     fxjj + fxij)
315     f(2,atom2) = f(2,atom2) + 0.5d0*SSD_v0*(-dsdr*drdy*w - dspdr*drdy*wp + &
316     fyjj + fyij)
317     f(3,atom2) = f(3,atom2) + 0.5d0*SSD_v0*(-dsdr*drdz*w - dspdr*drdz*wp + &
318     fzjj + fzij)
319 gezelter 319 #endif
320 gezelter 322
321     if (doStress()) then
322     tau_Temp(1) = tau_Temp(1) + fxradial * d(1)
323     tau_Temp(2) = tau_Temp(2) + fxradial * d(2)
324     tau_Temp(3) = tau_Temp(3) + fxradial * d(3)
325     tau_Temp(4) = tau_Temp(4) + fyradial * d(1)
326     tau_Temp(5) = tau_Temp(5) + fyradial * d(2)
327     tau_Temp(6) = tau_Temp(6) + fyradial * d(3)
328     tau_Temp(7) = tau_Temp(7) + fzradial * d(1)
329     tau_Temp(8) = tau_Temp(8) + fzradial * d(2)
330     tau_Temp(9) = tau_Temp(9) + fzradial * d(3)
331     virial_Temp = virial_Temp + (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
332     endif
333    
334     end subroutine do_sticky_pair
335 gezelter 319
336 gezelter 322 !! calculates the switching functions and their derivatives for a given
337     subroutine calc_sw_fnc(r, s, sp, dsdr, dspdr)
338 gezelter 319
339 gezelter 322 real (kind=dp), intent(in) :: r
340     real (kind=dp), intent(inout) :: s, sp, dsdr, dspdr
341 gezelter 319
342 gezelter 322 ! distances must be in angstroms
343 gezelter 319
344 gezelter 322 if (r.lt.SSD_rl) then
345     s = 1.0d0
346     sp = 1.0d0
347     dsdr = 0.0d0
348     dspdr = 0.0d0
349     elseif (r.gt.SSD_rup) then
350     s = 0.0d0
351     sp = 0.0d0
352     dsdr = 0.0d0
353     dspdr = 0.0d0
354     else
355     sp = ((SSD_rup + 2.0d0*r - 3.0d0*SSD_rl) * (SSD_rup-r)**2) / &
356     ((SSD_rup - SSD_rl)**3)
357     dspdr = 6.0d0*(r-SSD_rup)*(r-SSD_rl)/((SSD_rup - SSD_rl)**3)
358 gezelter 319
359 gezelter 322 if (r.gt.SSD_ru) then
360     s = 0.0d0
361     dsdr = 0.0d0
362 gezelter 319 else
363 gezelter 322 s = ((SSD_ru + 2.0d0*r - 3.0d0*SSD_rl) * (SSD_ru-r)**2) / &
364     ((SSD_ru - SSD_rl)**3)
365     dsdr = 6.0d0*(r-SSD_ru)*(r-SSD_rl)/((SSD_ru - SSD_rl)**3)
366 gezelter 319 endif
367 gezelter 322 endif
368 gezelter 319
369 gezelter 322 return
370     end subroutine calc_sw_fnc
371     end module sticky_pair