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 635 by gezelter, Thu Jul 17 20:32:24 2003 UTC vs.
Revision 1160 by gezelter, Tue May 11 21:31:15 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.11 2003-07-17 20:32:24 gezelter Exp $, $Date: 2003-07-17 20:32:24 $, $Name: not supported by cvs2svn $, $Revision: 1.11 $
12 > !! @version $Id: calc_sticky_pair.F90,v 1.18 2004-05-11 21:31:14 gezelter Exp $, $Date: 2004-05-11 21:31:14 $, $Name: not supported by cvs2svn $, $Revision: 1.18 $
13  
14   module sticky_pair
15  
# Line 32 | Line 32 | module sticky_pair
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 62 | Line 63 | contains
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, &
77 >  subroutine do_sticky_pair(atom1, atom2, d, rij, r2, sw, vpair, pot, A,f, t, &
78         do_pot, do_stress)
79 <
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 82 | 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
92 >    real (kind=dp) :: pot, vpair, sw
93 >    real (kind=dp), dimension(9,nLocal) :: A
94 >    real (kind=dp), dimension(3,nLocal) :: f
95 >    real (kind=dp), dimension(3,nLocal) :: t
96      logical, intent(in) :: do_pot, do_stress
97  
98      real (kind=dp) :: xi, yi, zi, xj, yj, zj, xi2, yi2, zi2, xj2, yj2, zj2
# Line 104 | Line 111 | contains
111      real (kind=dp) :: rijtest, rjitest
112      real (kind=dp) :: radcomxi, radcomyi, radcomzi
113      real (kind=dp) :: radcomxj, radcomyj, radcomzj
114 +    integer :: id1, id2
115  
108 –
116      if (.not.sticky_initialized) then
117         write(*,*) 'Sticky forces not initialized!'
118         return
119      endif
120  
114 –    if ( rij .LE. SSD_rup ) then
121  
122 +    if ( rij .LE. SSD_rbig ) then
123 +
124         r3 = r2*rij
125         r5 = r3*r2
126  
# Line 172 | Line 180 | contains
180         wjp = zjf*zjf*zjs*zjs - SSD_w0
181         wp = wip + wjp
182  
183 +       vpair = vpair + 0.5d0*(SSD_v0*s*w + SSD_v0p*sp*wp)
184         if (do_pot) then
185   #ifdef IS_MPI
186 <          pot_row(atom1) = pot_row(atom1) + 0.25d0*(SSD_v0*s*w + SSD_v0p*sp*wp)
187 <          pot_col(atom2) = pot_col(atom2) + 0.25d0*(SSD_v0*s*w + SSD_v0p*sp*wp)
186 >          pot_row(atom1) = pot_row(atom1) + 0.25d0*(SSD_v0*s*w + SSD_v0p*sp*wp)*sw
187 >          pot_col(atom2) = pot_col(atom2) + 0.25d0*(SSD_v0*s*w + SSD_v0p*sp*wp)*sw
188   #else
189 <          pot = pot + 0.5d0*(SSD_v0*s*w + SSD_v0p*sp*wp)
189 >          pot = pot + 0.5d0*(SSD_v0*s*w + SSD_v0p*sp*wp)*sw
190   #endif  
191         endif
192  
# Line 219 | Line 228 | contains
228         ! do the torques first since they are easy:
229         ! remember that these are still in the body fixed axes
230  
231 <       txi = 0.5d0*(SSD_v0*s*dwidux + SSD_v0p*sp*dwipdux)
232 <       tyi = 0.5d0*(SSD_v0*s*dwiduy + SSD_v0p*sp*dwipduy)
233 <       tzi = 0.5d0*(SSD_v0*s*dwiduz + SSD_v0p*sp*dwipduz)
231 >       txi = 0.5d0*(SSD_v0*s*dwidux + SSD_v0p*sp*dwipdux)*sw
232 >       tyi = 0.5d0*(SSD_v0*s*dwiduy + SSD_v0p*sp*dwipduy)*sw
233 >       tzi = 0.5d0*(SSD_v0*s*dwiduz + SSD_v0p*sp*dwipduz)*sw
234  
235 <       txj = 0.5d0*(SSD_v0*s*dwjdux + SSD_v0p*sp*dwjpdux)
236 <       tyj = 0.5d0*(SSD_v0*s*dwjduy + SSD_v0p*sp*dwjpduy)
237 <       tzj = 0.5d0*(SSD_v0*s*dwjduz + SSD_v0p*sp*dwjpduz)
235 >       txj = 0.5d0*(SSD_v0*s*dwjdux + SSD_v0p*sp*dwjpdux)*sw
236 >       tyj = 0.5d0*(SSD_v0*s*dwjduy + SSD_v0p*sp*dwjpduy)*sw
237 >       tzj = 0.5d0*(SSD_v0*s*dwjduz + SSD_v0p*sp*dwjpduz)*sw
238  
239         ! go back to lab frame using transpose of rotation matrix:
240  
# Line 256 | Line 265 | contains
265  
266         ! first rotate the i terms back into the lab frame:
267  
268 <       radcomxi = SSD_v0*s*dwidx+SSD_v0p*sp*dwipdx
269 <       radcomyi = SSD_v0*s*dwidy+SSD_v0p*sp*dwipdy
270 <       radcomzi = SSD_v0*s*dwidz+SSD_v0p*sp*dwipdz
268 >       radcomxi = (SSD_v0*s*dwidx+SSD_v0p*sp*dwipdx)*sw
269 >       radcomyi = (SSD_v0*s*dwidy+SSD_v0p*sp*dwipdy)*sw
270 >       radcomzi = (SSD_v0*s*dwidz+SSD_v0p*sp*dwipdz)*sw
271  
272 <       radcomxj = SSD_v0*s*dwjdx+SSD_v0p*sp*dwjpdx
273 <       radcomyj = SSD_v0*s*dwjdy+SSD_v0p*sp*dwjpdy
274 <       radcomzj = SSD_v0*s*dwjdz+SSD_v0p*sp*dwjpdz
272 >       radcomxj = (SSD_v0*s*dwjdx+SSD_v0p*sp*dwjpdx)*sw
273 >       radcomyj = (SSD_v0*s*dwjdy+SSD_v0p*sp*dwjpdy)*sw
274 >       radcomzj = (SSD_v0*s*dwjdz+SSD_v0p*sp*dwjpdz)*sw
275  
276   #ifdef IS_MPI    
277         fxii = a_Row(1,atom1)*(radcomxi) + &
# Line 339 | Line 348 | contains
348   #endif
349  
350         if (do_stress) then          
342 –          if (molMembershipList(atom1) .ne. molMembershipList(atom2)) then
351  
352 + #ifdef IS_MPI
353 +          id1 = tagRow(atom1)
354 +          id2 = tagColumn(atom2)
355 + #else
356 +          id1 = atom1
357 +          id2 = atom2
358 + #endif
359 +
360 +          if (molMembershipList(id1) .ne. molMembershipList(id2)) then
361 +
362               ! because the d vector is the rj - ri vector, and
363               ! because fxradial, fyradial, and fzradial are the
364               ! (positive) force on atom i (negative on atom j) we need

Diff Legend

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