--- trunk/OOPSE/libmdtools/calc_sticky_pair.F90 2003/07/17 19:38:23 634 +++ trunk/OOPSE/libmdtools/calc_sticky_pair.F90 2003/07/17 20:32:24 635 @@ -9,7 +9,7 @@ !! @author Matthew Meineke !! @author Christopher Fennel !! @author J. Daniel Gezelter -!! @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 $ +!! @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 $ module sticky_pair @@ -27,8 +27,10 @@ module sticky_pair logical, save :: sticky_initialized = .false. real( kind = dp ), save :: SSD_w0 = 0.0_dp real( kind = dp ), save :: SSD_v0 = 0.0_dp + real( kind = dp ), save :: SSD_v0p = 0.0_dp real( kind = dp ), save :: SSD_rl = 0.0_dp real( kind = dp ), save :: SSD_ru = 0.0_dp + real( kind = dp ), save :: SSD_rlp = 0.0_dp real( kind = dp ), save :: SSD_rup = 0.0_dp public :: check_sticky_FF @@ -44,16 +46,22 @@ contains return end subroutine check_sticky_FF - subroutine set_sticky_params(sticky_w0, sticky_v0) - real( kind = dp ), intent(in) :: sticky_w0, sticky_v0 + subroutine set_sticky_params(sticky_w0, sticky_v0, sticky_v0p, & + sticky_rl, sticky_ru, sticky_rlp, sticky_rup) + + real( kind = dp ), intent(in) :: sticky_w0, sticky_v0, sticky_v0p + real( kind = dp ), intent(in) :: sticky_rl, sticky_ru + real( kind = dp ), intent(in) :: sticky_rlp, sticky_rup ! we could pass all 5 parameters if we felt like it... SSD_w0 = sticky_w0 SSD_v0 = sticky_v0 - SSD_rl = 2.75_DP - SSD_ru = 3.35_DP - SSD_rup = 4.0_DP + SSD_v0p = sticky_v0p + SSD_rl = sticky_rl + SSD_ru = sticky_ru + SSD_rlp = sticky_rlp + SSD_rup = sticky_rup sticky_initialized = .true. return @@ -166,10 +174,10 @@ contains if (do_pot) then #ifdef IS_MPI - pot_row(atom1) = pot_row(atom1) + 0.25d0*SSD_v0*(s*w + sp*wp) - pot_col(atom2) = pot_col(atom2) + 0.25d0*SSD_v0*(s*w + sp*wp) + pot_row(atom1) = pot_row(atom1) + 0.25d0*(SSD_v0*s*w + SSD_v0p*sp*wp) + pot_col(atom2) = pot_col(atom2) + 0.25d0*(SSD_v0*s*w + SSD_v0p*sp*wp) #else - pot = pot + 0.5d0*SSD_v0*(s*w + sp*wp) + pot = pot + 0.5d0*(SSD_v0*s*w + SSD_v0p*sp*wp) #endif endif @@ -211,13 +219,13 @@ contains ! do the torques first since they are easy: ! remember that these are still in the body fixed axes - txi = 0.5d0*SSD_v0*(s*dwidux + sp*dwipdux) - tyi = 0.5d0*SSD_v0*(s*dwiduy + sp*dwipduy) - tzi = 0.5d0*SSD_v0*(s*dwiduz + sp*dwipduz) + txi = 0.5d0*(SSD_v0*s*dwidux + SSD_v0p*sp*dwipdux) + tyi = 0.5d0*(SSD_v0*s*dwiduy + SSD_v0p*sp*dwipduy) + tzi = 0.5d0*(SSD_v0*s*dwiduz + SSD_v0p*sp*dwipduz) - txj = 0.5d0*SSD_v0*(s*dwjdux + sp*dwjpdux) - tyj = 0.5d0*SSD_v0*(s*dwjduy + sp*dwjpduy) - tzj = 0.5d0*SSD_v0*(s*dwjduz + sp*dwjpduz) + txj = 0.5d0*(SSD_v0*s*dwjdux + SSD_v0p*sp*dwjpdux) + tyj = 0.5d0*(SSD_v0*s*dwjduy + SSD_v0p*sp*dwjpduy) + tzj = 0.5d0*(SSD_v0*s*dwjduz + SSD_v0p*sp*dwjpduz) ! go back to lab frame using transpose of rotation matrix: @@ -248,13 +256,13 @@ contains ! first rotate the i terms back into the lab frame: - radcomxi = s*dwidx+sp*dwipdx - radcomyi = s*dwidy+sp*dwipdy - radcomzi = s*dwidz+sp*dwipdz + radcomxi = SSD_v0*s*dwidx+SSD_v0p*sp*dwipdx + radcomyi = SSD_v0*s*dwidy+SSD_v0p*sp*dwipdy + radcomzi = SSD_v0*s*dwidz+SSD_v0p*sp*dwipdz - radcomxj = s*dwjdx+sp*dwjpdx - radcomyj = s*dwjdy+sp*dwjpdy - radcomzj = s*dwjdz+sp*dwjpdz + radcomxj = SSD_v0*s*dwjdx+SSD_v0p*sp*dwjpdx + radcomyj = SSD_v0*s*dwjdy+SSD_v0p*sp*dwjpdy + radcomzj = SSD_v0*s*dwjdz+SSD_v0p*sp*dwjpdz #ifdef IS_MPI fxii = a_Row(1,atom1)*(radcomxi) + & @@ -308,9 +316,9 @@ contains ! now assemble these with the radial-only terms: - fxradial = 0.5d0*SSD_v0*(dsdr*drdx*w + dspdr*drdx*wp + fxii + fxji) - fyradial = 0.5d0*SSD_v0*(dsdr*drdy*w + dspdr*drdy*wp + fyii + fyji) - fzradial = 0.5d0*SSD_v0*(dsdr*drdz*w + dspdr*drdz*wp + fzii + fzji) + fxradial = 0.5d0*(SSD_v0*dsdr*drdx*w + SSD_v0p*dspdr*drdx*wp + fxii + fxji) + fyradial = 0.5d0*(SSD_v0*dsdr*drdy*w + SSD_v0p*dspdr*drdy*wp + fyii + fyji) + fzradial = 0.5d0*(SSD_v0*dsdr*drdz*w + SSD_v0p*dspdr*drdz*wp + fzii + fzji) #ifdef IS_MPI f_Row(1,atom1) = f_Row(1,atom1) + fxradial @@ -357,35 +365,34 @@ contains !! calculates the switching functions and their derivatives for a given subroutine calc_sw_fnc(r, s, sp, dsdr, dspdr) - + real (kind=dp), intent(in) :: r real (kind=dp), intent(inout) :: s, sp, dsdr, dspdr - + ! distances must be in angstroms if (r.lt.SSD_rl) then s = 1.0d0 - sp = 1.0d0 dsdr = 0.0d0 + elseif (r.gt.SSD_ru) then + s = 0.0d0 + dsdr = 0.0d0 + else + s = ((SSD_ru + 2.0d0*r - 3.0d0*SSD_rl) * (SSD_ru-r)**2) / & + ((SSD_ru - SSD_rl)**3) + dsdr = 6.0d0*(r-SSD_ru)*(r-SSD_rl)/((SSD_ru - SSD_rl)**3) + endif + + if (r.lt.SSD_rlp) then + sp = 1.0d0 dspdr = 0.0d0 elseif (r.gt.SSD_rup) then - s = 0.0d0 sp = 0.0d0 - dsdr = 0.0d0 dspdr = 0.0d0 else - sp = ((SSD_rup + 2.0d0*r - 3.0d0*SSD_rl) * (SSD_rup-r)**2) / & - ((SSD_rup - SSD_rl)**3) - dspdr = 6.0d0*(r-SSD_rup)*(r-SSD_rl)/((SSD_rup - SSD_rl)**3) - - if (r.gt.SSD_ru) then - s = 0.0d0 - dsdr = 0.0d0 - else - s = ((SSD_ru + 2.0d0*r - 3.0d0*SSD_rl) * (SSD_ru-r)**2) / & - ((SSD_ru - SSD_rl)**3) - dsdr = 6.0d0*(r-SSD_ru)*(r-SSD_rl)/((SSD_ru - SSD_rl)**3) - endif + sp = ((SSD_rup + 2.0d0*r - 3.0d0*SSD_rlp) * (SSD_rup-r)**2) / & + ((SSD_rup - SSD_rlp)**3) + dspdr = 6.0d0*(r-SSD_rup)*(r-SSD_rlp)/((SSD_rup - SSD_rlp)**3) endif return