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.28 2007-04-20 18:15:47 chrisfen Exp $, $Date: 2007-04-20 18:15:47 $, $Name: not supported by cvs2svn $, $Revision: 1.28 $


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 :: 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, 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, shiftedFrc
    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
    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 = epsilon * (myPot - myPotC)
    vpair = vpair + pot_temp
    dudr = sw * 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

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