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.13 2005-08-09 22:33:48 gezelter Exp $, $Date: 2005-08-09 22:33:48 $, $Name: not supported by cvs2svn $, $Revision: 1.13 $


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.
  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

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