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.

    
  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
!         write(*,*) 'd2 = ', d2,l2, d1,l2 
!          GBMixingMap(i,j)%sigma0 = sqrt(d1*d1 + d2*d2)
          GBMixingMap(i,j)%sigma0 = 0.5_dp*(d1 + 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.
    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   
    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, 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

#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

    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 = 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)
    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
!!$    
         
   
#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:	997
Committed:	Sat Jul 1 22:27:39 2006 UTC (19 years ago) by kdaily
File size:	16609 byte(s)
Log Message:	fixed forces and torques for GB
#	User	Rev	Content
1	gezelter	246	!!
2			!! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			!!
4			!! The University of Notre Dame grants you ("Licensee") a
5			!! non-exclusive, royalty free, license to use, modify and
6			!! redistribute this software in source and binary code form, provided
7			!! that the following conditions are met:
8			!!
9			!! 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	676	module gayberne
44	gezelter	115	use force_globals
45			use definitions
46			use simulation
47	kdaily	668	use atype_module
48			use vector_class
49	gezelter	981	use linearalgebra
50	kdaily	668	use status
51			use lj
52	gezelter	982	use fForceOptions
53	gezelter	115	#ifdef IS_MPI
54			use mpiSimulation
55			#endif
56	kdaily	529
57	gezelter	115	implicit none
58
59	kdaily	668	private
60	gezelter	115
61	gezelter	676	#define __FORTRAN90
62			#include "UseTheForce/DarkSide/fInteractionMap.h"
63	gezelter	115
64	gezelter	981	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	676	public :: newGBtype
71	gezelter	981	public :: complete_GB_FF
72	gezelter	115	public :: do_gb_pair
73	chrisfen	578	public :: getGayBerneCut
74	gezelter	676	public :: destroyGBtypes
75	gezelter	115
76	gezelter	676	type :: GBtype
77			integer :: atid
78	gezelter	981	real(kind = dp ) :: d
79			real(kind = dp ) :: l
80	gezelter	676	real(kind = dp ) :: eps
81			real(kind = dp ) :: eps_ratio
82	gezelter	981	real(kind = dp ) :: dw
83			logical :: isLJ
84	gezelter	676	end type GBtype
85	gezelter	981
86	gezelter	676	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	981
93	gezelter	676	type(GBList), save :: GBMap
94	gezelter	981
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	115	contains
110	gezelter	981
111			subroutine newGBtype(c_ident, d, l, eps, eps_ratio, dw, status)
112	gezelter	676
113			integer, intent(in) :: c_ident
114	gezelter	981	real( kind = dp ), intent(in) :: d, l, eps, eps_ratio, dw
115	gezelter	676	integer, intent(out) :: status
116	gezelter	981
117			integer :: nGBTypes, nLJTypes, ntypes, myATID
118	gezelter	676	integer, pointer :: MatchList(:) => null()
119			integer :: current, i
120	kdaily	668	status = 0
121	gezelter	981
122	gezelter	676	if (.not.associated(GBMap%GBtypes)) then
123	gezelter	981
124	gezelter	676	call getMatchingElementList(atypes, "is_GayBerne", .true., &
125			nGBtypes, MatchList)
126
127	gezelter	981	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	676	ntypes = getSize(atypes)
135
136	gezelter	981	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	676
180	gezelter	981	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	982
202	gezelter	981
203			end subroutine complete_GB_FF
204	kdaily	668
205	gezelter	981	subroutine createGBMixingMap()
206			integer :: nGBtypes, i, j
207			real (kind = dp) :: d1, l1, e1, er1, dw1
208			real (kind = dp) :: d2, l2, e2, er2, dw2
209			real (kind = dp) :: er, ermu, xp, ap2
210	kdaily	668
211	gezelter	981	if (GBMap%currentGBtype == 0) then
212			call handleError("GB", "No members in GBMap")
213			return
214			end if
215
216			nGBtypes = GBMap%nGBtypes
217
218			if (.not. allocated(GBMixingMap)) then
219			allocate(GBMixingMap(nGBtypes, nGBtypes))
220	gezelter	676	endif
221	kdaily	668
222	gezelter	981	do i = 1, nGBtypes
223	kdaily	668
224	gezelter	981	d1 = GBMap%GBtypes(i)%d
225			l1 = GBMap%GBtypes(i)%l
226			e1 = GBMap%GBtypes(i)%eps
227			er1 = GBMap%GBtypes(i)%eps_ratio
228			dw1 = GBMap%GBtypes(i)%dw
229	gezelter	676
230	gezelter	981	do j = i, nGBtypes
231	gezelter	676
232	gezelter	981	d2 = GBMap%GBtypes(j)%d
233			l2 = GBMap%GBtypes(j)%l
234			e2 = GBMap%GBtypes(j)%eps
235			er2 = GBMap%GBtypes(j)%eps_ratio
236			dw2 = GBMap%GBtypes(j)%dw
237	kdaily	997	! write(,) 'd2 = ', d2,l2, d1,l2
238			! GBMixingMap(i,j)%sigma0 = sqrt(d1d1 + d2d2)
239			GBMixingMap(i,j)%sigma0 = 0.5_dp*(d1 + d2)
240	gezelter	981	GBMixingMap(i,j)%xa2 = (l1l1 - d1d1)/(l1l1 + d2d2)
241			GBMixingMap(i,j)%xai2 = (l2l2 - d2d2)/(l2l2 + d1d1)
242			GBMixingMap(i,j)%x2 = (l1l1 - d1d1) * (l2l2 - d2d2) / &
243			((l2l2 + d1d1) * (l1l1 + d2d2))
244
245			! assumed LB mixing rules for now:
246
247			GBMixingMap(i,j)%dw = 0.5_dp * (dw1 + dw2)
248			GBMixingMap(i,j)%eps0 = sqrt(e1 * e2)
249
250			er = sqrt(er1 * er2)
251			ermu = er**(1.0_dp / mu)
252			xp = (1.0_dp - ermu) / (1.0_dp + ermu)
253			ap2 = 1.0_dp / (1.0_dp + ermu)
254
255			GBMixingMap(i,j)%xp2 = xp*xp
256			GBMixingMap(i,j)%xpap2 = xp*ap2
257			GBMixingMap(i,j)%xpapi2 = xp/ap2
258
259			if (i.ne.j) then
260			GBMixingMap(j,i)%sigma0 = GBMixingMap(i,j)%sigma0
261			GBMixingMap(j,i)%dw = GBMixingMap(i,j)%dw
262			GBMixingMap(j,i)%eps0 = GBMixingMap(i,j)%eps0
263			GBMixingMap(j,i)%x2 = GBMixingMap(i,j)%x2
264			GBMixingMap(j,i)%xa2 = GBMixingMap(i,j)%xa2
265			GBMixingMap(j,i)%xai2 = GBMixingMap(i,j)%xai2
266			GBMixingMap(j,i)%xp2 = GBMixingMap(i,j)%xp2
267			GBMixingMap(j,i)%xpap2 = GBMixingMap(i,j)%xpap2
268			GBMixingMap(j,i)%xpapi2 = GBMixingMap(i,j)%xpapi2
269			endif
270			enddo
271			enddo
272			haveMixingMap = .true.
273	gezelter	983	mu = getGayBerneMu()
274			nu = getGayBerneNu()
275	gezelter	981	end subroutine createGBMixingMap
276	gezelter	676
277	gezelter	981
278	chrisfen	580	!! gay berne cutoff should be a parameter in globals, this is a temporary
279			!! work around - this should be fixed when gay berne is up and running
280	gezelter	676
281	chrisfen	578	function getGayBerneCut(atomID) result(cutValue)
282	gezelter	676	integer, intent(in) :: atomID
283			integer :: gbt1
284	gezelter	981	real(kind=dp) :: cutValue, l, d
285	gezelter	115
286	gezelter	676	if (GBMap%currentGBtype == 0) then
287			call handleError("GB", "No members in GBMap")
288			return
289			end if
290
291			gbt1 = GBMap%atidToGBtype(atomID)
292	gezelter	981	l = GBMap%GBtypes(gbt1)%l
293			d = GBMap%GBtypes(gbt1)%d
294			cutValue = 2.5_dp*max(l,d)
295	gezelter	676
296	chrisfen	578	end function getGayBerneCut
297
298	gezelter	115	subroutine do_gb_pair(atom1, atom2, d, r, r2, sw, vpair, fpair, &
299	gezelter	981	pot, Amat, f, t, do_pot)
300	kdaily	529
301	gezelter	115	integer, intent(in) :: atom1, atom2
302	gezelter	676	integer :: atid1, atid2, gbt1, gbt2, id1, id2
303	gezelter	115	real (kind=dp), intent(inout) :: r, r2
304			real (kind=dp), dimension(3), intent(in) :: d
305			real (kind=dp), dimension(3), intent(inout) :: fpair
306			real (kind=dp) :: pot, sw, vpair
307	gezelter	981	real (kind=dp), dimension(9,nLocal) :: Amat
308	gezelter	115	real (kind=dp), dimension(3,nLocal) :: f
309			real (kind=dp), dimension(3,nLocal) :: t
310			logical, intent(in) :: do_pot
311	gezelter	981	real (kind = dp), dimension(3) :: ul1, ul2, rxu1, rxu2, uxu, rhat
312	gezelter	115
313	gezelter	981	real (kind = dp) :: sigma0, dw, eps0, x2, xa2, xai2, xp2, xpap2, xpapi2
314	kdaily	997	real (kind = dp) :: e1, e2, eps, sigma, s3, s03, au2, bu2, au, bu, a, b, g, g2
315	gezelter	981	real (kind = dp) :: U, BigR, R3, R6, R7, R12, R13, H, Hp, fx, fy, fz
316			real (kind = dp) :: dUdr, dUda, dUdb, dUdg, pref1, pref2
317			logical :: i_is_lj, j_is_lj
318
319			if (.not.haveMixingMap) then
320			call createGBMixingMap()
321			endif
322
323	gezelter	676	#ifdef IS_MPI
324			atid1 = atid_Row(atom1)
325			atid2 = atid_Col(atom2)
326			#else
327			atid1 = atid(atom1)
328			atid2 = atid(atom2)
329			#endif
330	gezelter	115
331	gezelter	676	gbt1 = GBMap%atidToGBtype(atid1)
332	gezelter	981	gbt2 = GBMap%atidToGBtype(atid2)
333	gezelter	676
334	gezelter	981	i_is_LJ = GBMap%GBTypes(gbt1)%isLJ
335			j_is_LJ = GBMap%GBTypes(gbt2)%isLJ
336	gezelter	676
337	gezelter	981	sigma0 = GBMixingMap(gbt1, gbt2)%sigma0
338			dw = GBMixingMap(gbt1, gbt2)%dw
339			eps0 = GBMixingMap(gbt1, gbt2)%eps0
340			x2 = GBMixingMap(gbt1, gbt2)%x2
341			xa2 = GBMixingMap(gbt1, gbt2)%xa2
342			xai2 = GBMixingMap(gbt1, gbt2)%xai2
343			xp2 = GBMixingMap(gbt1, gbt2)%xp2
344			xpap2 = GBMixingMap(gbt1, gbt2)%xpap2
345			xpapi2 = GBMixingMap(gbt1, gbt2)%xpapi2
346
347	gezelter	115	#ifdef IS_MPI
348	gezelter	686	ul1(1) = A_Row(7,atom1)
349			ul1(2) = A_Row(8,atom1)
350	gezelter	246	ul1(3) = A_Row(9,atom1)
351	gezelter	115
352	gezelter	686	ul2(1) = A_Col(7,atom2)
353			ul2(2) = A_Col(8,atom2)
354	gezelter	246	ul2(3) = A_Col(9,atom2)
355	gezelter	115	#else
356	gezelter	981	ul1(1) = Amat(7,atom1)
357			ul1(2) = Amat(8,atom1)
358			ul1(3) = Amat(9,atom1)
359	gezelter	115
360	gezelter	981	ul2(1) = Amat(7,atom2)
361			ul2(2) = Amat(8,atom2)
362			ul2(3) = Amat(9,atom2)
363	gezelter	115	#endif
364	kdaily	529
365	gezelter	981	if (i_is_LJ) then
366			a = 0.0_dp
367			ul1 = 0.0_dp
368			else
369			a = d(1)ul1(1) + d(2)ul1(2) + d(3)*ul1(3)
370			endif
371	gezelter	115
372	gezelter	981	if (j_is_LJ) then
373			b = 0.0_dp
374			ul2 = 0.0_dp
375			else
376			b = d(1)ul2(1) + d(2)ul2(2) + d(3)*ul2(3)
377			endif
378	gezelter	115
379	gezelter	981	if (i_is_LJ.or.j_is_LJ) then
380			g = 0.0_dp
381			else
382			g = ul1(1)ul2(1) + ul1(2)ul2(2) + ul1(3)*ul2(3)
383			endif
384			au = a / r
385			bu = b / r
386	kdaily	997
387			au2 = au*au
388			bu2 = bu*bu
389	gezelter	981	g2 = g*g
390	gezelter	115
391	kdaily	997
392			H = (xa2 * au2 + xai2 * bu2 - 2.0_dpx2aubug) / (1.0_dp - x2*g2)
393			Hp = (xpap2au2 + xpapi2bu2 - 2.0_dpxp2aubug) / (1.0_dp - xp2*g2)
394	gezelter	981	sigma = sigma0 / sqrt(1.0_dp - H)
395			e1 = 1.0_dp / sqrt(1.0_dp - x2*g2)
396			e2 = 1.0_dp - Hp
397			eps = eps0 * (e1*nu) (e2**mu)
398			BigR = dwsigma0 / (r - sigma + dwsigma0)
399
400			R3 = BigRBigRBigR
401			R6 = R3*R3
402			R7 = R6 * BigR
403			R12 = R6*R6
404			R13 = R6*R7
405	gezelter	115
406	gezelter	981	U = 4.0_dp * eps * (R12 - R6)
407	gezelter	115
408	gezelter	981	s3 = sigmasigmasigma
409			s03 = sigma0sigma0sigma0
410	kdaily	529
411	gezelter	981	pref1 = - 8.0_dp * eps * mu * (R12 - R6) / (e2 * r)
412	gezelter	983
413	kdaily	997
414	gezelter	981	pref2 = 8.0_dp * eps * s3 * (6.0_dpR13 - 3.0_dpR7) / (dwrs03)
415	gezelter	507
416	kdaily	997	dUdr = - (pref1 * Hp + pref2 * (sigma0sigma0r/s3 + H))
417	kdaily	529
418	gezelter	981	dUda = pref1 * (xpap2au - xp2bug) / (1.0_dp - xp2 g2) &
419			+ pref2 * (xa2 * au - x2 bug) / (1.0_dp - x2 * g2)
420	kdaily	997
421	gezelter	981	dUdb = pref1 * (xpapi2bu - xp2aug) / (1.0_dp - xp2 g2) &
422			+ pref2 * (xai2 * bu - x2 aug) / (1.0_dp - x2 * g2)
423	gezelter	507
424	gezelter	981	dUdg = 4.0_dp * eps * nu * (R12 - R6) * x2 * g / (1.0_dp - x2*g2) &
425			+ 8.0_dp * eps * mu * (R12 - R6) * (xp2aubu - Hpxp2g) / &
426			(1.0_dp - xp2 * g2) / e2 &
427	gezelter	983	+ 8.0_dp * eps * s3 * (3.0_dp * R7 - 6.0_dp * R13) * &
428	gezelter	981	(x2 * au * bu - H * x2 * g) / (1.0_dp - x2 * g2) / (dw * s03)
429	kdaily	997
430	gezelter	981
431			rhat = d / r
432	gezelter	507
433	gezelter	983	fx = dUdr * rhat(1) + dUda * ul1(1) + dUdb * ul2(1)
434			fy = dUdr * rhat(2) + dUda * ul1(2) + dUdb * ul2(2)
435			fz = dUdr * rhat(3) + dUda * ul1(3) + dUdb * ul2(3)
436	gezelter	115
437	gezelter	981	rxu1 = cross_product(d, ul1)
438			rxu2 = cross_product(d, ul2)
439			uxu = cross_product(ul1, ul2)
440	kdaily	997
441			!!$ write(,) 'pref = ' , pref1, pref2
442	gezelter	983	!!$ write(,) 'rxu1 = ' , rxu1(1), rxu1(2), rxu1(3)
443			!!$ write(,) 'rxu2 = ' , rxu2(1), rxu2(2), rxu2(3)
444			!!$ write(,) 'uxu = ' , uxu(1), uxu(2), uxu(3)
445	kdaily	997	!!$ write(,) 'dUda = ', dUda, dudb, dudg, dudr
446			!!$ write(,) 'H = ', H,hp,sigma, e1, e2, BigR
447			!!$ write(,) 'chi = ', xa2, xai2, x2
448			!!$ write(,) 'chip = ', xpap2, xpapi2, xp2
449			!!$ write(,) 'eps = ', eps0, e1, e2, eps
450			!!$ write(,) 'U =', U, pref1, pref2
451			!!$ write(,) 'f =', fx, fy, fz
452			!!$ write(,) 'au =', au, bu, g
453			!!$
454
455
456	gezelter	115	#ifdef IS_MPI
457	gezelter	981	f_Row(1,atom1) = f_Row(1,atom1) + fx
458			f_Row(2,atom1) = f_Row(2,atom1) + fy
459			f_Row(3,atom1) = f_Row(3,atom1) + fz
460	kdaily	529
461	gezelter	981	f_Col(1,atom2) = f_Col(1,atom2) - fx
462			f_Col(2,atom2) = f_Col(2,atom2) - fy
463			f_Col(3,atom2) = f_Col(3,atom2) - fz
464	kdaily	529
465	gezelter	981	t_Row(1,atom1) = t_Row(1,atom1) + dUdarxu1(1) - dUdguxu(1)
466			t_Row(2,atom1) = t_Row(2,atom1) + dUdarxu1(2) - dUdguxu(2)
467			t_Row(3,atom1) = t_Row(3,atom1) + dUdarxu1(3) - dUdguxu(3)
468
469			t_Col(1,atom2) = t_Col(1,atom2) + dUdbrxu2(1) + dUdguxu(1)
470			t_Col(2,atom2) = t_Col(2,atom2) + dUdbrxu2(2) + dUdguxu(2)
471			t_Col(3,atom2) = t_Col(3,atom2) + dUdbrxu2(3) + dUdguxu(3)
472	gezelter	115	#else
473	gezelter	981	f(1,atom1) = f(1,atom1) + fx
474			f(2,atom1) = f(2,atom1) + fy
475			f(3,atom1) = f(3,atom1) + fz
476	kdaily	529
477	gezelter	981	f(1,atom2) = f(1,atom2) - fx
478			f(2,atom2) = f(2,atom2) - fy
479			f(3,atom2) = f(3,atom2) - fz
480	kdaily	529
481	gezelter	981	t(1,atom1) = t(1,atom1) + dUdarxu1(1) - dUdguxu(1)
482			t(2,atom1) = t(2,atom1) + dUdarxu1(2) - dUdguxu(2)
483			t(3,atom1) = t(3,atom1) + dUdarxu1(3) - dUdguxu(3)
484
485			t(1,atom2) = t(1,atom2) + dUdbrxu2(1) + dUdguxu(1)
486			t(2,atom2) = t(2,atom2) + dUdbrxu2(2) + dUdguxu(2)
487			t(3,atom2) = t(3,atom2) + dUdbrxu2(3) + dUdguxu(3)
488	gezelter	115	#endif
489	kdaily	668
490	gezelter	115	if (do_pot) then
491			#ifdef IS_MPI
492	gezelter	981	pot_row(VDW_POT,atom1) = pot_row(VDW_POT,atom1) + 0.5d0Usw
493			pot_col(VDW_POT,atom2) = pot_col(VDW_POT,atom2) + 0.5d0Usw
494	gezelter	115	#else
495	gezelter	981	pot = pot + U*sw
496	gezelter	115	#endif
497			endif
498	gezelter	676
499	gezelter	981	vpair = vpair + U*sw
500	gezelter	115	#ifdef IS_MPI
501			id1 = AtomRowToGlobal(atom1)
502			id2 = AtomColToGlobal(atom2)
503			#else
504			id1 = atom1
505			id2 = atom2
506			#endif
507	kdaily	529
508	gezelter	115	if (molMembershipList(id1) .ne. molMembershipList(id2)) then
509	kdaily	529
510	gezelter	981	fpair(1) = fpair(1) + fx
511			fpair(2) = fpair(2) + fy
512			fpair(3) = fpair(3) + fz
513	kdaily	529
514	gezelter	115	endif
515	kdaily	529
516	gezelter	115	return
517			end subroutine do_gb_pair
518	gezelter	981
519	gezelter	676	subroutine destroyGBTypes()
520
521			GBMap%nGBtypes = 0
522			GBMap%currentGBtype = 0
523	kdaily	668
524	gezelter	676	if (associated(GBMap%GBtypes)) then
525			deallocate(GBMap%GBtypes)
526			GBMap%GBtypes => null()
527	kdaily	668	end if
528
529	gezelter	676	if (associated(GBMap%atidToGBtype)) then
530			deallocate(GBMap%atidToGBtype)
531			GBMap%atidToGBtype => null()
532			end if
533	kdaily	668
534	gezelter	981	haveMixingMap = .false.
535
536	gezelter	676	end subroutine destroyGBTypes
537	kdaily	668
538	gezelter	676	end module gayberne
539	kdaily	668