--- trunk/OOPSE-3.0/src/UseTheForce/DarkSide/sticky.F90 2005/05/03 22:55:29 2219 +++ trunk/OOPSE-3.0/src/UseTheForce/DarkSide/sticky.F90 2005/05/05 14:47:35 2220 @@ -50,7 +50,7 @@ !! @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 @@ -70,6 +70,7 @@ module sticky public :: newStickyType public :: do_sticky_pair public :: destroyStickyTypes + public :: do_sticky_power_pair type :: StickyList @@ -510,4 +511,347 @@ end module sticky 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