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.31 2009-10-23 18:41:08 gezelter Exp $, $Date: 2009-10-23 18:41:08 $, $Name: not supported by cvs2svn $, $Revision: 1.31 $


module lj
  use definitions
  use atype_module
  use vector_class
  use simulation
  use status
  use fForceOptions
  use force_globals

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

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

  real(kind=DP), save :: defaultCutoff = 0.0_DP
  logical, save :: defaultShiftPot = .false.
  logical, save :: defaultShiftFrc = .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.
     logical       :: shiftedFrc
  end type MixParameters

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

  public :: newLJtype
  public :: setLJDefaultCutoff
  public :: getSigma
  public :: getEpsilon
  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, shiftedFrc)
    real(kind=dp), intent(in) :: thisRcut
    logical, intent(in) :: shiftedPot
    logical, intent(in) :: shiftedFrc
    defaultCutoff = thisRcut
    defaultShiftPot = shiftedPot
    defaultShiftFrc = shiftedFrc
    haveDefaultCutoff = .true.

    !we only want to build LJ Mixing map if LJ is being used.
    if(LJMap%nLJTypes /= 0) then
       call createMixingMap()
    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 = sqrt(s1 * s2)
          else
             MixingMap(i,j)%sigma = 0.5_dp * (s1 + s2)
          endif
          
          MixingMap(i,j)%epsilon = sqrt(e1 * e2)

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

          if (haveDefaultCutoff) then
             MixingMap(i,j)%shiftedPot = defaultShiftPot
             MixingMap(i,j)%shiftedFrc = defaultShiftFrc
          else
             MixingMap(i,j)%shiftedPot = defaultShiftPot
             MixingMap(i,j)%shiftedFrc = defaultShiftFrc
          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
             MixingMap(j,i)%shiftedFrc = MixingMap(i,j)%shiftedFrc
          endif

       enddo
    enddo
    
    haveMixingMap = .true.

  end subroutine createMixingMap
          
  subroutine do_lj_pair(atom1, atom2, atid1, atid2, d, rij, r2, rcut, sw, vdwMult, &
       vpair, fpair, pot, f1, do_pot)
    
    integer, intent(in) ::  atom1, atom2, atid1, atid2
    integer :: ljt1, ljt2
    real( kind = dp ), intent(in) :: rij, r2, rcut, vdwMult
    real( kind = dp ) :: pot, sw, vpair
    real( kind = dp ), intent(inout), dimension(3) :: f1
    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, shiftedFrc
    integer :: id1, id2, localError

    if (.not.haveMixingMap) then
       call createMixingMap()
    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
    shiftedFrc = MixingMap(ljt1,ljt2)%shiftedFrc

    ros = rij * sigmai

    if (isSoftCore) then

       call getSoftFunc(ros, myPot, myDeriv)

       if (shiftedPot) then
          rcos = rcut * sigmai
          call getSoftFunc(rcos, myPotC, myDerivC)
          myDerivC = 0.0_dp
       elseif (shiftedFrc) then
          rcos = rcut * sigmai
          call getSoftFunc(rcos, myPotC, myDerivC)
          myPotC = myPotC + myDerivC * (rij - rcut) * sigmai
       else
          myPotC = 0.0_dp
          myDerivC = 0.0_dp
       endif
              
    else

       call getLJfunc(ros, myPot, myDeriv)

       if (shiftedPot) then
          rcos = rcut * sigmai
          call getLJfunc(rcos, myPotC, myDerivC) 
          myDerivC = 0.0_dp
       elseif (shiftedFrc) then
          rcos = rcut * sigmai
          call getLJfunc(rcos, myPotC, myDerivC)
          myPotC = myPotC + myDerivC * (rij - rcut) * sigmai
       else
          myPotC = 0.0_dp
          myDerivC = 0.0_dp
       endif
       
    endif

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

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

    fx = dudr * drdx
    fy = dudr * drdy
    fz = dudr * drdz
    
    pot = pot + sw*pot_temp

    f1(1) = fx
    f1(2) = fy
    f1(3) = fz

    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

  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

    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)
    
    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
    
    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 
    
    return
  end subroutine getSoftFunc

end module lj
Revision:	3559
Committed:	Fri Oct 23 18:41:09 2009 UTC (15 years, 10 months ago) by gezelter
File size:	12485 byte(s)
Log Message:	removing MPI responsibilities from the lowest level force routines. This is in preparation for migrating these routines (LJ, electrostatic, eam, suttonchen, gay-berne, sticky) to C++
#	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	3559	!! @version $Id: LJ.F90,v 1.31 2009-10-23 18:41:08 gezelter Exp $, $Date: 2009-10-23 18:41:08 $, $Name: not supported by cvs2svn $, $Revision: 1.31 $
47	gezelter	1608
48	gezelter	1930
49	gezelter	1608	module lj
50	gezelter	2756	use definitions
51	gezelter	1608	use atype_module
52			use vector_class
53			use simulation
54	gezelter	1628	use status
55	chuckv	2497	use fForceOptions
56	gezelter	1608	use force_globals
57
58			implicit none
59			PRIVATE
60	chuckv	2355	#define __FORTRAN90
61			#include "UseTheForce/DarkSide/fInteractionMap.h"
62	gezelter	2204
63	gezelter	2271	logical, save :: useGeometricDistanceMixing = .false.
64			logical, save :: haveMixingMap = .false.
65
66			real(kind=DP), save :: defaultCutoff = 0.0_DP
67	chrisfen	3129	logical, save :: defaultShiftPot = .false.
68			logical, save :: defaultShiftFrc = .false.
69	gezelter	2271	logical, save :: haveDefaultCutoff = .false.
70
71			type, private :: LJtype
72			integer :: atid
73	gezelter	1628	real(kind=dp) :: sigma
74			real(kind=dp) :: epsilon
75	gezelter	2271	logical :: isSoftCore = .false.
76			end type LJtype
77	gezelter	2204
78	gezelter	2271	type, private :: LJList
79	chuckv	2319	integer :: Nljtypes = 0
80	gezelter	2271	integer :: currentLJtype = 0
81			type(LJtype), pointer :: LJtypes(:) => null()
82			integer, pointer :: atidToLJtype(:) => null()
83			end type LJList
84	gezelter	2204
85	gezelter	2271	type(LJList), save :: LJMap
86	gezelter	2204
87	gezelter	1628	type :: MixParameters
88			real(kind=DP) :: sigma
89			real(kind=DP) :: epsilon
90	gezelter	2717	real(kind=dp) :: sigmai
91	gezelter	2271	real(kind=dp) :: rCut
92			logical :: rCutWasSet = .false.
93			logical :: shiftedPot
94			logical :: isSoftCore = .false.
95	chrisfen	3129	logical :: shiftedFrc
96	gezelter	1628	end type MixParameters
97	gezelter	2204
98	gezelter	1628	type(MixParameters), dimension(:,:), allocatable :: MixingMap
99	gezelter	2204
100	gezelter	2271	public :: newLJtype
101			public :: setLJDefaultCutoff
102			public :: getSigma
103			public :: getEpsilon
104	gezelter	1628	public :: do_lj_pair
105	gezelter	2271	public :: destroyLJtypes
106	gezelter	2204
107	gezelter	1608	contains
108
109	gezelter	2271	subroutine newLJtype(c_ident, sigma, epsilon, isSoftCore, status)
110	gezelter	1930	integer,intent(in) :: c_ident
111	gezelter	1628	real(kind=dp),intent(in) :: sigma
112			real(kind=dp),intent(in) :: epsilon
113	gezelter	2271	integer, intent(in) :: isSoftCore
114	chuckv	1624	integer,intent(out) :: status
115	gezelter	2271	integer :: nLJTypes, ntypes, myATID
116	gezelter	1930	integer, pointer :: MatchList(:) => null()
117	gezelter	2271	integer :: current
118	gezelter	1628
119	chuckv	1624	status = 0
120	gezelter	2271	! check to see if this is the first time into this routine...
121			if (.not.associated(LJMap%LJtypes)) then
122	gezelter	2204
123	gezelter	2271	call getMatchingElementList(atypes, "is_LennardJones", .true., &
124			nLJTypes, MatchList)
125
126			LJMap%nLJtypes = nLJTypes
127
128			allocate(LJMap%LJtypes(nLJTypes))
129
130			ntypes = getSize(atypes)
131
132			allocate(LJMap%atidToLJtype(ntypes))
133			end if
134
135			LJMap%currentLJtype = LJMap%currentLJtype + 1
136			current = LJMap%currentLJtype
137
138	gezelter	1930	myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
139	gezelter	3432
140	gezelter	2271	LJMap%atidToLJtype(myATID) = current
141			LJMap%LJtypes(current)%atid = myATID
142			LJMap%LJtypes(current)%sigma = sigma
143			LJMap%LJtypes(current)%epsilon = epsilon
144			if (isSoftCore .eq. 1) then
145			LJMap%LJtypes(current)%isSoftCore = .true.
146			else
147			LJMap%LJtypes(current)%isSoftCore = .false.
148			endif
149			end subroutine newLJtype
150	chuckv	1624
151	chrisfen	3129	subroutine setLJDefaultCutoff(thisRcut, shiftedPot, shiftedFrc)
152	gezelter	2271	real(kind=dp), intent(in) :: thisRcut
153			logical, intent(in) :: shiftedPot
154	chrisfen	3129	logical, intent(in) :: shiftedFrc
155	gezelter	2271	defaultCutoff = thisRcut
156	chrisfen	3129	defaultShiftPot = shiftedPot
157			defaultShiftFrc = shiftedFrc
158	gezelter	2271	haveDefaultCutoff = .true.
159	chrisfen	3129
160	chrisfen	2727	!we only want to build LJ Mixing map if LJ is being used.
161	chuckv	2319	if(LJMap%nLJTypes /= 0) then
162			call createMixingMap()
163			end if
164	gezelter	2717
165	gezelter	2271	end subroutine setLJDefaultCutoff
166	gezelter	2204
167	gezelter	1633	function getSigma(atid) result (s)
168			integer, intent(in) :: atid
169	gezelter	2271	integer :: ljt1
170	gezelter	1633	real(kind=dp) :: s
171	gezelter	2204
172	gezelter	2271	if (LJMap%currentLJtype == 0) then
173			call handleError("LJ", "No members in LJMap")
174	gezelter	1633	return
175			end if
176	gezelter	2204
177	gezelter	2271	ljt1 = LJMap%atidToLJtype(atid)
178			s = LJMap%LJtypes(ljt1)%sigma
179
180	gezelter	1633	end function getSigma
181
182			function getEpsilon(atid) result (e)
183			integer, intent(in) :: atid
184	gezelter	2271	integer :: ljt1
185	gezelter	1633	real(kind=dp) :: e
186	gezelter	2204
187	gezelter	2271	if (LJMap%currentLJtype == 0) then
188			call handleError("LJ", "No members in LJMap")
189	gezelter	1633	return
190			end if
191	gezelter	2204
192	gezelter	2271	ljt1 = LJMap%atidToLJtype(atid)
193			e = LJMap%LJtypes(ljt1)%epsilon
194
195	gezelter	1633	end function getEpsilon
196
197	gezelter	2271	subroutine createMixingMap()
198			integer :: nLJtypes, i, j
199			real ( kind = dp ) :: s1, s2, e1, e2
200			real ( kind = dp ) :: rcut6, tp6, tp12
201	gezelter	2289	logical :: isSoftCore1, isSoftCore2, doShift
202	gezelter	1628
203	gezelter	2271	if (LJMap%currentLJtype == 0) then
204			call handleError("LJ", "No members in LJMap")
205	gezelter	2086	return
206	gezelter	2271	end if
207	gezelter	2204
208	gezelter	2271	nLJtypes = LJMap%nLJtypes
209	gezelter	2204
210	gezelter	2086	if (.not. allocated(MixingMap)) then
211	gezelter	2271	allocate(MixingMap(nLJtypes, nLJtypes))
212	gezelter	2086	endif
213
214	chuckv	2497	useGeometricDistanceMixing = usesGeometricDistanceMixing()
215	gezelter	2271	do i = 1, nLJtypes
216	gezelter	2204
217	gezelter	2271	s1 = LJMap%LJtypes(i)%sigma
218			e1 = LJMap%LJtypes(i)%epsilon
219			isSoftCore1 = LJMap%LJtypes(i)%isSoftCore
220	gezelter	2204
221	gezelter	2271	do j = i, nLJtypes
222
223			s2 = LJMap%LJtypes(j)%sigma
224			e2 = LJMap%LJtypes(j)%epsilon
225			isSoftCore2 = LJMap%LJtypes(j)%isSoftCore
226
227			MixingMap(i,j)%isSoftCore = isSoftCore1 .or. isSoftCore2
228	gezelter	2204
229	gezelter	2271	! only the distance parameter uses different mixing policies
230			if (useGeometricDistanceMixing) then
231	gezelter	2756	MixingMap(i,j)%sigma = sqrt(s1 * s2)
232	gezelter	2271	else
233			MixingMap(i,j)%sigma = 0.5_dp * (s1 + s2)
234			endif
235
236	gezelter	2756	MixingMap(i,j)%epsilon = sqrt(e1 * e2)
237	gezelter	2204
238	gezelter	2717	MixingMap(i,j)%sigmai = 1.0_DP / (MixingMap(i,j)%sigma)
239	gezelter	2204
240	gezelter	2461	if (haveDefaultCutoff) then
241	chrisfen	3129	MixingMap(i,j)%shiftedPot = defaultShiftPot
242			MixingMap(i,j)%shiftedFrc = defaultShiftFrc
243	gezelter	2271	else
244	chrisfen	3129	MixingMap(i,j)%shiftedPot = defaultShiftPot
245			MixingMap(i,j)%shiftedFrc = defaultShiftFrc
246	gezelter	2461	endif
247	gezelter	2204
248	gezelter	2289	if (i.ne.j) then
249			MixingMap(j,i)%sigma = MixingMap(i,j)%sigma
250			MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
251	gezelter	2717	MixingMap(j,i)%sigmai = MixingMap(i,j)%sigmai
252	gezelter	2289	MixingMap(j,i)%rCut = MixingMap(i,j)%rCut
253			MixingMap(j,i)%rCutWasSet = MixingMap(i,j)%rCutWasSet
254			MixingMap(j,i)%shiftedPot = MixingMap(i,j)%shiftedPot
255			MixingMap(j,i)%isSoftCore = MixingMap(i,j)%isSoftCore
256	chrisfen	3129	MixingMap(j,i)%shiftedFrc = MixingMap(i,j)%shiftedFrc
257	gezelter	2289	endif
258
259	gezelter	2271	enddo
260			enddo
261
262	chrisfen	1754	haveMixingMap = .true.
263	chrisfen	3129
264	gezelter	1628	end subroutine createMixingMap
265	chrisfen	2727
266	gezelter	3559	subroutine do_lj_pair(atom1, atom2, atid1, atid2, d, rij, r2, rcut, sw, vdwMult, &
267			vpair, fpair, pot, f1, do_pot)
268	gezelter	2461
269	gezelter	3559	integer, intent(in) :: atom1, atom2, atid1, atid2
270			integer :: ljt1, ljt2
271	gezelter	3441	real( kind = dp ), intent(in) :: rij, r2, rcut, vdwMult
272	gezelter	1608	real( kind = dp ) :: pot, sw, vpair
273	gezelter	3559	real( kind = dp ), intent(inout), dimension(3) :: f1
274	gezelter	1608	real( kind = dp ), intent(in), dimension(3) :: d
275			real( kind = dp ), intent(inout), dimension(3) :: fpair
276			logical, intent(in) :: do_pot
277
278			! local Variables
279			real( kind = dp ) :: drdx, drdy, drdz
280			real( kind = dp ) :: fx, fy, fz
281	gezelter	2717	real( kind = dp ) :: myPot, myPotC, myDeriv, myDerivC, ros, rcos
282	gezelter	1608	real( kind = dp ) :: pot_temp, dudr
283	gezelter	2717	real( kind = dp ) :: sigmai
284	gezelter	1608	real( kind = dp ) :: epsilon
285	chrisfen	3129	logical :: isSoftCore, shiftedPot, shiftedFrc
286	gezelter	1628	integer :: id1, id2, localError
287	gezelter	1608
288	gezelter	1628	if (.not.haveMixingMap) then
289	gezelter	2271	call createMixingMap()
290	gezelter	1628	endif
291
292	gezelter	2271	ljt1 = LJMap%atidToLJtype(atid1)
293			ljt2 = LJMap%atidToLJtype(atid2)
294
295	gezelter	2717	sigmai = MixingMap(ljt1,ljt2)%sigmai
296	gezelter	2271	epsilon = MixingMap(ljt1,ljt2)%epsilon
297			isSoftCore = MixingMap(ljt1,ljt2)%isSoftCore
298			shiftedPot = MixingMap(ljt1,ljt2)%shiftedPot
299	chrisfen	3129	shiftedFrc = MixingMap(ljt1,ljt2)%shiftedFrc
300	gezelter	2271
301	gezelter	2717	ros = rij * sigmai
302	gezelter	2204
303	gezelter	2717	if (isSoftCore) then
304	tim	2062
305	gezelter	2717	call getSoftFunc(ros, myPot, myDeriv)
306
307	gezelter	2271	if (shiftedPot) then
308	gezelter	2717	rcos = rcut * sigmai
309			call getSoftFunc(rcos, myPotC, myDerivC)
310	chrisfen	3129	myDerivC = 0.0_dp
311			elseif (shiftedFrc) then
312			rcos = rcut * sigmai
313			call getSoftFunc(rcos, myPotC, myDerivC)
314			myPotC = myPotC + myDerivC * (rij - rcut) * sigmai
315			else
316			myPotC = 0.0_dp
317			myDerivC = 0.0_dp
318	gezelter	2271	endif
319	gezelter	2717
320	tim	2062	else
321	gezelter	2717
322			call getLJfunc(ros, myPot, myDeriv)
323
324	gezelter	2271	if (shiftedPot) then
325	gezelter	2717	rcos = rcut * sigmai
326	gezelter	3126	call getLJfunc(rcos, myPotC, myDerivC)
327	chrisfen	3129	myDerivC = 0.0_dp
328			elseif (shiftedFrc) then
329			rcos = rcut * sigmai
330			call getLJfunc(rcos, myPotC, myDerivC)
331			myPotC = myPotC + myDerivC * (rij - rcut) * sigmai
332			else
333			myPotC = 0.0_dp
334			myDerivC = 0.0_dp
335	gezelter	2204	endif
336	gezelter	2271
337	gezelter	1608	endif
338	gezelter	2204
339	gezelter	3441	pot_temp = vdwMult * epsilon * (myPot - myPotC)
340	chrisfen	3129	vpair = vpair + pot_temp
341	gezelter	3441	dudr = sw * vdwMult * epsilon * (myDeriv - myDerivC) * sigmai
342	gezelter	2717
343	gezelter	1608	drdx = d(1) / rij
344			drdy = d(2) / rij
345			drdz = d(3) / rij
346	gezelter	2204
347	gezelter	1608	fx = dudr * drdx
348			fy = dudr * drdy
349			fz = dudr * drdz
350	gezelter	3559
351			pot = pot + sw*pot_temp
352	gezelter	2204
353	gezelter	3559	f1(1) = fx
354			f1(2) = fy
355			f1(3) = fz
356	gezelter	2204
357	gezelter	1608	return
358	gezelter	2204
359	gezelter	1608	end subroutine do_lj_pair
360	gezelter	2204
361	chuckv	2189	subroutine destroyLJTypes()
362	gezelter	2271
363			LJMap%nLJtypes = 0
364			LJMap%currentLJtype = 0
365
366			if (associated(LJMap%LJtypes)) then
367			deallocate(LJMap%LJtypes)
368			LJMap%LJtypes => null()
369			end if
370
371			if (associated(LJMap%atidToLJtype)) then
372			deallocate(LJMap%atidToLJtype)
373			LJMap%atidToLJtype => null()
374			end if
375
376	chuckv	2189	haveMixingMap = .false.
377	gezelter	2717
378	chuckv	2189	end subroutine destroyLJTypes
379
380	gezelter	2717	subroutine getLJfunc(r, myPot, myDeriv)
381
382			real(kind=dp), intent(in) :: r
383			real(kind=dp), intent(inout) :: myPot, myDeriv
384			real(kind=dp) :: ri, ri2, ri6, ri7, ri12, ri13
385			real(kind=dp) :: a, b, c, d, dx
386			integer :: j
387
388	chrisfen	2727	ri = 1.0_DP / r
389			ri2 = ri*ri
390			ri6 = ri2ri2ri2
391			ri7 = ri6*ri
392			ri12 = ri6*ri6
393			ri13 = ri12*ri
394
395			myPot = 4.0_DP * (ri12 - ri6)
396			myDeriv = 24.0_DP * (ri7 - 2.0_DP * ri13)
397
398	gezelter	2717	return
399			end subroutine getLJfunc
400
401			subroutine getSoftFunc(r, myPot, myDeriv)
402
403			real(kind=dp), intent(in) :: r
404			real(kind=dp), intent(inout) :: myPot, myDeriv
405			real(kind=dp) :: ri, ri2, ri6, ri7
406			real(kind=dp) :: a, b, c, d, dx
407			integer :: j
408
409	chrisfen	2727	ri = 1.0_DP / r
410			ri2 = ri*ri
411			ri6 = ri2ri2ri2
412			ri7 = ri6*ri
413			myPot = 4.0_DP * (ri6)
414			myDeriv = - 24.0_DP * ri7
415
416	gezelter	2717	return
417			end subroutine getSoftFunc
418
419	gezelter	1608	end module lj