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.27 2007-04-09 18:24:00 gezelter Exp $, $Date: 2007-04-09 18:24:00 $, $Name: not supported by cvs2svn $, $Revision: 1.27 $


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
    myDerivC = 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

    !! these are the shifted POTENTIAL variants.
    ! pot_temp = epsilon * (myPot - myPotC)
    ! dudr = sw * epsilon * myDeriv * sigmai

    !! these are the shifted FORCE variants.

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

    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:	1127
Committed:	Mon Apr 9 18:24:00 2007 UTC (18 years, 3 months ago) by gezelter
File size:	12900 byte(s)
Log Message:	more changes for atomic virials and Lennard-Jones force switching
#	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	1127	!! @version $Id: LJ.F90,v 1.27 2007-04-09 18:24:00 gezelter Exp $, $Date: 2007-04-09 18:24:00 $, $Name: not supported by cvs2svn $, $Revision: 1.27 $
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			logical, save :: defaultShift = .false.
71			logical, save :: haveDefaultCutoff = .false.
72
73			type, private :: LJtype
74			integer :: atid
75	gezelter	135	real(kind=dp) :: sigma
76			real(kind=dp) :: epsilon
77	gezelter	572	logical :: isSoftCore = .false.
78			end type LJtype
79	gezelter	507
80	gezelter	572	type, private :: LJList
81	chuckv	620	integer :: Nljtypes = 0
82	gezelter	572	integer :: currentLJtype = 0
83			type(LJtype), pointer :: LJtypes(:) => null()
84			integer, pointer :: atidToLJtype(:) => null()
85			end type LJList
86	gezelter	507
87	gezelter	572	type(LJList), save :: LJMap
88	gezelter	507
89	gezelter	135	type :: MixParameters
90			real(kind=DP) :: sigma
91			real(kind=DP) :: epsilon
92	gezelter	938	real(kind=dp) :: sigmai
93	gezelter	572	real(kind=dp) :: rCut
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 :: getSigma
104			public :: getEpsilon
105	gezelter	135	public :: do_lj_pair
106	gezelter	572	public :: destroyLJtypes
107	gezelter	507
108	gezelter	115	contains
109
110	gezelter	572	subroutine newLJtype(c_ident, sigma, epsilon, isSoftCore, status)
111	gezelter	246	integer,intent(in) :: c_ident
112	gezelter	135	real(kind=dp),intent(in) :: sigma
113			real(kind=dp),intent(in) :: epsilon
114	gezelter	572	integer, intent(in) :: isSoftCore
115	chuckv	131	integer,intent(out) :: status
116	gezelter	572	integer :: nLJTypes, ntypes, myATID
117	gezelter	246	integer, pointer :: MatchList(:) => null()
118	gezelter	572	integer :: current
119	gezelter	135
120	chuckv	131	status = 0
121	gezelter	572	! check to see if this is the first time into this routine...
122			if (.not.associated(LJMap%LJtypes)) then
123	gezelter	507
124	gezelter	572	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	246	myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
140	gezelter	572	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	131
151	gezelter	572	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	942	!we only want to build LJ Mixing map if LJ is being used.
158	chuckv	620	if(LJMap%nLJTypes /= 0) then
159			call createMixingMap()
160			end if
161	gezelter	938
162	gezelter	572	end subroutine setLJDefaultCutoff
163	gezelter	507
164	gezelter	140	function getSigma(atid) result (s)
165			integer, intent(in) :: atid
166	gezelter	572	integer :: ljt1
167	gezelter	140	real(kind=dp) :: s
168	gezelter	507
169	gezelter	572	if (LJMap%currentLJtype == 0) then
170			call handleError("LJ", "No members in LJMap")
171	gezelter	140	return
172			end if
173	gezelter	507
174	gezelter	572	ljt1 = LJMap%atidToLJtype(atid)
175			s = LJMap%LJtypes(ljt1)%sigma
176
177	gezelter	140	end function getSigma
178
179			function getEpsilon(atid) result (e)
180			integer, intent(in) :: atid
181	gezelter	572	integer :: ljt1
182	gezelter	140	real(kind=dp) :: e
183	gezelter	507
184	gezelter	572	if (LJMap%currentLJtype == 0) then
185			call handleError("LJ", "No members in LJMap")
186	gezelter	140	return
187			end if
188	gezelter	507
189	gezelter	572	ljt1 = LJMap%atidToLJtype(atid)
190			e = LJMap%LJtypes(ljt1)%epsilon
191
192	gezelter	140	end function getEpsilon
193
194	gezelter	572	subroutine createMixingMap()
195			integer :: nLJtypes, i, j
196			real ( kind = dp ) :: s1, s2, e1, e2
197			real ( kind = dp ) :: rcut6, tp6, tp12
198	gezelter	590	logical :: isSoftCore1, isSoftCore2, doShift
199	gezelter	135
200	gezelter	572	if (LJMap%currentLJtype == 0) then
201			call handleError("LJ", "No members in LJMap")
202	gezelter	402	return
203	gezelter	572	end if
204	gezelter	507
205	gezelter	572	nLJtypes = LJMap%nLJtypes
206	gezelter	507
207	gezelter	402	if (.not. allocated(MixingMap)) then
208	gezelter	572	allocate(MixingMap(nLJtypes, nLJtypes))
209	gezelter	402	endif
210
211	chuckv	798	useGeometricDistanceMixing = usesGeometricDistanceMixing()
212	gezelter	572	do i = 1, nLJtypes
213	gezelter	507
214	gezelter	572	s1 = LJMap%LJtypes(i)%sigma
215			e1 = LJMap%LJtypes(i)%epsilon
216			isSoftCore1 = LJMap%LJtypes(i)%isSoftCore
217	gezelter	507
218	gezelter	572	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	507
226	gezelter	572	! only the distance parameter uses different mixing policies
227			if (useGeometricDistanceMixing) then
228	gezelter	960	MixingMap(i,j)%sigma = sqrt(s1 * s2)
229	gezelter	572	else
230			MixingMap(i,j)%sigma = 0.5_dp * (s1 + s2)
231			endif
232
233	gezelter	960	MixingMap(i,j)%epsilon = sqrt(e1 * e2)
234	gezelter	507
235	gezelter	938	MixingMap(i,j)%sigmai = 1.0_DP / (MixingMap(i,j)%sigma)
236	gezelter	507
237	gezelter	762	if (haveDefaultCutoff) then
238			MixingMap(i,j)%shiftedPot = defaultShift
239	gezelter	572	else
240	gezelter	762	MixingMap(i,j)%shiftedPot = defaultShift
241			endif
242	gezelter	507
243	gezelter	590	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	938	MixingMap(j,i)%sigmai = MixingMap(i,j)%sigmai
247	gezelter	590	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	572	enddo
254			enddo
255
256	chrisfen	217	haveMixingMap = .true.
257	gezelter	572
258	gezelter	135	end subroutine createMixingMap
259	chrisfen	942
260	gezelter	762	subroutine do_lj_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
261	gezelter	115	pot, f, do_pot)
262	gezelter	762
263	gezelter	115	integer, intent(in) :: atom1, atom2
264	gezelter	572	integer :: atid1, atid2, ljt1, ljt2
265	gezelter	762	real( kind = dp ), intent(in) :: rij, r2, rcut
266	gezelter	115	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	938	real( kind = dp ) :: myPot, myPotC, myDeriv, myDerivC, ros, rcos
276	gezelter	115	real( kind = dp ) :: pot_temp, dudr
277	gezelter	938	real( kind = dp ) :: sigmai
278	gezelter	115	real( kind = dp ) :: epsilon
279	gezelter	572	logical :: isSoftCore, shiftedPot
280	gezelter	135	integer :: id1, id2, localError
281	gezelter	115
282	gezelter	135	if (.not.haveMixingMap) then
283	gezelter	572	call createMixingMap()
284	gezelter	135	endif
285
286	gezelter	115	! Look up the correct parameters in the mixing matrix
287			#ifdef IS_MPI
288	gezelter	572	atid1 = atid_Row(atom1)
289			atid2 = atid_Col(atom2)
290	gezelter	115	#else
291	gezelter	572	atid1 = atid(atom1)
292			atid2 = atid(atom2)
293	gezelter	115	#endif
294
295	gezelter	572	ljt1 = LJMap%atidToLJtype(atid1)
296			ljt2 = LJMap%atidToLJtype(atid2)
297
298	gezelter	938	sigmai = MixingMap(ljt1,ljt2)%sigmai
299	gezelter	572	epsilon = MixingMap(ljt1,ljt2)%epsilon
300			isSoftCore = MixingMap(ljt1,ljt2)%isSoftCore
301			shiftedPot = MixingMap(ljt1,ljt2)%shiftedPot
302
303	gezelter	938	ros = rij * sigmai
304			myPotC = 0.0_DP
305	gezelter	1127	myDerivC = 0.0_DP
306	gezelter	507
307	gezelter	938	if (isSoftCore) then
308	tim	378
309	gezelter	938	call getSoftFunc(ros, myPot, myDeriv)
310
311	gezelter	572	if (shiftedPot) then
312	gezelter	938	rcos = rcut * sigmai
313			call getSoftFunc(rcos, myPotC, myDerivC)
314	gezelter	572	endif
315	gezelter	938
316	tim	378	else
317	gezelter	938
318			call getLJfunc(ros, myPot, myDeriv)
319
320	gezelter	572	if (shiftedPot) then
321	gezelter	938	rcos = rcut * sigmai
322	gezelter	1126	call getLJfunc(rcos, myPotC, myDerivC)
323	gezelter	507	endif
324	gezelter	572
325	gezelter	115	endif
326	gezelter	507
327	gezelter	1127	!! these are the shifted POTENTIAL variants.
328			! pot_temp = epsilon * (myPot - myPotC)
329			! dudr = sw * epsilon * myDeriv * sigmai
330
331			!! these are the shifted FORCE variants.
332
333			pot_temp = epsilon * (myPot - myPotC - myDerivC * (rij - rcut) * sigmai)
334			dudr = sw * epsilon * (myDeriv - myDerivC) * sigmai
335
336	gezelter	938	vpair = vpair + pot_temp
337
338	gezelter	115	drdx = d(1) / rij
339			drdy = d(2) / rij
340			drdz = d(3) / rij
341	gezelter	507
342	gezelter	115	fx = dudr * drdx
343			fy = dudr * drdy
344			fz = dudr * drdz
345	gezelter	507
346	gezelter	115	#ifdef IS_MPI
347			if (do_pot) then
348	gezelter	662	pot_Row(VDW_POT,atom1) = pot_Row(VDW_POT,atom1) + swpot_temp0.5
349			pot_Col(VDW_POT,atom2) = pot_Col(VDW_POT,atom2) + swpot_temp0.5
350	gezelter	115	endif
351	gezelter	507
352	gezelter	115	f_Row(1,atom1) = f_Row(1,atom1) + fx
353			f_Row(2,atom1) = f_Row(2,atom1) + fy
354			f_Row(3,atom1) = f_Row(3,atom1) + fz
355	gezelter	507
356	gezelter	115	f_Col(1,atom2) = f_Col(1,atom2) - fx
357			f_Col(2,atom2) = f_Col(2,atom2) - fy
358			f_Col(3,atom2) = f_Col(3,atom2) - fz
359	gezelter	507
360	gezelter	115	#else
361			if (do_pot) pot = pot + sw*pot_temp
362
363			f(1,atom1) = f(1,atom1) + fx
364			f(2,atom1) = f(2,atom1) + fy
365			f(3,atom1) = f(3,atom1) + fz
366	gezelter	507
367	gezelter	115	f(1,atom2) = f(1,atom2) - fx
368			f(2,atom2) = f(2,atom2) - fy
369			f(3,atom2) = f(3,atom2) - fz
370			#endif
371	gezelter	507
372	gezelter	115	#ifdef IS_MPI
373			id1 = AtomRowToGlobal(atom1)
374			id2 = AtomColToGlobal(atom2)
375			#else
376			id1 = atom1
377			id2 = atom2
378			#endif
379
380			if (molMembershipList(id1) .ne. molMembershipList(id2)) then
381	gezelter	507
382	gezelter	115	fpair(1) = fpair(1) + fx
383			fpair(2) = fpair(2) + fy
384			fpair(3) = fpair(3) + fz
385
386			endif
387
388			return
389	gezelter	507
390	gezelter	115	end subroutine do_lj_pair
391	gezelter	507
392	chuckv	492	subroutine destroyLJTypes()
393	gezelter	572
394			LJMap%nLJtypes = 0
395			LJMap%currentLJtype = 0
396
397			if (associated(LJMap%LJtypes)) then
398			deallocate(LJMap%LJtypes)
399			LJMap%LJtypes => null()
400			end if
401
402			if (associated(LJMap%atidToLJtype)) then
403			deallocate(LJMap%atidToLJtype)
404			LJMap%atidToLJtype => null()
405			end if
406
407	chuckv	492	haveMixingMap = .false.
408	gezelter	938
409	chuckv	492	end subroutine destroyLJTypes
410
411	gezelter	938	subroutine getLJfunc(r, myPot, myDeriv)
412
413			real(kind=dp), intent(in) :: r
414			real(kind=dp), intent(inout) :: myPot, myDeriv
415			real(kind=dp) :: ri, ri2, ri6, ri7, ri12, ri13
416			real(kind=dp) :: a, b, c, d, dx
417			integer :: j
418
419	chrisfen	942	ri = 1.0_DP / r
420			ri2 = ri*ri
421			ri6 = ri2ri2ri2
422			ri7 = ri6*ri
423			ri12 = ri6*ri6
424			ri13 = ri12*ri
425
426			myPot = 4.0_DP * (ri12 - ri6)
427			myDeriv = 24.0_DP * (ri7 - 2.0_DP * ri13)
428
429	gezelter	938	return
430			end subroutine getLJfunc
431
432			subroutine getSoftFunc(r, myPot, myDeriv)
433
434			real(kind=dp), intent(in) :: r
435			real(kind=dp), intent(inout) :: myPot, myDeriv
436			real(kind=dp) :: ri, ri2, ri6, ri7
437			real(kind=dp) :: a, b, c, d, dx
438			integer :: j
439
440	chrisfen	942	ri = 1.0_DP / r
441			ri2 = ri*ri
442			ri6 = ri2ri2ri2
443			ri7 = ri6*ri
444			myPot = 4.0_DP * (ri6)
445			myDeriv = - 24.0_DP * ri7
446
447	gezelter	938	return
448			end subroutine getSoftFunc
449
450	gezelter	115	end module lj