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(atom1, atom2, atid1, atid2, d, r, r2, sw, vdwMult, vpair, fpair, &
       pot, A1, A2, f1, t1, t2, do_pot)
    
    integer, intent(in) :: atom1, atom2, 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
    logical, intent(in) :: do_pot
    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:	1390
Committed:	Wed Nov 25 20:02:06 2009 UTC (15 years, 8 months ago) by gezelter
File size:	14279 byte(s)
Log Message:	Almost all of the changes necessary to create OpenMD out of our old project (OOPSE-4)
#	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	1386	subroutine do_gb_pair(atom1, atom2, atid1, atid2, d, r, r2, sw, vdwMult, vpair, fpair, &
289			pot, A1, A2, f1, t1, t2, do_pot)
290	kdaily	529
291	gezelter	1386	integer, intent(in) :: atom1, atom2, atid1, atid2
292			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	115	logical, intent(in) :: do_pot
301	gezelter	981	real (kind = dp), dimension(3) :: ul1, ul2, rxu1, rxu2, uxu, rhat
302	gezelter	115
303	gezelter	981	real (kind = dp) :: sigma0, dw, eps0, x2, xa2, xai2, xp2, xpap2, xpapi2
304	kdaily	997	real (kind = dp) :: e1, e2, eps, sigma, s3, s03, au2, bu2, au, bu, a, b, g, g2
305	gezelter	981	real (kind = dp) :: U, BigR, R3, R6, R7, R12, R13, H, Hp, fx, fy, fz
306			real (kind = dp) :: dUdr, dUda, dUdb, dUdg, pref1, pref2
307			logical :: i_is_lj, j_is_lj
308
309			if (.not.haveMixingMap) then
310			call createGBMixingMap()
311			endif
312
313	gezelter	676	gbt1 = GBMap%atidToGBtype(atid1)
314	gezelter	981	gbt2 = GBMap%atidToGBtype(atid2)
315	gezelter	676
316	gezelter	981	i_is_LJ = GBMap%GBTypes(gbt1)%isLJ
317			j_is_LJ = GBMap%GBTypes(gbt2)%isLJ
318	gezelter	676
319	gezelter	981	sigma0 = GBMixingMap(gbt1, gbt2)%sigma0
320			dw = GBMixingMap(gbt1, gbt2)%dw
321			eps0 = GBMixingMap(gbt1, gbt2)%eps0
322			x2 = GBMixingMap(gbt1, gbt2)%x2
323			xa2 = GBMixingMap(gbt1, gbt2)%xa2
324			xai2 = GBMixingMap(gbt1, gbt2)%xai2
325			xp2 = GBMixingMap(gbt1, gbt2)%xp2
326			xpap2 = GBMixingMap(gbt1, gbt2)%xpap2
327			xpapi2 = GBMixingMap(gbt1, gbt2)%xpapi2
328
329	gezelter	1386	ul1(1) = A1(7)
330			ul1(2) = A1(8)
331			ul1(3) = A1(9)
332	gezelter	115
333	gezelter	1386	ul2(1) = A2(7)
334			ul2(2) = A2(8)
335			ul2(3) = A2(9)
336	kdaily	529
337	gezelter	981	if (i_is_LJ) then
338			a = 0.0_dp
339			ul1 = 0.0_dp
340			else
341			a = d(1)ul1(1) + d(2)ul1(2) + d(3)*ul1(3)
342			endif
343	gezelter	115
344	gezelter	981	if (j_is_LJ) then
345			b = 0.0_dp
346			ul2 = 0.0_dp
347			else
348			b = d(1)ul2(1) + d(2)ul2(2) + d(3)*ul2(3)
349			endif
350	gezelter	115
351	gezelter	981	if (i_is_LJ.or.j_is_LJ) then
352			g = 0.0_dp
353			else
354			g = ul1(1)ul2(1) + ul1(2)ul2(2) + ul1(3)*ul2(3)
355			endif
356	gezelter	1002
357	gezelter	981	au = a / r
358			bu = b / r
359	kdaily	997
360	gezelter	1002	au2 = au * au
361			bu2 = bu * bu
362			g2 = g * g
363	gezelter	115
364	kdaily	997	H = (xa2 * au2 + xai2 * bu2 - 2.0_dpx2aubug) / (1.0_dp - x2*g2)
365			Hp = (xpap2au2 + xpapi2bu2 - 2.0_dpxp2aubug) / (1.0_dp - xp2*g2)
366	gezelter	1002
367	gezelter	981	sigma = sigma0 / sqrt(1.0_dp - H)
368			e1 = 1.0_dp / sqrt(1.0_dp - x2*g2)
369			e2 = 1.0_dp - Hp
370			eps = eps0 * (e1*nu) (e2**mu)
371			BigR = dwsigma0 / (r - sigma + dwsigma0)
372
373			R3 = BigRBigRBigR
374			R6 = R3*R3
375			R7 = R6 * BigR
376			R12 = R6*R6
377			R13 = R6*R7
378	gezelter	115
379	gezelter	1286	U = vdwMult * 4.0_dp * eps * (R12 - R6)
380	gezelter	115
381	gezelter	981	s3 = sigmasigmasigma
382			s03 = sigma0sigma0sigma0
383	kdaily	529
384	gezelter	1286	pref1 = - vdwMult * 8.0_dp * eps * mu * (R12 - R6) / (e2 * r)
385	gezelter	983
386	gezelter	1286	pref2 = vdwMult * 8.0_dp * eps * s3 * (6.0_dpR13 - 3.0_dpR7) / (dwrs03)
387	gezelter	507
388	kdaily	997	dUdr = - (pref1 * Hp + pref2 * (sigma0sigma0r/s3 + H))
389	kdaily	529
390	gezelter	981	dUda = pref1 * (xpap2au - xp2bug) / (1.0_dp - xp2 g2) &
391			+ pref2 * (xa2 * au - x2 bug) / (1.0_dp - x2 * g2)
392	kdaily	997
393	gezelter	981	dUdb = pref1 * (xpapi2bu - xp2aug) / (1.0_dp - xp2 g2) &
394			+ pref2 * (xai2 * bu - x2 aug) / (1.0_dp - x2 * g2)
395	gezelter	507
396	gezelter	981	dUdg = 4.0_dp * eps * nu * (R12 - R6) * x2 * g / (1.0_dp - x2*g2) &
397			+ 8.0_dp * eps * mu * (R12 - R6) * (xp2aubu - Hpxp2g) / &
398			(1.0_dp - xp2 * g2) / e2 &
399	gezelter	983	+ 8.0_dp * eps * s3 * (3.0_dp * R7 - 6.0_dp * R13) * &
400	gezelter	981	(x2 * au * bu - H * x2 * g) / (1.0_dp - x2 * g2) / (dw * s03)
401	kdaily	997
402	gezelter	981
403			rhat = d / r
404	gezelter	507
405	gezelter	983	fx = dUdr * rhat(1) + dUda * ul1(1) + dUdb * ul2(1)
406			fy = dUdr * rhat(2) + dUda * ul1(2) + dUdb * ul2(2)
407			fz = dUdr * rhat(3) + dUda * ul1(3) + dUdb * ul2(3)
408	gezelter	115
409	gezelter	981	rxu1 = cross_product(d, ul1)
410			rxu2 = cross_product(d, ul2)
411			uxu = cross_product(ul1, ul2)
412	kdaily	997
413			!!$ write(,) 'pref = ' , pref1, pref2
414	gezelter	983	!!$ write(,) 'rxu1 = ' , rxu1(1), rxu1(2), rxu1(3)
415			!!$ write(,) 'rxu2 = ' , rxu2(1), rxu2(2), rxu2(3)
416			!!$ write(,) 'uxu = ' , uxu(1), uxu(2), uxu(3)
417	kdaily	997	!!$ write(,) 'dUda = ', dUda, dudb, dudg, dudr
418			!!$ write(,) 'H = ', H,hp,sigma, e1, e2, BigR
419			!!$ write(,) 'chi = ', xa2, xai2, x2
420			!!$ write(,) 'chip = ', xpap2, xpapi2, xp2
421			!!$ write(,) 'eps = ', eps0, e1, e2, eps
422			!!$ write(,) 'U =', U, pref1, pref2
423			!!$ write(,) 'f =', fx, fy, fz
424			!!$ write(,) 'au =', au, bu, g
425			!!$
426	gezelter	1386
427			pot = pot + U*sw
428	kdaily	529
429	gezelter	1386	f1(1) = f1(1) + fx
430			f1(2) = f1(2) + fy
431			f1(3) = f1(3) + fz
432
433			t1(1) = t1(1) + dUdarxu1(1) - dUdguxu(1)
434			t1(2) = t1(2) + dUdarxu1(2) - dUdguxu(2)
435			t1(3) = t1(3) + dUdarxu1(3) - dUdguxu(3)
436
437			t2(1) = t2(1) + dUdbrxu2(1) + dUdguxu(1)
438			t2(2) = t2(2) + dUdbrxu2(2) + dUdguxu(2)
439			t2(3) = t2(3) + dUdbrxu2(3) + dUdguxu(3)
440
441	gezelter	981	vpair = vpair + U*sw
442	kdaily	529
443	gezelter	115	return
444			end subroutine do_gb_pair
445	gezelter	981
446	gezelter	676	subroutine destroyGBTypes()
447
448			GBMap%nGBtypes = 0
449			GBMap%currentGBtype = 0
450	kdaily	668
451	gezelter	676	if (associated(GBMap%GBtypes)) then
452			deallocate(GBMap%GBtypes)
453			GBMap%GBtypes => null()
454	kdaily	668	end if
455
456	gezelter	676	if (associated(GBMap%atidToGBtype)) then
457			deallocate(GBMap%atidToGBtype)
458			GBMap%atidToGBtype => null()
459			end if
460	kdaily	668
461	gezelter	981	haveMixingMap = .false.
462
463	gezelter	676	end subroutine destroyGBTypes
464	kdaily	668
465	gezelter	676	end module gayberne
466	kdaily	668