UseTheForce/DarkSide/gb.F90

!!
!! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
!!
!! The University of Notre Dame grants you ("Licensee") a
!! non-exclusive, royalty free, license to use, modify and
!! redistribute this software in source and binary code form, provided
!! that the following conditions are met:
!!
!! 1. Redistributions of source code must retain the above copyright
!!    notice, this list of conditions and the following disclaimer.
!!
!! 2. Redistributions in binary form must reproduce the above copyright
!!    notice, this list of conditions and the following disclaimer in the
!!    documentation and/or other materials provided with the
!!    distribution.
!!
!! This software is provided "AS IS," without a warranty of any
!! kind. All express or implied conditions, representations and
!! warranties, including any implied warranty of merchantability,
!! fitness for a particular purpose or non-infringement, are hereby
!! excluded.  The University of Notre Dame and its licensors shall not
!! be liable for any damages suffered by licensee as a result of
!! using, modifying or distributing the software or its
!! derivatives. In no event will the University of Notre Dame or its
!! licensors be liable for any lost revenue, profit or data, or for
!! direct, indirect, special, consequential, incidental or punitive
!! damages, however caused and regardless of the theory of liability,
!! arising out of the use of or inability to use software, even if the
!! University of Notre Dame has been advised of the possibility of
!! such damages.
!!
!! SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
!! research, please cite the appropriate papers when you publish your
!! work.  Good starting points are:
!!                                                                      
!! [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
!! [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
!! [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
!! [4]  Vardeman & Gezelter, in progress (2009).
!!

!! gayberne is the Gay-Berne interaction for ellipsoidal particles.  The original 
!! paper (for identical uniaxial particles) is:
!!    J. G. Gay and B. J. Berne, J. Chem. Phys., 74, 3316-3319 (1981).
!! A more-general GB potential for dissimilar uniaxial particles:
!!    D. J. Cleaver, C. M. Care, M. P. Allen and M. P. Neal, Phys. Rev. E, 
!!    54, 559-567 (1996).
!! Further parameterizations can be found in:
!!    A. P. J. Emerson, G. R. Luckhurst and S. G. Whatling, Mol. Phys., 
!!    82, 113-124 (1994).
!! And a nice force expression:
!!    G. R. Luckhurst and R. A. Stephens, Liq. Cryst. 8, 451-464 (1990).
!! Even clearer force and torque expressions:
!!    P. A. Golubkov and P. Y. Ren, J. Chem. Phys., 125, 64103 (2006).
!! New expressions for cross interactions of strength parameters:
!!    J. Wu, X. Zhen, H. Shen, G. Li, and P. Ren, J. Chem. Phys., 
!!    135, 155104 (2011).
!!
!! In this version of the GB interaction, each uniaxial ellipsoidal type 
!! is described using a set of 6 parameters:
!!  d:  range parameter for side-by-side (S) and cross (X) configurations
!!  l:  range parameter for end-to-end (E) configuration
!!  epsilon_X:  well-depth parameter for cross (X) configuration
!!  epsilon_S:  well-depth parameter for side-by-side (S) configuration
!!  epsilon_E:  well depth parameter for end-to-end (E) configuration
!!  dw: "softness" of the potential
!! 
!! Additionally, there are two "universal" paramters to govern the overall 
!! importance of the purely orientational (nu) and the mixed
!! orientational / translational (mu) parts of strength of the interactions. 
!! These parameters have default or "canonical" values, but may be changed
!! as a force field option:
!! nu_: purely orientational part : defaults to 1
!! mu_: mixed orientational / translational part : defaults to 2

module gayberne
  use force_globals
  use definitions
  use simulation
  use atype_module
  use vector_class
  use linearalgebra
  use status
  use lj
  use fForceOptions
  
  implicit none

  private

#define __FORTRAN90
#include "UseTheForce/DarkSide/fInteractionMap.h"

  logical, save :: haveGBMap = .false.
  logical, save :: haveMixingMap = .false.
  real(kind=dp), save :: mu = 2.0_dp
  real(kind=dp), save :: nu = 1.0_dp


  public :: newGBtype
  public :: complete_GB_FF
  public :: do_gb_pair
  public :: getGayBerneCut
  public :: destroyGBtypes

  type :: GBtype
     integer          :: atid
     real(kind = dp ) :: d
     real(kind = dp ) :: l
     real(kind = dp ) :: epsX
     real(kind = dp ) :: epsS
     real(kind = dp ) :: epsE
     real(kind = dp ) :: dw
     logical          :: isLJ
  end type GBtype
  
  type, private :: GBList
     integer               :: nGBtypes = 0
     integer               :: currentGBtype = 0
     type(GBtype), pointer :: GBtypes(:)      => null()
     integer, pointer      :: atidToGBtype(:) => null()
  end type GBList
  
  type(GBList), save :: GBMap
  
  type :: GBMixParameters
     real(kind=DP) :: sigma0
     real(kind=DP) :: eps0
     real(kind=DP) :: dw
     real(kind=DP) :: x2
     real(kind=DP) :: xa2
     real(kind=DP) :: xai2
     real(kind=DP) :: xp2
     real(kind=DP) :: xpap2
     real(kind=DP) :: xpapi2
  end type GBMixParameters
  
  type(GBMixParameters), dimension(:,:), allocatable :: GBMixingMap
  
contains
  
  subroutine newGBtype(c_ident, d, l, epsX, epsS, epsE, dw, status)
    
    integer, intent(in) :: c_ident
    real( kind = dp ), intent(in) :: d, l, epsX, epsS, epsE, dw
    integer, intent(out) :: status
    
    integer :: nGBTypes, nLJTypes, ntypes, myATID
    integer, pointer :: MatchList(:) => null()
    integer :: current, i
    status = 0
    
    if (.not.associated(GBMap%GBtypes)) then
       
       call getMatchingElementList(atypes, "is_GayBerne", .true., &
            nGBtypes, MatchList)
       
       call getMatchingElementList(atypes, "is_LennardJones", .true., &
            nLJTypes, MatchList)
       
       GBMap%nGBtypes = nGBtypes + nLJTypes
       
       allocate(GBMap%GBtypes(nGBtypes + nLJTypes))
       
       ntypes = getSize(atypes)
       
       allocate(GBMap%atidToGBtype(ntypes))       
    endif
    
    GBMap%currentGBtype = GBMap%currentGBtype + 1
    current = GBMap%currentGBtype
    
    myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
    
    GBMap%atidToGBtype(myATID)       = current
    GBMap%GBtypes(current)%atid      = myATID
    GBMap%GBtypes(current)%d         = d
    GBMap%GBtypes(current)%l         = l
    GBMap%GBtypes(current)%epsX      = epsX
    GBMap%GBtypes(current)%epsS      = epsS
    GBMap%GBtypes(current)%epsE      = epsE
    GBMap%GBtypes(current)%dw        = dw
    GBMap%GBtypes(current)%isLJ      = .false.
    
    return
  end subroutine newGBtype
  
  subroutine complete_GB_FF(status)
    integer :: status
    integer :: i, j, l, m, lm, function_type
    real(kind=dp) :: thisDP, sigma
    integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, myATID, current
    logical :: thisProperty
    
    status = 0
    if (GBMap%currentGBtype == 0) then
       call handleError("complete_GB_FF", "No members in GBMap")
       status = -1
       return
    end if
    
    nAtypes = getSize(atypes)
    
    if (nAtypes == 0) then
       status = -1
       return
    end if
    
    ! atypes comes from c side
    do i = 1, nAtypes
       
       myATID = getFirstMatchingElement(atypes, 'c_ident', i)
       call getElementProperty(atypes, myATID, "is_LennardJones", thisProperty)
       
       if (thisProperty) then
          GBMap%currentGBtype = GBMap%currentGBtype + 1
          current = GBMap%currentGBtype
          
          GBMap%atidToGBtype(myATID) = current
          GBMap%GBtypes(current)%atid      = myATID       
          GBMap%GBtypes(current)%isLJ      = .true.          
          GBMap%GBtypes(current)%d         = getSigma(myATID) / sqrt(2.0_dp)
          GBMap%GBtypes(current)%l         = GBMap%GBtypes(current)%d
          GBMap%GBtypes(current)%epsX      = getEpsilon(myATID)
          GBMap%GBtypes(current)%epsS      = GBMap%GBtypes(current)%epsX
          GBMap%GBtypes(current)%epsE      = GBMap%GBtypes(current)%epsX
          GBMap%GBtypes(current)%dw        = 1.0_dp
          
       endif
       
    end do
    
    haveGBMap = .true.

    
  end subroutine complete_GB_FF

  subroutine createGBMixingMap()
    integer :: nGBtypes, i, j
    real (kind = dp) :: d1, l1, eX1, eS1, eE1, dw1
    real (kind = dp) :: d2, l2, eX2, eS2, eE2, dw2
    real (kind = dp) :: xp, ap2, mi

    if (GBMap%currentGBtype == 0) then
       call handleError("GB", "No members in GBMap")
       return
    end if
    
    nGBtypes = GBMap%nGBtypes

    if (.not. allocated(GBMixingMap)) then
       allocate(GBMixingMap(nGBtypes, nGBtypes))
    endif

    do i = 1, nGBtypes

       d1 = GBMap%GBtypes(i)%d
       l1 = GBMap%GBtypes(i)%l
       eX1 = GBMap%GBtypes(i)%epsX
       eS1 = GBMap%GBtypes(i)%epsS
       eE1 = GBMap%GBtypes(i)%epsE
       dw1 = GBMap%GBtypes(i)%dw

       do j = 1, nGBtypes

          d2 = GBMap%GBtypes(j)%d
          l2 = GBMap%GBtypes(j)%l
          eX2 = GBMap%GBtypes(j)%epsX
          eS2 = GBMap%GBtypes(j)%epsS
          eE2 = GBMap%GBtypes(j)%epsE
          dw2 = GBMap%GBtypes(j)%dw

!  Cleaver paper uses sqrt of squares to get sigma0 for
!  mixed interactions.
            
          GBMixingMap(i,j)%sigma0 = sqrt(d1*d1 + d2*d2)
          GBMixingMap(i,j)%xa2 = (l1*l1 - d1*d1)/(l1*l1 + d2*d2)
          GBMixingMap(i,j)%xai2 = (l2*l2 - d2*d2)/(l2*l2 + d1*d1)
          GBMixingMap(i,j)%x2 = (l1*l1 - d1*d1) * (l2*l2 - d2*d2) / &
               ((l2*l2 + d1*d1) * (l1*l1 + d2*d2))

          ! assumed LB mixing rules for now:

          GBMixingMap(i,j)%dw = 0.5_dp * (dw1 + dw2)
          GBMixingMap(i,j)%eps0 = sqrt(eX1 * eX2)

          mi = 1.0 / mu
      
          GBMixingMap(i,j)%xpap2  = ((eS1**mi) - (eE1**mi)) / &
               ((eS1**mi) + (eE2**mi))
          GBMixingMap(i,j)%xpapi2 = ((eS2**mi) - (eE2**mi)) / &
               ((eS2**mi) + (eE1**mi))
          GBMixingMap(i,j)%xp2    = ((eS1**mi) - (eE1**mi)) * &
               ((eS2**mi) - (eE2**mi)) / &
               ((eS2**mi) + (eE1**mi)) / ((eS1**mi) + (eE2**mi))

       enddo
    enddo
    haveMixingMap = .true.
    mu = getGayBerneMu()
    nu = getGayBerneNu()    
  end subroutine createGBMixingMap
  

  !! gay berne cutoff should be a parameter in globals, this is a temporary 
  !! work around - this should be fixed when gay berne is up and running

  function getGayBerneCut(atomID) result(cutValue)
    integer, intent(in) :: atomID 
    integer :: gbt1
    real(kind=dp) :: cutValue, l, d

    if (GBMap%currentGBtype == 0) then
       call handleError("GB", "No members in GBMap")
       return
    end if

    gbt1 = GBMap%atidToGBtype(atomID)
    l = GBMap%GBtypes(gbt1)%l
    d = GBMap%GBtypes(gbt1)%d   

    ! sigma is actually sqrt(2)*l  for prolate ellipsoids
    
    cutValue = 2.5_dp*sqrt(2.0_dp)*max(l,d)

  end function getGayBerneCut

  subroutine do_gb_pair(atid1, atid2, d, r, r2, sw, vdwMult, vpair, fpair, &
       pot, A1, A2, f1, t1, t2)
    
    integer, intent(in) :: atid1, atid2
    integer :: gbt1, gbt2, id1, id2
    real (kind=dp), intent(inout) :: r, r2, vdwMult
    real (kind=dp), dimension(3), intent(in) :: d
    real (kind=dp), dimension(3), intent(inout) :: fpair
    real (kind=dp) :: pot, sw, vpair
    real (kind=dp), dimension(9) :: A1, A2
    real (kind=dp), dimension(3) :: f1
    real (kind=dp), dimension(3) :: t1, t2
    real (kind = dp), dimension(3) :: ul1, ul2, rxu1, rxu2, uxu, rhat

    real (kind = dp) :: sigma0, dw, eps0, x2, xa2, xai2, xp2, xpap2, xpapi2
    real (kind = dp) :: e1, e2, eps, sigma, s3, s03, au2, bu2, au, bu, a, b, g, g2
    real (kind = dp) :: U, BigR, R3, R6, R7, R12, R13, H, Hp, fx, fy, fz
    real (kind = dp) :: dUdr, dUda, dUdb, dUdg, pref1, pref2
    logical :: i_is_lj, j_is_lj

    if (.not.haveMixingMap) then
       call createGBMixingMap()
    endif

    gbt1 = GBMap%atidToGBtype(atid1)
    gbt2 = GBMap%atidToGBtype(atid2)    

    i_is_LJ = GBMap%GBTypes(gbt1)%isLJ
    j_is_LJ = GBMap%GBTypes(gbt2)%isLJ

    sigma0 = GBMixingMap(gbt1, gbt2)%sigma0
    dw     = GBMixingMap(gbt1, gbt2)%dw    
    eps0   = GBMixingMap(gbt1, gbt2)%eps0  
    x2     = GBMixingMap(gbt1, gbt2)%x2    
    xa2    = GBMixingMap(gbt1, gbt2)%xa2   
    xai2   = GBMixingMap(gbt1, gbt2)%xai2  
    xp2    = GBMixingMap(gbt1, gbt2)%xp2   
    xpap2  = GBMixingMap(gbt1, gbt2)%xpap2 
    xpapi2 = GBMixingMap(gbt1, gbt2)%xpapi2
    
!!$    write(*,*) 'atypes = ',atid1, atid2
!!$    write(*,*) 'sigma0 = ',sigma0 
!!$    write(*,*) 'dw     = ',dw 
!!$    write(*,*) 'eps0   = ',eps0   
!!$    write(*,*) 'x2     = ',x2     
!!$    write(*,*) 'xa2    = ',xa2    
!!$    write(*,*) 'xai2   = ',xai2   
!!$    write(*,*) 'xp2    = ',xp2    
!!$    write(*,*) 'xpap2  = ',xpap2  
!!$    write(*,*) 'xpapi2 = ',xpapi2 


    ul1(1) = A1(7)
    ul1(2) = A1(8)
    ul1(3) = A1(9)

    ul2(1) = A2(7)
    ul2(2) = A2(8)
    ul2(3) = A2(9)
    
!!$    write(*,*) 'ul1 = ', ul1(1), ul1(2), ul1(3)
!!$    write(*,*) 'ul2 = ', ul2(1), ul2(2), ul2(3)

    if (i_is_LJ) then
       a = 0.0_dp
       ul1 = 0.0_dp
    else
       a = d(1)*ul1(1)   + d(2)*ul1(2)   + d(3)*ul1(3)
    endif

    if (j_is_LJ) then
       b = 0.0_dp
       ul2 = 0.0_dp
    else       
       b = d(1)*ul2(1)   + d(2)*ul2(2)   + d(3)*ul2(3)
    endif

    if (i_is_LJ.or.j_is_LJ) then
       g = 0.0_dp
    else
       g = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3)
    endif

    au = a / r
    bu = b / r

    au2 = au * au
    bu2 = bu * bu
    g2 = g * g

    H  = (xa2 * au2 + xai2 * bu2 - 2.0_dp*x2*au*bu*g)  / (1.0_dp - x2*g2)
    Hp = (xpap2*au2 + xpapi2*bu2 - 2.0_dp*xp2*au*bu*g) / (1.0_dp - xp2*g2)

!!$    write(*,*) 'au2 = ',au2
!!$    write(*,*) 'bu2 = ',bu2
!!$    write(*,*) 'g2 = ',g2
!!$    write(*,*) 'H = ',H
!!$    write(*,*) 'Hp = ',Hp

    sigma = sigma0 / sqrt(1.0_dp - H)
    e1 = 1.0_dp / sqrt(1.0_dp - x2*g2)
    e2 = 1.0_dp - Hp
    eps = eps0 * (e1**nu) * (e2**mu)
    BigR = dw*sigma0 / (r - sigma + dw*sigma0)
    
    R3 = BigR*BigR*BigR
    R6 = R3*R3
    R7 = R6 * BigR
    R12 = R6*R6
    R13 = R6*R7

    U = vdwMult * 4.0_dp * eps * (R12 - R6)

    s3 = sigma*sigma*sigma
    s03 = sigma0*sigma0*sigma0

!!$    write(*,*) 'vdwMult = ', vdwMult
!!$    write(*,*) 'eps = ', eps
!!$    write(*,*) 'mu = ', mu
!!$    write(*,*) 'R12 = ', R12
!!$    write(*,*) 'R6 = ', R6 
!!$    write(*,*) 'R13 = ', R13
!!$    write(*,*) 'R7 = ', R7 
!!$    write(*,*) 'e2 = ', e2 
!!$    write(*,*) 'rij = ', r
!!$    write(*,*) 's3 = ', s3 
!!$    write(*,*) 's03 = ', s03 
!!$    write(*,*) 'dw = ', dw 

    pref1 = - vdwMult * 8.0_dp * eps * mu * (R12 - R6) / (e2 * r)

    pref2 = vdwMult * 8.0_dp * eps * s3 * (6.0_dp*R13 - 3.0_dp*R7) / (dw*r*s03)

    dUdr = - (pref1 * Hp + pref2 * (sigma0*sigma0*r/s3 + H))
    
    dUda = pref1 * (xpap2*au - xp2*bu*g) / (1.0_dp - xp2 * g2) &
         + pref2 * (xa2 * au - x2 *bu*g) / (1.0_dp - x2 * g2)

    dUdb = pref1 * (xpapi2*bu - xp2*au*g) / (1.0_dp - xp2 * g2) &
         + pref2 * (xai2 * bu - x2 *au*g) / (1.0_dp - x2 * g2)

    dUdg = 4.0_dp * eps * nu * (R12 - R6) * x2 * g / (1.0_dp - x2*g2) &
         + 8.0_dp * eps * mu * (R12 - R6) * (xp2*au*bu - Hp*xp2*g) / &
         (1.0_dp - xp2 * g2) / e2 &
         + 8.0_dp * eps * s3 * (3.0_dp * R7 - 6.0_dp * R13) * &  
         (x2 * au * bu - H * x2 * g) / (1.0_dp - x2 * g2) / (dw * s03)
!!$    write(*,*) 'pref = ',pref1 ,  pref2
!!$    write(*,*) 'dU = ',dUdr ,  dUda, dUdb ,  dUdg

            
    rhat = d / r

    fx = dUdr * rhat(1) + dUda * ul1(1) + dUdb * ul2(1)
    fy = dUdr * rhat(2) + dUda * ul1(2) + dUdb * ul2(2)
    fz = dUdr * rhat(3) + dUda * ul1(3) + dUdb * ul2(3)    

    rxu1 = cross_product(d, ul1)
    rxu2 = cross_product(d, ul2)    
    uxu = cross_product(ul1, ul2)

    pot = pot + U*sw

    f1(1) = f1(1) + fx*sw
    f1(2) = f1(2) + fy*sw
    f1(3) = f1(3) + fz*sw

    t1(1) = t1(1) + (dUda*rxu1(1) - dUdg*uxu(1))*sw
    t1(2) = t1(2) + (dUda*rxu1(2) - dUdg*uxu(2))*sw
    t1(3) = t1(3) + (dUda*rxu1(3) - dUdg*uxu(3))*sw
                                                   
    t2(1) = t2(1) + (dUdb*rxu2(1) + dUdg*uxu(1))*sw
    t2(2) = t2(2) + (dUdb*rxu2(2) + dUdg*uxu(2))*sw
    t2(3) = t2(3) + (dUdb*rxu2(3) + dUdg*uxu(3))*sw

    vpair = vpair + U

!!$    write(*,*) 'f1 term = ', fx*sw, fy*sw, fz*sw
!!$    write(*,*) 't1 term = ', (dUda*rxu1(1) - dUdg*uxu(1))*sw, &
!!$         (dUda*rxu1(2) - dUdg*uxu(2))*sw, &
!!$         (dUda*rxu1(3) - dUdg*uxu(3))*sw                              
!!$    write(*,*) 't2 term = ', (dUdb*rxu2(1) + dUdg*uxu(1))*sw, &
!!$         (dUdb*rxu2(2) + dUdg*uxu(2))*sw, &
!!$         (dUdb*rxu2(3) + dUdg*uxu(3))*sw
!!$    write(*,*) 'vp term = ', U

    return
  end subroutine do_gb_pair
  
  subroutine destroyGBTypes()

    GBMap%nGBtypes = 0
    GBMap%currentGBtype = 0
    
    if (associated(GBMap%GBtypes)) then
       deallocate(GBMap%GBtypes)
       GBMap%GBtypes => null()
    end if
    
    if (associated(GBMap%atidToGBtype)) then
       deallocate(GBMap%atidToGBtype)
       GBMap%atidToGBtype => null()
    end if
    
    haveMixingMap = .false.
    
  end subroutine destroyGBTypes

end module gayberne
    
Revision:	1689
Committed:	Wed Mar 14 17:57:08 2012 UTC (13 years, 5 months ago) by gezelter
File size:	17400 byte(s)
Log Message:	Bug fixes for GB. Now using strength parameter mixing ideas from Wu et al. [J. Chem. Phys. 135, 155104 (2011)]. This helps get the dissimilar particle mixing behavior to be the same whichever order the two particles come in. This does require that the force field file to specify explicitly the values for epsilon in the cross (X), side-by-side (S), and end-to-end (E) configurations.
#	User	Rev	Content
1	gezelter	246	!!
2			!! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			!!
4			!! The University of Notre Dame grants you ("Licensee") a
5			!! non-exclusive, royalty free, license to use, modify and
6			!! redistribute this software in source and binary code form, provided
7			!! that the following conditions are met:
8			!!
9	gezelter	1390	!! 1. Redistributions of source code must retain the above copyright
10	gezelter	246	!! notice, this list of conditions and the following disclaimer.
11			!!
12	gezelter	1390	!! 2. Redistributions in binary form must reproduce the above copyright
13	gezelter	246	!! notice, this list of conditions and the following disclaimer in the
14			!! documentation and/or other materials provided with the
15			!! distribution.
16			!!
17			!! This software is provided "AS IS," without a warranty of any
18			!! kind. All express or implied conditions, representations and
19			!! warranties, including any implied warranty of merchantability,
20			!! fitness for a particular purpose or non-infringement, are hereby
21			!! excluded. The University of Notre Dame and its licensors shall not
22			!! be liable for any damages suffered by licensee as a result of
23			!! using, modifying or distributing the software or its
24			!! derivatives. In no event will the University of Notre Dame or its
25			!! licensors be liable for any lost revenue, profit or data, or for
26			!! direct, indirect, special, consequential, incidental or punitive
27			!! damages, however caused and regardless of the theory of liability,
28			!! arising out of the use of or inability to use software, even if the
29			!! University of Notre Dame has been advised of the possibility of
30			!! such damages.
31			!!
32	gezelter	1390	!! SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33			!! research, please cite the appropriate papers when you publish your
34			!! work. Good starting points are:
35			!!
36			!! [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37			!! [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38			!! [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39			!! [4] Vardeman & Gezelter, in progress (2009).
40			!!
41	gezelter	246
42	gezelter	1689	!! gayberne is the Gay-Berne interaction for ellipsoidal particles. The original
43			!! paper (for identical uniaxial particles) is:
44			!! J. G. Gay and B. J. Berne, J. Chem. Phys., 74, 3316-3319 (1981).
45			!! A more-general GB potential for dissimilar uniaxial particles:
46			!! D. J. Cleaver, C. M. Care, M. P. Allen and M. P. Neal, Phys. Rev. E,
47			!! 54, 559-567 (1996).
48			!! Further parameterizations can be found in:
49			!! A. P. J. Emerson, G. R. Luckhurst and S. G. Whatling, Mol. Phys.,
50			!! 82, 113-124 (1994).
51			!! And a nice force expression:
52			!! G. R. Luckhurst and R. A. Stephens, Liq. Cryst. 8, 451-464 (1990).
53			!! Even clearer force and torque expressions:
54			!! P. A. Golubkov and P. Y. Ren, J. Chem. Phys., 125, 64103 (2006).
55			!! New expressions for cross interactions of strength parameters:
56			!! J. Wu, X. Zhen, H. Shen, G. Li, and P. Ren, J. Chem. Phys.,
57			!! 135, 155104 (2011).
58			!!
59			!! In this version of the GB interaction, each uniaxial ellipsoidal type
60			!! is described using a set of 6 parameters:
61			!! d: range parameter for side-by-side (S) and cross (X) configurations
62			!! l: range parameter for end-to-end (E) configuration
63			!! epsilon_X: well-depth parameter for cross (X) configuration
64			!! epsilon_S: well-depth parameter for side-by-side (S) configuration
65			!! epsilon_E: well depth parameter for end-to-end (E) configuration
66			!! dw: "softness" of the potential
67			!!
68			!! Additionally, there are two "universal" paramters to govern the overall
69			!! importance of the purely orientational (nu) and the mixed
70			!! orientational / translational (mu) parts of strength of the interactions.
71			!! These parameters have default or "canonical" values, but may be changed
72			!! as a force field option:
73			!! nu_: purely orientational part : defaults to 1
74			!! mu_: mixed orientational / translational part : defaults to 2
75	gezelter	246
76	gezelter	676	module gayberne
77	gezelter	115	use force_globals
78			use definitions
79			use simulation
80	kdaily	668	use atype_module
81			use vector_class
82	gezelter	981	use linearalgebra
83	kdaily	668	use status
84			use lj
85	gezelter	982	use fForceOptions
86	kdaily	529
87	gezelter	115	implicit none
88
89	kdaily	668	private
90	gezelter	115
91	gezelter	676	#define __FORTRAN90
92			#include "UseTheForce/DarkSide/fInteractionMap.h"
93	gezelter	115
94	gezelter	981	logical, save :: haveGBMap = .false.
95			logical, save :: haveMixingMap = .false.
96			real(kind=dp), save :: mu = 2.0_dp
97			real(kind=dp), save :: nu = 1.0_dp
98
99
100	gezelter	676	public :: newGBtype
101	gezelter	981	public :: complete_GB_FF
102	gezelter	115	public :: do_gb_pair
103	chrisfen	578	public :: getGayBerneCut
104	gezelter	676	public :: destroyGBtypes
105	gezelter	115
106	gezelter	676	type :: GBtype
107			integer :: atid
108	gezelter	981	real(kind = dp ) :: d
109			real(kind = dp ) :: l
110	gezelter	1689	real(kind = dp ) :: epsX
111			real(kind = dp ) :: epsS
112			real(kind = dp ) :: epsE
113	gezelter	981	real(kind = dp ) :: dw
114			logical :: isLJ
115	gezelter	676	end type GBtype
116	gezelter	981
117	gezelter	676	type, private :: GBList
118			integer :: nGBtypes = 0
119			integer :: currentGBtype = 0
120			type(GBtype), pointer :: GBtypes(:) => null()
121			integer, pointer :: atidToGBtype(:) => null()
122			end type GBList
123	gezelter	981
124	gezelter	676	type(GBList), save :: GBMap
125	gezelter	981
126			type :: GBMixParameters
127			real(kind=DP) :: sigma0
128			real(kind=DP) :: eps0
129			real(kind=DP) :: dw
130			real(kind=DP) :: x2
131			real(kind=DP) :: xa2
132			real(kind=DP) :: xai2
133			real(kind=DP) :: xp2
134			real(kind=DP) :: xpap2
135			real(kind=DP) :: xpapi2
136			end type GBMixParameters
137
138			type(GBMixParameters), dimension(:,:), allocatable :: GBMixingMap
139
140	gezelter	115	contains
141	gezelter	981
142	gezelter	1689	subroutine newGBtype(c_ident, d, l, epsX, epsS, epsE, dw, status)
143	gezelter	676
144			integer, intent(in) :: c_ident
145	gezelter	1689	real( kind = dp ), intent(in) :: d, l, epsX, epsS, epsE, dw
146	gezelter	676	integer, intent(out) :: status
147	gezelter	981
148			integer :: nGBTypes, nLJTypes, ntypes, myATID
149	gezelter	676	integer, pointer :: MatchList(:) => null()
150			integer :: current, i
151	kdaily	668	status = 0
152	gezelter	981
153	gezelter	676	if (.not.associated(GBMap%GBtypes)) then
154	gezelter	981
155	gezelter	676	call getMatchingElementList(atypes, "is_GayBerne", .true., &
156			nGBtypes, MatchList)
157
158	gezelter	981	call getMatchingElementList(atypes, "is_LennardJones", .true., &
159			nLJTypes, MatchList)
160
161			GBMap%nGBtypes = nGBtypes + nLJTypes
162
163			allocate(GBMap%GBtypes(nGBtypes + nLJTypes))
164
165	gezelter	676	ntypes = getSize(atypes)
166
167	gezelter	981	allocate(GBMap%atidToGBtype(ntypes))
168			endif
169
170			GBMap%currentGBtype = GBMap%currentGBtype + 1
171			current = GBMap%currentGBtype
172
173			myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
174
175			GBMap%atidToGBtype(myATID) = current
176			GBMap%GBtypes(current)%atid = myATID
177			GBMap%GBtypes(current)%d = d
178			GBMap%GBtypes(current)%l = l
179	gezelter	1689	GBMap%GBtypes(current)%epsX = epsX
180			GBMap%GBtypes(current)%epsS = epsS
181			GBMap%GBtypes(current)%epsE = epsE
182	gezelter	981	GBMap%GBtypes(current)%dw = dw
183			GBMap%GBtypes(current)%isLJ = .false.
184
185			return
186			end subroutine newGBtype
187
188			subroutine complete_GB_FF(status)
189			integer :: status
190			integer :: i, j, l, m, lm, function_type
191			real(kind=dp) :: thisDP, sigma
192			integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, myATID, current
193			logical :: thisProperty
194
195			status = 0
196			if (GBMap%currentGBtype == 0) then
197			call handleError("complete_GB_FF", "No members in GBMap")
198			status = -1
199			return
200			end if
201
202			nAtypes = getSize(atypes)
203
204			if (nAtypes == 0) then
205			status = -1
206			return
207			end if
208
209			! atypes comes from c side
210			do i = 1, nAtypes
211	gezelter	676
212	gezelter	981	myATID = getFirstMatchingElement(atypes, 'c_ident', i)
213			call getElementProperty(atypes, myATID, "is_LennardJones", thisProperty)
214
215			if (thisProperty) then
216			GBMap%currentGBtype = GBMap%currentGBtype + 1
217			current = GBMap%currentGBtype
218
219			GBMap%atidToGBtype(myATID) = current
220			GBMap%GBtypes(current)%atid = myATID
221			GBMap%GBtypes(current)%isLJ = .true.
222	xsun	1006	GBMap%GBtypes(current)%d = getSigma(myATID) / sqrt(2.0_dp)
223	gezelter	981	GBMap%GBtypes(current)%l = GBMap%GBtypes(current)%d
224	gezelter	1689	GBMap%GBtypes(current)%epsX = getEpsilon(myATID)
225			GBMap%GBtypes(current)%epsS = GBMap%GBtypes(current)%epsX
226			GBMap%GBtypes(current)%epsE = GBMap%GBtypes(current)%epsX
227	gezelter	981	GBMap%GBtypes(current)%dw = 1.0_dp
228
229			endif
230
231			end do
232
233			haveGBMap = .true.
234	gezelter	982
235	gezelter	981
236			end subroutine complete_GB_FF
237	kdaily	668
238	gezelter	981	subroutine createGBMixingMap()
239			integer :: nGBtypes, i, j
240	gezelter	1689	real (kind = dp) :: d1, l1, eX1, eS1, eE1, dw1
241			real (kind = dp) :: d2, l2, eX2, eS2, eE2, dw2
242			real (kind = dp) :: xp, ap2, mi
243	kdaily	668
244	gezelter	981	if (GBMap%currentGBtype == 0) then
245			call handleError("GB", "No members in GBMap")
246			return
247			end if
248
249			nGBtypes = GBMap%nGBtypes
250
251			if (.not. allocated(GBMixingMap)) then
252			allocate(GBMixingMap(nGBtypes, nGBtypes))
253	gezelter	676	endif
254	kdaily	668
255	gezelter	981	do i = 1, nGBtypes
256	kdaily	668
257	gezelter	981	d1 = GBMap%GBtypes(i)%d
258			l1 = GBMap%GBtypes(i)%l
259	gezelter	1689	eX1 = GBMap%GBtypes(i)%epsX
260			eS1 = GBMap%GBtypes(i)%epsS
261			eE1 = GBMap%GBtypes(i)%epsE
262	gezelter	981	dw1 = GBMap%GBtypes(i)%dw
263	gezelter	676
264	xsun	1010	do j = 1, nGBtypes
265	gezelter	676
266	gezelter	981	d2 = GBMap%GBtypes(j)%d
267			l2 = GBMap%GBtypes(j)%l
268	gezelter	1689	eX2 = GBMap%GBtypes(j)%epsX
269			eS2 = GBMap%GBtypes(j)%epsS
270			eE2 = GBMap%GBtypes(j)%epsE
271	gezelter	981	dw2 = GBMap%GBtypes(j)%dw
272	xsun	1006
273			! Cleaver paper uses sqrt of squares to get sigma0 for
274			! mixed interactions.
275
276			GBMixingMap(i,j)%sigma0 = sqrt(d1d1 + d2d2)
277	gezelter	981	GBMixingMap(i,j)%xa2 = (l1l1 - d1d1)/(l1l1 + d2d2)
278			GBMixingMap(i,j)%xai2 = (l2l2 - d2d2)/(l2l2 + d1d1)
279			GBMixingMap(i,j)%x2 = (l1l1 - d1d1) * (l2l2 - d2d2) / &
280			((l2l2 + d1d1) * (l1l1 + d2d2))
281
282			! assumed LB mixing rules for now:
283
284			GBMixingMap(i,j)%dw = 0.5_dp * (dw1 + dw2)
285	gezelter	1689	GBMixingMap(i,j)%eps0 = sqrt(eX1 * eX2)
286	gezelter	981
287	gezelter	1689	mi = 1.0 / mu
288
289			GBMixingMap(i,j)%xpap2 = ((eS1mi) - (eE1mi)) / &
290			((eS1mi) + (eE2mi))
291			GBMixingMap(i,j)%xpapi2 = ((eS2mi) - (eE2mi)) / &
292			((eS2mi) + (eE1mi))
293			GBMixingMap(i,j)%xp2 = ((eS1mi) - (eE1mi)) * &
294			((eS2mi) - (eE2mi)) / &
295			((eS2mi) + (eE1mi)) / ((eS1mi) + (eE2mi))
296	gezelter	981
297			enddo
298			enddo
299			haveMixingMap = .true.
300	gezelter	983	mu = getGayBerneMu()
301			nu = getGayBerneNu()
302	gezelter	981	end subroutine createGBMixingMap
303	gezelter	676
304	gezelter	981
305	chrisfen	580	!! gay berne cutoff should be a parameter in globals, this is a temporary
306			!! work around - this should be fixed when gay berne is up and running
307	gezelter	676
308	chrisfen	578	function getGayBerneCut(atomID) result(cutValue)
309	gezelter	676	integer, intent(in) :: atomID
310			integer :: gbt1
311	gezelter	981	real(kind=dp) :: cutValue, l, d
312	gezelter	115
313	gezelter	676	if (GBMap%currentGBtype == 0) then
314			call handleError("GB", "No members in GBMap")
315			return
316			end if
317
318			gbt1 = GBMap%atidToGBtype(atomID)
319	gezelter	981	l = GBMap%GBtypes(gbt1)%l
320			d = GBMap%GBtypes(gbt1)%d
321	gezelter	676
322	xsun	1010	! sigma is actually sqrt(2)*l for prolate ellipsoids
323
324			cutValue = 2.5_dpsqrt(2.0_dp)max(l,d)
325
326	chrisfen	578	end function getGayBerneCut
327
328	gezelter	1464	subroutine do_gb_pair(atid1, atid2, d, r, r2, sw, vdwMult, vpair, fpair, &
329			pot, A1, A2, f1, t1, t2)
330	kdaily	529
331	gezelter	1464	integer, intent(in) :: atid1, atid2
332	gezelter	1386	integer :: gbt1, gbt2, id1, id2
333	gezelter	1286	real (kind=dp), intent(inout) :: r, r2, vdwMult
334	gezelter	115	real (kind=dp), dimension(3), intent(in) :: d
335			real (kind=dp), dimension(3), intent(inout) :: fpair
336			real (kind=dp) :: pot, sw, vpair
337	gezelter	1386	real (kind=dp), dimension(9) :: A1, A2
338			real (kind=dp), dimension(3) :: f1
339			real (kind=dp), dimension(3) :: t1, t2
340	gezelter	981	real (kind = dp), dimension(3) :: ul1, ul2, rxu1, rxu2, uxu, rhat
341	gezelter	115
342	gezelter	981	real (kind = dp) :: sigma0, dw, eps0, x2, xa2, xai2, xp2, xpap2, xpapi2
343	kdaily	997	real (kind = dp) :: e1, e2, eps, sigma, s3, s03, au2, bu2, au, bu, a, b, g, g2
344	gezelter	981	real (kind = dp) :: U, BigR, R3, R6, R7, R12, R13, H, Hp, fx, fy, fz
345			real (kind = dp) :: dUdr, dUda, dUdb, dUdg, pref1, pref2
346			logical :: i_is_lj, j_is_lj
347
348			if (.not.haveMixingMap) then
349			call createGBMixingMap()
350			endif
351
352	gezelter	676	gbt1 = GBMap%atidToGBtype(atid1)
353	gezelter	981	gbt2 = GBMap%atidToGBtype(atid2)
354	gezelter	676
355	gezelter	981	i_is_LJ = GBMap%GBTypes(gbt1)%isLJ
356			j_is_LJ = GBMap%GBTypes(gbt2)%isLJ
357	gezelter	676
358	gezelter	981	sigma0 = GBMixingMap(gbt1, gbt2)%sigma0
359			dw = GBMixingMap(gbt1, gbt2)%dw
360			eps0 = GBMixingMap(gbt1, gbt2)%eps0
361			x2 = GBMixingMap(gbt1, gbt2)%x2
362			xa2 = GBMixingMap(gbt1, gbt2)%xa2
363			xai2 = GBMixingMap(gbt1, gbt2)%xai2
364			xp2 = GBMixingMap(gbt1, gbt2)%xp2
365			xpap2 = GBMixingMap(gbt1, gbt2)%xpap2
366			xpapi2 = GBMixingMap(gbt1, gbt2)%xpapi2
367
368	gezelter	1689	!!$ write(,) 'atypes = ',atid1, atid2
369			!!$ write(,) 'sigma0 = ',sigma0
370			!!$ write(,) 'dw = ',dw
371			!!$ write(,) 'eps0 = ',eps0
372			!!$ write(,) 'x2 = ',x2
373			!!$ write(,) 'xa2 = ',xa2
374			!!$ write(,) 'xai2 = ',xai2
375			!!$ write(,) 'xp2 = ',xp2
376			!!$ write(,) 'xpap2 = ',xpap2
377			!!$ write(,) 'xpapi2 = ',xpapi2
378
379
380	gezelter	1386	ul1(1) = A1(7)
381			ul1(2) = A1(8)
382			ul1(3) = A1(9)
383	gezelter	115
384	gezelter	1386	ul2(1) = A2(7)
385			ul2(2) = A2(8)
386			ul2(3) = A2(9)
387	kdaily	529
388	gezelter	1689	!!$ write(,) 'ul1 = ', ul1(1), ul1(2), ul1(3)
389			!!$ write(,) 'ul2 = ', ul2(1), ul2(2), ul2(3)
390
391	gezelter	981	if (i_is_LJ) then
392			a = 0.0_dp
393			ul1 = 0.0_dp
394			else
395			a = d(1)ul1(1) + d(2)ul1(2) + d(3)*ul1(3)
396			endif
397	gezelter	115
398	gezelter	981	if (j_is_LJ) then
399			b = 0.0_dp
400			ul2 = 0.0_dp
401			else
402			b = d(1)ul2(1) + d(2)ul2(2) + d(3)*ul2(3)
403			endif
404	gezelter	115
405	gezelter	981	if (i_is_LJ.or.j_is_LJ) then
406			g = 0.0_dp
407			else
408			g = ul1(1)ul2(1) + ul1(2)ul2(2) + ul1(3)*ul2(3)
409			endif
410	gezelter	1002
411	gezelter	981	au = a / r
412			bu = b / r
413	kdaily	997
414	gezelter	1002	au2 = au * au
415			bu2 = bu * bu
416			g2 = g * g
417	gezelter	115
418	kdaily	997	H = (xa2 * au2 + xai2 * bu2 - 2.0_dpx2aubug) / (1.0_dp - x2*g2)
419			Hp = (xpap2au2 + xpapi2bu2 - 2.0_dpxp2aubug) / (1.0_dp - xp2*g2)
420	gezelter	1002
421	gezelter	1689	!!$ write(,) 'au2 = ',au2
422			!!$ write(,) 'bu2 = ',bu2
423			!!$ write(,) 'g2 = ',g2
424			!!$ write(,) 'H = ',H
425			!!$ write(,) 'Hp = ',Hp
426
427	gezelter	981	sigma = sigma0 / sqrt(1.0_dp - H)
428			e1 = 1.0_dp / sqrt(1.0_dp - x2*g2)
429			e2 = 1.0_dp - Hp
430			eps = eps0 * (e1*nu) (e2**mu)
431			BigR = dwsigma0 / (r - sigma + dwsigma0)
432
433			R3 = BigRBigRBigR
434			R6 = R3*R3
435			R7 = R6 * BigR
436			R12 = R6*R6
437			R13 = R6*R7
438	gezelter	115
439	gezelter	1286	U = vdwMult * 4.0_dp * eps * (R12 - R6)
440	gezelter	115
441	gezelter	981	s3 = sigmasigmasigma
442			s03 = sigma0sigma0sigma0
443	kdaily	529
444	gezelter	1689	!!$ write(,) 'vdwMult = ', vdwMult
445			!!$ write(,) 'eps = ', eps
446			!!$ write(,) 'mu = ', mu
447			!!$ write(,) 'R12 = ', R12
448			!!$ write(,) 'R6 = ', R6
449			!!$ write(,) 'R13 = ', R13
450			!!$ write(,) 'R7 = ', R7
451			!!$ write(,) 'e2 = ', e2
452			!!$ write(,) 'rij = ', r
453			!!$ write(,) 's3 = ', s3
454			!!$ write(,) 's03 = ', s03
455			!!$ write(,) 'dw = ', dw
456
457	gezelter	1286	pref1 = - vdwMult * 8.0_dp * eps * mu * (R12 - R6) / (e2 * r)
458	gezelter	983
459	gezelter	1286	pref2 = vdwMult * 8.0_dp * eps * s3 * (6.0_dpR13 - 3.0_dpR7) / (dwrs03)
460	gezelter	507
461	kdaily	997	dUdr = - (pref1 * Hp + pref2 * (sigma0sigma0r/s3 + H))
462	kdaily	529
463	gezelter	981	dUda = pref1 * (xpap2au - xp2bug) / (1.0_dp - xp2 g2) &
464			+ pref2 * (xa2 * au - x2 bug) / (1.0_dp - x2 * g2)
465	kdaily	997
466	gezelter	981	dUdb = pref1 * (xpapi2bu - xp2aug) / (1.0_dp - xp2 g2) &
467			+ pref2 * (xai2 * bu - x2 aug) / (1.0_dp - x2 * g2)
468	gezelter	507
469	gezelter	981	dUdg = 4.0_dp * eps * nu * (R12 - R6) * x2 * g / (1.0_dp - x2*g2) &
470			+ 8.0_dp * eps * mu * (R12 - R6) * (xp2aubu - Hpxp2g) / &
471			(1.0_dp - xp2 * g2) / e2 &
472	gezelter	983	+ 8.0_dp * eps * s3 * (3.0_dp * R7 - 6.0_dp * R13) * &
473	gezelter	981	(x2 * au * bu - H * x2 * g) / (1.0_dp - x2 * g2) / (dw * s03)
474	gezelter	1689	!!$ write(,) 'pref = ',pref1 , pref2
475			!!$ write(,) 'dU = ',dUdr , dUda, dUdb , dUdg
476	kdaily	997
477	gezelter	981
478			rhat = d / r
479	gezelter	507
480	gezelter	983	fx = dUdr * rhat(1) + dUda * ul1(1) + dUdb * ul2(1)
481			fy = dUdr * rhat(2) + dUda * ul1(2) + dUdb * ul2(2)
482			fz = dUdr * rhat(3) + dUda * ul1(3) + dUdb * ul2(3)
483	gezelter	115
484	gezelter	981	rxu1 = cross_product(d, ul1)
485			rxu2 = cross_product(d, ul2)
486			uxu = cross_product(ul1, ul2)
487	gezelter	1386
488			pot = pot + U*sw
489
490	gezelter	1687	f1(1) = f1(1) + fx*sw
491			f1(2) = f1(2) + fy*sw
492			f1(3) = f1(3) + fz*sw
493
494			t1(1) = t1(1) + (dUdarxu1(1) - dUdguxu(1))*sw
495			t1(2) = t1(2) + (dUdarxu1(2) - dUdguxu(2))*sw
496			t1(3) = t1(3) + (dUdarxu1(3) - dUdguxu(3))*sw
497	gezelter	1689
498	gezelter	1687	t2(1) = t2(1) + (dUdbrxu2(1) + dUdguxu(1))*sw
499			t2(2) = t2(2) + (dUdbrxu2(2) + dUdguxu(2))*sw
500			t2(3) = t2(3) + (dUdbrxu2(3) + dUdguxu(3))*sw
501
502			vpair = vpair + U
503	gezelter	1689
504			!!$ write(,) 'f1 term = ', fxsw, fysw, fz*sw
505			!!$ write(,) 't1 term = ', (dUdarxu1(1) - dUdguxu(1))*sw, &
506			!!$ (dUdarxu1(2) - dUdguxu(2))*sw, &
507			!!$ (dUdarxu1(3) - dUdguxu(3))*sw
508			!!$ write(,) 't2 term = ', (dUdbrxu2(1) + dUdguxu(1))*sw, &
509			!!$ (dUdbrxu2(2) + dUdguxu(2))*sw, &
510			!!$ (dUdbrxu2(3) + dUdguxu(3))*sw
511			!!$ write(,) 'vp term = ', U
512
513	gezelter	115	return
514			end subroutine do_gb_pair
515	gezelter	981
516	gezelter	676	subroutine destroyGBTypes()
517
518			GBMap%nGBtypes = 0
519			GBMap%currentGBtype = 0
520	kdaily	668
521	gezelter	676	if (associated(GBMap%GBtypes)) then
522			deallocate(GBMap%GBtypes)
523			GBMap%GBtypes => null()
524	kdaily	668	end if
525
526	gezelter	676	if (associated(GBMap%atidToGBtype)) then
527			deallocate(GBMap%atidToGBtype)
528			GBMap%atidToGBtype => null()
529			end if
530	kdaily	668
531	gezelter	981	haveMixingMap = .false.
532
533	gezelter	676	end subroutine destroyGBTypes
534	kdaily	668
535	gezelter	676	end module gayberne
536	kdaily	668