[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 2219 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 2220 by chrisfen, Thu May 5 14:47:35 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.8 2005-05-05 14:47:35 chrisfen Exp $, $Date: 2005-05-05 14:47:35 $, $Name: not supported by cvs2svn $, $Revision: 1.8 $
+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) :: r3, r5, r6, s, sp, dsdr, dspdr
-+
+    real (kind=dp) :: wi, wj, w, wip, wjp, wp, wi2, wj2
-+
+    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
-+
-+
+    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
-+
-+
+       r3 = r2*rij
-+
+       r5 = r3*r2
-+
-+
+       drdx = d(1) / rij
-+
+       drdy = d(2) / rij
-+
+       drdz = d(3) / rij
-+
-+
+#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
-+
-+
+       xj2 = xj*xj
-+
+       yj2 = yj*yj
-+
+       zj2 = zj*zj
-+
-+
+       call calc_sw_fnc(rij, rl, ru, rlp, rup, s, sp, dsdr, dspdr)
-+
-+
+       wi = 2.0d0*(xi2-yi2)*zi / r3
-+
+       wj = 2.0d0*(xj2-yj2)*zj / r3
-+
+       !rootwi = sqrt(abs(wi))
-+
+       !rootwj = sqrt(abs(wj))
-+
+       wi2 = wi*wi
-+
+       wj2 = wj*wj
-+
-+
-+
+       w = wi*wi2+wj*wj2
-+
-+
+       zif = zi/rij - 0.6d0
-+
+       zis = zi/rij + 0.8d0
-+
-+
+       zjf = zj/rij - 0.6d0
-+
+       zjs = zj/rij + 0.8d0
-+
-+
+       wip = zif*zif*zis*zis - w0
-+
+       wjp = zjf*zjf*zjs*zjs - w0
-+
+       wp = wip + 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 = 1.5d0*rootwi*( 4.0d0*xi*zi/r3 - 6.0d0*xi*zi*(xi2-yi2)/r5 )
-+
+!       dwidy = 1.5d0*rootwi*( -4.0d0*yi*zi/r3 - 6.0d0*yi*zi*(xi2-yi2)/r5 )
-+
+!       dwidz = 1.5d0*rootwi*( 2.0d0*(xi2-yi2)/r3 - 6.0d0*zi2*(xi2-yi2)/r5 )
-+
-+
+!       dwjdx = 1.5d0*rootwj*( 4.0d0*xj*zj/r3  - 6.0d0*xj*zj*(xj2-yj2)/r5 )
-+
+!       dwjdy = 1.5d0*rootwj*( -4.0d0*yj*zj/r3  - 6.0d0*yj*zj*(xj2-yj2)/r5 )
-+
+!       dwjdz = 1.5d0*rootwj*( 2.0d0*(xj2-yj2)/r3  - 6.0d0*zj2*(xj2-yj2)/r5 )
-+
-+
+       dwidx = 3.0d0*wi2*( 4.0d0*xi*zi/r3 - 6.0d0*xi*zi*(xi2-yi2)/r5 )
-+
+       dwidy = 3.0d0*wi2*( -4.0d0*yi*zi/r3 - 6.0d0*yi*zi*(xi2-yi2)/r5 )
-+
+       dwidz = 3.0d0*wi2*( 2.0d0*(xi2-yi2)/r3 - 6.0d0*zi2*(xi2-yi2)/r5 )
-+
-+
+       dwjdx = 3.0d0*wj2*( 4.0d0*xj*zj/r3  - 6.0d0*xj*zj*(xj2-yj2)/r5 )
-+
+       dwjdy = 3.0d0*wj2*( -4.0d0*yj*zj/r3  - 6.0d0*yj*zj*(xj2-yj2)/r5 )
-+
+       dwjdz = 3.0d0*wj2*( 2.0d0*(xj2-yj2)/r3  - 6.0d0*zj2*(xj2-yj2)/r5 )
-+
-+
+       uglyi = zif*zif*zis + zif*zis*zis
-+
+       uglyj = zjf*zjf*zjs + zjf*zjs*zjs
-+
-+
+       dwipdx = -2.0d0*xi*zi*uglyi/r3
-+
+       dwipdy = -2.0d0*yi*zi*uglyi/r3
-+
+       dwipdz = 2.0d0*(1.0d0/rij - zi2/r3)*uglyi
-+
-+
+       dwjpdx = -2.0d0*xj*zj*uglyj/r3
-+
+       dwjpdy = -2.0d0*yj*zj*uglyj/r3
-+
+       dwjpdz = 2.0d0*(1.0d0/rij - zj2/r3)*uglyj
-+
-+
+!       dwidux = 1.5d0*rootwi*( 4.0d0*(yi*zi2 + 0.5d0*yi*(xi2-yi2))/r3 )
-+
+!       dwiduy = 1.5d0*rootwi*( 4.0d0*(xi*zi2 - 0.5d0*xi*(xi2-yi2))/r3 )
-+
+!       dwiduz = 1.5d0*rootwi*( -8.0d0*xi*yi*zi/r3 )
-+
-+
+!       dwjdux = 1.5d0*rootwj*( 4.0d0*(yj*zj2 + 0.5d0*yj*(xj2-yj2))/r3 )
-+
+!       dwjduy = 1.5d0*rootwj*( 4.0d0*(xj*zj2 - 0.5d0*xj*(xj2-yj2))/r3 )
-+
+!       dwjduz = 1.5d0*rootwj*( -8.0d0*xj*yj*zj/r3 )
-+
-+
+       dwidux = 3.0d0*wi2*( 4.0d0*(yi*zi2 + 0.5d0*yi*(xi2-yi2))/r3 )
-+
+       dwiduy = 3.0d0*wi2*( 4.0d0*(xi*zi2 - 0.5d0*xi*(xi2-yi2))/r3 )
-+
+       dwiduz = 3.0d0*wi2*( -8.0d0*xi*yi*zi/r3 )
-+
-+
+       dwjdux = 3.0d0*wj2*( 4.0d0*(yj*zj2 + 0.5d0*yj*(xj2-yj2))/r3 )
-+
+       dwjduy = 3.0d0*wj2*( 4.0d0*(xj*zj2 - 0.5d0*xj*(xj2-yj2))/r3 )
-+
+       dwjduz = 3.0d0*wj2*( -8.0d0*xj*yj*zj/r3 )
-+
-+
+       dwipdux =  2.0d0*yi*uglyi/rij
-+
+       dwipduy = -2.0d0*xi*uglyi/rij
-+
+       dwipduz =  0.0d0
-+
-+
+       dwjpdux =  2.0d0*yj*uglyj/rij
-+
+       dwjpduy = -2.0d0*xj*uglyj/rij
-+
+       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*drdx*w + v0p*dspdr*drdx*wp + fxii + fxji)
-+
+       fyradial = 0.5d0*(v0*dsdr*drdy*w + v0p*dspdr*drdy*wp + fyii + fyji)
-+
+       fzradial = 0.5d0*(v0*dsdr*drdz*w + v0p*dspdr*drdz*wp + 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 2219 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs. Revision 2220 by chrisfen, Thu May 5 14:47:35 2005 UTC

Diff Legend

Comparing trunk/OOPSE-3.0/src/UseTheForce/DarkSide/sticky.F90 (file contents):
Revision 2219 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 2220 by chrisfen, Thu May 5 14:47:35 2005 UTC