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
#ifdef IS_MPI
  use mpiSimulation
#endif
  
  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)
          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.

    mu = getGayBerneMu()
    nu = getGayBerneNu()

    
  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 = i, 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

          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
          
          if (i.ne.j) then
             GBMixingMap(j,i)%sigma0 = GBMixingMap(i,j)%sigma0
             GBMixingMap(j,i)%dw     = GBMixingMap(i,j)%dw    
             GBMixingMap(j,i)%eps0   = GBMixingMap(i,j)%eps0  
             GBMixingMap(j,i)%x2     = GBMixingMap(i,j)%x2    
             GBMixingMap(j,i)%xa2    = GBMixingMap(i,j)%xa2   
             GBMixingMap(j,i)%xai2   = GBMixingMap(i,j)%xai2  
             GBMixingMap(j,i)%xp2    = GBMixingMap(i,j)%xp2   
             GBMixingMap(j,i)%xpap2  = GBMixingMap(i,j)%xpap2 
             GBMixingMap(j,i)%xpapi2 = GBMixingMap(i,j)%xpapi2
          endif
       enddo
    enddo
    haveMixingMap = .true.
    
  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   
    cutValue = 2.5_dp*max(l,d)

  end function getGayBerneCut

  subroutine do_gb_pair(atom1, atom2, d, r, r2, sw, vpair, fpair, &
       pot, Amat, f, t, do_pot)
    
    integer, intent(in) :: atom1, atom2
    integer :: atid1, atid2, gbt1, gbt2, id1, id2
    real (kind=dp), intent(inout) :: r, r2
    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,nLocal) :: Amat
    real (kind=dp), dimension(3,nLocal) :: f
    real (kind=dp), dimension(3,nLocal) :: t
    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, 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

#ifdef IS_MPI
    atid1 = atid_Row(atom1)
    atid2 = atid_Col(atom2)
#else
    atid1 = atid(atom1)
    atid2 = atid(atom2)
#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
    
#ifdef IS_MPI
    ul1(1) = A_Row(7,atom1)
    ul1(2) = A_Row(8,atom1)
    ul1(3) = A_Row(9,atom1)

    ul2(1) = A_Col(7,atom2)
    ul2(2) = A_Col(8,atom2)
    ul2(3) = A_Col(9,atom2)
#else
    ul1(1) = Amat(7,atom1)
    ul1(2) = Amat(8,atom1)
    ul1(3) = Amat(9,atom1)

    ul2(1) = Amat(7,atom2)
    ul2(2) = Amat(8,atom2)
    ul2(3) = Amat(9,atom2)
#endif
    
    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
    g2 = g*g

    H  = (xa2 * au + xai2 * bu - 2.0_dp*x2*au*bu*g)  / (1.0_dp - x2*g2)
    Hp = (xpap2*au + xpapi2*bu - 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 = 4.0_dp * eps * (R12 - R6)

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

    pref1 = - 8.0_dp * eps * mu * (R12 - R6) / (e2 * r)
    pref2 = 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)
    fx = -dUdr * rhat(3) - dUda * ul1(3) - dUdb * ul2(3)

    rxu1 = cross_product(d, ul1)
    rxu2 = cross_product(d, ul2)    
    uxu = cross_product(ul1, ul2)
           
#ifdef IS_MPI
    f_Row(1,atom1) = f_Row(1,atom1) + fx
    f_Row(2,atom1) = f_Row(2,atom1) + fy
    f_Row(3,atom1) = f_Row(3,atom1) + fz
    
    f_Col(1,atom2) = f_Col(1,atom2) - fx
    f_Col(2,atom2) = f_Col(2,atom2) - fy
    f_Col(3,atom2) = f_Col(3,atom2) - fz
    
    t_Row(1,atom1) = t_Row(1,atom1) + dUda*rxu1(1) - dUdg*uxu(1)
    t_Row(2,atom1) = t_Row(2,atom1) + dUda*rxu1(2) - dUdg*uxu(2)
    t_Row(3,atom1) = t_Row(3,atom1) + dUda*rxu1(3) - dUdg*uxu(3)
                                                                
    t_Col(1,atom2) = t_Col(1,atom2) + dUdb*rxu2(1) + dUdg*uxu(1)
    t_Col(2,atom2) = t_Col(2,atom2) + dUdb*rxu2(2) + dUdg*uxu(2)
    t_Col(3,atom2) = t_Col(3,atom2) + dUdb*rxu2(3) + dUdg*uxu(3)
#else
    f(1,atom1) = f(1,atom1) + fx
    f(2,atom1) = f(2,atom1) + fy
    f(3,atom1) = f(3,atom1) + fz
    
    f(1,atom2) = f(1,atom2) - fx
    f(2,atom2) = f(2,atom2) - fy
    f(3,atom2) = f(3,atom2) - fz
    
    t(1,atom1) = t(1,atom1) +  dUda*rxu1(1) - dUdg*uxu(1)
    t(2,atom1) = t(2,atom1) +  dUda*rxu1(2) - dUdg*uxu(2)
    t(3,atom1) = t(3,atom1) +  dUda*rxu1(3) - dUdg*uxu(3)
                                                         
    t(1,atom2) = t(1,atom2) +  dUdb*rxu2(1) + dUdg*uxu(1)
    t(2,atom2) = t(2,atom2) +  dUdb*rxu2(2) + dUdg*uxu(2)
    t(3,atom2) = t(3,atom2) +  dUdb*rxu2(3) + dUdg*uxu(3)
#endif
   
    if (do_pot) then
#ifdef IS_MPI 
       pot_row(VDW_POT,atom1) = pot_row(VDW_POT,atom1) + 0.5d0*U*sw
       pot_col(VDW_POT,atom2) = pot_col(VDW_POT,atom2) + 0.5d0*U*sw
#else
       pot = pot + U*sw
#endif
    endif
    
    vpair = vpair + U*sw
#ifdef IS_MPI
    id1 = AtomRowToGlobal(atom1)
    id2 = AtomColToGlobal(atom2)
#else
    id1 = atom1
    id2 = atom2
#endif
    
    if (molMembershipList(id1) .ne. molMembershipList(id2)) then
       
       fpair(1) = fpair(1) + fx
       fpair(2) = fpair(2) + fy
       fpair(3) = fpair(3) + fz
       
    endif
    
    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:	2788
Committed:	Mon Jun 5 18:44:05 2006 UTC (18 years, 2 months ago) by gezelter
File size:	15902 byte(s)
Log Message:	Added methods to access GB_mu and GB_nu from the force field options
#	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	gezelter	1608	#ifdef IS_MPI
54			use mpiSimulation
55			#endif
56	kdaily	2226
57	gezelter	1608	implicit none
58
59	kdaily	2367	private
60	gezelter	1608
61	gezelter	2375	#define __FORTRAN90
62			#include "UseTheForce/DarkSide/fInteractionMap.h"
63	gezelter	1608
64	gezelter	2787	logical, save :: haveGBMap = .false.
65			logical, save :: haveMixingMap = .false.
66			real(kind=dp), save :: mu = 2.0_dp
67			real(kind=dp), save :: nu = 1.0_dp
68
69
70	gezelter	2375	public :: newGBtype
71	gezelter	2787	public :: complete_GB_FF
72	gezelter	1608	public :: do_gb_pair
73	chrisfen	2277	public :: getGayBerneCut
74	gezelter	2375	public :: destroyGBtypes
75	gezelter	1608
76	gezelter	2375	type :: GBtype
77			integer :: atid
78	gezelter	2787	real(kind = dp ) :: d
79			real(kind = dp ) :: l
80	gezelter	2375	real(kind = dp ) :: eps
81			real(kind = dp ) :: eps_ratio
82	gezelter	2787	real(kind = dp ) :: dw
83			logical :: isLJ
84	gezelter	2375	end type GBtype
85	gezelter	2787
86	gezelter	2375	type, private :: GBList
87			integer :: nGBtypes = 0
88			integer :: currentGBtype = 0
89			type(GBtype), pointer :: GBtypes(:) => null()
90			integer, pointer :: atidToGBtype(:) => null()
91			end type GBList
92	gezelter	2787
93	gezelter	2375	type(GBList), save :: GBMap
94	gezelter	2787
95			type :: GBMixParameters
96			real(kind=DP) :: sigma0
97			real(kind=DP) :: eps0
98			real(kind=DP) :: dw
99			real(kind=DP) :: x2
100			real(kind=DP) :: xa2
101			real(kind=DP) :: xai2
102			real(kind=DP) :: xp2
103			real(kind=DP) :: xpap2
104			real(kind=DP) :: xpapi2
105			end type GBMixParameters
106
107			type(GBMixParameters), dimension(:,:), allocatable :: GBMixingMap
108
109	gezelter	1608	contains
110	gezelter	2787
111			subroutine newGBtype(c_ident, d, l, eps, eps_ratio, dw, status)
112	gezelter	2375
113			integer, intent(in) :: c_ident
114	gezelter	2787	real( kind = dp ), intent(in) :: d, l, eps, eps_ratio, dw
115	gezelter	2375	integer, intent(out) :: status
116	gezelter	2787
117			integer :: nGBTypes, nLJTypes, ntypes, myATID
118	gezelter	2375	integer, pointer :: MatchList(:) => null()
119			integer :: current, i
120	kdaily	2367	status = 0
121	gezelter	2787
122	gezelter	2375	if (.not.associated(GBMap%GBtypes)) then
123	gezelter	2787
124	gezelter	2375	call getMatchingElementList(atypes, "is_GayBerne", .true., &
125			nGBtypes, MatchList)
126
127	gezelter	2787	call getMatchingElementList(atypes, "is_LennardJones", .true., &
128			nLJTypes, MatchList)
129
130			GBMap%nGBtypes = nGBtypes + nLJTypes
131
132			allocate(GBMap%GBtypes(nGBtypes + nLJTypes))
133
134	gezelter	2375	ntypes = getSize(atypes)
135
136	gezelter	2787	allocate(GBMap%atidToGBtype(ntypes))
137			endif
138
139			GBMap%currentGBtype = GBMap%currentGBtype + 1
140			current = GBMap%currentGBtype
141
142			myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
143
144			GBMap%atidToGBtype(myATID) = current
145			GBMap%GBtypes(current)%atid = myATID
146			GBMap%GBtypes(current)%d = d
147			GBMap%GBtypes(current)%l = l
148			GBMap%GBtypes(current)%eps = eps
149			GBMap%GBtypes(current)%eps_ratio = eps_ratio
150			GBMap%GBtypes(current)%dw = dw
151			GBMap%GBtypes(current)%isLJ = .false.
152
153			return
154			end subroutine newGBtype
155
156			subroutine complete_GB_FF(status)
157			integer :: status
158			integer :: i, j, l, m, lm, function_type
159			real(kind=dp) :: thisDP, sigma
160			integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, myATID, current
161			logical :: thisProperty
162
163			status = 0
164			if (GBMap%currentGBtype == 0) then
165			call handleError("complete_GB_FF", "No members in GBMap")
166			status = -1
167			return
168			end if
169
170			nAtypes = getSize(atypes)
171
172			if (nAtypes == 0) then
173			status = -1
174			return
175			end if
176
177			! atypes comes from c side
178			do i = 1, nAtypes
179	gezelter	2375
180	gezelter	2787	myATID = getFirstMatchingElement(atypes, 'c_ident', i)
181			call getElementProperty(atypes, myATID, "is_LennardJones", thisProperty)
182
183			if (thisProperty) then
184			GBMap%currentGBtype = GBMap%currentGBtype + 1
185			current = GBMap%currentGBtype
186
187			GBMap%atidToGBtype(myATID) = current
188			GBMap%GBtypes(current)%atid = myATID
189			GBMap%GBtypes(current)%isLJ = .true.
190			GBMap%GBtypes(current)%d = getSigma(myATID)
191			GBMap%GBtypes(current)%l = GBMap%GBtypes(current)%d
192			GBMap%GBtypes(current)%eps = getEpsilon(myATID)
193			GBMap%GBtypes(current)%eps_ratio = 1.0_dp
194			GBMap%GBtypes(current)%dw = 1.0_dp
195
196			endif
197
198			end do
199
200			haveGBMap = .true.
201	gezelter	2788
202			mu = getGayBerneMu()
203			nu = getGayBerneNu()
204
205	gezelter	2787
206			end subroutine complete_GB_FF
207	kdaily	2367
208	gezelter	2787	subroutine createGBMixingMap()
209			integer :: nGBtypes, i, j
210			real (kind = dp) :: d1, l1, e1, er1, dw1
211			real (kind = dp) :: d2, l2, e2, er2, dw2
212			real (kind = dp) :: er, ermu, xp, ap2
213	kdaily	2367
214	gezelter	2787	if (GBMap%currentGBtype == 0) then
215			call handleError("GB", "No members in GBMap")
216			return
217			end if
218
219			nGBtypes = GBMap%nGBtypes
220
221			if (.not. allocated(GBMixingMap)) then
222			allocate(GBMixingMap(nGBtypes, nGBtypes))
223	gezelter	2375	endif
224	kdaily	2367
225	gezelter	2787	do i = 1, nGBtypes
226	kdaily	2367
227	gezelter	2787	d1 = GBMap%GBtypes(i)%d
228			l1 = GBMap%GBtypes(i)%l
229			e1 = GBMap%GBtypes(i)%eps
230			er1 = GBMap%GBtypes(i)%eps_ratio
231			dw1 = GBMap%GBtypes(i)%dw
232	gezelter	2375
233	gezelter	2787	do j = i, nGBtypes
234	gezelter	2375
235	gezelter	2787	d2 = GBMap%GBtypes(j)%d
236			l2 = GBMap%GBtypes(j)%l
237			e2 = GBMap%GBtypes(j)%eps
238			er2 = GBMap%GBtypes(j)%eps_ratio
239			dw2 = GBMap%GBtypes(j)%dw
240
241			GBMixingMap(i,j)%sigma0 = sqrt(d1d1 + d2d2)
242			GBMixingMap(i,j)%xa2 = (l1l1 - d1d1)/(l1l1 + d2d2)
243			GBMixingMap(i,j)%xai2 = (l2l2 - d2d2)/(l2l2 + d1d1)
244			GBMixingMap(i,j)%x2 = (l1l1 - d1d1) * (l2l2 - d2d2) / &
245			((l2l2 + d1d1) * (l1l1 + d2d2))
246
247			! assumed LB mixing rules for now:
248
249			GBMixingMap(i,j)%dw = 0.5_dp * (dw1 + dw2)
250			GBMixingMap(i,j)%eps0 = sqrt(e1 * e2)
251
252			er = sqrt(er1 * er2)
253			ermu = er**(1.0_dp / mu)
254			xp = (1.0_dp - ermu) / (1.0_dp + ermu)
255			ap2 = 1.0_dp / (1.0_dp + ermu)
256
257			GBMixingMap(i,j)%xp2 = xp*xp
258			GBMixingMap(i,j)%xpap2 = xp*ap2
259			GBMixingMap(i,j)%xpapi2 = xp/ap2
260
261			if (i.ne.j) then
262			GBMixingMap(j,i)%sigma0 = GBMixingMap(i,j)%sigma0
263			GBMixingMap(j,i)%dw = GBMixingMap(i,j)%dw
264			GBMixingMap(j,i)%eps0 = GBMixingMap(i,j)%eps0
265			GBMixingMap(j,i)%x2 = GBMixingMap(i,j)%x2
266			GBMixingMap(j,i)%xa2 = GBMixingMap(i,j)%xa2
267			GBMixingMap(j,i)%xai2 = GBMixingMap(i,j)%xai2
268			GBMixingMap(j,i)%xp2 = GBMixingMap(i,j)%xp2
269			GBMixingMap(j,i)%xpap2 = GBMixingMap(i,j)%xpap2
270			GBMixingMap(j,i)%xpapi2 = GBMixingMap(i,j)%xpapi2
271			endif
272			enddo
273			enddo
274			haveMixingMap = .true.
275
276			end subroutine createGBMixingMap
277	gezelter	2375
278	gezelter	2787
279	chrisfen	2279	!! gay berne cutoff should be a parameter in globals, this is a temporary
280			!! work around - this should be fixed when gay berne is up and running
281	gezelter	2375
282	chrisfen	2277	function getGayBerneCut(atomID) result(cutValue)
283	gezelter	2375	integer, intent(in) :: atomID
284			integer :: gbt1
285	gezelter	2787	real(kind=dp) :: cutValue, l, d
286	gezelter	1608
287	gezelter	2375	if (GBMap%currentGBtype == 0) then
288			call handleError("GB", "No members in GBMap")
289			return
290			end if
291
292			gbt1 = GBMap%atidToGBtype(atomID)
293	gezelter	2787	l = GBMap%GBtypes(gbt1)%l
294			d = GBMap%GBtypes(gbt1)%d
295			cutValue = 2.5_dp*max(l,d)
296	gezelter	2375
297	chrisfen	2277	end function getGayBerneCut
298
299	gezelter	1608	subroutine do_gb_pair(atom1, atom2, d, r, r2, sw, vpair, fpair, &
300	gezelter	2787	pot, Amat, f, t, do_pot)
301	kdaily	2226
302	gezelter	1608	integer, intent(in) :: atom1, atom2
303	gezelter	2375	integer :: atid1, atid2, gbt1, gbt2, id1, id2
304	gezelter	1608	real (kind=dp), intent(inout) :: r, r2
305			real (kind=dp), dimension(3), intent(in) :: d
306			real (kind=dp), dimension(3), intent(inout) :: fpair
307			real (kind=dp) :: pot, sw, vpair
308	gezelter	2787	real (kind=dp), dimension(9,nLocal) :: Amat
309	gezelter	1608	real (kind=dp), dimension(3,nLocal) :: f
310			real (kind=dp), dimension(3,nLocal) :: t
311			logical, intent(in) :: do_pot
312	gezelter	2787	real (kind = dp), dimension(3) :: ul1, ul2, rxu1, rxu2, uxu, rhat
313	gezelter	1608
314	gezelter	2787	real (kind = dp) :: sigma0, dw, eps0, x2, xa2, xai2, xp2, xpap2, xpapi2
315			real (kind = dp) :: e1, e2, eps, sigma, s3, s03, au, bu, a, b, g, g2
316			real (kind = dp) :: U, BigR, R3, R6, R7, R12, R13, H, Hp, fx, fy, fz
317			real (kind = dp) :: dUdr, dUda, dUdb, dUdg, pref1, pref2
318			logical :: i_is_lj, j_is_lj
319
320			if (.not.haveMixingMap) then
321			call createGBMixingMap()
322			endif
323
324	gezelter	2375	#ifdef IS_MPI
325			atid1 = atid_Row(atom1)
326			atid2 = atid_Col(atom2)
327			#else
328			atid1 = atid(atom1)
329			atid2 = atid(atom2)
330			#endif
331	gezelter	1608
332	gezelter	2375	gbt1 = GBMap%atidToGBtype(atid1)
333	gezelter	2787	gbt2 = GBMap%atidToGBtype(atid2)
334	gezelter	2375
335	gezelter	2787	i_is_LJ = GBMap%GBTypes(gbt1)%isLJ
336			j_is_LJ = GBMap%GBTypes(gbt2)%isLJ
337	gezelter	2375
338	gezelter	2787	sigma0 = GBMixingMap(gbt1, gbt2)%sigma0
339			dw = GBMixingMap(gbt1, gbt2)%dw
340			eps0 = GBMixingMap(gbt1, gbt2)%eps0
341			x2 = GBMixingMap(gbt1, gbt2)%x2
342			xa2 = GBMixingMap(gbt1, gbt2)%xa2
343			xai2 = GBMixingMap(gbt1, gbt2)%xai2
344			xp2 = GBMixingMap(gbt1, gbt2)%xp2
345			xpap2 = GBMixingMap(gbt1, gbt2)%xpap2
346			xpapi2 = GBMixingMap(gbt1, gbt2)%xpapi2
347
348	gezelter	1608	#ifdef IS_MPI
349	gezelter	2385	ul1(1) = A_Row(7,atom1)
350			ul1(2) = A_Row(8,atom1)
351	gezelter	1930	ul1(3) = A_Row(9,atom1)
352	gezelter	1608
353	gezelter	2385	ul2(1) = A_Col(7,atom2)
354			ul2(2) = A_Col(8,atom2)
355	gezelter	1930	ul2(3) = A_Col(9,atom2)
356	gezelter	1608	#else
357	gezelter	2787	ul1(1) = Amat(7,atom1)
358			ul1(2) = Amat(8,atom1)
359			ul1(3) = Amat(9,atom1)
360	gezelter	1608
361	gezelter	2787	ul2(1) = Amat(7,atom2)
362			ul2(2) = Amat(8,atom2)
363			ul2(3) = Amat(9,atom2)
364	gezelter	1608	#endif
365	kdaily	2226
366	gezelter	2787	if (i_is_LJ) then
367			a = 0.0_dp
368			ul1 = 0.0_dp
369			else
370			a = d(1)ul1(1) + d(2)ul1(2) + d(3)*ul1(3)
371			endif
372	gezelter	1608
373	gezelter	2787	if (j_is_LJ) then
374			b = 0.0_dp
375			ul2 = 0.0_dp
376			else
377			b = d(1)ul2(1) + d(2)ul2(2) + d(3)*ul2(3)
378			endif
379	gezelter	1608
380	gezelter	2787	if (i_is_LJ.or.j_is_LJ) then
381			g = 0.0_dp
382			else
383			g = ul1(1)ul2(1) + ul1(2)ul2(2) + ul1(3)*ul2(3)
384			endif
385	gezelter	1608
386	gezelter	2787	au = a / r
387			bu = b / r
388			g2 = g*g
389	gezelter	1608
390	gezelter	2787	H = (xa2 * au + xai2 * bu - 2.0_dpx2aubug) / (1.0_dp - x2*g2)
391			Hp = (xpap2au + xpapi2bu - 2.0_dpxp2aubug) / (1.0_dp - xp2*g2)
392			sigma = sigma0 / sqrt(1.0_dp - H)
393			e1 = 1.0_dp / sqrt(1.0_dp - x2*g2)
394			e2 = 1.0_dp - Hp
395			eps = eps0 * (e1*nu) (e2**mu)
396			BigR = dwsigma0 / (r - sigma + dwsigma0)
397
398			R3 = BigRBigRBigR
399			R6 = R3*R3
400			R7 = R6 * BigR
401			R12 = R6*R6
402			R13 = R6*R7
403	gezelter	1608
404	gezelter	2787	U = 4.0_dp * eps * (R12 - R6)
405	gezelter	1608
406	gezelter	2787	s3 = sigmasigmasigma
407			s03 = sigma0sigma0sigma0
408	kdaily	2226
409	gezelter	2787	pref1 = - 8.0_dp * eps * mu * (R12 - R6) / (e2 * r)
410			pref2 = 8.0_dp * eps * s3 * (6.0_dpR13 - 3.0_dpR7) / (dwrs03)
411	gezelter	2204
412	gezelter	2787	dUdr = - (pref1 * Hp + pref2 * (sigma0sigma0r/s3 - H))
413	kdaily	2226
414	gezelter	2787	dUda = pref1 * (xpap2au - xp2bug) / (1.0_dp - xp2 g2) &
415			+ pref2 * (xa2 * au - x2 bug) / (1.0_dp - x2 * g2)
416	kdaily	2226
417	gezelter	2787	dUdb = pref1 * (xpapi2bu - xp2aug) / (1.0_dp - xp2 g2) &
418			+ pref2 * (xai2 * bu - x2 aug) / (1.0_dp - x2 * g2)
419	gezelter	2204
420	gezelter	2787	dUdg = 4.0_dp * eps * nu * (R12 - R6) * x2 * g / (1.0_dp - x2*g2) &
421			+ 8.0_dp * eps * mu * (R12 - R6) * (xp2aubu - Hpxp2g) / &
422			(1.0_dp - xp2 * g2) / e2 &
423			+ 8.0_dp * eps * s3 * (3.0_dp * R7 - 6.0_dp * R13) * &
424			(x2 * au * bu - H * x2 * g) / (1.0_dp - x2 * g2) / (dw * s03)
425
426			rhat = d / r
427	gezelter	2204
428	gezelter	2787	fx = -dUdr * rhat(1) - dUda * ul1(1) - dUdb * ul2(1)
429			fy = -dUdr * rhat(2) - dUda * ul1(2) - dUdb * ul2(2)
430			fx = -dUdr * rhat(3) - dUda * ul1(3) - dUdb * ul2(3)
431	gezelter	1608
432	gezelter	2787	rxu1 = cross_product(d, ul1)
433			rxu2 = cross_product(d, ul2)
434			uxu = cross_product(ul1, ul2)
435
436	gezelter	1608	#ifdef IS_MPI
437	gezelter	2787	f_Row(1,atom1) = f_Row(1,atom1) + fx
438			f_Row(2,atom1) = f_Row(2,atom1) + fy
439			f_Row(3,atom1) = f_Row(3,atom1) + fz
440	kdaily	2226
441	gezelter	2787	f_Col(1,atom2) = f_Col(1,atom2) - fx
442			f_Col(2,atom2) = f_Col(2,atom2) - fy
443			f_Col(3,atom2) = f_Col(3,atom2) - fz
444	kdaily	2226
445	gezelter	2787	t_Row(1,atom1) = t_Row(1,atom1) + dUdarxu1(1) - dUdguxu(1)
446			t_Row(2,atom1) = t_Row(2,atom1) + dUdarxu1(2) - dUdguxu(2)
447			t_Row(3,atom1) = t_Row(3,atom1) + dUdarxu1(3) - dUdguxu(3)
448
449			t_Col(1,atom2) = t_Col(1,atom2) + dUdbrxu2(1) + dUdguxu(1)
450			t_Col(2,atom2) = t_Col(2,atom2) + dUdbrxu2(2) + dUdguxu(2)
451			t_Col(3,atom2) = t_Col(3,atom2) + dUdbrxu2(3) + dUdguxu(3)
452	gezelter	1608	#else
453	gezelter	2787	f(1,atom1) = f(1,atom1) + fx
454			f(2,atom1) = f(2,atom1) + fy
455			f(3,atom1) = f(3,atom1) + fz
456	kdaily	2226
457	gezelter	2787	f(1,atom2) = f(1,atom2) - fx
458			f(2,atom2) = f(2,atom2) - fy
459			f(3,atom2) = f(3,atom2) - fz
460	kdaily	2226
461	gezelter	2787	t(1,atom1) = t(1,atom1) + dUdarxu1(1) - dUdguxu(1)
462			t(2,atom1) = t(2,atom1) + dUdarxu1(2) - dUdguxu(2)
463			t(3,atom1) = t(3,atom1) + dUdarxu1(3) - dUdguxu(3)
464
465			t(1,atom2) = t(1,atom2) + dUdbrxu2(1) + dUdguxu(1)
466			t(2,atom2) = t(2,atom2) + dUdbrxu2(2) + dUdguxu(2)
467			t(3,atom2) = t(3,atom2) + dUdbrxu2(3) + dUdguxu(3)
468	gezelter	1608	#endif
469	kdaily	2367
470	gezelter	1608	if (do_pot) then
471			#ifdef IS_MPI
472	gezelter	2787	pot_row(VDW_POT,atom1) = pot_row(VDW_POT,atom1) + 0.5d0Usw
473			pot_col(VDW_POT,atom2) = pot_col(VDW_POT,atom2) + 0.5d0Usw
474	gezelter	1608	#else
475	gezelter	2787	pot = pot + U*sw
476	gezelter	1608	#endif
477			endif
478	gezelter	2375
479	gezelter	2787	vpair = vpair + U*sw
480	gezelter	1608	#ifdef IS_MPI
481			id1 = AtomRowToGlobal(atom1)
482			id2 = AtomColToGlobal(atom2)
483			#else
484			id1 = atom1
485			id2 = atom2
486			#endif
487	kdaily	2226
488	gezelter	1608	if (molMembershipList(id1) .ne. molMembershipList(id2)) then
489	kdaily	2226
490	gezelter	2787	fpair(1) = fpair(1) + fx
491			fpair(2) = fpair(2) + fy
492			fpair(3) = fpair(3) + fz
493	kdaily	2226
494	gezelter	1608	endif
495	kdaily	2226
496	gezelter	1608	return
497			end subroutine do_gb_pair
498	gezelter	2787
499	gezelter	2375	subroutine destroyGBTypes()
500
501			GBMap%nGBtypes = 0
502			GBMap%currentGBtype = 0
503	kdaily	2367
504	gezelter	2375	if (associated(GBMap%GBtypes)) then
505			deallocate(GBMap%GBtypes)
506			GBMap%GBtypes => null()
507	kdaily	2367	end if
508
509	gezelter	2375	if (associated(GBMap%atidToGBtype)) then
510			deallocate(GBMap%atidToGBtype)
511			GBMap%atidToGBtype => null()
512			end if
513	kdaily	2367
514	gezelter	2787	haveMixingMap = .false.
515
516	gezelter	2375	end subroutine destroyGBTypes
517	kdaily	2367
518	gezelter	2375	end module gayberne
519	kdaily	2367