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