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.15 2005-09-21 17:20:14 chrisfen Exp $, $Date: 2005-09-21 17:20:14 $, $Name: not supported by cvs2svn $, $Revision: 1.15 $


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

  implicit none
  PRIVATE

  integer, parameter :: DP = selected_real_kind(15)

  logical, save :: useGeometricDistanceMixing = .false.
  logical, save :: haveMixingMap = .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) :: sigma6
     real(kind=dp) :: rCut
     real(kind=dp) :: delta
     logical       :: rCutWasSet = .false.
     logical       :: shiftedPot
     logical       :: isSoftCore = .false.
  end type MixParameters

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

  public :: newLJtype
  public :: setLJDefaultCutoff
  public :: setLJUniformCutoff
  public :: setLJCutoffByTypes
  public :: getSigma
  public :: getEpsilon
  public :: useGeometricMixing
  public :: do_lj_pair
  public :: destroyLJtypes

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.
    call createMixingMap()
  end subroutine setLJDefaultCutoff

  subroutine setLJUniformCutoff(thisRcut, shiftedPot)
    real(kind=dp), intent(in) :: thisRcut
    logical, intent(in) :: shiftedPot
    integer :: nLJtypes, i, j

    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

    do i = 1, nLJtypes
       do j = 1, nLJtypes
          MixingMap(i,j)%rCut = thisRcut
          MixingMap(i,j)%shiftedPot = shiftedPot
          MixingMap(i,j)%rCutWasSet = .true.
       enddo
    enddo
    call createMixingMap()
  end subroutine setLJUniformCutoff

  subroutine setLJCutoffByTypes(atid1, atid2, thisRcut, shiftedPot)
    integer, intent(in) :: atid1, atid2
    real(kind=dp), intent(in) :: thisRcut
    logical, intent(in) :: shiftedPot
    integer :: nLJtypes, ljt1, ljt2

    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

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

    MixingMap(ljt1,ljt2)%rCut = thisRcut
    MixingMap(ljt1,ljt2)%shiftedPot = shiftedPot
    MixingMap(ljt1,ljt2)%rCutWasSet = .true.
    MixingMap(ljt2,ljt1)%rCut = thisRcut
    MixingMap(ljt2,ljt1)%shiftedPot = shiftedPot
    MixingMap(ljt2,ljt1)%rCutWasSet = .true.

    call createMixingMap()
  end subroutine setLJCutoffByTypes

  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 useGeometricMixing() 
    useGeometricDistanceMixing = .true.
    haveMixingMap = .false.
    return
  end subroutine useGeometricMixing

  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

    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)%sigma6 = (MixingMap(i,j)%sigma)**6

          if (MixingMap(i,j)%rCutWasSet) then
             rcut6 = (MixingMap(i,j)%rcut)**6
          else
             if (haveDefaultCutoff) then
                rcut6 = defaultCutoff**6
                doShift = defaultShift
             else
                call handleError("LJ", "No specified or default cutoff value!")
             endif
          endif
          
          tp6    = MixingMap(i,j)%sigma6/rcut6
          tp12    = tp6**2          
          MixingMap(i,j)%delta =-4.0_DP*MixingMap(i,j)%epsilon*(tp12 - tp6)
          MixingMap(i,j)%shiftedPot = doShift

          if (i.ne.j) then
             MixingMap(j,i)%sigma      = MixingMap(i,j)%sigma
             MixingMap(j,i)%epsilon    = MixingMap(i,j)%epsilon
             MixingMap(j,i)%sigma6     = MixingMap(i,j)%sigma6
             MixingMap(j,i)%rCut       = MixingMap(i,j)%rCut
             MixingMap(j,i)%delta      = MixingMap(i,j)%delta
             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 do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
       pot, f, do_pot)

    integer, intent(in) ::  atom1, atom2
    integer :: atid1, atid2, ljt1, ljt2
    real( kind = dp ), intent(in) :: rij, r2
    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 ) :: pot_temp, dudr
    real( kind = dp ) :: sigma6
    real( kind = dp ) :: epsilon
    real( kind = dp ) :: r6
    real( kind = dp ) :: t6
    real( kind = dp ) :: t12
    real( kind = dp ) :: delta
    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)

    sigma6     = MixingMap(ljt1,ljt2)%sigma6
    epsilon    = MixingMap(ljt1,ljt2)%epsilon
    delta      = MixingMap(ljt1,ljt2)%delta
    isSoftCore = MixingMap(ljt1,ljt2)%isSoftCore
    shiftedPot = MixingMap(ljt1,ljt2)%shiftedPot

    r6 = r2 * r2 * r2

    t6  = sigma6/ r6
    t12 = t6 * t6     

    if (isSoftCore) then
       
       pot_temp = 4.0E0_DP * epsilon * t6
       if (shiftedPot) then
          pot_temp = pot_temp + delta
       endif
       
       vpair = vpair + pot_temp
       
       dudr = -sw * 24.0E0_DP * epsilon * t6 / rij
       
    else
       pot_temp = 4.0E0_DP * epsilon * (t12 - t6) 
       if (shiftedPot) then
          pot_temp = pot_temp + delta
       endif
       
       vpair = vpair + pot_temp
       
       dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
    endif

    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(atom1) = pot_Row(atom1) + sw*pot_temp*0.5
       pot_Col(atom2) = pot_Col(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.
  end subroutine destroyLJTypes

end module lj
Revision:	2317
Committed:	Wed Sep 21 17:20:14 2005 UTC (18 years, 11 months ago) by chrisfen
File size:	13871 byte(s)
Log Message:	Fixed a defaultCutoff bug (HEMES!)
#	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	chrisfen	2317	!! @version $Id: LJ.F90,v 1.15 2005-09-21 17:20:14 chrisfen Exp $, $Date: 2005-09-21 17:20:14 $, $Name: not supported by cvs2svn $, $Revision: 1.15 $
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	gezelter	1608	#ifdef IS_MPI
55			use mpiSimulation
56			#endif
57			use force_globals
58
59			implicit none
60			PRIVATE
61	gezelter	2204
62	chuckv	1624	integer, parameter :: DP = selected_real_kind(15)
63	gezelter	2204
64	gezelter	2271	logical, save :: useGeometricDistanceMixing = .false.
65			logical, save :: haveMixingMap = .false.
66
67			real(kind=DP), save :: defaultCutoff = 0.0_DP
68			logical, save :: defaultShift = .false.
69			logical, save :: haveDefaultCutoff = .false.
70
71
72			type, private :: LJtype
73			integer :: atid
74	gezelter	1628	real(kind=dp) :: sigma
75			real(kind=dp) :: epsilon
76	gezelter	2271	logical :: isSoftCore = .false.
77			end type LJtype
78	gezelter	2204
79	gezelter	2271	type, private :: LJList
80			integer :: nLJtypes = 0
81			integer :: currentLJtype = 0
82			type(LJtype), pointer :: LJtypes(:) => null()
83			integer, pointer :: atidToLJtype(:) => null()
84			end type LJList
85	gezelter	2204
86	gezelter	2271	type(LJList), save :: LJMap
87	gezelter	2204
88	gezelter	1628	type :: MixParameters
89			real(kind=DP) :: sigma
90			real(kind=DP) :: epsilon
91	gezelter	2271	real(kind=dp) :: sigma6
92			real(kind=dp) :: rCut
93			real(kind=dp) :: delta
94			logical :: rCutWasSet = .false.
95			logical :: shiftedPot
96			logical :: isSoftCore = .false.
97	gezelter	1628	end type MixParameters
98	gezelter	2204
99	gezelter	1628	type(MixParameters), dimension(:,:), allocatable :: MixingMap
100	gezelter	2204
101	gezelter	2271	public :: newLJtype
102			public :: setLJDefaultCutoff
103			public :: setLJUniformCutoff
104			public :: setLJCutoffByTypes
105			public :: getSigma
106			public :: getEpsilon
107	gezelter	1628	public :: useGeometricMixing
108			public :: do_lj_pair
109	gezelter	2271	public :: destroyLJtypes
110	gezelter	2204
111	gezelter	1608	contains
112
113	gezelter	2271	subroutine newLJtype(c_ident, sigma, epsilon, isSoftCore, status)
114	gezelter	1930	integer,intent(in) :: c_ident
115	gezelter	1628	real(kind=dp),intent(in) :: sigma
116			real(kind=dp),intent(in) :: epsilon
117	gezelter	2271	integer, intent(in) :: isSoftCore
118	chuckv	1624	integer,intent(out) :: status
119	gezelter	2271	integer :: nLJTypes, ntypes, myATID
120	gezelter	1930	integer, pointer :: MatchList(:) => null()
121	gezelter	2271	integer :: current
122	gezelter	1628
123	chuckv	1624	status = 0
124	gezelter	2271	! check to see if this is the first time into this routine...
125			if (.not.associated(LJMap%LJtypes)) then
126	gezelter	2204
127	gezelter	2271	call getMatchingElementList(atypes, "is_LennardJones", .true., &
128			nLJTypes, MatchList)
129
130			LJMap%nLJtypes = nLJTypes
131
132			allocate(LJMap%LJtypes(nLJTypes))
133
134			ntypes = getSize(atypes)
135
136			allocate(LJMap%atidToLJtype(ntypes))
137			end if
138
139			LJMap%currentLJtype = LJMap%currentLJtype + 1
140			current = LJMap%currentLJtype
141
142	gezelter	1930	myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
143	gezelter	2271	LJMap%atidToLJtype(myATID) = current
144			LJMap%LJtypes(current)%atid = myATID
145			LJMap%LJtypes(current)%sigma = sigma
146			LJMap%LJtypes(current)%epsilon = epsilon
147			if (isSoftCore .eq. 1) then
148			LJMap%LJtypes(current)%isSoftCore = .true.
149			else
150			LJMap%LJtypes(current)%isSoftCore = .false.
151			endif
152			end subroutine newLJtype
153	chuckv	1624
154	gezelter	2271	subroutine setLJDefaultCutoff(thisRcut, shiftedPot)
155			real(kind=dp), intent(in) :: thisRcut
156			logical, intent(in) :: shiftedPot
157			defaultCutoff = thisRcut
158			defaultShift = shiftedPot
159			haveDefaultCutoff = .true.
160	chrisfen	2317	call createMixingMap()
161	gezelter	2271	end subroutine setLJDefaultCutoff
162	gezelter	2204
163	gezelter	2271	subroutine setLJUniformCutoff(thisRcut, shiftedPot)
164			real(kind=dp), intent(in) :: thisRcut
165			logical, intent(in) :: shiftedPot
166			integer :: nLJtypes, i, j
167	gezelter	2204
168	gezelter	2271	if (LJMap%currentLJtype == 0) then
169			call handleError("LJ", "No members in LJMap")
170			return
171	chuckv	1624	end if
172
173	gezelter	2271	nLJtypes = LJMap%nLJtypes
174			if (.not. allocated(MixingMap)) then
175			allocate(MixingMap(nLJtypes, nLJtypes))
176	gezelter	1628	endif
177	gezelter	2204
178	gezelter	2271	do i = 1, nLJtypes
179			do j = 1, nLJtypes
180			MixingMap(i,j)%rCut = thisRcut
181			MixingMap(i,j)%shiftedPot = shiftedPot
182			MixingMap(i,j)%rCutWasSet = .true.
183			enddo
184			enddo
185			call createMixingMap()
186			end subroutine setLJUniformCutoff
187	gezelter	1628
188	gezelter	2271	subroutine setLJCutoffByTypes(atid1, atid2, thisRcut, shiftedPot)
189			integer, intent(in) :: atid1, atid2
190			real(kind=dp), intent(in) :: thisRcut
191			logical, intent(in) :: shiftedPot
192			integer :: nLJtypes, ljt1, ljt2
193
194			if (LJMap%currentLJtype == 0) then
195			call handleError("LJ", "No members in LJMap")
196			return
197			end if
198
199			nLJtypes = LJMap%nLJtypes
200			if (.not. allocated(MixingMap)) then
201			allocate(MixingMap(nLJtypes, nLJtypes))
202	tim	2062	endif
203	gezelter	1633
204	gezelter	2271	ljt1 = LJMap%atidToLJtype(atid1)
205			ljt2 = LJMap%atidToLJtype(atid2)
206
207			MixingMap(ljt1,ljt2)%rCut = thisRcut
208			MixingMap(ljt1,ljt2)%shiftedPot = shiftedPot
209			MixingMap(ljt1,ljt2)%rCutWasSet = .true.
210			MixingMap(ljt2,ljt1)%rCut = thisRcut
211			MixingMap(ljt2,ljt1)%shiftedPot = shiftedPot
212			MixingMap(ljt2,ljt1)%rCutWasSet = .true.
213
214			call createMixingMap()
215			end subroutine setLJCutoffByTypes
216
217	gezelter	1633	function getSigma(atid) result (s)
218			integer, intent(in) :: atid
219	gezelter	2271	integer :: ljt1
220	gezelter	1633	real(kind=dp) :: s
221	gezelter	2204
222	gezelter	2271	if (LJMap%currentLJtype == 0) then
223			call handleError("LJ", "No members in LJMap")
224	gezelter	1633	return
225			end if
226	gezelter	2204
227	gezelter	2271	ljt1 = LJMap%atidToLJtype(atid)
228			s = LJMap%LJtypes(ljt1)%sigma
229
230	gezelter	1633	end function getSigma
231
232			function getEpsilon(atid) result (e)
233			integer, intent(in) :: atid
234	gezelter	2271	integer :: ljt1
235	gezelter	1633	real(kind=dp) :: e
236	gezelter	2204
237	gezelter	2271	if (LJMap%currentLJtype == 0) then
238			call handleError("LJ", "No members in LJMap")
239	gezelter	1633	return
240			end if
241	gezelter	2204
242	gezelter	2271	ljt1 = LJMap%atidToLJtype(atid)
243			e = LJMap%LJtypes(ljt1)%epsilon
244
245	gezelter	1633	end function getEpsilon
246
247	gezelter	1628	subroutine useGeometricMixing()
248			useGeometricDistanceMixing = .true.
249			haveMixingMap = .false.
250			return
251			end subroutine useGeometricMixing
252	gezelter	2204
253	gezelter	2271	subroutine createMixingMap()
254			integer :: nLJtypes, i, j
255			real ( kind = dp ) :: s1, s2, e1, e2
256			real ( kind = dp ) :: rcut6, tp6, tp12
257	gezelter	2289	logical :: isSoftCore1, isSoftCore2, doShift
258	gezelter	1628
259	gezelter	2271	if (LJMap%currentLJtype == 0) then
260			call handleError("LJ", "No members in LJMap")
261	gezelter	2086	return
262	gezelter	2271	end if
263	gezelter	2204
264	gezelter	2271	nLJtypes = LJMap%nLJtypes
265	gezelter	2204
266	gezelter	2086	if (.not. allocated(MixingMap)) then
267	gezelter	2271	allocate(MixingMap(nLJtypes, nLJtypes))
268	gezelter	2086	endif
269
270	gezelter	2271	do i = 1, nLJtypes
271	gezelter	2204
272	gezelter	2271	s1 = LJMap%LJtypes(i)%sigma
273			e1 = LJMap%LJtypes(i)%epsilon
274			isSoftCore1 = LJMap%LJtypes(i)%isSoftCore
275	gezelter	2204
276	gezelter	2271	do j = i, nLJtypes
277
278			s2 = LJMap%LJtypes(j)%sigma
279			e2 = LJMap%LJtypes(j)%epsilon
280			isSoftCore2 = LJMap%LJtypes(j)%isSoftCore
281
282			MixingMap(i,j)%isSoftCore = isSoftCore1 .or. isSoftCore2
283	gezelter	2204
284	gezelter	2271	! only the distance parameter uses different mixing policies
285			if (useGeometricDistanceMixing) then
286			MixingMap(i,j)%sigma = dsqrt(s1 * s2)
287			else
288			MixingMap(i,j)%sigma = 0.5_dp * (s1 + s2)
289			endif
290
291			MixingMap(i,j)%epsilon = dsqrt(e1 * e2)
292	gezelter	2204
293	gezelter	2271	MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
294	gezelter	2204
295	gezelter	2271	if (MixingMap(i,j)%rCutWasSet) then
296			rcut6 = (MixingMap(i,j)%rcut)**6
297			else
298			if (haveDefaultCutoff) then
299			rcut6 = defaultCutoff**6
300	gezelter	2289	doShift = defaultShift
301	gezelter	2271	else
302			call handleError("LJ", "No specified or default cutoff value!")
303	gezelter	2086	endif
304	gezelter	2271	endif
305
306			tp6 = MixingMap(i,j)%sigma6/rcut6
307			tp12 = tp6**2
308	gezelter	2289	MixingMap(i,j)%delta =-4.0_DPMixingMap(i,j)%epsilon(tp12 - tp6)
309			MixingMap(i,j)%shiftedPot = doShift
310	gezelter	2204
311	gezelter	2289	if (i.ne.j) then
312			MixingMap(j,i)%sigma = MixingMap(i,j)%sigma
313			MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
314			MixingMap(j,i)%sigma6 = MixingMap(i,j)%sigma6
315			MixingMap(j,i)%rCut = MixingMap(i,j)%rCut
316			MixingMap(j,i)%delta = MixingMap(i,j)%delta
317			MixingMap(j,i)%rCutWasSet = MixingMap(i,j)%rCutWasSet
318			MixingMap(j,i)%shiftedPot = MixingMap(i,j)%shiftedPot
319			MixingMap(j,i)%isSoftCore = MixingMap(i,j)%isSoftCore
320			endif
321
322	gezelter	2271	enddo
323			enddo
324
325	chrisfen	1754	haveMixingMap = .true.
326	gezelter	2271
327	gezelter	1628	end subroutine createMixingMap
328	gezelter	2271
329	gezelter	1608	subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
330			pot, f, do_pot)
331
332			integer, intent(in) :: atom1, atom2
333	gezelter	2271	integer :: atid1, atid2, ljt1, ljt2
334	gezelter	1608	real( kind = dp ), intent(in) :: rij, r2
335			real( kind = dp ) :: pot, sw, vpair
336			real( kind = dp ), dimension(3,nLocal) :: f
337			real( kind = dp ), intent(in), dimension(3) :: d
338			real( kind = dp ), intent(inout), dimension(3) :: fpair
339			logical, intent(in) :: do_pot
340
341			! local Variables
342			real( kind = dp ) :: drdx, drdy, drdz
343			real( kind = dp ) :: fx, fy, fz
344			real( kind = dp ) :: pot_temp, dudr
345			real( kind = dp ) :: sigma6
346			real( kind = dp ) :: epsilon
347			real( kind = dp ) :: r6
348			real( kind = dp ) :: t6
349			real( kind = dp ) :: t12
350			real( kind = dp ) :: delta
351	gezelter	2271	logical :: isSoftCore, shiftedPot
352	gezelter	1628	integer :: id1, id2, localError
353	gezelter	1608
354	gezelter	1628	if (.not.haveMixingMap) then
355	gezelter	2271	call createMixingMap()
356	gezelter	1628	endif
357
358	gezelter	1608	! Look up the correct parameters in the mixing matrix
359			#ifdef IS_MPI
360	gezelter	2271	atid1 = atid_Row(atom1)
361			atid2 = atid_Col(atom2)
362	gezelter	1608	#else
363	gezelter	2271	atid1 = atid(atom1)
364			atid2 = atid(atom2)
365	gezelter	1608	#endif
366
367	gezelter	2271	ljt1 = LJMap%atidToLJtype(atid1)
368			ljt2 = LJMap%atidToLJtype(atid2)
369
370			sigma6 = MixingMap(ljt1,ljt2)%sigma6
371			epsilon = MixingMap(ljt1,ljt2)%epsilon
372			delta = MixingMap(ljt1,ljt2)%delta
373			isSoftCore = MixingMap(ljt1,ljt2)%isSoftCore
374			shiftedPot = MixingMap(ljt1,ljt2)%shiftedPot
375
376	gezelter	1608	r6 = r2 * r2 * r2
377	gezelter	2204
378	gezelter	1608	t6 = sigma6/ r6
379			t12 = t6 * t6
380	tim	2062
381	gezelter	2271	if (isSoftCore) then
382
383			pot_temp = 4.0E0_DP * epsilon * t6
384			if (shiftedPot) then
385			pot_temp = pot_temp + delta
386			endif
387
388			vpair = vpair + pot_temp
389
390			dudr = -sw * 24.0E0_DP * epsilon * t6 / rij
391
392	tim	2062	else
393	gezelter	2204	pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
394	gezelter	2271	if (shiftedPot) then
395	gezelter	2204	pot_temp = pot_temp + delta
396			endif
397	gezelter	2271
398	gezelter	2204	vpair = vpair + pot_temp
399	gezelter	2271
400	gezelter	2204	dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
401	gezelter	1608	endif
402	gezelter	2204
403	gezelter	1608	drdx = d(1) / rij
404			drdy = d(2) / rij
405			drdz = d(3) / rij
406	gezelter	2204
407	gezelter	1608	fx = dudr * drdx
408			fy = dudr * drdy
409			fz = dudr * drdz
410	gezelter	2204
411	gezelter	1608	#ifdef IS_MPI
412			if (do_pot) then
413			pot_Row(atom1) = pot_Row(atom1) + swpot_temp0.5
414			pot_Col(atom2) = pot_Col(atom2) + swpot_temp0.5
415			endif
416	gezelter	2204
417	gezelter	1608	f_Row(1,atom1) = f_Row(1,atom1) + fx
418			f_Row(2,atom1) = f_Row(2,atom1) + fy
419			f_Row(3,atom1) = f_Row(3,atom1) + fz
420	gezelter	2204
421	gezelter	1608	f_Col(1,atom2) = f_Col(1,atom2) - fx
422			f_Col(2,atom2) = f_Col(2,atom2) - fy
423			f_Col(3,atom2) = f_Col(3,atom2) - fz
424	gezelter	2204
425	gezelter	1608	#else
426			if (do_pot) pot = pot + sw*pot_temp
427
428			f(1,atom1) = f(1,atom1) + fx
429			f(2,atom1) = f(2,atom1) + fy
430			f(3,atom1) = f(3,atom1) + fz
431	gezelter	2204
432	gezelter	1608	f(1,atom2) = f(1,atom2) - fx
433			f(2,atom2) = f(2,atom2) - fy
434			f(3,atom2) = f(3,atom2) - fz
435			#endif
436	gezelter	2204
437	gezelter	1608	#ifdef IS_MPI
438			id1 = AtomRowToGlobal(atom1)
439			id2 = AtomColToGlobal(atom2)
440			#else
441			id1 = atom1
442			id2 = atom2
443			#endif
444
445			if (molMembershipList(id1) .ne. molMembershipList(id2)) then
446	gezelter	2204
447	gezelter	1608	fpair(1) = fpair(1) + fx
448			fpair(2) = fpair(2) + fy
449			fpair(3) = fpair(3) + fz
450
451			endif
452
453			return
454	gezelter	2204
455	gezelter	1608	end subroutine do_lj_pair
456	gezelter	2204
457	chuckv	2189	subroutine destroyLJTypes()
458	gezelter	2271
459			LJMap%nLJtypes = 0
460			LJMap%currentLJtype = 0
461
462			if (associated(LJMap%LJtypes)) then
463			deallocate(LJMap%LJtypes)
464			LJMap%LJtypes => null()
465			end if
466
467			if (associated(LJMap%atidToLJtype)) then
468			deallocate(LJMap%atidToLJtype)
469			LJMap%atidToLJtype => null()
470			end if
471
472	chuckv	2189	haveMixingMap = .false.
473			end subroutine destroyLJTypes
474
475	gezelter	1608	end module lj