mdtools/libmdCode/calc_sticky_pair.F90

!! This Module Calculates forces due to SSD potential and VDW interactions
!! [Chandra and Ichiye, J. Chem. Phys. 111, 2701 (1999)].

!! This module contains the Public procedures:


!! Corresponds to the force field defined in ssd_FF.cpp 
!! @author Charles F. Vardeman II
!! @author Matthew Meineke
!! @author Christopher Fennel
!! @author J. Daniel Gezelter
!! @version $Id: calc_sticky_pair.F90,v 1.5 2003-03-12 19:31:55 gezelter Exp $, $Date: 2003-03-12 19:31:55 $, $Name: not supported by cvs2svn $, $Revision: 1.5 $

module sticky_pair

  use simulation
  use definitions
  use forceGlobals
#ifdef IS_MPI
  use mpiSimulation
#endif

  implicit none

  PRIVATE

  logical, save :: sticky_initialized = .false.
  real( kind = dp ), save :: SSD_w0 
  real( kind = dp ), save :: SSD_v0
  real( kind = dp ), save :: SSD_rl
  real( kind = dp ), save :: SSD_ru
  real( kind = dp ), save :: SSD_rup

  public :: check_sticky_FF
  public :: set_sticky_params
  public :: do_sticky_pair

contains

  subroutine check_sticky_FF(status)
    integer :: status 
    status = -1
    if (sticky_initialized) status = 0
    return
  end subroutine check_sticky_FF

  subroutine set_sticky_params(sticky_w0, sticky_v0)
    real( kind = dp ), intent(in) :: sticky_w0, sticky_v0
    
    ! 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
   
    sticky_initialized = .true.
    return
  end subroutine set_sticky_params

  subroutine do_sticky_pair(atom1, atom2, d, rij, r2, A, pot, f, t, &
       do_pot, do_stress)
    
    !! This routine does only the sticky portion of the SSD potential
    !! [Chandra and Ichiye, J. Chem. Phys. 111, 2701 (1999)].
    !! The Lennard-Jones and dipolar interaction must be handled separately.
    
    !! 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) :: pot
    real (kind=dp), dimension(9,getNlocal()) :: A
    real (kind=dp), dimension(3,getNlocal()) :: f
    real (kind=dp), dimension(3,getNlocal()) :: t
    logical, intent(in) :: do_pot, do_stress

    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
    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
       
    if (.not.sticky_initialized) then
       write(*,*) 'Sticky forces not initialized!'
       return
    endif

    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, s, sp, dsdr, dspdr)
    
    wi = 2.0d0*(xi2-yi2)*zi / r3
    wj = 2.0d0*(xj2-yj2)*zj / r3
    w = wi+wj
    
    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 - SSD_w0
    wjp = zjf*zjf*zjs*zjs - SSD_w0
    wp = wip + wjp


    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)
#else
       pot = pot + 0.5d0*SSD_v0*(s*w + sp*wp)
#endif  
    endif
    
    dwidx =   4.0d0*xi*zi/r3  - 6.0d0*xi*zi*(xi2-yi2)/r5
    dwidy = - 4.0d0*yi*zi/r3  - 6.0d0*yi*zi*(xi2-yi2)/r5
    dwidz =   2.0d0*(xi2-yi2)/r3  - 6.0d0*zi2*(xi2-yi2)/r5
    
    dwjdx =   4.0d0*xj*zj/r3  - 6.0d0*xj*zj*(xj2-yj2)/r5
    dwjdy = - 4.0d0*yj*zj/r3  - 6.0d0*yj*zj*(xj2-yj2)/r5
    dwjdz =   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 = 4.0d0*(yi*zi2 + 0.5d0*yi*(xi2-yi2))/r3
    dwiduy = 4.0d0*(xi*zi2 - 0.5d0*xi*(xi2-yi2))/r3
    dwiduz = - 8.0d0*xi*yi*zi/r3
    
    dwjdux = 4.0d0*(yj*zj2 + 0.5d0*yj*(xj2-yj2))/r3
    dwjduy = 4.0d0*(xj*zj2 - 0.5d0*xj*(xj2-yj2))/r3
    dwjduz = - 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*SSD_v0*(s*dwidux + sp*dwipdux)
    tyi = 0.5d0*SSD_v0*(s*dwiduy + sp*dwipduy)
    tzi = 0.5d0*SSD_v0*(s*dwiduz + 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)
    
    ! 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:

#ifdef IS_MPI    
    fxii = a_Row(1,atom1)*(s*dwidx+sp*dwipdx) + &
         a_Row(4,atom1)*(s*dwidy+sp*dwipdy) + &
         a_Row(7,atom1)*(s*dwidz+sp*dwipdz)
    fyii = a_Row(2,atom1)*(s*dwidx+sp*dwipdx) + &
         a_Row(5,atom1)*(s*dwidy+sp*dwipdy) + &
         a_Row(8,atom1)*(s*dwidz+sp*dwipdz)
    fzii = a_Row(3,atom1)*(s*dwidx+sp*dwipdx) + &
         a_Row(6,atom1)*(s*dwidy+sp*dwipdy) + &
         a_Row(9,atom1)*(s*dwidz+sp*dwipdz)

    fxjj = a_Col(1,atom2)*(s*dwjdx+sp*dwjpdx) + &
         a_Col(4,atom2)*(s*dwjdy+sp*dwjpdy) + &
         a_Col(7,atom2)*(s*dwjdz+sp*dwjpdz)
    fyjj = a_Col(2,atom2)*(s*dwjdx+sp*dwjpdx) + &
         a_Col(5,atom2)*(s*dwjdy+sp*dwjpdy) + &
         a_Col(8,atom2)*(s*dwjdz+sp*dwjpdz)
    fzjj = a_Col(3,atom2)*(s*dwjdx+sp*dwjpdx)+ &
         a_Col(6,atom2)*(s*dwjdy+sp*dwjpdy) + &
         a_Col(9,atom2)*(s*dwjdz+sp*dwjpdz)
#else
    fxii = a(1,atom1)*(s*dwidx+sp*dwipdx) + &
         a(4,atom1)*(s*dwidy+sp*dwipdy) + &
         a(7,atom1)*(s*dwidz+sp*dwipdz)
    fyii = a(2,atom1)*(s*dwidx+sp*dwipdx) + &
         a(5,atom1)*(s*dwidy+sp*dwipdy) + &
         a(8,atom1)*(s*dwidz+sp*dwipdz)
    fzii = a(3,atom1)*(s*dwidx+sp*dwipdx) + &
         a(6,atom1)*(s*dwidy+sp*dwipdy) + &
         a(9,atom1)*(s*dwidz+sp*dwipdz)

    fxjj = a(1,atom2)*(s*dwjdx+sp*dwjpdx) + &
         a(4,atom2)*(s*dwjdy+sp*dwjpdy) + &
         a(7,atom2)*(s*dwjdz+sp*dwjpdz)
    fyjj = a(2,atom2)*(s*dwjdx+sp*dwjpdx) + &
         a(5,atom2)*(s*dwjdy+sp*dwjpdy) + &
         a(8,atom2)*(s*dwjdz+sp*dwjpdz)
    fzjj = a(3,atom2)*(s*dwjdx+sp*dwjpdx)+ &
         a(6,atom2)*(s*dwjdy+sp*dwjpdy) + &
         a(9,atom2)*(s*dwjdz+sp*dwjpdz)
#endif
    
    fxij = -fxii
    fyij = -fyii
    fzij = -fzii
        
    fxji = -fxjj
    fyji = -fyjj
    fzji = -fzjj
    
    ! 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)
       
#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) + 0.5d0*SSD_v0*(-dsdr*drdx*w - &
         dspdr*drdx*wp + fxjj + fxij)
    f_Col(2,atom2) = f_Col(2,atom2) + 0.5d0*SSD_v0*(-dsdr*drdy*w - &
         dspdr*drdy*wp + fyjj + fyij)
    f_Col(3,atom2) = f_Col(3,atom2) + 0.5d0*SSD_v0*(-dsdr*drdz*w - &
         dspdr*drdz*wp + fzjj + fzij)
#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) + 0.5d0*SSD_v0*(-dsdr*drdx*w - dspdr*drdx*wp + &
         fxjj + fxij)
    f(2,atom2) = f(2,atom2) + 0.5d0*SSD_v0*(-dsdr*drdy*w - dspdr*drdy*wp + &
         fyjj + fyij)
    f(3,atom2) = f(3,atom2) + 0.5d0*SSD_v0*(-dsdr*drdz*w - dspdr*drdz*wp + &
         fzjj + fzij)
#endif
    
    if (do_stress) then          
       tau_Temp(1) = tau_Temp(1) + fxradial * d(1)
       tau_Temp(2) = tau_Temp(2) + fxradial * d(2)
       tau_Temp(3) = tau_Temp(3) + fxradial * d(3)
       tau_Temp(4) = tau_Temp(4) + fyradial * d(1)
       tau_Temp(5) = tau_Temp(5) + fyradial * d(2)
       tau_Temp(6) = tau_Temp(6) + fyradial * d(3)
       tau_Temp(7) = tau_Temp(7) + fzradial * d(1)
       tau_Temp(8) = tau_Temp(8) + fzradial * d(2)
       tau_Temp(9) = tau_Temp(9) + fzradial * d(3)
       virial_Temp = virial_Temp + (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
    endif
   
  end subroutine do_sticky_pair

  !! 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
       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
    endif
    
    return
  end subroutine calc_sw_fnc
end module sticky_pair
Revision:	326
Committed:	Wed Mar 12 19:31:55 2003 UTC (21 years, 5 months ago) by gezelter
File size:	12766 byte(s)
Log Message:	fixes for major rewrite
#	User	Rev	Content
1	gezelter	319	!! This Module Calculates forces due to SSD potential and VDW interactions
2			!! [Chandra and Ichiye, J. Chem. Phys. 111, 2701 (1999)].
3
4			!! This module contains the Public procedures:
5
6
7			!! Corresponds to the force field defined in ssd_FF.cpp
8			!! @author Charles F. Vardeman II
9			!! @author Matthew Meineke
10	gezelter	322	!! @author Christopher Fennel
11			!! @author J. Daniel Gezelter
12	gezelter	326	!! @version $Id: calc_sticky_pair.F90,v 1.5 2003-03-12 19:31:55 gezelter Exp $, $Date: 2003-03-12 19:31:55 $, $Name: not supported by cvs2svn $, $Revision: 1.5 $
13	gezelter	319
14	gezelter	322	module sticky_pair
15	gezelter	319
16			use simulation
17			use definitions
18	gezelter	322	use forceGlobals
19	gezelter	319	#ifdef IS_MPI
20			use mpiSimulation
21			#endif
22	gezelter	322
23	gezelter	319	implicit none
24	gezelter	322
25	gezelter	319	PRIVATE
26
27	gezelter	322	logical, save :: sticky_initialized = .false.
28	gezelter	323	real( kind = dp ), save :: SSD_w0
29	gezelter	322	real( kind = dp ), save :: SSD_v0
30			real( kind = dp ), save :: SSD_rl
31			real( kind = dp ), save :: SSD_ru
32			real( kind = dp ), save :: SSD_rup
33	gezelter	319
34	gezelter	323	public :: check_sticky_FF
35			public :: set_sticky_params
36	gezelter	322	public :: do_sticky_pair
37	gezelter	319
38			contains
39
40	gezelter	323	subroutine check_sticky_FF(status)
41	gezelter	324	integer :: status
42			status = -1
43	gezelter	323	if (sticky_initialized) status = 0
44			return
45			end subroutine check_sticky_FF
46
47			subroutine set_sticky_params(sticky_w0, sticky_v0)
48	gezelter	322	real( kind = dp ), intent(in) :: sticky_w0, sticky_v0
49
50			! we could pass all 5 parameters if we felt like it...
51	gezelter	323
52	gezelter	322	SSD_w0 = sticky_w0
53			SSD_v0 = sticky_v0
54			SSD_rl = 2.75_DP
55			SSD_ru = 3.35_DP
56			SSD_rup = 4.0_DP
57	gezelter	323
58	gezelter	322	sticky_initialized = .true.
59			return
60	gezelter	323	end subroutine set_sticky_params
61	gezelter	319
62	gezelter	326	subroutine do_sticky_pair(atom1, atom2, d, rij, r2, A, pot, f, t, &
63			do_pot, do_stress)
64	gezelter	322
65			!! This routine does only the sticky portion of the SSD potential
66			!! [Chandra and Ichiye, J. Chem. Phys. 111, 2701 (1999)].
67			!! The Lennard-Jones and dipolar interaction must be handled separately.
68
69			!! We assume that the rotation matrices have already been calculated
70			!! and placed in the A array.
71
72			!! i and j are pointers to the two SSD atoms
73	gezelter	319
74	gezelter	322	integer, intent(in) :: atom1, atom2
75	gezelter	323	real (kind=dp), intent(inout) :: rij, r2
76	gezelter	322	real (kind=dp), dimension(3), intent(in) :: d
77	gezelter	323	real (kind=dp) :: pot
78	gezelter	322	real (kind=dp), dimension(9,getNlocal()) :: A
79			real (kind=dp), dimension(3,getNlocal()) :: f
80			real (kind=dp), dimension(3,getNlocal()) :: t
81	gezelter	326	logical, intent(in) :: do_pot, do_stress
82	gezelter	319
83	gezelter	322	real (kind=dp) :: xi, yi, zi, xj, yj, zj, xi2, yi2, zi2, xj2, yj2, zj2
84	gezelter	323	real (kind=dp) :: r3, r5, r6, s, sp, dsdr, dspdr
85	gezelter	322	real (kind=dp) :: wi, wj, w, wip, wjp, wp
86			real (kind=dp) :: dwidx, dwidy, dwidz, dwjdx, dwjdy, dwjdz
87			real (kind=dp) :: dwipdx, dwipdy, dwipdz, dwjpdx, dwjpdy, dwjpdz
88			real (kind=dp) :: dwidux, dwiduy, dwiduz, dwjdux, dwjduy, dwjduz
89			real (kind=dp) :: dwipdux, dwipduy, dwipduz, dwjpdux, dwjpduy, dwjpduz
90			real (kind=dp) :: zif, zis, zjf, zjs, uglyi, uglyj
91			real (kind=dp) :: drdx, drdy, drdz
92			real (kind=dp) :: txi, tyi, tzi, txj, tyj, tzj
93			real (kind=dp) :: fxii, fyii, fzii, fxjj, fyjj, fzjj
94			real (kind=dp) :: fxij, fyij, fzij, fxji, fyji, fzji
95			real (kind=dp) :: fxradial, fyradial, fzradial
96
97			if (.not.sticky_initialized) then
98			write(,) 'Sticky forces not initialized!'
99	gezelter	319	return
100			endif
101
102	gezelter	322	r3 = r2*rij
103			r5 = r3*r2
104
105			drdx = d(1) / rij
106			drdy = d(2) / rij
107			drdz = d(3) / rij
108
109	gezelter	319	#ifdef IS_MPI
110	gezelter	322	! rotate the inter-particle separation into the two different
111			! body-fixed coordinate systems:
112	gezelter	319
113	gezelter	322	xi = A_row(1,atom1)d(1) + A_row(2,atom1)d(2) + A_row(3,atom1)*d(3)
114			yi = A_row(4,atom1)d(1) + A_row(5,atom1)d(2) + A_row(6,atom1)*d(3)
115			zi = A_row(7,atom1)d(1) + A_row(8,atom1)d(2) + A_row(9,atom1)*d(3)
116	gezelter	319
117	gezelter	322	! negative sign because this is the vector from j to i:
118
119			xj = -(A_Col(1,atom2)d(1) + A_Col(2,atom2)d(2) + A_Col(3,atom2)*d(3))
120			yj = -(A_Col(4,atom2)d(1) + A_Col(5,atom2)d(2) + A_Col(6,atom2)*d(3))
121			zj = -(A_Col(7,atom2)d(1) + A_Col(8,atom2)d(2) + A_Col(9,atom2)*d(3))
122	gezelter	319	#else
123	gezelter	322	! rotate the inter-particle separation into the two different
124			! body-fixed coordinate systems:
125	gezelter	319
126	gezelter	322	xi = a(1,atom1)d(1) + a(2,atom1)d(2) + a(3,atom1)*d(3)
127			yi = a(4,atom1)d(1) + a(5,atom1)d(2) + a(6,atom1)*d(3)
128			zi = a(7,atom1)d(1) + a(8,atom1)d(2) + a(9,atom1)*d(3)
129
130			! negative sign because this is the vector from j to i:
131
132			xj = -(a(1,atom2)d(1) + a(2,atom2)d(2) + a(3,atom2)*d(3))
133			yj = -(a(4,atom2)d(1) + a(5,atom2)d(2) + a(6,atom2)*d(3))
134			zj = -(a(7,atom2)d(1) + a(8,atom2)d(2) + a(9,atom2)*d(3))
135	gezelter	319	#endif
136
137	gezelter	322	xi2 = xi*xi
138			yi2 = yi*yi
139			zi2 = zi*zi
140	gezelter	319
141	gezelter	322	xj2 = xj*xj
142			yj2 = yj*yj
143			zj2 = zj*zj
144	gezelter	319
145	gezelter	322	call calc_sw_fnc(rij, s, sp, dsdr, dspdr)
146
147			wi = 2.0d0(xi2-yi2)zi / r3
148			wj = 2.0d0(xj2-yj2)zj / r3
149			w = wi+wj
150
151			zif = zi/rij - 0.6d0
152			zis = zi/rij + 0.8d0
153
154			zjf = zj/rij - 0.6d0
155			zjs = zj/rij + 0.8d0
156
157			wip = zifzifzis*zis - SSD_w0
158			wjp = zjfzjfzjs*zjs - SSD_w0
159			wp = wip + wjp
160	gezelter	319
161	gezelter	326
162			if (do_pot) then
163	gezelter	322	#ifdef IS_MPI
164	gezelter	326	pot_row(atom1) = pot_row(atom1) + 0.25d0SSD_v0(sw + spwp)
165			pot_col(atom2) = pot_col(atom2) + 0.25d0SSD_v0(sw + spwp)
166	gezelter	322	#else
167	gezelter	326	pot = pot + 0.5d0SSD_v0(sw + spwp)
168			#endif
169			endif
170	gezelter	319
171	gezelter	322	dwidx = 4.0d0xizi/r3 - 6.0d0xizi*(xi2-yi2)/r5
172			dwidy = - 4.0d0yizi/r3 - 6.0d0yizi*(xi2-yi2)/r5
173			dwidz = 2.0d0(xi2-yi2)/r3 - 6.0d0zi2*(xi2-yi2)/r5
174
175			dwjdx = 4.0d0xjzj/r3 - 6.0d0xjzj*(xj2-yj2)/r5
176			dwjdy = - 4.0d0yjzj/r3 - 6.0d0yjzj*(xj2-yj2)/r5
177			dwjdz = 2.0d0(xj2-yj2)/r3 - 6.0d0zj2*(xj2-yj2)/r5
178
179			uglyi = zifzifzis + zifziszis
180			uglyj = zjfzjfzjs + zjfzjszjs
181
182			dwipdx = -2.0d0xizi*uglyi/r3
183			dwipdy = -2.0d0yizi*uglyi/r3
184			dwipdz = 2.0d0(1.0d0/rij - zi2/r3)uglyi
185
186			dwjpdx = -2.0d0xjzj*uglyj/r3
187			dwjpdy = -2.0d0yjzj*uglyj/r3
188			dwjpdz = 2.0d0(1.0d0/rij - zj2/r3)uglyj
189
190			dwidux = 4.0d0(yizi2 + 0.5d0yi(xi2-yi2))/r3
191			dwiduy = 4.0d0(xizi2 - 0.5d0xi(xi2-yi2))/r3
192			dwiduz = - 8.0d0xiyi*zi/r3
193
194			dwjdux = 4.0d0(yjzj2 + 0.5d0yj(xj2-yj2))/r3
195			dwjduy = 4.0d0(xjzj2 - 0.5d0xj(xj2-yj2))/r3
196			dwjduz = - 8.0d0xjyj*zj/r3
197
198			dwipdux = 2.0d0yiuglyi/rij
199			dwipduy = -2.0d0xiuglyi/rij
200			dwipduz = 0.0d0
201
202			dwjpdux = 2.0d0yjuglyj/rij
203			dwjpduy = -2.0d0xjuglyj/rij
204			dwjpduz = 0.0d0
205
206			! do the torques first since they are easy:
207			! remember that these are still in the body fixed axes
208
209			txi = 0.5d0SSD_v0(sdwidux + spdwipdux)
210			tyi = 0.5d0SSD_v0(sdwiduy + spdwipduy)
211			tzi = 0.5d0SSD_v0(sdwiduz + spdwipduz)
212
213			txj = 0.5d0SSD_v0(sdwjdux + spdwjpdux)
214			tyj = 0.5d0SSD_v0(sdwjduy + spdwjpduy)
215			tzj = 0.5d0SSD_v0(sdwjduz + spdwjpduz)
216
217			! go back to lab frame using transpose of rotation matrix:
218	gezelter	319
219			#ifdef IS_MPI
220	gezelter	322	t_Row(1,atom1) = t_Row(1,atom1) + a_Row(1,atom1)*txi + &
221			a_Row(4,atom1)tyi + a_Row(7,atom1)tzi
222			t_Row(2,atom1) = t_Row(2,atom1) + a_Row(2,atom1)*txi + &
223			a_Row(5,atom1)tyi + a_Row(8,atom1)tzi
224			t_Row(3,atom1) = t_Row(3,atom1) + a_Row(3,atom1)*txi + &
225			a_Row(6,atom1)tyi + a_Row(9,atom1)tzi
226	gezelter	319
227	gezelter	322	t_Col(1,atom2) = t_Col(1,atom2) + a_Col(1,atom2)*txj + &
228			a_Col(4,atom2)tyj + a_Col(7,atom2)tzj
229			t_Col(2,atom2) = t_Col(2,atom2) + a_Col(2,atom2)*txj + &
230			a_Col(5,atom2)tyj + a_Col(8,atom2)tzj
231			t_Col(3,atom2) = t_Col(3,atom2) + a_Col(3,atom2)*txj + &
232			a_Col(6,atom2)tyj + a_Col(9,atom2)tzj
233	gezelter	319	#else
234	gezelter	322	t(1,atom1) = t(1,atom1) + a(1,atom1)txi + a(4,atom1)tyi + a(7,atom1)*tzi
235			t(2,atom1) = t(2,atom1) + a(2,atom1)txi + a(5,atom1)tyi + a(8,atom1)*tzi
236			t(3,atom1) = t(3,atom1) + a(3,atom1)txi + a(6,atom1)tyi + a(9,atom1)*tzi
237
238			t(1,atom2) = t(1,atom2) + a(1,atom2)txj + a(4,atom2)tyj + a(7,atom2)*tzj
239			t(2,atom2) = t(2,atom2) + a(2,atom2)txj + a(5,atom2)tyj + a(8,atom2)*tzj
240			t(3,atom2) = t(3,atom2) + a(3,atom2)txj + a(6,atom2)tyj + a(9,atom2)*tzj
241			#endif
242			! Now, on to the forces:
243
244			! first rotate the i terms back into the lab frame:
245	gezelter	319
246	gezelter	322	#ifdef IS_MPI
247			fxii = a_Row(1,atom1)(sdwidx+sp*dwipdx) + &
248			a_Row(4,atom1)(sdwidy+sp*dwipdy) + &
249			a_Row(7,atom1)(sdwidz+sp*dwipdz)
250			fyii = a_Row(2,atom1)(sdwidx+sp*dwipdx) + &
251			a_Row(5,atom1)(sdwidy+sp*dwipdy) + &
252			a_Row(8,atom1)(sdwidz+sp*dwipdz)
253			fzii = a_Row(3,atom1)(sdwidx+sp*dwipdx) + &
254			a_Row(6,atom1)(sdwidy+sp*dwipdy) + &
255			a_Row(9,atom1)(sdwidz+sp*dwipdz)
256	gezelter	319
257	gezelter	322	fxjj = a_Col(1,atom2)(sdwjdx+sp*dwjpdx) + &
258			a_Col(4,atom2)(sdwjdy+sp*dwjpdy) + &
259			a_Col(7,atom2)(sdwjdz+sp*dwjpdz)
260			fyjj = a_Col(2,atom2)(sdwjdx+sp*dwjpdx) + &
261			a_Col(5,atom2)(sdwjdy+sp*dwjpdy) + &
262			a_Col(8,atom2)(sdwjdz+sp*dwjpdz)
263			fzjj = a_Col(3,atom2)(sdwjdx+sp*dwjpdx)+ &
264			a_Col(6,atom2)(sdwjdy+sp*dwjpdy) + &
265			a_Col(9,atom2)(sdwjdz+sp*dwjpdz)
266	gezelter	319	#else
267	gezelter	322	fxii = a(1,atom1)(sdwidx+sp*dwipdx) + &
268			a(4,atom1)(sdwidy+sp*dwipdy) + &
269			a(7,atom1)(sdwidz+sp*dwipdz)
270			fyii = a(2,atom1)(sdwidx+sp*dwipdx) + &
271			a(5,atom1)(sdwidy+sp*dwipdy) + &
272			a(8,atom1)(sdwidz+sp*dwipdz)
273			fzii = a(3,atom1)(sdwidx+sp*dwipdx) + &
274			a(6,atom1)(sdwidy+sp*dwipdy) + &
275			a(9,atom1)(sdwidz+sp*dwipdz)
276	gezelter	319
277	gezelter	322	fxjj = a(1,atom2)(sdwjdx+sp*dwjpdx) + &
278			a(4,atom2)(sdwjdy+sp*dwjpdy) + &
279			a(7,atom2)(sdwjdz+sp*dwjpdz)
280			fyjj = a(2,atom2)(sdwjdx+sp*dwjpdx) + &
281			a(5,atom2)(sdwjdy+sp*dwjpdy) + &
282			a(8,atom2)(sdwjdz+sp*dwjpdz)
283			fzjj = a(3,atom2)(sdwjdx+sp*dwjpdx)+ &
284			a(6,atom2)(sdwjdy+sp*dwjpdy) + &
285			a(9,atom2)(sdwjdz+sp*dwjpdz)
286	gezelter	319	#endif
287	gezelter	322
288			fxij = -fxii
289			fyij = -fyii
290			fzij = -fzii
291
292			fxji = -fxjj
293			fyji = -fyjj
294			fzji = -fzjj
295
296			! now assemble these with the radial-only terms:
297	gezelter	319
298	gezelter	322	fxradial = 0.5d0SSD_v0(dsdrdrdxw + dspdrdrdxwp + fxii + fxji)
299			fyradial = 0.5d0SSD_v0(dsdrdrdyw + dspdrdrdywp + fyii + fyji)
300			fzradial = 0.5d0SSD_v0(dsdrdrdzw + dspdrdrdzwp + fzii + fzji)
301	gezelter	319
302			#ifdef IS_MPI
303	gezelter	322	f_Row(1,atom1) = f_Row(1,atom1) + fxradial
304			f_Row(2,atom1) = f_Row(2,atom1) + fyradial
305			f_Row(3,atom1) = f_Row(3,atom1) + fzradial
306
307			f_Col(1,atom2) = f_Col(1,atom2) + 0.5d0SSD_v0(-dsdrdrdxw - &
308			dspdrdrdxwp + fxjj + fxij)
309			f_Col(2,atom2) = f_Col(2,atom2) + 0.5d0SSD_v0(-dsdrdrdyw - &
310			dspdrdrdywp + fyjj + fyij)
311			f_Col(3,atom2) = f_Col(3,atom2) + 0.5d0SSD_v0(-dsdrdrdzw - &
312			dspdrdrdzwp + fzjj + fzij)
313	gezelter	319	#else
314	gezelter	322	f(1,atom1) = f(1,atom1) + fxradial
315			f(2,atom1) = f(2,atom1) + fyradial
316			f(3,atom1) = f(3,atom1) + fzradial
317
318			f(1,atom2) = f(1,atom2) + 0.5d0SSD_v0(-dsdrdrdxw - dspdrdrdxwp + &
319			fxjj + fxij)
320			f(2,atom2) = f(2,atom2) + 0.5d0SSD_v0(-dsdrdrdyw - dspdrdrdywp + &
321			fyjj + fyij)
322			f(3,atom2) = f(3,atom2) + 0.5d0SSD_v0(-dsdrdrdzw - dspdrdrdzwp + &
323			fzjj + fzij)
324	gezelter	319	#endif
325	gezelter	322
326	gezelter	326	if (do_stress) then
327	gezelter	322	tau_Temp(1) = tau_Temp(1) + fxradial * d(1)
328			tau_Temp(2) = tau_Temp(2) + fxradial * d(2)
329			tau_Temp(3) = tau_Temp(3) + fxradial * d(3)
330			tau_Temp(4) = tau_Temp(4) + fyradial * d(1)
331			tau_Temp(5) = tau_Temp(5) + fyradial * d(2)
332			tau_Temp(6) = tau_Temp(6) + fyradial * d(3)
333			tau_Temp(7) = tau_Temp(7) + fzradial * d(1)
334			tau_Temp(8) = tau_Temp(8) + fzradial * d(2)
335			tau_Temp(9) = tau_Temp(9) + fzradial * d(3)
336			virial_Temp = virial_Temp + (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
337			endif
338
339			end subroutine do_sticky_pair
340	gezelter	319
341	gezelter	322	!! calculates the switching functions and their derivatives for a given
342			subroutine calc_sw_fnc(r, s, sp, dsdr, dspdr)
343	gezelter	319
344	gezelter	322	real (kind=dp), intent(in) :: r
345			real (kind=dp), intent(inout) :: s, sp, dsdr, dspdr
346	gezelter	319
347	gezelter	322	! distances must be in angstroms
348	gezelter	319
349	gezelter	322	if (r.lt.SSD_rl) then
350			s = 1.0d0
351			sp = 1.0d0
352			dsdr = 0.0d0
353			dspdr = 0.0d0
354			elseif (r.gt.SSD_rup) then
355			s = 0.0d0
356			sp = 0.0d0
357			dsdr = 0.0d0
358			dspdr = 0.0d0
359			else
360			sp = ((SSD_rup + 2.0d0r - 3.0d0SSD_rl) * (SSD_rup-r)**2) / &
361			((SSD_rup - SSD_rl)**3)
362			dspdr = 6.0d0(r-SSD_rup)(r-SSD_rl)/((SSD_rup - SSD_rl)**3)
363	gezelter	319
364	gezelter	322	if (r.gt.SSD_ru) then
365			s = 0.0d0
366			dsdr = 0.0d0
367	gezelter	319	else
368	gezelter	322	s = ((SSD_ru + 2.0d0r - 3.0d0SSD_rl) * (SSD_ru-r)**2) / &
369			((SSD_ru - SSD_rl)**3)
370			dsdr = 6.0d0(r-SSD_ru)(r-SSD_rl)/((SSD_ru - SSD_rl)**3)
371	gezelter	319	endif
372	gezelter	322	endif
373	gezelter	319
374	gezelter	322	return
375			end subroutine calc_sw_fnc
376			end module sticky_pair