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.25 2006-05-17 15:37:14 gezelter Exp $, $Date: 2006-05-17 15:37:14 $, $Name: not supported by cvs2svn $, $Revision: 1.25 $


module lj
  use definitions
  use atype_module
  use vector_class
  use simulation
  use status
  use fForceOptions
#ifdef IS_MPI
  use mpiSimulation
#endif
  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 :: 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

  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)
    real(kind=dp), intent(in) :: thisRcut
    logical, intent(in) :: shiftedPot
    defaultCutoff = thisRcut
    defaultShift = shiftedPot
    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 = 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 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

    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.

  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:	2756
Committed:	Wed May 17 15:37:15 2006 UTC (18 years, 4 months ago) by gezelter
File size:	12639 byte(s)
Log Message:	Getting fortran side prepped for single precision...
#	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	2756	!! @version $Id: LJ.F90,v 1.25 2006-05-17 15:37:14 gezelter Exp $, $Date: 2006-05-17 15:37:14 $, $Name: not supported by cvs2svn $, $Revision: 1.25 $
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	#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	gezelter	2271	logical, save :: useGeometricDistanceMixing = .false.
67			logical, save :: haveMixingMap = .false.
68
69			real(kind=DP), save :: defaultCutoff = 0.0_DP
70			logical, save :: defaultShift = .false.
71			logical, save :: haveDefaultCutoff = .false.
72
73			type, private :: LJtype
74			integer :: atid
75	gezelter	1628	real(kind=dp) :: sigma
76			real(kind=dp) :: epsilon
77	gezelter	2271	logical :: isSoftCore = .false.
78			end type LJtype
79	gezelter	2204
80	gezelter	2271	type, private :: LJList
81	chuckv	2319	integer :: Nljtypes = 0
82	gezelter	2271	integer :: currentLJtype = 0
83			type(LJtype), pointer :: LJtypes(:) => null()
84			integer, pointer :: atidToLJtype(:) => null()
85			end type LJList
86	gezelter	2204
87	gezelter	2271	type(LJList), save :: LJMap
88	gezelter	2204
89	gezelter	1628	type :: MixParameters
90			real(kind=DP) :: sigma
91			real(kind=DP) :: epsilon
92	gezelter	2717	real(kind=dp) :: sigmai
93	gezelter	2271	real(kind=dp) :: rCut
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 :: getSigma
104			public :: getEpsilon
105	gezelter	1628	public :: do_lj_pair
106	gezelter	2271	public :: destroyLJtypes
107	gezelter	2204
108	gezelter	1608	contains
109
110	gezelter	2271	subroutine newLJtype(c_ident, sigma, epsilon, isSoftCore, status)
111	gezelter	1930	integer,intent(in) :: c_ident
112	gezelter	1628	real(kind=dp),intent(in) :: sigma
113			real(kind=dp),intent(in) :: epsilon
114	gezelter	2271	integer, intent(in) :: isSoftCore
115	chuckv	1624	integer,intent(out) :: status
116	gezelter	2271	integer :: nLJTypes, ntypes, myATID
117	gezelter	1930	integer, pointer :: MatchList(:) => null()
118	gezelter	2271	integer :: current
119	gezelter	1628
120	chuckv	1624	status = 0
121	gezelter	2271	! check to see if this is the first time into this routine...
122			if (.not.associated(LJMap%LJtypes)) then
123	gezelter	2204
124	gezelter	2271	call getMatchingElementList(atypes, "is_LennardJones", .true., &
125			nLJTypes, MatchList)
126
127			LJMap%nLJtypes = nLJTypes
128
129			allocate(LJMap%LJtypes(nLJTypes))
130
131			ntypes = getSize(atypes)
132
133			allocate(LJMap%atidToLJtype(ntypes))
134			end if
135
136			LJMap%currentLJtype = LJMap%currentLJtype + 1
137			current = LJMap%currentLJtype
138
139	gezelter	1930	myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
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	gezelter	2271	subroutine setLJDefaultCutoff(thisRcut, shiftedPot)
152			real(kind=dp), intent(in) :: thisRcut
153			logical, intent(in) :: shiftedPot
154			defaultCutoff = thisRcut
155			defaultShift = shiftedPot
156			haveDefaultCutoff = .true.
157	chrisfen	2727	!we only want to build LJ Mixing map if LJ is being used.
158	chuckv	2319	if(LJMap%nLJTypes /= 0) then
159			call createMixingMap()
160			end if
161	gezelter	2717
162	gezelter	2271	end subroutine setLJDefaultCutoff
163	gezelter	2204
164	gezelter	1633	function getSigma(atid) result (s)
165			integer, intent(in) :: atid
166	gezelter	2271	integer :: ljt1
167	gezelter	1633	real(kind=dp) :: s
168	gezelter	2204
169	gezelter	2271	if (LJMap%currentLJtype == 0) then
170			call handleError("LJ", "No members in LJMap")
171	gezelter	1633	return
172			end if
173	gezelter	2204
174	gezelter	2271	ljt1 = LJMap%atidToLJtype(atid)
175			s = LJMap%LJtypes(ljt1)%sigma
176
177	gezelter	1633	end function getSigma
178
179			function getEpsilon(atid) result (e)
180			integer, intent(in) :: atid
181	gezelter	2271	integer :: ljt1
182	gezelter	1633	real(kind=dp) :: e
183	gezelter	2204
184	gezelter	2271	if (LJMap%currentLJtype == 0) then
185			call handleError("LJ", "No members in LJMap")
186	gezelter	1633	return
187			end if
188	gezelter	2204
189	gezelter	2271	ljt1 = LJMap%atidToLJtype(atid)
190			e = LJMap%LJtypes(ljt1)%epsilon
191
192	gezelter	1633	end function getEpsilon
193
194	gezelter	2271	subroutine createMixingMap()
195			integer :: nLJtypes, i, j
196			real ( kind = dp ) :: s1, s2, e1, e2
197			real ( kind = dp ) :: rcut6, tp6, tp12
198	gezelter	2289	logical :: isSoftCore1, isSoftCore2, doShift
199	gezelter	1628
200	gezelter	2271	if (LJMap%currentLJtype == 0) then
201			call handleError("LJ", "No members in LJMap")
202	gezelter	2086	return
203	gezelter	2271	end if
204	gezelter	2204
205	gezelter	2271	nLJtypes = LJMap%nLJtypes
206	gezelter	2204
207	gezelter	2086	if (.not. allocated(MixingMap)) then
208	gezelter	2271	allocate(MixingMap(nLJtypes, nLJtypes))
209	gezelter	2086	endif
210
211	chuckv	2497	useGeometricDistanceMixing = usesGeometricDistanceMixing()
212	gezelter	2271	do i = 1, nLJtypes
213	gezelter	2204
214	gezelter	2271	s1 = LJMap%LJtypes(i)%sigma
215			e1 = LJMap%LJtypes(i)%epsilon
216			isSoftCore1 = LJMap%LJtypes(i)%isSoftCore
217	gezelter	2204
218	gezelter	2271	do j = i, nLJtypes
219
220			s2 = LJMap%LJtypes(j)%sigma
221			e2 = LJMap%LJtypes(j)%epsilon
222			isSoftCore2 = LJMap%LJtypes(j)%isSoftCore
223
224			MixingMap(i,j)%isSoftCore = isSoftCore1 .or. isSoftCore2
225	gezelter	2204
226	gezelter	2271	! only the distance parameter uses different mixing policies
227			if (useGeometricDistanceMixing) then
228	gezelter	2756	MixingMap(i,j)%sigma = sqrt(s1 * s2)
229	gezelter	2271	else
230			MixingMap(i,j)%sigma = 0.5_dp * (s1 + s2)
231			endif
232
233	gezelter	2756	MixingMap(i,j)%epsilon = sqrt(e1 * e2)
234	gezelter	2204
235	gezelter	2717	MixingMap(i,j)%sigmai = 1.0_DP / (MixingMap(i,j)%sigma)
236	gezelter	2204
237	gezelter	2461	if (haveDefaultCutoff) then
238			MixingMap(i,j)%shiftedPot = defaultShift
239	gezelter	2271	else
240	gezelter	2461	MixingMap(i,j)%shiftedPot = defaultShift
241			endif
242	gezelter	2204
243	gezelter	2289	if (i.ne.j) then
244			MixingMap(j,i)%sigma = MixingMap(i,j)%sigma
245			MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
246	gezelter	2717	MixingMap(j,i)%sigmai = MixingMap(i,j)%sigmai
247	gezelter	2289	MixingMap(j,i)%rCut = MixingMap(i,j)%rCut
248			MixingMap(j,i)%rCutWasSet = MixingMap(i,j)%rCutWasSet
249			MixingMap(j,i)%shiftedPot = MixingMap(i,j)%shiftedPot
250			MixingMap(j,i)%isSoftCore = MixingMap(i,j)%isSoftCore
251			endif
252
253	gezelter	2271	enddo
254			enddo
255
256	chrisfen	1754	haveMixingMap = .true.
257	gezelter	2271
258	gezelter	1628	end subroutine createMixingMap
259	chrisfen	2727
260	gezelter	2461	subroutine do_lj_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
261	gezelter	1608	pot, f, do_pot)
262	gezelter	2461
263	gezelter	1608	integer, intent(in) :: atom1, atom2
264	gezelter	2271	integer :: atid1, atid2, ljt1, ljt2
265	gezelter	2461	real( kind = dp ), intent(in) :: rij, r2, rcut
266	gezelter	1608	real( kind = dp ) :: pot, sw, vpair
267			real( kind = dp ), dimension(3,nLocal) :: f
268			real( kind = dp ), intent(in), dimension(3) :: d
269			real( kind = dp ), intent(inout), dimension(3) :: fpair
270			logical, intent(in) :: do_pot
271
272			! local Variables
273			real( kind = dp ) :: drdx, drdy, drdz
274			real( kind = dp ) :: fx, fy, fz
275	gezelter	2717	real( kind = dp ) :: myPot, myPotC, myDeriv, myDerivC, ros, rcos
276	gezelter	1608	real( kind = dp ) :: pot_temp, dudr
277	gezelter	2717	real( kind = dp ) :: sigmai
278	gezelter	1608	real( kind = dp ) :: epsilon
279	gezelter	2271	logical :: isSoftCore, shiftedPot
280	gezelter	1628	integer :: id1, id2, localError
281	gezelter	1608
282	gezelter	1628	if (.not.haveMixingMap) then
283	gezelter	2271	call createMixingMap()
284	gezelter	1628	endif
285
286	gezelter	1608	! Look up the correct parameters in the mixing matrix
287			#ifdef IS_MPI
288	gezelter	2271	atid1 = atid_Row(atom1)
289			atid2 = atid_Col(atom2)
290	gezelter	1608	#else
291	gezelter	2271	atid1 = atid(atom1)
292			atid2 = atid(atom2)
293	gezelter	1608	#endif
294
295	gezelter	2271	ljt1 = LJMap%atidToLJtype(atid1)
296			ljt2 = LJMap%atidToLJtype(atid2)
297
298	gezelter	2717	sigmai = MixingMap(ljt1,ljt2)%sigmai
299	gezelter	2271	epsilon = MixingMap(ljt1,ljt2)%epsilon
300			isSoftCore = MixingMap(ljt1,ljt2)%isSoftCore
301			shiftedPot = MixingMap(ljt1,ljt2)%shiftedPot
302
303	gezelter	2717	ros = rij * sigmai
304			myPotC = 0.0_DP
305	gezelter	2204
306	gezelter	2717	if (isSoftCore) then
307	tim	2062
308	gezelter	2717	call getSoftFunc(ros, myPot, myDeriv)
309
310	gezelter	2271	if (shiftedPot) then
311	gezelter	2717	rcos = rcut * sigmai
312			call getSoftFunc(rcos, myPotC, myDerivC)
313	gezelter	2271	endif
314	gezelter	2717
315	tim	2062	else
316	gezelter	2717
317			call getLJfunc(ros, myPot, myDeriv)
318
319	gezelter	2271	if (shiftedPot) then
320	gezelter	2717	rcos = rcut * sigmai
321			call getLJfunc(rcos, myPotC, myDerivC)
322	gezelter	2204	endif
323	gezelter	2271
324	gezelter	1608	endif
325	gezelter	2204
326	gezelter	2717	pot_temp = epsilon * (myPot - myPotC)
327			vpair = vpair + pot_temp
328			dudr = sw * epsilon * myDeriv * sigmai
329
330	gezelter	1608	drdx = d(1) / rij
331			drdy = d(2) / rij
332			drdz = d(3) / rij
333	gezelter	2204
334	gezelter	1608	fx = dudr * drdx
335			fy = dudr * drdy
336			fz = dudr * drdz
337	gezelter	2204
338	gezelter	1608	#ifdef IS_MPI
339			if (do_pot) then
340	gezelter	2361	pot_Row(VDW_POT,atom1) = pot_Row(VDW_POT,atom1) + swpot_temp0.5
341			pot_Col(VDW_POT,atom2) = pot_Col(VDW_POT,atom2) + swpot_temp0.5
342	gezelter	1608	endif
343	gezelter	2204
344	gezelter	1608	f_Row(1,atom1) = f_Row(1,atom1) + fx
345			f_Row(2,atom1) = f_Row(2,atom1) + fy
346			f_Row(3,atom1) = f_Row(3,atom1) + fz
347	gezelter	2204
348	gezelter	1608	f_Col(1,atom2) = f_Col(1,atom2) - fx
349			f_Col(2,atom2) = f_Col(2,atom2) - fy
350			f_Col(3,atom2) = f_Col(3,atom2) - fz
351	gezelter	2204
352	gezelter	1608	#else
353			if (do_pot) pot = pot + sw*pot_temp
354
355			f(1,atom1) = f(1,atom1) + fx
356			f(2,atom1) = f(2,atom1) + fy
357			f(3,atom1) = f(3,atom1) + fz
358	gezelter	2204
359	gezelter	1608	f(1,atom2) = f(1,atom2) - fx
360			f(2,atom2) = f(2,atom2) - fy
361			f(3,atom2) = f(3,atom2) - fz
362			#endif
363	gezelter	2204
364	gezelter	1608	#ifdef IS_MPI
365			id1 = AtomRowToGlobal(atom1)
366			id2 = AtomColToGlobal(atom2)
367			#else
368			id1 = atom1
369			id2 = atom2
370			#endif
371
372			if (molMembershipList(id1) .ne. molMembershipList(id2)) then
373	gezelter	2204
374	gezelter	1608	fpair(1) = fpair(1) + fx
375			fpair(2) = fpair(2) + fy
376			fpair(3) = fpair(3) + fz
377
378			endif
379
380			return
381	gezelter	2204
382	gezelter	1608	end subroutine do_lj_pair
383	gezelter	2204
384	chuckv	2189	subroutine destroyLJTypes()
385	gezelter	2271
386			LJMap%nLJtypes = 0
387			LJMap%currentLJtype = 0
388
389			if (associated(LJMap%LJtypes)) then
390			deallocate(LJMap%LJtypes)
391			LJMap%LJtypes => null()
392			end if
393
394			if (associated(LJMap%atidToLJtype)) then
395			deallocate(LJMap%atidToLJtype)
396			LJMap%atidToLJtype => null()
397			end if
398
399	chuckv	2189	haveMixingMap = .false.
400	gezelter	2717
401	chuckv	2189	end subroutine destroyLJTypes
402
403	gezelter	2717	subroutine getLJfunc(r, myPot, myDeriv)
404
405			real(kind=dp), intent(in) :: r
406			real(kind=dp), intent(inout) :: myPot, myDeriv
407			real(kind=dp) :: ri, ri2, ri6, ri7, ri12, ri13
408			real(kind=dp) :: a, b, c, d, dx
409			integer :: j
410
411	chrisfen	2727	ri = 1.0_DP / r
412			ri2 = ri*ri
413			ri6 = ri2ri2ri2
414			ri7 = ri6*ri
415			ri12 = ri6*ri6
416			ri13 = ri12*ri
417
418			myPot = 4.0_DP * (ri12 - ri6)
419			myDeriv = 24.0_DP * (ri7 - 2.0_DP * ri13)
420
421	gezelter	2717	return
422			end subroutine getLJfunc
423
424			subroutine getSoftFunc(r, myPot, myDeriv)
425
426			real(kind=dp), intent(in) :: r
427			real(kind=dp), intent(inout) :: myPot, myDeriv
428			real(kind=dp) :: ri, ri2, ri6, ri7
429			real(kind=dp) :: a, b, c, d, dx
430			integer :: j
431
432	chrisfen	2727	ri = 1.0_DP / r
433			ri2 = ri*ri
434			ri6 = ri2ri2ri2
435			ri7 = ri6*ri
436			myPot = 4.0_DP * (ri6)
437			myDeriv = - 24.0_DP * ri7
438
439	gezelter	2717	return
440			end subroutine getSoftFunc
441
442	gezelter	1608	end module lj