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).
!!


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 ) :: eps
     real(kind = dp ) :: eps_ratio 
     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, eps, eps_ratio, dw, status)
    
    integer, intent(in) :: c_ident
    real( kind = dp ), intent(in) :: d, l, eps, eps_ratio, 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)%eps       = eps
    GBMap%GBtypes(current)%eps_ratio = eps_ratio
    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)%eps       = getEpsilon(myATID)
          GBMap%GBtypes(current)%eps_ratio = 1.0_dp
          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, e1, er1, dw1
    real (kind = dp) :: d2, l2, e2, er2, dw2
    real (kind = dp) :: er, ermu, xp, ap2

    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
       e1 = GBMap%GBtypes(i)%eps
       er1 = GBMap%GBtypes(i)%eps_ratio
       dw1 = GBMap%GBtypes(i)%dw

       do j = 1, nGBtypes

          d2 = GBMap%GBtypes(j)%d
          l2 = GBMap%GBtypes(j)%l
          e2 = GBMap%GBtypes(j)%eps
          er2 = GBMap%GBtypes(j)%eps_ratio
          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(e1 * e2)

          er = sqrt(er1 * er2)
          ermu = er**(1.0_dp / mu)
          xp = (1.0_dp - ermu) / (1.0_dp + ermu)
          ap2 = 1.0_dp / (1.0_dp + ermu)

          GBMixingMap(i,j)%xp2 = xp*xp
          GBMixingMap(i,j)%xpap2 = xp*ap2
          GBMixingMap(i,j)%xpapi2 = xp/ap2
       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
    
    ul1(1) = A1(7)
    ul1(2) = A1(8)
    ul1(3) = A1(9)

    ul2(1) = A2(7)
    ul2(2) = A2(8)
    ul2(3) = A2(9)
    
    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)

    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

    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)

            
    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)
    
!!$    write(*,*) 'pref = ' , pref1, pref2
!!$    write(*,*) 'rxu1 = ' , rxu1(1), rxu1(2), rxu1(3)
!!$    write(*,*) 'rxu2 = ' , rxu2(1), rxu2(2), rxu2(3)
!!$    write(*,*) 'uxu = ' , uxu(1), uxu(2), uxu(3)
!!$    write(*,*) 'dUda = ', dUda, dudb, dudg, dudr
!!$    write(*,*) 'H = ', H,hp,sigma, e1, e2, BigR
!!$    write(*,*) 'chi = ', xa2, xai2, x2 
!!$    write(*,*) 'chip = ', xpap2, xpapi2, xp2
!!$    write(*,*) 'eps = ', eps0, e1, e2, eps
!!$    write(*,*) 'U =', U, pref1, pref2
!!$    write(*,*) 'f =', fx, fy, fz
!!$    write(*,*) 'au =', au, bu, g
!!$    

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

    t1(1) = t1(1) + dUda*rxu1(1) - dUdg*uxu(1)
    t1(2) = t1(2) + dUda*rxu1(2) - dUdg*uxu(2)
    t1(3) = t1(3) + dUda*rxu1(3) - dUdg*uxu(3)
                                                    
    t2(1) = t2(1) + dUdb*rxu2(1) + dUdg*uxu(1)
    t2(2) = t2(2) + dUdb*rxu2(2) + dUdg*uxu(2)
    t2(3) = t2(3) + dUdb*rxu2(3) + dUdg*uxu(3)
       
    vpair = vpair + U*sw
    
    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:	1464
Committed:	Fri Jul 9 19:29:05 2010 UTC (15 years ago) by gezelter
File size:	14209 byte(s)
Log Message:	removing cruft (atom numbers, do_pot, do_stress) from many modules and force managers
#	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
43	gezelter	676	module gayberne
44	gezelter	115	use force_globals
45			use definitions
46			use simulation
47	kdaily	668	use atype_module
48			use vector_class
49	gezelter	981	use linearalgebra
50	kdaily	668	use status
51			use lj
52	gezelter	982	use fForceOptions
53	kdaily	529
54	gezelter	115	implicit none
55
56	kdaily	668	private
57	gezelter	115
58	gezelter	676	#define __FORTRAN90
59			#include "UseTheForce/DarkSide/fInteractionMap.h"
60	gezelter	115
61	gezelter	981	logical, save :: haveGBMap = .false.
62			logical, save :: haveMixingMap = .false.
63			real(kind=dp), save :: mu = 2.0_dp
64			real(kind=dp), save :: nu = 1.0_dp
65
66
67	gezelter	676	public :: newGBtype
68	gezelter	981	public :: complete_GB_FF
69	gezelter	115	public :: do_gb_pair
70	chrisfen	578	public :: getGayBerneCut
71	gezelter	676	public :: destroyGBtypes
72	gezelter	115
73	gezelter	676	type :: GBtype
74			integer :: atid
75	gezelter	981	real(kind = dp ) :: d
76			real(kind = dp ) :: l
77	gezelter	676	real(kind = dp ) :: eps
78			real(kind = dp ) :: eps_ratio
79	gezelter	981	real(kind = dp ) :: dw
80			logical :: isLJ
81	gezelter	676	end type GBtype
82	gezelter	981
83	gezelter	676	type, private :: GBList
84			integer :: nGBtypes = 0
85			integer :: currentGBtype = 0
86			type(GBtype), pointer :: GBtypes(:) => null()
87			integer, pointer :: atidToGBtype(:) => null()
88			end type GBList
89	gezelter	981
90	gezelter	676	type(GBList), save :: GBMap
91	gezelter	981
92			type :: GBMixParameters
93			real(kind=DP) :: sigma0
94			real(kind=DP) :: eps0
95			real(kind=DP) :: dw
96			real(kind=DP) :: x2
97			real(kind=DP) :: xa2
98			real(kind=DP) :: xai2
99			real(kind=DP) :: xp2
100			real(kind=DP) :: xpap2
101			real(kind=DP) :: xpapi2
102			end type GBMixParameters
103
104			type(GBMixParameters), dimension(:,:), allocatable :: GBMixingMap
105
106	gezelter	115	contains
107	gezelter	981
108			subroutine newGBtype(c_ident, d, l, eps, eps_ratio, dw, status)
109	gezelter	676
110			integer, intent(in) :: c_ident
111	gezelter	981	real( kind = dp ), intent(in) :: d, l, eps, eps_ratio, dw
112	gezelter	676	integer, intent(out) :: status
113	gezelter	981
114			integer :: nGBTypes, nLJTypes, ntypes, myATID
115	gezelter	676	integer, pointer :: MatchList(:) => null()
116			integer :: current, i
117	kdaily	668	status = 0
118	gezelter	981
119	gezelter	676	if (.not.associated(GBMap%GBtypes)) then
120	gezelter	981
121	gezelter	676	call getMatchingElementList(atypes, "is_GayBerne", .true., &
122			nGBtypes, MatchList)
123
124	gezelter	981	call getMatchingElementList(atypes, "is_LennardJones", .true., &
125			nLJTypes, MatchList)
126
127			GBMap%nGBtypes = nGBtypes + nLJTypes
128
129			allocate(GBMap%GBtypes(nGBtypes + nLJTypes))
130
131	gezelter	676	ntypes = getSize(atypes)
132
133	gezelter	981	allocate(GBMap%atidToGBtype(ntypes))
134			endif
135
136			GBMap%currentGBtype = GBMap%currentGBtype + 1
137			current = GBMap%currentGBtype
138
139			myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
140
141			GBMap%atidToGBtype(myATID) = current
142			GBMap%GBtypes(current)%atid = myATID
143			GBMap%GBtypes(current)%d = d
144			GBMap%GBtypes(current)%l = l
145			GBMap%GBtypes(current)%eps = eps
146			GBMap%GBtypes(current)%eps_ratio = eps_ratio
147			GBMap%GBtypes(current)%dw = dw
148			GBMap%GBtypes(current)%isLJ = .false.
149
150			return
151			end subroutine newGBtype
152
153			subroutine complete_GB_FF(status)
154			integer :: status
155			integer :: i, j, l, m, lm, function_type
156			real(kind=dp) :: thisDP, sigma
157			integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, myATID, current
158			logical :: thisProperty
159
160			status = 0
161			if (GBMap%currentGBtype == 0) then
162			call handleError("complete_GB_FF", "No members in GBMap")
163			status = -1
164			return
165			end if
166
167			nAtypes = getSize(atypes)
168
169			if (nAtypes == 0) then
170			status = -1
171			return
172			end if
173
174			! atypes comes from c side
175			do i = 1, nAtypes
176	gezelter	676
177	gezelter	981	myATID = getFirstMatchingElement(atypes, 'c_ident', i)
178			call getElementProperty(atypes, myATID, "is_LennardJones", thisProperty)
179
180			if (thisProperty) then
181			GBMap%currentGBtype = GBMap%currentGBtype + 1
182			current = GBMap%currentGBtype
183
184			GBMap%atidToGBtype(myATID) = current
185			GBMap%GBtypes(current)%atid = myATID
186			GBMap%GBtypes(current)%isLJ = .true.
187	xsun	1006	GBMap%GBtypes(current)%d = getSigma(myATID) / sqrt(2.0_dp)
188	gezelter	981	GBMap%GBtypes(current)%l = GBMap%GBtypes(current)%d
189			GBMap%GBtypes(current)%eps = getEpsilon(myATID)
190			GBMap%GBtypes(current)%eps_ratio = 1.0_dp
191			GBMap%GBtypes(current)%dw = 1.0_dp
192
193			endif
194
195			end do
196
197			haveGBMap = .true.
198	gezelter	982
199	gezelter	981
200			end subroutine complete_GB_FF
201	kdaily	668
202	gezelter	981	subroutine createGBMixingMap()
203			integer :: nGBtypes, i, j
204			real (kind = dp) :: d1, l1, e1, er1, dw1
205			real (kind = dp) :: d2, l2, e2, er2, dw2
206			real (kind = dp) :: er, ermu, xp, ap2
207	kdaily	668
208	gezelter	981	if (GBMap%currentGBtype == 0) then
209			call handleError("GB", "No members in GBMap")
210			return
211			end if
212
213			nGBtypes = GBMap%nGBtypes
214
215			if (.not. allocated(GBMixingMap)) then
216			allocate(GBMixingMap(nGBtypes, nGBtypes))
217	gezelter	676	endif
218	kdaily	668
219	gezelter	981	do i = 1, nGBtypes
220	kdaily	668
221	gezelter	981	d1 = GBMap%GBtypes(i)%d
222			l1 = GBMap%GBtypes(i)%l
223			e1 = GBMap%GBtypes(i)%eps
224			er1 = GBMap%GBtypes(i)%eps_ratio
225			dw1 = GBMap%GBtypes(i)%dw
226	gezelter	676
227	xsun	1010	do j = 1, nGBtypes
228	gezelter	676
229	gezelter	981	d2 = GBMap%GBtypes(j)%d
230			l2 = GBMap%GBtypes(j)%l
231			e2 = GBMap%GBtypes(j)%eps
232			er2 = GBMap%GBtypes(j)%eps_ratio
233			dw2 = GBMap%GBtypes(j)%dw
234	xsun	1006
235			! Cleaver paper uses sqrt of squares to get sigma0 for
236			! mixed interactions.
237
238			GBMixingMap(i,j)%sigma0 = sqrt(d1d1 + d2d2)
239	gezelter	981	GBMixingMap(i,j)%xa2 = (l1l1 - d1d1)/(l1l1 + d2d2)
240			GBMixingMap(i,j)%xai2 = (l2l2 - d2d2)/(l2l2 + d1d1)
241			GBMixingMap(i,j)%x2 = (l1l1 - d1d1) * (l2l2 - d2d2) / &
242			((l2l2 + d1d1) * (l1l1 + d2d2))
243
244			! assumed LB mixing rules for now:
245
246			GBMixingMap(i,j)%dw = 0.5_dp * (dw1 + dw2)
247			GBMixingMap(i,j)%eps0 = sqrt(e1 * e2)
248
249			er = sqrt(er1 * er2)
250			ermu = er**(1.0_dp / mu)
251			xp = (1.0_dp - ermu) / (1.0_dp + ermu)
252			ap2 = 1.0_dp / (1.0_dp + ermu)
253
254			GBMixingMap(i,j)%xp2 = xp*xp
255			GBMixingMap(i,j)%xpap2 = xp*ap2
256			GBMixingMap(i,j)%xpapi2 = xp/ap2
257			enddo
258			enddo
259			haveMixingMap = .true.
260	gezelter	983	mu = getGayBerneMu()
261			nu = getGayBerneNu()
262	gezelter	981	end subroutine createGBMixingMap
263	gezelter	676
264	gezelter	981
265	chrisfen	580	!! gay berne cutoff should be a parameter in globals, this is a temporary
266			!! work around - this should be fixed when gay berne is up and running
267	gezelter	676
268	chrisfen	578	function getGayBerneCut(atomID) result(cutValue)
269	gezelter	676	integer, intent(in) :: atomID
270			integer :: gbt1
271	gezelter	981	real(kind=dp) :: cutValue, l, d
272	gezelter	115
273	gezelter	676	if (GBMap%currentGBtype == 0) then
274			call handleError("GB", "No members in GBMap")
275			return
276			end if
277
278			gbt1 = GBMap%atidToGBtype(atomID)
279	gezelter	981	l = GBMap%GBtypes(gbt1)%l
280			d = GBMap%GBtypes(gbt1)%d
281	gezelter	676
282	xsun	1010	! sigma is actually sqrt(2)*l for prolate ellipsoids
283
284			cutValue = 2.5_dpsqrt(2.0_dp)max(l,d)
285
286	chrisfen	578	end function getGayBerneCut
287
288	gezelter	1464	subroutine do_gb_pair(atid1, atid2, d, r, r2, sw, vdwMult, vpair, fpair, &
289			pot, A1, A2, f1, t1, t2)
290	kdaily	529
291	gezelter	1464	integer, intent(in) :: atid1, atid2
292	gezelter	1386	integer :: gbt1, gbt2, id1, id2
293	gezelter	1286	real (kind=dp), intent(inout) :: r, r2, vdwMult
294	gezelter	115	real (kind=dp), dimension(3), intent(in) :: d
295			real (kind=dp), dimension(3), intent(inout) :: fpair
296			real (kind=dp) :: pot, sw, vpair
297	gezelter	1386	real (kind=dp), dimension(9) :: A1, A2
298			real (kind=dp), dimension(3) :: f1
299			real (kind=dp), dimension(3) :: t1, t2
300	gezelter	981	real (kind = dp), dimension(3) :: ul1, ul2, rxu1, rxu2, uxu, rhat
301	gezelter	115
302	gezelter	981	real (kind = dp) :: sigma0, dw, eps0, x2, xa2, xai2, xp2, xpap2, xpapi2
303	kdaily	997	real (kind = dp) :: e1, e2, eps, sigma, s3, s03, au2, bu2, au, bu, a, b, g, g2
304	gezelter	981	real (kind = dp) :: U, BigR, R3, R6, R7, R12, R13, H, Hp, fx, fy, fz
305			real (kind = dp) :: dUdr, dUda, dUdb, dUdg, pref1, pref2
306			logical :: i_is_lj, j_is_lj
307
308			if (.not.haveMixingMap) then
309			call createGBMixingMap()
310			endif
311
312	gezelter	676	gbt1 = GBMap%atidToGBtype(atid1)
313	gezelter	981	gbt2 = GBMap%atidToGBtype(atid2)
314	gezelter	676
315	gezelter	981	i_is_LJ = GBMap%GBTypes(gbt1)%isLJ
316			j_is_LJ = GBMap%GBTypes(gbt2)%isLJ
317	gezelter	676
318	gezelter	981	sigma0 = GBMixingMap(gbt1, gbt2)%sigma0
319			dw = GBMixingMap(gbt1, gbt2)%dw
320			eps0 = GBMixingMap(gbt1, gbt2)%eps0
321			x2 = GBMixingMap(gbt1, gbt2)%x2
322			xa2 = GBMixingMap(gbt1, gbt2)%xa2
323			xai2 = GBMixingMap(gbt1, gbt2)%xai2
324			xp2 = GBMixingMap(gbt1, gbt2)%xp2
325			xpap2 = GBMixingMap(gbt1, gbt2)%xpap2
326			xpapi2 = GBMixingMap(gbt1, gbt2)%xpapi2
327
328	gezelter	1386	ul1(1) = A1(7)
329			ul1(2) = A1(8)
330			ul1(3) = A1(9)
331	gezelter	115
332	gezelter	1386	ul2(1) = A2(7)
333			ul2(2) = A2(8)
334			ul2(3) = A2(9)
335	kdaily	529
336	gezelter	981	if (i_is_LJ) then
337			a = 0.0_dp
338			ul1 = 0.0_dp
339			else
340			a = d(1)ul1(1) + d(2)ul1(2) + d(3)*ul1(3)
341			endif
342	gezelter	115
343	gezelter	981	if (j_is_LJ) then
344			b = 0.0_dp
345			ul2 = 0.0_dp
346			else
347			b = d(1)ul2(1) + d(2)ul2(2) + d(3)*ul2(3)
348			endif
349	gezelter	115
350	gezelter	981	if (i_is_LJ.or.j_is_LJ) then
351			g = 0.0_dp
352			else
353			g = ul1(1)ul2(1) + ul1(2)ul2(2) + ul1(3)*ul2(3)
354			endif
355	gezelter	1002
356	gezelter	981	au = a / r
357			bu = b / r
358	kdaily	997
359	gezelter	1002	au2 = au * au
360			bu2 = bu * bu
361			g2 = g * g
362	gezelter	115
363	kdaily	997	H = (xa2 * au2 + xai2 * bu2 - 2.0_dpx2aubug) / (1.0_dp - x2*g2)
364			Hp = (xpap2au2 + xpapi2bu2 - 2.0_dpxp2aubug) / (1.0_dp - xp2*g2)
365	gezelter	1002
366	gezelter	981	sigma = sigma0 / sqrt(1.0_dp - H)
367			e1 = 1.0_dp / sqrt(1.0_dp - x2*g2)
368			e2 = 1.0_dp - Hp
369			eps = eps0 * (e1*nu) (e2**mu)
370			BigR = dwsigma0 / (r - sigma + dwsigma0)
371
372			R3 = BigRBigRBigR
373			R6 = R3*R3
374			R7 = R6 * BigR
375			R12 = R6*R6
376			R13 = R6*R7
377	gezelter	115
378	gezelter	1286	U = vdwMult * 4.0_dp * eps * (R12 - R6)
379	gezelter	115
380	gezelter	981	s3 = sigmasigmasigma
381			s03 = sigma0sigma0sigma0
382	kdaily	529
383	gezelter	1286	pref1 = - vdwMult * 8.0_dp * eps * mu * (R12 - R6) / (e2 * r)
384	gezelter	983
385	gezelter	1286	pref2 = vdwMult * 8.0_dp * eps * s3 * (6.0_dpR13 - 3.0_dpR7) / (dwrs03)
386	gezelter	507
387	kdaily	997	dUdr = - (pref1 * Hp + pref2 * (sigma0sigma0r/s3 + H))
388	kdaily	529
389	gezelter	981	dUda = pref1 * (xpap2au - xp2bug) / (1.0_dp - xp2 g2) &
390			+ pref2 * (xa2 * au - x2 bug) / (1.0_dp - x2 * g2)
391	kdaily	997
392	gezelter	981	dUdb = pref1 * (xpapi2bu - xp2aug) / (1.0_dp - xp2 g2) &
393			+ pref2 * (xai2 * bu - x2 aug) / (1.0_dp - x2 * g2)
394	gezelter	507
395	gezelter	981	dUdg = 4.0_dp * eps * nu * (R12 - R6) * x2 * g / (1.0_dp - x2*g2) &
396			+ 8.0_dp * eps * mu * (R12 - R6) * (xp2aubu - Hpxp2g) / &
397			(1.0_dp - xp2 * g2) / e2 &
398	gezelter	983	+ 8.0_dp * eps * s3 * (3.0_dp * R7 - 6.0_dp * R13) * &
399	gezelter	981	(x2 * au * bu - H * x2 * g) / (1.0_dp - x2 * g2) / (dw * s03)
400	kdaily	997
401	gezelter	981
402			rhat = d / r
403	gezelter	507
404	gezelter	983	fx = dUdr * rhat(1) + dUda * ul1(1) + dUdb * ul2(1)
405			fy = dUdr * rhat(2) + dUda * ul1(2) + dUdb * ul2(2)
406			fz = dUdr * rhat(3) + dUda * ul1(3) + dUdb * ul2(3)
407	gezelter	115
408	gezelter	981	rxu1 = cross_product(d, ul1)
409			rxu2 = cross_product(d, ul2)
410			uxu = cross_product(ul1, ul2)
411	kdaily	997
412			!!$ write(,) 'pref = ' , pref1, pref2
413	gezelter	983	!!$ write(,) 'rxu1 = ' , rxu1(1), rxu1(2), rxu1(3)
414			!!$ write(,) 'rxu2 = ' , rxu2(1), rxu2(2), rxu2(3)
415			!!$ write(,) 'uxu = ' , uxu(1), uxu(2), uxu(3)
416	kdaily	997	!!$ write(,) 'dUda = ', dUda, dudb, dudg, dudr
417			!!$ write(,) 'H = ', H,hp,sigma, e1, e2, BigR
418			!!$ write(,) 'chi = ', xa2, xai2, x2
419			!!$ write(,) 'chip = ', xpap2, xpapi2, xp2
420			!!$ write(,) 'eps = ', eps0, e1, e2, eps
421			!!$ write(,) 'U =', U, pref1, pref2
422			!!$ write(,) 'f =', fx, fy, fz
423			!!$ write(,) 'au =', au, bu, g
424			!!$
425	gezelter	1386
426			pot = pot + U*sw
427	kdaily	529
428	gezelter	1386	f1(1) = f1(1) + fx
429			f1(2) = f1(2) + fy
430			f1(3) = f1(3) + fz
431
432			t1(1) = t1(1) + dUdarxu1(1) - dUdguxu(1)
433			t1(2) = t1(2) + dUdarxu1(2) - dUdguxu(2)
434			t1(3) = t1(3) + dUdarxu1(3) - dUdguxu(3)
435
436			t2(1) = t2(1) + dUdbrxu2(1) + dUdguxu(1)
437			t2(2) = t2(2) + dUdbrxu2(2) + dUdguxu(2)
438			t2(3) = t2(3) + dUdbrxu2(3) + dUdguxu(3)
439
440	gezelter	981	vpair = vpair + U*sw
441	kdaily	529
442	gezelter	115	return
443			end subroutine do_gb_pair
444	gezelter	981
445	gezelter	676	subroutine destroyGBTypes()
446
447			GBMap%nGBtypes = 0
448			GBMap%currentGBtype = 0
449	kdaily	668
450	gezelter	676	if (associated(GBMap%GBtypes)) then
451			deallocate(GBMap%GBtypes)
452			GBMap%GBtypes => null()
453	kdaily	668	end if
454
455	gezelter	676	if (associated(GBMap%atidToGBtype)) then
456			deallocate(GBMap%atidToGBtype)
457			GBMap%atidToGBtype => null()
458			end if
459	kdaily	668
460	gezelter	981	haveMixingMap = .false.
461
462	gezelter	676	end subroutine destroyGBTypes
463	kdaily	668
464	gezelter	676	end module gayberne
465	kdaily	668