ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/calc_sticky_pair.F90
Revision: 730
Committed: Wed Aug 27 16:25:11 2003 UTC (20 years, 10 months ago) by gezelter
File size: 14039 byte(s)
Log Message:
More fixes for stress tensor parallel bug.

File Contents

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