UseTheForce/DarkSide/LJ.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.
!!


!! Calculates Long Range forces Lennard-Jones interactions.
!! @author Charles F. Vardeman II
!! @author Matthew Meineke
!! @version $Id: LJ.F90,v 1.21 2006-04-17 21:49:12 gezelter Exp $, $Date: 2006-04-17 21:49:12 $, $Name: not supported by cvs2svn $, $Revision: 1.21 $


module lj
  use atype_module
  use vector_class
  use simulation
  use status
  use fForceOptions
  use interpolation
#ifdef IS_MPI
  use mpiSimulation
#endif
  use force_globals

  implicit none
  PRIVATE
#define __FORTRAN90
#include "UseTheForce/DarkSide/fInteractionMap.h"

  integer, parameter :: DP = selected_real_kind(15)

  logical, save :: useGeometricDistanceMixing = .false.
  logical, save :: haveMixingMap = .false.
  logical, save :: useSplines = .false.

  real(kind=DP), save :: defaultCutoff = 0.0_DP
  logical, save :: defaultShift = .false.
  logical, save :: haveDefaultCutoff = .false.

  type, private :: LJtype
     integer       :: atid
     real(kind=dp) :: sigma
     real(kind=dp) :: epsilon
     logical       :: isSoftCore = .false.
  end type LJtype

  type, private :: LJList
     integer               :: Nljtypes = 0
     integer               :: currentLJtype = 0
     type(LJtype), pointer :: LJtypes(:)      => null()
     integer, pointer      :: atidToLJtype(:) => null()
  end type LJList

  type(LJList), save :: LJMap

  type :: MixParameters
     real(kind=DP) :: sigma
     real(kind=DP) :: epsilon
     real(kind=dp) :: sigmai
     real(kind=dp) :: rCut
     logical       :: rCutWasSet = .false.
     logical       :: shiftedPot
     logical       :: isSoftCore = .false.
  end type MixParameters

  type(MixParameters), dimension(:,:), allocatable :: MixingMap

  type(cubicSpline), save :: vLJspline
  type(cubicSpline), save :: vLJpspline
  type(cubicSpline), save :: vSoftSpline
  type(cubicSpline), save :: vSoftpSpline

  public :: newLJtype
  public :: setLJDefaultCutoff
  public :: getSigma
  public :: getEpsilon
  public :: do_lj_pair
  public :: destroyLJtypes
  public :: setLJsplineRmax

contains

  subroutine newLJtype(c_ident, sigma, epsilon, isSoftCore, status)
    integer,intent(in) :: c_ident
    real(kind=dp),intent(in) :: sigma
    real(kind=dp),intent(in) :: epsilon
    integer, intent(in) :: isSoftCore
    integer,intent(out) :: status
    integer :: nLJTypes, ntypes, myATID
    integer, pointer :: MatchList(:) => null()
    integer :: current

    status = 0
    ! check to see if this is the first time into this routine...
    if (.not.associated(LJMap%LJtypes)) then

       call getMatchingElementList(atypes, "is_LennardJones", .true., &
            nLJTypes, MatchList)
       
       LJMap%nLJtypes =  nLJTypes

       allocate(LJMap%LJtypes(nLJTypes))

       ntypes = getSize(atypes)

       allocate(LJMap%atidToLJtype(ntypes))
    end if

    LJMap%currentLJtype = LJMap%currentLJtype + 1
    current = LJMap%currentLJtype

    myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
    LJMap%atidToLJtype(myATID)        = current
    LJMap%LJtypes(current)%atid       = myATID
    LJMap%LJtypes(current)%sigma      = sigma
    LJMap%LJtypes(current)%epsilon    = epsilon 
    if (isSoftCore .eq. 1) then
       LJMap%LJtypes(current)%isSoftCore = .true.
    else
       LJMap%LJtypes(current)%isSoftCore = .false.
    endif
  end subroutine newLJtype

  subroutine setLJDefaultCutoff(thisRcut, shiftedPot)
    real(kind=dp), intent(in) :: thisRcut
    logical, intent(in) :: shiftedPot
    defaultCutoff = thisRcut
    defaultShift = shiftedPot
    haveDefaultCutoff = .true.
    !we only want to build LJ Mixing map and spline if LJ is being used.
    if(LJMap%nLJTypes /= 0) then
       call createMixingMap()
       call setLJsplineRmax(defaultCutoff)
    end if

  end subroutine setLJDefaultCutoff

  function getSigma(atid) result (s)
    integer, intent(in) :: atid
    integer :: ljt1
    real(kind=dp) :: s

    if (LJMap%currentLJtype == 0) then
       call handleError("LJ", "No members in LJMap")
       return
    end if

    ljt1 = LJMap%atidToLJtype(atid)
    s = LJMap%LJtypes(ljt1)%sigma

  end function getSigma

  function getEpsilon(atid) result (e)
    integer, intent(in) :: atid
    integer :: ljt1
    real(kind=dp) :: e

    if (LJMap%currentLJtype == 0) then
       call handleError("LJ", "No members in LJMap")
       return
    end if

    ljt1 = LJMap%atidToLJtype(atid)
    e = LJMap%LJtypes(ljt1)%epsilon

  end function getEpsilon

  subroutine createMixingMap()
    integer :: nLJtypes, i, j
    real ( kind = dp ) :: s1, s2, e1, e2
    real ( kind = dp ) :: rcut6, tp6, tp12
    logical :: isSoftCore1, isSoftCore2, doShift

    if (LJMap%currentLJtype == 0) then
       call handleError("LJ", "No members in LJMap")
       return
    end if

    nLJtypes = LJMap%nLJtypes

    if (.not. allocated(MixingMap)) then
       allocate(MixingMap(nLJtypes, nLJtypes))
    endif

    useGeometricDistanceMixing = usesGeometricDistanceMixing()
    do i = 1, nLJtypes

       s1 = LJMap%LJtypes(i)%sigma
       e1 = LJMap%LJtypes(i)%epsilon
       isSoftCore1 = LJMap%LJtypes(i)%isSoftCore

       do j = i, nLJtypes
          
          s2 = LJMap%LJtypes(j)%sigma
          e2 = LJMap%LJtypes(j)%epsilon
          isSoftCore2 = LJMap%LJtypes(j)%isSoftCore
          
          MixingMap(i,j)%isSoftCore = isSoftCore1 .or. isSoftCore2

          ! only the distance parameter uses different mixing policies
          if (useGeometricDistanceMixing) then
             MixingMap(i,j)%sigma = dsqrt(s1 * s2)
          else
             MixingMap(i,j)%sigma = 0.5_dp * (s1 + s2)
          endif
          
          MixingMap(i,j)%epsilon = dsqrt(e1 * e2)

          MixingMap(i,j)%sigmai = 1.0_DP  / (MixingMap(i,j)%sigma)

          if (haveDefaultCutoff) then
             MixingMap(i,j)%shiftedPot = defaultShift
          else
             MixingMap(i,j)%shiftedPot = defaultShift
          endif          

          if (i.ne.j) then
             MixingMap(j,i)%sigma      = MixingMap(i,j)%sigma
             MixingMap(j,i)%epsilon    = MixingMap(i,j)%epsilon
             MixingMap(j,i)%sigmai     = MixingMap(i,j)%sigmai
             MixingMap(j,i)%rCut       = MixingMap(i,j)%rCut
             MixingMap(j,i)%rCutWasSet = MixingMap(i,j)%rCutWasSet
             MixingMap(j,i)%shiftedPot = MixingMap(i,j)%shiftedPot
             MixingMap(j,i)%isSoftCore = MixingMap(i,j)%isSoftCore
          endif

       enddo
    enddo
    
    haveMixingMap = .true.
    
  end subroutine createMixingMap
  
  subroutine setLJsplineRmax(largestRcut)
    real( kind = dp ), intent(in) :: largestRcut
    real( kind = dp ) :: s, bigS, smallS, rmax, rmin
    integer :: np, i

    if (LJMap%nLJtypes .ne. 0) then

       !
       ! find the largest and smallest values of sigma that we'll need
       !
       bigS = 0.0_DP
       smallS = 1.0e9
       do i = 1, LJMap%nLJtypes
          s = LJMap%LJtypes(i)%sigma
          if (s .gt. bigS) bigS = s
          if (s .lt. smallS) smallS = s
       enddo
       
       !
       ! give ourselves a 20% margin just in case
       !
       rmax = 1.2 * largestRcut / smallS    
       !
       ! assume atoms will never get closer than 1 angstrom
       !
       rmin = 1 / bigS
       !
       ! assume 500 points is enough
       !
       np = 500
       
       write(*,*) 'calling setupSplines with rmin = ', rmin, ' rmax = ', rmax, &
            ' np = ', np
       
       call setupSplines(rmin, rmax, np)

    endif
    return
  end subroutine setLJsplineRmax
        
  subroutine do_lj_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
       pot, f, do_pot)
    
    integer, intent(in) ::  atom1, atom2
    integer :: atid1, atid2, ljt1, ljt2
    real( kind = dp ), intent(in) :: rij, r2, rcut
    real( kind = dp ) :: pot, sw, vpair
    real( kind = dp ), dimension(3,nLocal) :: f    
    real( kind = dp ), intent(in), dimension(3) :: d
    real( kind = dp ), intent(inout), dimension(3) :: fpair
    logical, intent(in) :: do_pot

    ! local Variables
    real( kind = dp ) :: drdx, drdy, drdz
    real( kind = dp ) :: fx, fy, fz
    real( kind = dp ) :: myPot, myPotC, myDeriv, myDerivC, ros, rcos
    real( kind = dp ) :: pot_temp, dudr
    real( kind = dp ) :: sigmai
    real( kind = dp ) :: epsilon
    logical :: isSoftCore, shiftedPot
    integer :: id1, id2, localError

    if (.not.haveMixingMap) then
       call createMixingMap()
    endif

    ! Look up the correct parameters in the mixing matrix
#ifdef IS_MPI
    atid1 = atid_Row(atom1)
    atid2 = atid_Col(atom2)
#else
    atid1 = atid(atom1)
    atid2 = atid(atom2)
#endif

    ljt1 = LJMap%atidToLJtype(atid1)
    ljt2 = LJMap%atidToLJtype(atid2)

    sigmai     = MixingMap(ljt1,ljt2)%sigmai
    epsilon    = MixingMap(ljt1,ljt2)%epsilon
    isSoftCore = MixingMap(ljt1,ljt2)%isSoftCore
    shiftedPot = MixingMap(ljt1,ljt2)%shiftedPot

    ros = rij * sigmai
    myPotC = 0.0_DP

    if (isSoftCore) then

       call getSoftFunc(ros, myPot, myDeriv)

       if (shiftedPot) then
          rcos = rcut * sigmai
          call getSoftFunc(rcos, myPotC, myDerivC)
       endif
              
    else

       call getLJfunc(ros, myPot, myDeriv)

       if (shiftedPot) then
          rcos = rcut * sigmai
          call getLJfunc(rcos, myPotC, myDerivC)
       endif
       
    endif

    !write(*,*) rij, ros, rcos, myPot, myDeriv, myPotC

    pot_temp = epsilon * (myPot - myPotC)
    vpair = vpair + pot_temp
    dudr = sw * epsilon * myDeriv * sigmai

    drdx = d(1) / rij
    drdy = d(2) / rij
    drdz = d(3) / rij

    fx = dudr * drdx
    fy = dudr * drdy
    fz = dudr * drdz

#ifdef IS_MPI
    if (do_pot) then
       pot_Row(VDW_POT,atom1) = pot_Row(VDW_POT,atom1) + sw*pot_temp*0.5
       pot_Col(VDW_POT,atom2) = pot_Col(VDW_POT,atom2) + sw*pot_temp*0.5
    endif

    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       

#else
    if (do_pot) pot = pot + sw*pot_temp

    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
#endif

#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_lj_pair

  subroutine destroyLJTypes()

    LJMap%nLJtypes = 0
    LJMap%currentLJtype = 0
    
    if (associated(LJMap%LJtypes)) then
       deallocate(LJMap%LJtypes)
       LJMap%LJtypes => null()
    end if
    
    if (associated(LJMap%atidToLJtype)) then
       deallocate(LJMap%atidToLJtype)
       LJMap%atidToLJtype => null()
    end if
    
    haveMixingMap = .false.

    call deleteSpline(vLJspline)
    call deleteSpline(vLJpspline) 
    call deleteSpline(vSoftSpline)
    call deleteSpline(vSoftpSpline)

  end subroutine destroyLJTypes

  subroutine getLJfunc(r, myPot, myDeriv)

    real(kind=dp), intent(in) :: r
    real(kind=dp), intent(inout) :: myPot, myDeriv
    real(kind=dp) :: ri, ri2, ri6, ri7, ri12, ri13
    real(kind=dp) :: a, b, c, d, dx
    integer :: j

    if (useSplines) then
       j = MAX(1, MIN(vLJspline%np, idint((r-vLJspline%x(1)) * vLJspline%dx_i) + 1))
       
       dx = r - vLJspline%x(j)
       
       a = vLJspline%c(1,j)
       b = vLJspline%c(2,j)
       c = vLJspline%c(3,j)
       d = vLJspline%c(4,j)
       
       myPot = c + dx * d
       myPot = b + dx * myPot
       myPot = a + dx * myPot

       a = vLJpspline%c(1,j)
       b = vLJpspline%c(2,j)
       c = vLJpspline%c(3,j)
       d = vLJpspline%c(4,j)
       
       myDeriv = c + dx * d
       myDeriv = b + dx * myDeriv  
       myDeriv = a + dx * myDeriv
       
    else
       ri = 1.0_DP / r
       ri2 = ri*ri
       ri6 = ri2*ri2*ri2
       ri7 = ri6*ri
       ri12 = ri6*ri6
       ri13 = ri12*ri
       
       myPot = 4.0_DP * (ri12 - ri6)
       myDeriv = 24.0_DP * (ri7 - 2.0_DP * ri13)
    endif

    return
  end subroutine getLJfunc

  subroutine getSoftFunc(r, myPot, myDeriv)
    
    real(kind=dp), intent(in) :: r
    real(kind=dp), intent(inout) :: myPot, myDeriv
    real(kind=dp) :: ri, ri2, ri6, ri7
    real(kind=dp) :: a, b, c, d, dx
    integer :: j
    
    if (useSplines) then
       j = MAX(1, MIN(vSoftSpline%np, idint((r-vSoftSpline%x(1)) * vSoftSpline%dx_i) + 1))
       
       dx = r - vSoftSpline%x(j)
       
       a = vSoftSpline%c(1,j)
       b = vSoftSpline%c(2,j)
       c = vSoftSpline%c(3,j)
       d = vSoftSpline%c(4,j)
       
       myPot = c + dx * d
       myPot = b + dx * myPot
       myPot = a + dx * myPot

       a = vSoftPspline%c(1,j)
       b = vSoftPspline%c(2,j)
       c = vSoftPspline%c(3,j)
       d = vSoftPspline%c(4,j)
       
       myDeriv = c + dx * d
       myDeriv = b + dx * myDeriv  
       myDeriv = a + dx * myDeriv
       
    else
       ri = 1.0_DP / r    
       ri2 = ri*ri
       ri6 = ri2*ri2*ri2
       ri7 = ri6*ri
       myPot = 4.0_DP * (ri6)
       myDeriv = - 24.0_DP * ri7 
    endif

    return
  end subroutine getSoftFunc

  subroutine setupSplines(rmin, rmax, np)
    real( kind = dp ), intent(in) :: rmin, rmax
    integer, intent(in) :: np
    real( kind = dp ) :: rvals(np), vLJ(np), vLJp(np), vSoft(np), vSoftp(np)
    real( kind = dp ) :: dr, r, ri, ri2, ri6, ri7, ri12, ri13
    real( kind = dp ) :: vljpp1, vljppn, vsoftpp1, vsoftppn
    integer :: i

    dr = (rmax-rmin) / float(np-1)
    
    do i = 1, np
       r = rmin + dble(i-1)*dr
       ri = 1.0_DP / r
       ri2 = ri*ri
       ri6 = ri2*ri2*ri2
       ri7 = ri6*ri
       ri12 = ri6*ri6
       ri13 = ri12*ri

       rvals(i) = r
       vLJ(i) = 4.0_DP * (ri12 - ri6)
       vLJp(i) = 24.0_DP * (ri7 - 2.0_DP * ri13)

       vSoft(i) = 4.0_DP * (ri6)
       vSoftp(i) = - 24.0_DP * ri7       
    enddo

    vljpp1 = 624.0_DP / (rmin)**(14)  - 168.0_DP / (rmin)**(8)
    vljppn = 624.0_DP / (rmax)**(14)  - 168.0_DP / (rmax)**(8)

    vsoftpp1 = 168.0_DP / (rmin)**(8)
    vsoftppn = 168.0_DP / (rmax)**(8)

    call newSpline(vLJspline, rvals, vLJ, vLJp(1), vLJp(np), .true.)
    call newSpline(vLJpspline, rvals, vLJp, vljpp1, vljppn, .true.)
    call newSpline(vSoftSpline, rvals, vSoft, vSoftp(1), vSoftp(np), .true.)
    call newSpline(vSoftpSpline, rvals, vSoftp, vsoftpp1, vsoftppn, .true.)

    return
  end subroutine setupSplines
end module lj
Revision:	2717
Committed:	Mon Apr 17 21:49:12 2006 UTC (18 years, 5 months ago) by gezelter
File size:	16756 byte(s)
Log Message:	Many performance improvements
#	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	1608	!! Calculates Long Range forces Lennard-Jones interactions.
44			!! @author Charles F. Vardeman II
45			!! @author Matthew Meineke
46	gezelter	2717	!! @version $Id: LJ.F90,v 1.21 2006-04-17 21:49:12 gezelter Exp $, $Date: 2006-04-17 21:49:12 $, $Name: not supported by cvs2svn $, $Revision: 1.21 $
47	gezelter	1608
48	gezelter	1930
49	gezelter	1608	module lj
50			use atype_module
51			use vector_class
52			use simulation
53	gezelter	1628	use status
54	chuckv	2497	use fForceOptions
55	gezelter	2717	use interpolation
56	gezelter	1608	#ifdef IS_MPI
57			use mpiSimulation
58			#endif
59			use force_globals
60
61			implicit none
62			PRIVATE
63	chuckv	2355	#define __FORTRAN90
64			#include "UseTheForce/DarkSide/fInteractionMap.h"
65	gezelter	2204
66	chuckv	1624	integer, parameter :: DP = selected_real_kind(15)
67	gezelter	2204
68	gezelter	2271	logical, save :: useGeometricDistanceMixing = .false.
69			logical, save :: haveMixingMap = .false.
70	gezelter	2717	logical, save :: useSplines = .false.
71	gezelter	2271
72			real(kind=DP), save :: defaultCutoff = 0.0_DP
73			logical, save :: defaultShift = .false.
74			logical, save :: haveDefaultCutoff = .false.
75
76			type, private :: LJtype
77			integer :: atid
78	gezelter	1628	real(kind=dp) :: sigma
79			real(kind=dp) :: epsilon
80	gezelter	2271	logical :: isSoftCore = .false.
81			end type LJtype
82	gezelter	2204
83	gezelter	2271	type, private :: LJList
84	chuckv	2319	integer :: Nljtypes = 0
85	gezelter	2271	integer :: currentLJtype = 0
86			type(LJtype), pointer :: LJtypes(:) => null()
87			integer, pointer :: atidToLJtype(:) => null()
88			end type LJList
89	gezelter	2204
90	gezelter	2271	type(LJList), save :: LJMap
91	gezelter	2204
92	gezelter	1628	type :: MixParameters
93			real(kind=DP) :: sigma
94			real(kind=DP) :: epsilon
95	gezelter	2717	real(kind=dp) :: sigmai
96	gezelter	2271	real(kind=dp) :: rCut
97			logical :: rCutWasSet = .false.
98			logical :: shiftedPot
99			logical :: isSoftCore = .false.
100	gezelter	1628	end type MixParameters
101	gezelter	2204
102	gezelter	1628	type(MixParameters), dimension(:,:), allocatable :: MixingMap
103	gezelter	2204
104	gezelter	2717	type(cubicSpline), save :: vLJspline
105			type(cubicSpline), save :: vLJpspline
106			type(cubicSpline), save :: vSoftSpline
107			type(cubicSpline), save :: vSoftpSpline
108
109	gezelter	2271	public :: newLJtype
110			public :: setLJDefaultCutoff
111			public :: getSigma
112			public :: getEpsilon
113	gezelter	1628	public :: do_lj_pair
114	gezelter	2271	public :: destroyLJtypes
115	gezelter	2717	public :: setLJsplineRmax
116	gezelter	2204
117	gezelter	1608	contains
118
119	gezelter	2271	subroutine newLJtype(c_ident, sigma, epsilon, isSoftCore, status)
120	gezelter	1930	integer,intent(in) :: c_ident
121	gezelter	1628	real(kind=dp),intent(in) :: sigma
122			real(kind=dp),intent(in) :: epsilon
123	gezelter	2271	integer, intent(in) :: isSoftCore
124	chuckv	1624	integer,intent(out) :: status
125	gezelter	2271	integer :: nLJTypes, ntypes, myATID
126	gezelter	1930	integer, pointer :: MatchList(:) => null()
127	gezelter	2271	integer :: current
128	gezelter	1628
129	chuckv	1624	status = 0
130	gezelter	2271	! check to see if this is the first time into this routine...
131			if (.not.associated(LJMap%LJtypes)) then
132	gezelter	2204
133	gezelter	2271	call getMatchingElementList(atypes, "is_LennardJones", .true., &
134			nLJTypes, MatchList)
135
136			LJMap%nLJtypes = nLJTypes
137
138			allocate(LJMap%LJtypes(nLJTypes))
139
140			ntypes = getSize(atypes)
141
142			allocate(LJMap%atidToLJtype(ntypes))
143			end if
144
145			LJMap%currentLJtype = LJMap%currentLJtype + 1
146			current = LJMap%currentLJtype
147
148	gezelter	1930	myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
149	gezelter	2271	LJMap%atidToLJtype(myATID) = current
150			LJMap%LJtypes(current)%atid = myATID
151			LJMap%LJtypes(current)%sigma = sigma
152			LJMap%LJtypes(current)%epsilon = epsilon
153			if (isSoftCore .eq. 1) then
154			LJMap%LJtypes(current)%isSoftCore = .true.
155			else
156			LJMap%LJtypes(current)%isSoftCore = .false.
157			endif
158			end subroutine newLJtype
159	chuckv	1624
160	gezelter	2271	subroutine setLJDefaultCutoff(thisRcut, shiftedPot)
161			real(kind=dp), intent(in) :: thisRcut
162			logical, intent(in) :: shiftedPot
163			defaultCutoff = thisRcut
164			defaultShift = shiftedPot
165			haveDefaultCutoff = .true.
166	gezelter	2717	!we only want to build LJ Mixing map and spline if LJ is being used.
167	chuckv	2319	if(LJMap%nLJTypes /= 0) then
168			call createMixingMap()
169	gezelter	2717	call setLJsplineRmax(defaultCutoff)
170	chuckv	2319	end if
171	gezelter	2717
172	gezelter	2271	end subroutine setLJDefaultCutoff
173	gezelter	2204
174	gezelter	1633	function getSigma(atid) result (s)
175			integer, intent(in) :: atid
176	gezelter	2271	integer :: ljt1
177	gezelter	1633	real(kind=dp) :: s
178	gezelter	2204
179	gezelter	2271	if (LJMap%currentLJtype == 0) then
180			call handleError("LJ", "No members in LJMap")
181	gezelter	1633	return
182			end if
183	gezelter	2204
184	gezelter	2271	ljt1 = LJMap%atidToLJtype(atid)
185			s = LJMap%LJtypes(ljt1)%sigma
186
187	gezelter	1633	end function getSigma
188
189			function getEpsilon(atid) result (e)
190			integer, intent(in) :: atid
191	gezelter	2271	integer :: ljt1
192	gezelter	1633	real(kind=dp) :: e
193	gezelter	2204
194	gezelter	2271	if (LJMap%currentLJtype == 0) then
195			call handleError("LJ", "No members in LJMap")
196	gezelter	1633	return
197			end if
198	gezelter	2204
199	gezelter	2271	ljt1 = LJMap%atidToLJtype(atid)
200			e = LJMap%LJtypes(ljt1)%epsilon
201
202	gezelter	1633	end function getEpsilon
203
204	gezelter	2271	subroutine createMixingMap()
205			integer :: nLJtypes, i, j
206			real ( kind = dp ) :: s1, s2, e1, e2
207			real ( kind = dp ) :: rcut6, tp6, tp12
208	gezelter	2289	logical :: isSoftCore1, isSoftCore2, doShift
209	gezelter	1628
210	gezelter	2271	if (LJMap%currentLJtype == 0) then
211			call handleError("LJ", "No members in LJMap")
212	gezelter	2086	return
213	gezelter	2271	end if
214	gezelter	2204
215	gezelter	2271	nLJtypes = LJMap%nLJtypes
216	gezelter	2204
217	gezelter	2086	if (.not. allocated(MixingMap)) then
218	gezelter	2271	allocate(MixingMap(nLJtypes, nLJtypes))
219	gezelter	2086	endif
220
221	chuckv	2497	useGeometricDistanceMixing = usesGeometricDistanceMixing()
222	gezelter	2271	do i = 1, nLJtypes
223	gezelter	2204
224	gezelter	2271	s1 = LJMap%LJtypes(i)%sigma
225			e1 = LJMap%LJtypes(i)%epsilon
226			isSoftCore1 = LJMap%LJtypes(i)%isSoftCore
227	gezelter	2204
228	gezelter	2271	do j = i, nLJtypes
229
230			s2 = LJMap%LJtypes(j)%sigma
231			e2 = LJMap%LJtypes(j)%epsilon
232			isSoftCore2 = LJMap%LJtypes(j)%isSoftCore
233
234			MixingMap(i,j)%isSoftCore = isSoftCore1 .or. isSoftCore2
235	gezelter	2204
236	gezelter	2271	! only the distance parameter uses different mixing policies
237			if (useGeometricDistanceMixing) then
238			MixingMap(i,j)%sigma = dsqrt(s1 * s2)
239			else
240			MixingMap(i,j)%sigma = 0.5_dp * (s1 + s2)
241			endif
242
243			MixingMap(i,j)%epsilon = dsqrt(e1 * e2)
244	gezelter	2204
245	gezelter	2717	MixingMap(i,j)%sigmai = 1.0_DP / (MixingMap(i,j)%sigma)
246	gezelter	2204
247	gezelter	2461	if (haveDefaultCutoff) then
248			MixingMap(i,j)%shiftedPot = defaultShift
249	gezelter	2271	else
250	gezelter	2461	MixingMap(i,j)%shiftedPot = defaultShift
251			endif
252	gezelter	2204
253	gezelter	2289	if (i.ne.j) then
254			MixingMap(j,i)%sigma = MixingMap(i,j)%sigma
255			MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
256	gezelter	2717	MixingMap(j,i)%sigmai = MixingMap(i,j)%sigmai
257	gezelter	2289	MixingMap(j,i)%rCut = MixingMap(i,j)%rCut
258			MixingMap(j,i)%rCutWasSet = MixingMap(i,j)%rCutWasSet
259			MixingMap(j,i)%shiftedPot = MixingMap(i,j)%shiftedPot
260			MixingMap(j,i)%isSoftCore = MixingMap(i,j)%isSoftCore
261			endif
262
263	gezelter	2271	enddo
264			enddo
265
266	chrisfen	1754	haveMixingMap = .true.
267	gezelter	2271
268	gezelter	1628	end subroutine createMixingMap
269	gezelter	2271
270	gezelter	2717	subroutine setLJsplineRmax(largestRcut)
271			real( kind = dp ), intent(in) :: largestRcut
272			real( kind = dp ) :: s, bigS, smallS, rmax, rmin
273			integer :: np, i
274
275			if (LJMap%nLJtypes .ne. 0) then
276
277			!
278			! find the largest and smallest values of sigma that we'll need
279			!
280			bigS = 0.0_DP
281			smallS = 1.0e9
282			do i = 1, LJMap%nLJtypes
283			s = LJMap%LJtypes(i)%sigma
284			if (s .gt. bigS) bigS = s
285			if (s .lt. smallS) smallS = s
286			enddo
287
288			!
289			! give ourselves a 20% margin just in case
290			!
291			rmax = 1.2 * largestRcut / smallS
292			!
293			! assume atoms will never get closer than 1 angstrom
294			!
295			rmin = 1 / bigS
296			!
297			! assume 500 points is enough
298			!
299			np = 500
300
301			write(,) 'calling setupSplines with rmin = ', rmin, ' rmax = ', rmax, &
302			' np = ', np
303
304			call setupSplines(rmin, rmax, np)
305
306			endif
307			return
308			end subroutine setLJsplineRmax
309
310	gezelter	2461	subroutine do_lj_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
311	gezelter	1608	pot, f, do_pot)
312	gezelter	2461
313	gezelter	1608	integer, intent(in) :: atom1, atom2
314	gezelter	2271	integer :: atid1, atid2, ljt1, ljt2
315	gezelter	2461	real( kind = dp ), intent(in) :: rij, r2, rcut
316	gezelter	1608	real( kind = dp ) :: pot, sw, vpair
317			real( kind = dp ), dimension(3,nLocal) :: f
318			real( kind = dp ), intent(in), dimension(3) :: d
319			real( kind = dp ), intent(inout), dimension(3) :: fpair
320			logical, intent(in) :: do_pot
321
322			! local Variables
323			real( kind = dp ) :: drdx, drdy, drdz
324			real( kind = dp ) :: fx, fy, fz
325	gezelter	2717	real( kind = dp ) :: myPot, myPotC, myDeriv, myDerivC, ros, rcos
326	gezelter	1608	real( kind = dp ) :: pot_temp, dudr
327	gezelter	2717	real( kind = dp ) :: sigmai
328	gezelter	1608	real( kind = dp ) :: epsilon
329	gezelter	2271	logical :: isSoftCore, shiftedPot
330	gezelter	1628	integer :: id1, id2, localError
331	gezelter	1608
332	gezelter	1628	if (.not.haveMixingMap) then
333	gezelter	2271	call createMixingMap()
334	gezelter	1628	endif
335
336	gezelter	1608	! Look up the correct parameters in the mixing matrix
337			#ifdef IS_MPI
338	gezelter	2271	atid1 = atid_Row(atom1)
339			atid2 = atid_Col(atom2)
340	gezelter	1608	#else
341	gezelter	2271	atid1 = atid(atom1)
342			atid2 = atid(atom2)
343	gezelter	1608	#endif
344
345	gezelter	2271	ljt1 = LJMap%atidToLJtype(atid1)
346			ljt2 = LJMap%atidToLJtype(atid2)
347
348	gezelter	2717	sigmai = MixingMap(ljt1,ljt2)%sigmai
349	gezelter	2271	epsilon = MixingMap(ljt1,ljt2)%epsilon
350			isSoftCore = MixingMap(ljt1,ljt2)%isSoftCore
351			shiftedPot = MixingMap(ljt1,ljt2)%shiftedPot
352
353	gezelter	2717	ros = rij * sigmai
354			myPotC = 0.0_DP
355	gezelter	2204
356	gezelter	2717	if (isSoftCore) then
357	tim	2062
358	gezelter	2717	call getSoftFunc(ros, myPot, myDeriv)
359
360	gezelter	2271	if (shiftedPot) then
361	gezelter	2717	rcos = rcut * sigmai
362			call getSoftFunc(rcos, myPotC, myDerivC)
363	gezelter	2271	endif
364	gezelter	2717
365	tim	2062	else
366	gezelter	2717
367			call getLJfunc(ros, myPot, myDeriv)
368
369	gezelter	2271	if (shiftedPot) then
370	gezelter	2717	rcos = rcut * sigmai
371			call getLJfunc(rcos, myPotC, myDerivC)
372	gezelter	2204	endif
373	gezelter	2271
374	gezelter	1608	endif
375	gezelter	2204
376	gezelter	2717	!write(,) rij, ros, rcos, myPot, myDeriv, myPotC
377
378			pot_temp = epsilon * (myPot - myPotC)
379			vpair = vpair + pot_temp
380			dudr = sw * epsilon * myDeriv * sigmai
381
382	gezelter	1608	drdx = d(1) / rij
383			drdy = d(2) / rij
384			drdz = d(3) / rij
385	gezelter	2204
386	gezelter	1608	fx = dudr * drdx
387			fy = dudr * drdy
388			fz = dudr * drdz
389	gezelter	2204
390	gezelter	1608	#ifdef IS_MPI
391			if (do_pot) then
392	gezelter	2361	pot_Row(VDW_POT,atom1) = pot_Row(VDW_POT,atom1) + swpot_temp0.5
393			pot_Col(VDW_POT,atom2) = pot_Col(VDW_POT,atom2) + swpot_temp0.5
394	gezelter	1608	endif
395	gezelter	2204
396	gezelter	1608	f_Row(1,atom1) = f_Row(1,atom1) + fx
397			f_Row(2,atom1) = f_Row(2,atom1) + fy
398			f_Row(3,atom1) = f_Row(3,atom1) + fz
399	gezelter	2204
400	gezelter	1608	f_Col(1,atom2) = f_Col(1,atom2) - fx
401			f_Col(2,atom2) = f_Col(2,atom2) - fy
402			f_Col(3,atom2) = f_Col(3,atom2) - fz
403	gezelter	2204
404	gezelter	1608	#else
405			if (do_pot) pot = pot + sw*pot_temp
406
407			f(1,atom1) = f(1,atom1) + fx
408			f(2,atom1) = f(2,atom1) + fy
409			f(3,atom1) = f(3,atom1) + fz
410	gezelter	2204
411	gezelter	1608	f(1,atom2) = f(1,atom2) - fx
412			f(2,atom2) = f(2,atom2) - fy
413			f(3,atom2) = f(3,atom2) - fz
414			#endif
415	gezelter	2204
416	gezelter	1608	#ifdef IS_MPI
417			id1 = AtomRowToGlobal(atom1)
418			id2 = AtomColToGlobal(atom2)
419			#else
420			id1 = atom1
421			id2 = atom2
422			#endif
423
424			if (molMembershipList(id1) .ne. molMembershipList(id2)) then
425	gezelter	2204
426	gezelter	1608	fpair(1) = fpair(1) + fx
427			fpair(2) = fpair(2) + fy
428			fpair(3) = fpair(3) + fz
429
430			endif
431
432			return
433	gezelter	2204
434	gezelter	1608	end subroutine do_lj_pair
435	gezelter	2204
436	chuckv	2189	subroutine destroyLJTypes()
437	gezelter	2271
438			LJMap%nLJtypes = 0
439			LJMap%currentLJtype = 0
440
441			if (associated(LJMap%LJtypes)) then
442			deallocate(LJMap%LJtypes)
443			LJMap%LJtypes => null()
444			end if
445
446			if (associated(LJMap%atidToLJtype)) then
447			deallocate(LJMap%atidToLJtype)
448			LJMap%atidToLJtype => null()
449			end if
450
451	chuckv	2189	haveMixingMap = .false.
452	gezelter	2717
453			call deleteSpline(vLJspline)
454			call deleteSpline(vLJpspline)
455			call deleteSpline(vSoftSpline)
456			call deleteSpline(vSoftpSpline)
457
458	chuckv	2189	end subroutine destroyLJTypes
459
460	gezelter	2717	subroutine getLJfunc(r, myPot, myDeriv)
461
462			real(kind=dp), intent(in) :: r
463			real(kind=dp), intent(inout) :: myPot, myDeriv
464			real(kind=dp) :: ri, ri2, ri6, ri7, ri12, ri13
465			real(kind=dp) :: a, b, c, d, dx
466			integer :: j
467
468			if (useSplines) then
469			j = MAX(1, MIN(vLJspline%np, idint((r-vLJspline%x(1)) * vLJspline%dx_i) + 1))
470
471			dx = r - vLJspline%x(j)
472
473			a = vLJspline%c(1,j)
474			b = vLJspline%c(2,j)
475			c = vLJspline%c(3,j)
476			d = vLJspline%c(4,j)
477
478			myPot = c + dx * d
479			myPot = b + dx * myPot
480			myPot = a + dx * myPot
481
482			a = vLJpspline%c(1,j)
483			b = vLJpspline%c(2,j)
484			c = vLJpspline%c(3,j)
485			d = vLJpspline%c(4,j)
486
487			myDeriv = c + dx * d
488			myDeriv = b + dx * myDeriv
489			myDeriv = a + dx * myDeriv
490
491			else
492			ri = 1.0_DP / r
493			ri2 = ri*ri
494			ri6 = ri2ri2ri2
495			ri7 = ri6*ri
496			ri12 = ri6*ri6
497			ri13 = ri12*ri
498
499			myPot = 4.0_DP * (ri12 - ri6)
500			myDeriv = 24.0_DP * (ri7 - 2.0_DP * ri13)
501			endif
502
503			return
504			end subroutine getLJfunc
505
506			subroutine getSoftFunc(r, myPot, myDeriv)
507
508			real(kind=dp), intent(in) :: r
509			real(kind=dp), intent(inout) :: myPot, myDeriv
510			real(kind=dp) :: ri, ri2, ri6, ri7
511			real(kind=dp) :: a, b, c, d, dx
512			integer :: j
513
514			if (useSplines) then
515			j = MAX(1, MIN(vSoftSpline%np, idint((r-vSoftSpline%x(1)) * vSoftSpline%dx_i) + 1))
516
517			dx = r - vSoftSpline%x(j)
518
519			a = vSoftSpline%c(1,j)
520			b = vSoftSpline%c(2,j)
521			c = vSoftSpline%c(3,j)
522			d = vSoftSpline%c(4,j)
523
524			myPot = c + dx * d
525			myPot = b + dx * myPot
526			myPot = a + dx * myPot
527
528			a = vSoftPspline%c(1,j)
529			b = vSoftPspline%c(2,j)
530			c = vSoftPspline%c(3,j)
531			d = vSoftPspline%c(4,j)
532
533			myDeriv = c + dx * d
534			myDeriv = b + dx * myDeriv
535			myDeriv = a + dx * myDeriv
536
537			else
538			ri = 1.0_DP / r
539			ri2 = ri*ri
540			ri6 = ri2ri2ri2
541			ri7 = ri6*ri
542			myPot = 4.0_DP * (ri6)
543			myDeriv = - 24.0_DP * ri7
544			endif
545
546			return
547			end subroutine getSoftFunc
548
549			subroutine setupSplines(rmin, rmax, np)
550			real( kind = dp ), intent(in) :: rmin, rmax
551			integer, intent(in) :: np
552			real( kind = dp ) :: rvals(np), vLJ(np), vLJp(np), vSoft(np), vSoftp(np)
553			real( kind = dp ) :: dr, r, ri, ri2, ri6, ri7, ri12, ri13
554			real( kind = dp ) :: vljpp1, vljppn, vsoftpp1, vsoftppn
555			integer :: i
556
557			dr = (rmax-rmin) / float(np-1)
558
559			do i = 1, np
560			r = rmin + dble(i-1)*dr
561			ri = 1.0_DP / r
562			ri2 = ri*ri
563			ri6 = ri2ri2ri2
564			ri7 = ri6*ri
565			ri12 = ri6*ri6
566			ri13 = ri12*ri
567
568			rvals(i) = r
569			vLJ(i) = 4.0_DP * (ri12 - ri6)
570			vLJp(i) = 24.0_DP * (ri7 - 2.0_DP * ri13)
571
572			vSoft(i) = 4.0_DP * (ri6)
573			vSoftp(i) = - 24.0_DP * ri7
574			enddo
575
576			vljpp1 = 624.0_DP / (rmin)(14) - 168.0_DP / (rmin)(8)
577			vljppn = 624.0_DP / (rmax)(14) - 168.0_DP / (rmax)(8)
578
579			vsoftpp1 = 168.0_DP / (rmin)**(8)
580			vsoftppn = 168.0_DP / (rmax)**(8)
581
582			call newSpline(vLJspline, rvals, vLJ, vLJp(1), vLJp(np), .true.)
583			call newSpline(vLJpspline, rvals, vLJp, vljpp1, vljppn, .true.)
584			call newSpline(vSoftSpline, rvals, vSoft, vSoftp(1), vSoftp(np), .true.)
585			call newSpline(vSoftpSpline, rvals, vSoftp, vsoftpp1, vsoftppn, .true.)
586
587			return
588			end subroutine setupSplines
589	gezelter	1608	end module lj