ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/calc_sticky_pair.F90
(Generate patch)

Comparing trunk/OOPSE/libmdtools/calc_sticky_pair.F90 (file contents):
Revision 611 by gezelter, Tue Jul 15 17:10:50 2003 UTC vs.
Revision 1192 by gezelter, Mon May 24 21:03:30 2004 UTC

# Line 9 | Line 9
9   !! @author Matthew Meineke
10   !! @author Christopher Fennel
11   !! @author J. Daniel Gezelter
12 < !! @version $Id: calc_sticky_pair.F90,v 1.10 2003-07-15 17:10:50 gezelter Exp $, $Date: 2003-07-15 17:10:50 $, $Name: not supported by cvs2svn $, $Revision: 1.10 $
12 > !! @version $Id: calc_sticky_pair.F90,v 1.19 2004-05-24 21:03:25 gezelter Exp $, $Date: 2004-05-24 21:03:25 $, $Name: not supported by cvs2svn $, $Revision: 1.19 $
13  
14   module sticky_pair
15  
# Line 27 | Line 27 | module sticky_pair
27    logical, save :: sticky_initialized = .false.
28    real( kind = dp ), save :: SSD_w0 = 0.0_dp
29    real( kind = dp ), save :: SSD_v0 = 0.0_dp
30 +  real( kind = dp ), save :: SSD_v0p = 0.0_dp
31    real( kind = dp ), save :: SSD_rl = 0.0_dp
32    real( kind = dp ), save :: SSD_ru = 0.0_dp
33 +  real( kind = dp ), save :: SSD_rlp = 0.0_dp
34    real( kind = dp ), save :: SSD_rup = 0.0_dp
35 +  real( kind = dp ), save :: SSD_rbig = 0.0_dp
36  
37    public :: check_sticky_FF
38    public :: set_sticky_params
# Line 44 | Line 47 | contains
47      return
48    end subroutine check_sticky_FF
49  
50 <  subroutine set_sticky_params(sticky_w0, sticky_v0)
51 <    real( kind = dp ), intent(in) :: sticky_w0, sticky_v0
50 >  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      
57      ! we could pass all 5 parameters if we felt like it...
58      
59      SSD_w0 = sticky_w0
60      SSD_v0 = sticky_v0
61 <    SSD_rl = 2.75_DP
62 <    SSD_ru = 3.35_DP
63 <    SSD_rup = 4.0_DP
61 >    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 >
67 >    if (SSD_ru .gt. SSD_rup) then
68 >       SSD_rbig = SSD_ru
69 >    else
70 >       SSD_rbig = SSD_rup
71 >    endif
72    
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 <
77 >  subroutine do_sticky_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
78 >       pot, A, f, t, do_pot)
79 >    
80      !! 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 <
83 >    
84      !! We assume that the rotation matrices have already been calculated
85      !! and placed in the A array.
86  
# Line 74 | Line 89 | contains
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 <    real (kind=dp), dimension(9,getNlocal()) :: A
94 <    real (kind=dp), dimension(3,getNlocal()) :: f
95 <    real (kind=dp), dimension(3,getNlocal()) :: t
96 <    logical, intent(in) :: do_pot, do_stress
92 >    real (kind=dp), dimension(3), intent(inout) :: fpair
93 >    real (kind=dp) :: pot, vpair, sw
94 >    real (kind=dp), dimension(9,nLocal) :: A
95 >    real (kind=dp), dimension(3,nLocal) :: f
96 >    real (kind=dp), dimension(3,nLocal) :: t
97 >    logical, intent(in) :: do_pot
98  
99      real (kind=dp) :: xi, yi, zi, xj, yj, zj, xi2, yi2, zi2, xj2, yj2, zj2
100      real (kind=dp) :: r3, r5, r6, s, sp, dsdr, dspdr
# Line 96 | Line 112 | contains
112      real (kind=dp) :: rijtest, rjitest
113      real (kind=dp) :: radcomxi, radcomyi, radcomzi
114      real (kind=dp) :: radcomxj, radcomyj, radcomzj
115 +    integer :: id1, id2
116  
100 –
117      if (.not.sticky_initialized) then
118         write(*,*) 'Sticky forces not initialized!'
119         return
120      endif
121  
106 –    if ( rij .LE. SSD_rup ) then
122  
123 +    if ( rij .LE. SSD_rbig ) then
124 +
125         r3 = r2*rij
126         r5 = r3*r2
127  
# Line 164 | Line 181 | contains
181         wjp = zjf*zjf*zjs*zjs - SSD_w0
182         wp = wip + wjp
183  
184 +       vpair = vpair + 0.5d0*(SSD_v0*s*w + SSD_v0p*sp*wp)
185         if (do_pot) then
186   #ifdef IS_MPI
187 <          pot_row(atom1) = pot_row(atom1) + 0.25d0*SSD_v0*(s*w + sp*wp)
188 <          pot_col(atom2) = pot_col(atom2) + 0.25d0*SSD_v0*(s*w + sp*wp)
187 >          pot_row(atom1) = pot_row(atom1) + 0.25d0*(SSD_v0*s*w + SSD_v0p*sp*wp)*sw
188 >          pot_col(atom2) = pot_col(atom2) + 0.25d0*(SSD_v0*s*w + SSD_v0p*sp*wp)*sw
189   #else
190 <          pot = pot + 0.5d0*SSD_v0*(s*w + sp*wp)
190 >          pot = pot + 0.5d0*(SSD_v0*s*w + SSD_v0p*sp*wp)*sw
191   #endif  
192         endif
193  
# Line 211 | Line 229 | contains
229         ! do the torques first since they are easy:
230         ! remember that these are still in the body fixed axes
231  
232 <       txi = 0.5d0*SSD_v0*(s*dwidux + sp*dwipdux)
233 <       tyi = 0.5d0*SSD_v0*(s*dwiduy + sp*dwipduy)
234 <       tzi = 0.5d0*SSD_v0*(s*dwiduz + sp*dwipduz)
232 >       txi = 0.5d0*(SSD_v0*s*dwidux + SSD_v0p*sp*dwipdux)*sw
233 >       tyi = 0.5d0*(SSD_v0*s*dwiduy + SSD_v0p*sp*dwipduy)*sw
234 >       tzi = 0.5d0*(SSD_v0*s*dwiduz + SSD_v0p*sp*dwipduz)*sw
235  
236 <       txj = 0.5d0*SSD_v0*(s*dwjdux + sp*dwjpdux)
237 <       tyj = 0.5d0*SSD_v0*(s*dwjduy + sp*dwjpduy)
238 <       tzj = 0.5d0*SSD_v0*(s*dwjduz + sp*dwjpduz)
236 >       txj = 0.5d0*(SSD_v0*s*dwjdux + SSD_v0p*sp*dwjpdux)*sw
237 >       tyj = 0.5d0*(SSD_v0*s*dwjduy + SSD_v0p*sp*dwjpduy)*sw
238 >       tzj = 0.5d0*(SSD_v0*s*dwjduz + SSD_v0p*sp*dwjpduz)*sw
239  
240         ! go back to lab frame using transpose of rotation matrix:
241  
# Line 248 | Line 266 | contains
266  
267         ! first rotate the i terms back into the lab frame:
268  
269 <       radcomxi = s*dwidx+sp*dwipdx
270 <       radcomyi = s*dwidy+sp*dwipdy
271 <       radcomzi = s*dwidz+sp*dwipdz
269 >       radcomxi = (SSD_v0*s*dwidx+SSD_v0p*sp*dwipdx)*sw
270 >       radcomyi = (SSD_v0*s*dwidy+SSD_v0p*sp*dwipdy)*sw
271 >       radcomzi = (SSD_v0*s*dwidz+SSD_v0p*sp*dwipdz)*sw
272  
273 <       radcomxj = s*dwjdx+sp*dwjpdx
274 <       radcomyj = s*dwjdy+sp*dwjpdy
275 <       radcomzj = s*dwjdz+sp*dwjpdz
273 >       radcomxj = (SSD_v0*s*dwjdx+SSD_v0p*sp*dwjpdx)*sw
274 >       radcomyj = (SSD_v0*s*dwjdy+SSD_v0p*sp*dwjpdy)*sw
275 >       radcomzj = (SSD_v0*s*dwjdz+SSD_v0p*sp*dwjpdz)*sw
276  
277   #ifdef IS_MPI    
278         fxii = a_Row(1,atom1)*(radcomxi) + &
# Line 308 | Line 326 | contains
326  
327         ! now assemble these with the radial-only terms:
328  
329 <       fxradial = 0.5d0*SSD_v0*(dsdr*drdx*w + dspdr*drdx*wp + fxii + fxji)
330 <       fyradial = 0.5d0*SSD_v0*(dsdr*drdy*w + dspdr*drdy*wp + fyii + fyji)
331 <       fzradial = 0.5d0*SSD_v0*(dsdr*drdz*w + dspdr*drdz*wp + fzii + fzji)
329 >       fxradial = 0.5d0*(SSD_v0*dsdr*drdx*w + SSD_v0p*dspdr*drdx*wp + fxii + fxji)
330 >       fyradial = 0.5d0*(SSD_v0*dsdr*drdy*w + SSD_v0p*dspdr*drdy*wp + fyii + fyji)
331 >       fzradial = 0.5d0*(SSD_v0*dsdr*drdz*w + SSD_v0p*dspdr*drdz*wp + fzii + fzji)
332  
333   #ifdef IS_MPI
334         f_Row(1,atom1) = f_Row(1,atom1) + fxradial
# Line 330 | Line 348 | contains
348         f(3,atom2) = f(3,atom2) - fzradial
349   #endif
350  
351 <       if (do_stress) then          
352 <          if (molMembershipList(atom1) .ne. molMembershipList(atom2)) then
353 <
354 <             ! because the d vector is the rj - ri vector, and
355 <             ! because fxradial, fyradial, and fzradial are the
356 <             ! (positive) force on atom i (negative on atom j) we need
357 <             ! a negative sign here:
358 <
359 <             tau_Temp(1) = tau_Temp(1) - d(1) * fxradial
360 <             tau_Temp(2) = tau_Temp(2) - d(1) * fyradial
361 <             tau_Temp(3) = tau_Temp(3) - d(1) * fzradial
362 <             tau_Temp(4) = tau_Temp(4) - d(2) * fxradial
363 <             tau_Temp(5) = tau_Temp(5) - d(2) * fyradial
364 <             tau_Temp(6) = tau_Temp(6) - d(2) * fzradial
347 <             tau_Temp(7) = tau_Temp(7) - d(3) * fxradial
348 <             tau_Temp(8) = tau_Temp(8) - d(3) * fyradial
349 <             tau_Temp(9) = tau_Temp(9) - d(3) * fzradial
350 <
351 <             virial_Temp = virial_Temp + (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
352 <          endif
351 > #ifdef IS_MPI
352 >       id1 = tagRow(atom1)
353 >       id2 = tagColumn(atom2)
354 > #else
355 >       id1 = atom1
356 >       id2 = atom2
357 > #endif
358 >      
359 >       if (molMembershipList(id1) .ne. molMembershipList(id2)) then
360 >          
361 >          fpair(1) = fpair(1) + fxradial
362 >          fpair(2) = fpair(2) + fyradial
363 >          fpair(3) = fpair(3) + fzradial
364 >          
365         endif
366      endif
355 –
367    end subroutine do_sticky_pair
368  
369    !! calculates the switching functions and their derivatives for a given
370    subroutine calc_sw_fnc(r, s, sp, dsdr, dspdr)
371 <          
371 >    
372      real (kind=dp), intent(in) :: r
373      real (kind=dp), intent(inout) :: s, sp, dsdr, dspdr
374 <
374 >    
375      ! distances must be in angstroms
376      
377      if (r.lt.SSD_rl) then
378         s = 1.0d0
368 –       sp = 1.0d0
379         dsdr = 0.0d0
380 +    elseif (r.gt.SSD_ru) then
381 +       s = 0.0d0
382 +       dsdr = 0.0d0
383 +    else
384 +       s = ((SSD_ru + 2.0d0*r - 3.0d0*SSD_rl) * (SSD_ru-r)**2) / &
385 +            ((SSD_ru - SSD_rl)**3)
386 +       dsdr = 6.0d0*(r-SSD_ru)*(r-SSD_rl)/((SSD_ru - SSD_rl)**3)
387 +    endif
388 +
389 +    if (r.lt.SSD_rlp) then
390 +       sp = 1.0d0      
391         dspdr = 0.0d0
392      elseif (r.gt.SSD_rup) then
372 –       s = 0.0d0
393         sp = 0.0d0
374 –       dsdr = 0.0d0
394         dspdr = 0.0d0
395      else
396 <       sp = ((SSD_rup + 2.0d0*r - 3.0d0*SSD_rl) * (SSD_rup-r)**2) / &
397 <            ((SSD_rup - SSD_rl)**3)
398 <       dspdr = 6.0d0*(r-SSD_rup)*(r-SSD_rl)/((SSD_rup - SSD_rl)**3)
380 <      
381 <       if (r.gt.SSD_ru) then
382 <          s = 0.0d0
383 <          dsdr = 0.0d0
384 <       else
385 <          s = ((SSD_ru + 2.0d0*r - 3.0d0*SSD_rl) * (SSD_ru-r)**2) / &
386 <               ((SSD_ru - SSD_rl)**3)
387 <          dsdr = 6.0d0*(r-SSD_ru)*(r-SSD_rl)/((SSD_ru - SSD_rl)**3)
388 <       endif
396 >       sp = ((SSD_rup + 2.0d0*r - 3.0d0*SSD_rlp) * (SSD_rup-r)**2) / &
397 >            ((SSD_rup - SSD_rlp)**3)
398 >       dspdr = 6.0d0*(r-SSD_rup)*(r-SSD_rlp)/((SSD_rup - SSD_rlp)**3)      
399      endif
400      
401      return

Diff Legend

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