[ViewVC] Diff of: group/trunk/OOPSE-3.0/src/UseTheForce/DarkSide/sticky.F90

Comparing trunk/OOPSE-3.0/src/UseTheForce/DarkSide/sticky.F90 (file contents):
Revision 2204 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 2229 by chrisfen, Tue May 17 22:35:01 2005 UTC

+!! @author Matthew Meineke
+!! @author Christopher Fennell
+!! @author J. Daniel Gezelter
-<
+!! @version $Id: sticky.F90,v 1.7 2005-04-15 22:03:49 gezelter Exp $, $Date: 2005-04-15 22:03:49 $, $Name: not supported by cvs2svn $, $Revision: 1.7 $
->
+!! @version $Id: sticky.F90,v 1.10 2005-05-17 22:35:01 chrisfen Exp $, $Date: 2005-05-17 22:35:01 $, $Name: not supported by cvs2svn $, $Revision: 1.10 $
+module sticky
+  public :: newStickyType
+  public :: do_sticky_pair
+  public :: destroyStickyTypes
-+
+  public :: do_sticky_power_pair
+  type :: StickyList
+  subroutine destroyStickyTypes()
+    if(allocated(StickyMap)) deallocate(StickyMap)
+  end subroutine destroyStickyTypes
-+
-+
+    subroutine do_sticky_power_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
-+
+       pot, A, f, t, do_pot)
-+
+    !! We assume that the rotation matrices have already been calculated
-+
+    !! and placed in the A array.
-+
-+
+    !! i and j are pointers to the two SSD atoms
-+
-+
+    integer, intent(in) :: atom1, atom2
-+
+    real (kind=dp), intent(inout) :: rij, r2
-+
+    real (kind=dp), dimension(3), intent(in) :: d
-+
+    real (kind=dp), dimension(3), intent(inout) :: fpair
-+
+    real (kind=dp) :: pot, vpair, sw
-+
+    real (kind=dp), dimension(9,nLocal) :: A
-+
+    real (kind=dp), dimension(3,nLocal) :: f
-+
+    real (kind=dp), dimension(3,nLocal) :: t
-+
+    logical, intent(in) :: do_pot
-+
-+
+    real (kind=dp) :: xi, yi, zi, xj, yj, zj, xi2, yi2, zi2, xj2, yj2, zj2
-+
+    real (kind=dp) :: xihat, yihat, zihat, xjhat, yjhat, zjhat
-+
+    real (kind=dp) :: rI, rI2, rI3, rI4, rI5, rI6, rI7, s, sp, dsdr, dspdr
-+
+    real (kind=dp) :: wi, wj, w, wip, wjp, wp, wi2, wj2, wip3, wjp3
-+
+    real (kind=dp) :: dwidx, dwidy, dwidz, dwjdx, dwjdy, dwjdz
-+
+    real (kind=dp) :: dwipdx, dwipdy, dwipdz, dwjpdx, dwjpdy, dwjpdz
-+
+    real (kind=dp) :: dwidux, dwiduy, dwiduz, dwjdux, dwjduy, dwjduz
-+
+    real (kind=dp) :: dwipdux, dwipduy, dwipduz, dwjpdux, dwjpduy, dwjpduz
-+
+    real (kind=dp) :: zif, zis, zjf, zjs, uglyi, uglyj
-+
+    real (kind=dp) :: drdx, drdy, drdz
-+
+    real (kind=dp) :: txi, tyi, tzi, txj, tyj, tzj
-+
+    real (kind=dp) :: fxii, fyii, fzii, fxjj, fyjj, fzjj
-+
+    real (kind=dp) :: fxij, fyij, fzij, fxji, fyji, fzji
-+
+    real (kind=dp) :: fxradial, fyradial, fzradial
-+
+    real (kind=dp) :: rijtest, rjitest
-+
+    real (kind=dp) :: radcomxi, radcomyi, radcomzi
-+
+    real (kind=dp) :: radcomxj, radcomyj, radcomzj
-+
+    integer :: id1, id2
-+
+    integer :: me1, me2
-+
+    real (kind=dp) :: w0, v0, v0p, rl, ru, rlp, rup, rbig
-+
+    real (kind=dp) :: zi3, zi4, zi5, zj3, zj4, zj5
-+
+    real (kind=dp) :: frac1, frac2, prodVal
-+
+    real (kind=dp) :: prei1, prei2, prei, prej1, prej2, prej
-+
+    real (kind=dp) :: walt, walti, waltj, dwaltidx, dwaltidy, dwaltidz
-+
+    real (kind=dp) :: dwaltjdx, dwaltjdy, dwaltjdz
-+
+    real (kind=dp) :: dwaltidux, dwaltiduy, dwaltiduz
-+
+    real (kind=dp) :: dwaltjdux, dwaltjduy, dwaltjduz
-+
+    real (kind=dp) :: doSw1idx, doSw1idy, doSw1idz, doSw1jdx, doSw1jdy, doSw1jdz
-+
+    real (kind=dp) :: doSw1idux, doSw1iduy, doSw1iduz
-+
+    real (kind=dp) :: doSw1jdux, doSw1jduy, doSw1jduz
-+
+    real (kind=dp) :: doSw2idx, doSw2idy, doSw2idz, doSw2jdx, doSw2jdy, doSw2jdz
-+
+    real (kind=dp) :: doSw2idux, doSw2iduy, doSw2iduz
-+
+    real (kind=dp) :: doSw2jdux, doSw2jduy, doSw2jduz
-+
-+
+    if (.not.allocated(StickyMap)) then
-+
+       call handleError("sticky", "no StickyMap was present before first call of do_sticky_power_pair!")
-+
+       return
-+
+    end if
-+
-+
+#ifdef IS_MPI
-+
+    me1 = atid_Row(atom1)
-+
+    me2 = atid_Col(atom2)
-+
+#else
-+
+    me1 = atid(atom1)
-+
+    me2 = atid(atom2)
-+
+#endif
-+
-+
+    if (me1.eq.me2) then
-+
+       w0  = StickyMap(me1)%w0
-+
+       v0  = StickyMap(me1)%v0
-+
+       v0p = StickyMap(me1)%v0p
-+
+       rl  = StickyMap(me1)%rl
-+
+       ru  = StickyMap(me1)%ru
-+
+       rlp = StickyMap(me1)%rlp
-+
+       rup = StickyMap(me1)%rup
-+
+       rbig = StickyMap(me1)%rbig
-+
+    else
-+
+       ! This is silly, but if you want 2 sticky types in your
-+
+       ! simulation, we'll let you do it with the Lorentz-
-+
+       ! Berthelot mixing rules.
-+
+       ! (Warning: you'll be SLLLLLLLLLLLLLLLOOOOOOOOOOWWWWWWWWWWW)
-+
+       rl   = 0.5_dp * ( StickyMap(me1)%rl + StickyMap(me2)%rl )
-+
+       ru   = 0.5_dp * ( StickyMap(me1)%ru + StickyMap(me2)%ru )
-+
+       rlp  = 0.5_dp * ( StickyMap(me1)%rlp + StickyMap(me2)%rlp )
-+
+       rup  = 0.5_dp * ( StickyMap(me1)%rup + StickyMap(me2)%rup )
-+
+       rbig = max(ru, rup)
-+
+       w0  = sqrt( StickyMap(me1)%w0   * StickyMap(me2)%w0  )
-+
+       v0  = sqrt( StickyMap(me1)%v0   * StickyMap(me2)%v0  )
-+
+       v0p = sqrt( StickyMap(me1)%v0p  * StickyMap(me2)%v0p )
-+
+    endif
-+
-+
+    if ( rij .LE. rbig ) then
-+
-+
+       rI = 1.0d0/rij
-+
+       rI2 = rI*rI
-+
+       rI3 = rI2*rI
-+
+       rI4 = rI2*rI2
-+
+       rI5 = rI3*rI2
-+
+       rI6 = rI3*rI3
-+
+       rI7 = rI4*rI3
-+
-+
+       drdx = d(1) * rI
-+
+       drdy = d(2) * rI
-+
+       drdz = d(3) * rI
-+
-+
+#ifdef IS_MPI
-+
+       ! rotate the inter-particle separation into the two different
-+
+       ! body-fixed coordinate systems:
-+
-+
+       xi = A_row(1,atom1)*d(1) + A_row(2,atom1)*d(2) + A_row(3,atom1)*d(3)
-+
+       yi = A_row(4,atom1)*d(1) + A_row(5,atom1)*d(2) + A_row(6,atom1)*d(3)
-+
+       zi = A_row(7,atom1)*d(1) + A_row(8,atom1)*d(2) + A_row(9,atom1)*d(3)
-+
-+
+       ! negative sign because this is the vector from j to i:
-+
-+
+       xj = -(A_Col(1,atom2)*d(1) + A_Col(2,atom2)*d(2) + A_Col(3,atom2)*d(3))
-+
+       yj = -(A_Col(4,atom2)*d(1) + A_Col(5,atom2)*d(2) + A_Col(6,atom2)*d(3))
-+
+       zj = -(A_Col(7,atom2)*d(1) + A_Col(8,atom2)*d(2) + A_Col(9,atom2)*d(3))
-+
+#else
-+
+       ! rotate the inter-particle separation into the two different
-+
+       ! body-fixed coordinate systems:
-+
-+
+       xi = a(1,atom1)*d(1) + a(2,atom1)*d(2) + a(3,atom1)*d(3)
-+
+       yi = a(4,atom1)*d(1) + a(5,atom1)*d(2) + a(6,atom1)*d(3)
-+
+       zi = a(7,atom1)*d(1) + a(8,atom1)*d(2) + a(9,atom1)*d(3)
-+
-+
+       ! negative sign because this is the vector from j to i:
-+
-+
+       xj = -(a(1,atom2)*d(1) + a(2,atom2)*d(2) + a(3,atom2)*d(3))
-+
+       yj = -(a(4,atom2)*d(1) + a(5,atom2)*d(2) + a(6,atom2)*d(3))
-+
+       zj = -(a(7,atom2)*d(1) + a(8,atom2)*d(2) + a(9,atom2)*d(3))
-+
+#endif
-+
-+
+       xi2 = xi*xi
-+
+       yi2 = yi*yi
-+
+       zi2 = zi*zi
-+
+       zi3 = zi2*zi
-+
+       zi4 = zi2*zi2
-+
+       zi5 = zi4*zi
-+
+       xihat = xi*rI
-+
+       yihat = yi*rI
-+
+       zihat = zi*rI
-+
-+
+       xj2 = xj*xj
-+
+       yj2 = yj*yj
-+
+       zj2 = zj*zj
-+
+       zj3 = zj2*zj
-+
+       zj4 = zj2*zj2
-+
+       zj5 = zj4*zj
-+
+       xjhat = xj*rI
-+
+       yjhat = yj*rI
-+
+       zjhat = zj*rI
-+
-+
+       call calc_sw_fnc(rij, rl, ru, rlp, rup, s, sp, dsdr, dspdr)
-+
-+
+       frac1 = 0.5d0
-+
+       frac2 = 0.5d0
-+
-+
+       wi = 2.0d0*(xi2-yi2)*zi*rI3
-+
+       wj = 2.0d0*(xj2-yj2)*zj*rI3
-+
-+
+!       prodVal = zihat*zjhat
-+
+!       if (prodVal .ge. 0.0d0) then
-+
+!         wi = 0.0d0
-+
+!         wj = 0.0d0
-+
+!       endif
-+
-+
+       wi2 = wi*wi
-+
+       wj2 = wj*wj
-+
-+
+       w = frac1*wi*wi2 + frac2*wi + wj*wj2
-+
-+
+       zif = zihat - 0.6d0
-+
+       zis = zihat + 0.8d0
-+
-+
+       zjf = zjhat - 0.6d0
-+
+       zjs = zjhat + 0.8d0
-+
-+
+       wip = zif*zif*zis*zis - w0
-+
+       wjp = zjf*zjf*zjs*zjs - w0
-+
+       wp = wip + wjp
-+
-+
+       !wip = zihat - 0.2d0
-+
+       !wjp = zjhat - 0.2d0
-+
+       !wip3 = wip*wip*wip
-+
+       !wjp3 = wjp*wjp*wjp
-+
-+
+       !wp = wip3*wip + wjp3*wjp
-+
-+
+       vpair = vpair + 0.5d0*(v0*s*w + v0p*sp*wp)
-+
-+
+       if (do_pot) then
-+
+#ifdef IS_MPI
-+
+         pot_row(atom1) = pot_row(atom1) + 0.25d0*(v0*s*w + v0p*sp*wp)*sw
-+
+         pot_col(atom2) = pot_col(atom2) + 0.25d0*(v0*s*w + v0p*sp*wp)*sw
-+
+#else
-+
+         pot = pot + 0.5d0*(v0*s*w + v0p*sp*wp)*sw
-+
+#endif
-+
+       endif
-+
-+
+       dwidx = ( 4.0d0*xi*zi*rI3 - 6.0d0*xi*zi*(xi2-yi2)*rI5 )
-+
+       dwidy = ( -4.0d0*yi*zi*rI3 - 6.0d0*yi*zi*(xi2-yi2)*rI5 )
-+
+       dwidz = ( 2.0d0*(xi2-yi2)*rI3 - 6.0d0*zi2*(xi2-yi2)*rI5 )
-+
-+
+       dwidx = frac1*3.0d0*wi2*dwidx + frac2*dwidx
-+
+       dwidy = frac1*3.0d0*wi2*dwidy + frac2*dwidy
-+
+       dwidz = frac1*3.0d0*wi2*dwidz + frac2*dwidz
-+
-+
+       dwjdx = ( 4.0d0*xj*zj*rI3  - 6.0d0*xj*zj*(xj2-yj2)*rI5 )
-+
+       dwjdy = ( -4.0d0*yj*zj*rI3  - 6.0d0*yj*zj*(xj2-yj2)*rI5 )
-+
+       dwjdz = ( 2.0d0*(xj2-yj2)*rI3  - 6.0d0*zj2*(xj2-yj2)*rI5 )
-+
-+
+       dwjdx = frac1*3.0d0*wj2*dwjdx + frac2*dwjdx
-+
+       dwjdy = frac1*3.0d0*wj2*dwjdy + frac2*dwjdy
-+
+       dwjdz = frac1*3.0d0*wj2*dwjdz + frac2*dwjdz
-+
-+
+       uglyi = zif*zif*zis + zif*zis*zis
-+
+       uglyj = zjf*zjf*zjs + zjf*zjs*zjs
-+
-+
+       dwipdx = -2.0d0*xi*zi*uglyi*rI3
-+
+       dwipdy = -2.0d0*yi*zi*uglyi*rI3
-+
+       dwipdz = 2.0d0*(rI - zi2*rI3)*uglyi
-+
-+
+       dwjpdx = -2.0d0*xj*zj*uglyj*rI3
-+
+       dwjpdy = -2.0d0*yj*zj*uglyj*rI3
-+
+       dwjpdz = 2.0d0*(rI - zj2*rI3)*uglyj
-+
-+
+       !dwipdx = -4.0d0*wip3*zi*xihat
-+
+       !dwipdy = -4.0d0*wip3*zi*yihat
-+
+       !dwipdz = -4.0d0*wip3*(zi2 - 1.0d0)*rI
-+
-+
+       !dwjpdx = -4.0d0*wjp3*zj*xjhat
-+
+       !dwjpdy = -4.0d0*wjp3*zj*yjhat
-+
+       !dwjpdz = -4.0d0*wjp3*(zj2 - 1.0d0)*rI
-+
-+
+       !dwipdx = 0.0d0
-+
+       !dwipdy = 0.0d0
-+
+       !dwipdz = 0.0d0
-+
-+
+       !dwjpdx = 0.0d0
-+
+       !dwjpdy = 0.0d0
-+
+       !dwjpdz = 0.0d0
-+
-+
+       dwidux = ( 4.0d0*(yi*zi2 + 0.5d0*yi*(xi2-yi2))*rI3 )
-+
+       dwiduy = ( 4.0d0*(xi*zi2 - 0.5d0*xi*(xi2-yi2))*rI3 )
-+
+       dwiduz = ( -8.0d0*xi*yi*zi*rI3 )
-+
-+
+       dwidux = frac1*3.0d0*wi2*dwidux + frac2*dwidux
-+
+       dwiduy = frac1*3.0d0*wi2*dwiduy + frac2*dwiduy
-+
+       dwiduz = frac1*3.0d0*wi2*dwiduz + frac2*dwiduz
-+
-+
+       dwjdux = ( 4.0d0*(yj*zj2 + 0.5d0*yj*(xj2-yj2))*rI3 )
-+
+       dwjduy = ( 4.0d0*(xj*zj2 - 0.5d0*xj*(xj2-yj2))*rI3 )
-+
+       dwjduz = ( -8.0d0*xj*yj*zj*rI3 )
-+
-+
+       dwjdux = frac1*3.0d0*wj2*dwjdux + frac2*dwjdux
-+
+       dwjduy = frac1*3.0d0*wj2*dwjduy + frac2*dwjduy
-+
+       dwjduz = frac1*3.0d0*wj2*dwjduz + frac2*dwjduz
-+
-+
+       dwipdux =  2.0d0*yi*uglyi*rI
-+
+       dwipduy = -2.0d0*xi*uglyi*rI
-+
+       dwipduz =  0.0d0
-+
-+
+       dwjpdux =  2.0d0*yj*uglyj*rI
-+
+       dwjpduy = -2.0d0*xj*uglyj*rI
-+
+       dwjpduz =  0.0d0
-+
-+
+       !dwipdux =  4.0d0*wip3*yihat
-+
+       !dwipduy = -4.0d0*wip3*xihat
-+
+       !dwipduz =  0.0d0
-+
-+
+       !dwjpdux =  4.0d0*wjp3*yjhat
-+
+       !dwjpduy = -4.0d0*wjp3*xjhat
-+
+       !dwjpduz =  0.0d0
-+
-+
+       !dwipdux = 0.0d0
-+
+       !dwipduy = 0.0d0
-+
+       !dwipduz = 0.0d0
-+
-+
+       !dwjpdux = 0.0d0
-+
+       !dwjpduy = 0.0d0
-+
+       !dwjpduz = 0.0d0
-+
-+
+       ! do the torques first since they are easy:
-+
+       ! remember that these are still in the body fixed axes
-+
-+
+       txi = 0.5d0*(v0*s*dwidux + v0p*sp*dwipdux)*sw
-+
+       tyi = 0.5d0*(v0*s*dwiduy + v0p*sp*dwipduy)*sw
-+
+       tzi = 0.5d0*(v0*s*dwiduz + v0p*sp*dwipduz)*sw
-+
-+
+       txj = 0.5d0*(v0*s*dwjdux + v0p*sp*dwjpdux)*sw
-+
+       tyj = 0.5d0*(v0*s*dwjduy + v0p*sp*dwjpduy)*sw
-+
+       tzj = 0.5d0*(v0*s*dwjduz + v0p*sp*dwjpduz)*sw
-+
-+
+       ! go back to lab frame using transpose of rotation matrix:
-+
-+
+#ifdef IS_MPI
-+
+       t_Row(1,atom1) = t_Row(1,atom1) + a_Row(1,atom1)*txi + &
-+
+            a_Row(4,atom1)*tyi + a_Row(7,atom1)*tzi
-+
+       t_Row(2,atom1) = t_Row(2,atom1) + a_Row(2,atom1)*txi + &
-+
+            a_Row(5,atom1)*tyi + a_Row(8,atom1)*tzi
-+
+       t_Row(3,atom1) = t_Row(3,atom1) + a_Row(3,atom1)*txi + &
-+
+            a_Row(6,atom1)*tyi + a_Row(9,atom1)*tzi
-+
-+
+       t_Col(1,atom2) = t_Col(1,atom2) + a_Col(1,atom2)*txj + &
-+
+            a_Col(4,atom2)*tyj + a_Col(7,atom2)*tzj
-+
+       t_Col(2,atom2) = t_Col(2,atom2) + a_Col(2,atom2)*txj + &
-+
+            a_Col(5,atom2)*tyj + a_Col(8,atom2)*tzj
-+
+       t_Col(3,atom2) = t_Col(3,atom2) + a_Col(3,atom2)*txj + &
-+
+            a_Col(6,atom2)*tyj + a_Col(9,atom2)*tzj
-+
+#else
-+
+       t(1,atom1) = t(1,atom1) + a(1,atom1)*txi + a(4,atom1)*tyi + a(7,atom1)*tzi
-+
+       t(2,atom1) = t(2,atom1) + a(2,atom1)*txi + a(5,atom1)*tyi + a(8,atom1)*tzi
-+
+       t(3,atom1) = t(3,atom1) + a(3,atom1)*txi + a(6,atom1)*tyi + a(9,atom1)*tzi
-+
-+
+       t(1,atom2) = t(1,atom2) + a(1,atom2)*txj + a(4,atom2)*tyj + a(7,atom2)*tzj
-+
+       t(2,atom2) = t(2,atom2) + a(2,atom2)*txj + a(5,atom2)*tyj + a(8,atom2)*tzj
-+
+       t(3,atom2) = t(3,atom2) + a(3,atom2)*txj + a(6,atom2)*tyj + a(9,atom2)*tzj
-+
+#endif
-+
+       ! Now, on to the forces:
-+
-+
+       ! first rotate the i terms back into the lab frame:
-+
-+
+       radcomxi = (v0*s*dwidx+v0p*sp*dwipdx)*sw
-+
+       radcomyi = (v0*s*dwidy+v0p*sp*dwipdy)*sw
-+
+       radcomzi = (v0*s*dwidz+v0p*sp*dwipdz)*sw
-+
-+
+       radcomxj = (v0*s*dwjdx+v0p*sp*dwjpdx)*sw
-+
+       radcomyj = (v0*s*dwjdy+v0p*sp*dwjpdy)*sw
-+
+       radcomzj = (v0*s*dwjdz+v0p*sp*dwjpdz)*sw
-+
-+
+#ifdef IS_MPI
-+
+       fxii = a_Row(1,atom1)*(radcomxi) + &
-+
+            a_Row(4,atom1)*(radcomyi) + &
-+
+            a_Row(7,atom1)*(radcomzi)
-+
+       fyii = a_Row(2,atom1)*(radcomxi) + &
-+
+            a_Row(5,atom1)*(radcomyi) + &
-+
+            a_Row(8,atom1)*(radcomzi)
-+
+       fzii = a_Row(3,atom1)*(radcomxi) + &
-+
+            a_Row(6,atom1)*(radcomyi) + &
-+
+            a_Row(9,atom1)*(radcomzi)
-+
-+
+       fxjj = a_Col(1,atom2)*(radcomxj) + &
-+
+            a_Col(4,atom2)*(radcomyj) + &
-+
+            a_Col(7,atom2)*(radcomzj)
-+
+       fyjj = a_Col(2,atom2)*(radcomxj) + &
-+
+            a_Col(5,atom2)*(radcomyj) + &
-+
+            a_Col(8,atom2)*(radcomzj)
-+
+       fzjj = a_Col(3,atom2)*(radcomxj)+ &
-+
+            a_Col(6,atom2)*(radcomyj) + &
-+
+            a_Col(9,atom2)*(radcomzj)
-+
+#else
-+
+       fxii = a(1,atom1)*(radcomxi) + &
-+
+            a(4,atom1)*(radcomyi) + &
-+
+            a(7,atom1)*(radcomzi)
-+
+       fyii = a(2,atom1)*(radcomxi) + &
-+
+            a(5,atom1)*(radcomyi) + &
-+
+            a(8,atom1)*(radcomzi)
-+
+       fzii = a(3,atom1)*(radcomxi) + &
-+
+            a(6,atom1)*(radcomyi) + &
-+
+            a(9,atom1)*(radcomzi)
-+
-+
+       fxjj = a(1,atom2)*(radcomxj) + &
-+
+            a(4,atom2)*(radcomyj) + &
-+
+            a(7,atom2)*(radcomzj)
-+
+       fyjj = a(2,atom2)*(radcomxj) + &
-+
+            a(5,atom2)*(radcomyj) + &
-+
+            a(8,atom2)*(radcomzj)
-+
+       fzjj = a(3,atom2)*(radcomxj)+ &
-+
+            a(6,atom2)*(radcomyj) + &
-+
+            a(9,atom2)*(radcomzj)
-+
+#endif
-+
-+
+       fxij = -fxii
-+
+       fyij = -fyii
-+
+       fzij = -fzii
-+
-+
+       fxji = -fxjj
-+
+       fyji = -fyjj
-+
+       fzji = -fzjj
-+
-+
+       ! now assemble these with the radial-only terms:
-+
-+
+       fxradial = 0.5d0*((v0*dsdr*w + v0p*dspdr*wp)*drdx + fxii + fxji)
-+
+       fyradial = 0.5d0*((v0*dsdr*w + v0p*dspdr*wp)*drdy + fyii + fyji)
-+
+       fzradial = 0.5d0*((v0*dsdr*w + v0p*dspdr*wp)*drdz + fzii + fzji)
-+
-+
+#ifdef IS_MPI
-+
+       f_Row(1,atom1) = f_Row(1,atom1) + fxradial
-+
+       f_Row(2,atom1) = f_Row(2,atom1) + fyradial
-+
+       f_Row(3,atom1) = f_Row(3,atom1) + fzradial
-+
-+
+       f_Col(1,atom2) = f_Col(1,atom2) - fxradial
-+
+       f_Col(2,atom2) = f_Col(2,atom2) - fyradial
-+
+       f_Col(3,atom2) = f_Col(3,atom2) - fzradial
-+
+#else
-+
+       f(1,atom1) = f(1,atom1) + fxradial
-+
+       f(2,atom1) = f(2,atom1) + fyradial
-+
+       f(3,atom1) = f(3,atom1) + fzradial
-+
-+
+       f(1,atom2) = f(1,atom2) - fxradial
-+
+       f(2,atom2) = f(2,atom2) - fyradial
-+
+       f(3,atom2) = f(3,atom2) - fzradial
-+
+#endif
-+
-+
+#ifdef IS_MPI
-+
+       id1 = AtomRowToGlobal(atom1)
-+
+       id2 = AtomColToGlobal(atom2)
-+
+#else
-+
+       id1 = atom1
-+
+       id2 = atom2
-+
+#endif
-+
-+
+       if (molMembershipList(id1) .ne. molMembershipList(id2)) then
-+
-+
+          fpair(1) = fpair(1) + fxradial
-+
+          fpair(2) = fpair(2) + fyradial
-+
+          fpair(3) = fpair(3) + fzradial
-+
-+
+       endif
-+
+    endif
-+
+  end subroutine do_sticky_power_pair
-+
+end module sticky

Diff Legend

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

Comparing trunk/OOPSE-3.0/src/UseTheForce/DarkSide/sticky.F90 (file contents): Revision 2204 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs. Revision 2229 by chrisfen, Tue May 17 22:35:01 2005 UTC

Diff Legend

Comparing trunk/OOPSE-3.0/src/UseTheForce/DarkSide/sticky.F90 (file contents):
Revision 2204 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 2229 by chrisfen, Tue May 17 22:35:01 2005 UTC