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. Acknowledgement of the program authors must be made in any
!!    publication of scientific results based in part on use of the
!!    program.  An acceptable form of acknowledgement is citation of
!!    the article in which the program was described (Matthew
!!    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
!!    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
!!    Parallel Simulation Engine for Molecular Dynamics,"
!!    J. Comput. Chem. 26, pp. 252-271 (2005))
!!
!! 2. Redistributions of source code must retain the above copyright
!!    notice, this list of conditions and the following disclaimer.
!!
!! 3. 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.
!!


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:	3559
Committed:	Fri Oct 23 18:41:09 2009 UTC (14 years, 10 months ago) by gezelter
File size:	14262 byte(s)
Log Message:	removing MPI responsibilities from the lowest level force routines. This is in preparation for migrating these routines (LJ, electrostatic, eam, suttonchen, gay-berne, sticky) to C++
#	User	Rev	Content
1	gezelter	1930	!!
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			!! 1. Acknowledgement of the program authors must be made in any
10			!! publication of scientific results based in part on use of the
11			!! program. An acceptable form of acknowledgement is citation of
12			!! the article in which the program was described (Matthew
13			!! A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14			!! J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15			!! Parallel Simulation Engine for Molecular Dynamics,"
16			!! J. Comput. Chem. 26, pp. 252-271 (2005))
17			!!
18			!! 2. Redistributions of source code must retain the above copyright
19			!! notice, this list of conditions and the following disclaimer.
20			!!
21			!! 3. Redistributions in binary form must reproduce the above copyright
22			!! notice, this list of conditions and the following disclaimer in the
23			!! documentation and/or other materials provided with the
24			!! distribution.
25			!!
26			!! This software is provided "AS IS," without a warranty of any
27			!! kind. All express or implied conditions, representations and
28			!! warranties, including any implied warranty of merchantability,
29			!! fitness for a particular purpose or non-infringement, are hereby
30			!! excluded. The University of Notre Dame and its licensors shall not
31			!! be liable for any damages suffered by licensee as a result of
32			!! using, modifying or distributing the software or its
33			!! derivatives. In no event will the University of Notre Dame or its
34			!! licensors be liable for any lost revenue, profit or data, or for
35			!! direct, indirect, special, consequential, incidental or punitive
36			!! damages, however caused and regardless of the theory of liability,
37			!! arising out of the use of or inability to use software, even if the
38			!! University of Notre Dame has been advised of the possibility of
39			!! such damages.
40			!!
41
42
43	gezelter	2375	module gayberne
44	gezelter	1608	use force_globals
45			use definitions
46			use simulation
47	kdaily	2367	use atype_module
48			use vector_class
49	gezelter	2787	use linearalgebra
50	kdaily	2367	use status
51			use lj
52	gezelter	2788	use fForceOptions
53	kdaily	2226
54	gezelter	1608	implicit none
55
56	kdaily	2367	private
57	gezelter	1608
58	gezelter	2375	#define __FORTRAN90
59			#include "UseTheForce/DarkSide/fInteractionMap.h"
60	gezelter	1608
61	gezelter	2787	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	2375	public :: newGBtype
68	gezelter	2787	public :: complete_GB_FF
69	gezelter	1608	public :: do_gb_pair
70	chrisfen	2277	public :: getGayBerneCut
71	gezelter	2375	public :: destroyGBtypes
72	gezelter	1608
73	gezelter	2375	type :: GBtype
74			integer :: atid
75	gezelter	2787	real(kind = dp ) :: d
76			real(kind = dp ) :: l
77	gezelter	2375	real(kind = dp ) :: eps
78			real(kind = dp ) :: eps_ratio
79	gezelter	2787	real(kind = dp ) :: dw
80			logical :: isLJ
81	gezelter	2375	end type GBtype
82	gezelter	2787
83	gezelter	2375	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	2787
90	gezelter	2375	type(GBList), save :: GBMap
91	gezelter	2787
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	1608	contains
107	gezelter	2787
108			subroutine newGBtype(c_ident, d, l, eps, eps_ratio, dw, status)
109	gezelter	2375
110			integer, intent(in) :: c_ident
111	gezelter	2787	real( kind = dp ), intent(in) :: d, l, eps, eps_ratio, dw
112	gezelter	2375	integer, intent(out) :: status
113	gezelter	2787
114			integer :: nGBTypes, nLJTypes, ntypes, myATID
115	gezelter	2375	integer, pointer :: MatchList(:) => null()
116			integer :: current, i
117	kdaily	2367	status = 0
118	gezelter	2787
119	gezelter	2375	if (.not.associated(GBMap%GBtypes)) then
120	gezelter	2787
121	gezelter	2375	call getMatchingElementList(atypes, "is_GayBerne", .true., &
122			nGBtypes, MatchList)
123
124	gezelter	2787	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	2375	ntypes = getSize(atypes)
132
133	gezelter	2787	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	2375
177	gezelter	2787	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	2940	GBMap%GBtypes(current)%d = getSigma(myATID) / sqrt(2.0_dp)
188	gezelter	2787	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	2788
199	gezelter	2787
200			end subroutine complete_GB_FF
201	kdaily	2367
202	gezelter	2787	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	2367
208	gezelter	2787	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	2375	endif
218	kdaily	2367
219	gezelter	2787	do i = 1, nGBtypes
220	kdaily	2367
221	gezelter	2787	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	2375
227	xsun	2958	do j = 1, nGBtypes
228	gezelter	2375
229	gezelter	2787	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	2940
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	2787	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	2802	mu = getGayBerneMu()
261			nu = getGayBerneNu()
262	gezelter	2787	end subroutine createGBMixingMap
263	gezelter	2375
264	gezelter	2787
265	chrisfen	2279	!! 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	2375
268	chrisfen	2277	function getGayBerneCut(atomID) result(cutValue)
269	gezelter	2375	integer, intent(in) :: atomID
270			integer :: gbt1
271	gezelter	2787	real(kind=dp) :: cutValue, l, d
272	gezelter	1608
273	gezelter	2375	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	2787	l = GBMap%GBtypes(gbt1)%l
280			d = GBMap%GBtypes(gbt1)%d
281	gezelter	2375
282	xsun	2958	! sigma is actually sqrt(2)*l for prolate ellipsoids
283
284			cutValue = 2.5_dpsqrt(2.0_dp)max(l,d)
285
286	chrisfen	2277	end function getGayBerneCut
287
288	gezelter	3559	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	2226
291	gezelter	3559	integer, intent(in) :: atom1, atom2, atid1, atid2
292			integer :: gbt1, gbt2, id1, id2
293	gezelter	3441	real (kind=dp), intent(inout) :: r, r2, vdwMult
294	gezelter	1608	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	3559	real (kind=dp), dimension(9) :: A1, A2
298			real (kind=dp), dimension(3) :: f1
299			real (kind=dp), dimension(3) :: t1, t2
300	gezelter	1608	logical, intent(in) :: do_pot
301	gezelter	2787	real (kind = dp), dimension(3) :: ul1, ul2, rxu1, rxu2, uxu, rhat
302	gezelter	1608
303	gezelter	2787	real (kind = dp) :: sigma0, dw, eps0, x2, xa2, xai2, xp2, xpap2, xpapi2
304	kdaily	2915	real (kind = dp) :: e1, e2, eps, sigma, s3, s03, au2, bu2, au, bu, a, b, g, g2
305	gezelter	2787	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	2375	gbt1 = GBMap%atidToGBtype(atid1)
314	gezelter	2787	gbt2 = GBMap%atidToGBtype(atid2)
315	gezelter	2375
316	gezelter	2787	i_is_LJ = GBMap%GBTypes(gbt1)%isLJ
317			j_is_LJ = GBMap%GBTypes(gbt2)%isLJ
318	gezelter	2375
319	gezelter	2787	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	3559	ul1(1) = A1(7)
330			ul1(2) = A1(8)
331			ul1(3) = A1(9)
332	gezelter	1608
333	gezelter	3559	ul2(1) = A2(7)
334			ul2(2) = A2(8)
335			ul2(3) = A2(9)
336	kdaily	2226
337	gezelter	2787	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	1608
344	gezelter	2787	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	1608
351	gezelter	2787	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	2926
357	gezelter	2787	au = a / r
358			bu = b / r
359	kdaily	2915
360	gezelter	2926	au2 = au * au
361			bu2 = bu * bu
362			g2 = g * g
363	gezelter	1608
364	kdaily	2915	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	2926
367	gezelter	2787	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	1608
379	gezelter	3441	U = vdwMult * 4.0_dp * eps * (R12 - R6)
380	gezelter	1608
381	gezelter	2787	s3 = sigmasigmasigma
382			s03 = sigma0sigma0sigma0
383	kdaily	2226
384	gezelter	3441	pref1 = - vdwMult * 8.0_dp * eps * mu * (R12 - R6) / (e2 * r)
385	gezelter	2802
386	gezelter	3441	pref2 = vdwMult * 8.0_dp * eps * s3 * (6.0_dpR13 - 3.0_dpR7) / (dwrs03)
387	gezelter	2204
388	kdaily	2915	dUdr = - (pref1 * Hp + pref2 * (sigma0sigma0r/s3 + H))
389	kdaily	2226
390	gezelter	2787	dUda = pref1 * (xpap2au - xp2bug) / (1.0_dp - xp2 g2) &
391			+ pref2 * (xa2 * au - x2 bug) / (1.0_dp - x2 * g2)
392	kdaily	2915
393	gezelter	2787	dUdb = pref1 * (xpapi2bu - xp2aug) / (1.0_dp - xp2 g2) &
394			+ pref2 * (xai2 * bu - x2 aug) / (1.0_dp - x2 * g2)
395	gezelter	2204
396	gezelter	2787	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	2802	+ 8.0_dp * eps * s3 * (3.0_dp * R7 - 6.0_dp * R13) * &
400	gezelter	2787	(x2 * au * bu - H * x2 * g) / (1.0_dp - x2 * g2) / (dw * s03)
401	kdaily	2915
402	gezelter	2787
403			rhat = d / r
404	gezelter	2204
405	gezelter	2802	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	1608
409	gezelter	2787	rxu1 = cross_product(d, ul1)
410			rxu2 = cross_product(d, ul2)
411			uxu = cross_product(ul1, ul2)
412	kdaily	2915
413			!!$ write(,) 'pref = ' , pref1, pref2
414	gezelter	2802	!!$ 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	2915	!!$ 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	3559
427			pot = pot + U*sw
428	kdaily	2226
429	gezelter	3559	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	2787	vpair = vpair + U*sw
442	kdaily	2226
443	gezelter	1608	return
444			end subroutine do_gb_pair
445	gezelter	2787
446	gezelter	2375	subroutine destroyGBTypes()
447
448			GBMap%nGBtypes = 0
449			GBMap%currentGBtype = 0
450	kdaily	2367
451	gezelter	2375	if (associated(GBMap%GBtypes)) then
452			deallocate(GBMap%GBtypes)
453			GBMap%GBtypes => null()
454	kdaily	2367	end if
455
456	gezelter	2375	if (associated(GBMap%atidToGBtype)) then
457			deallocate(GBMap%atidToGBtype)
458			GBMap%atidToGBtype => null()
459			end if
460	kdaily	2367
461	gezelter	2787	haveMixingMap = .false.
462
463	gezelter	2375	end subroutine destroyGBTypes
464	kdaily	2367
465	gezelter	2375	end module gayberne
466	kdaily	2367