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:	960
Committed:	Wed May 17 15:37:15 2006 UTC (19 years ago) by gezelter
File size:	12639 byte(s)
Log Message:	Getting fortran side prepped for single precision...
#	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	960	!! @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	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	507
306	gezelter	938	if (isSoftCore) then
307	tim	378
308	gezelter	938	call getSoftFunc(ros, myPot, myDeriv)
309
310	gezelter	572	if (shiftedPot) then
311	gezelter	938	rcos = rcut * sigmai
312			call getSoftFunc(rcos, myPotC, myDerivC)
313	gezelter	572	endif
314	gezelter	938
315	tim	378	else
316	gezelter	938
317			call getLJfunc(ros, myPot, myDeriv)
318
319	gezelter	572	if (shiftedPot) then
320	gezelter	938	rcos = rcut * sigmai
321			call getLJfunc(rcos, myPotC, myDerivC)
322	gezelter	507	endif
323	gezelter	572
324	gezelter	115	endif
325	gezelter	507
326	gezelter	938	pot_temp = epsilon * (myPot - myPotC)
327			vpair = vpair + pot_temp
328			dudr = sw * epsilon * myDeriv * sigmai
329
330	gezelter	115	drdx = d(1) / rij
331			drdy = d(2) / rij
332			drdz = d(3) / rij
333	gezelter	507
334	gezelter	115	fx = dudr * drdx
335			fy = dudr * drdy
336			fz = dudr * drdz
337	gezelter	507
338	gezelter	115	#ifdef IS_MPI
339			if (do_pot) then
340	gezelter	662	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	115	endif
343	gezelter	507
344	gezelter	115	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	507
348	gezelter	115	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	507
352	gezelter	115	#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	507
359	gezelter	115	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	507
364	gezelter	115	#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	507
374	gezelter	115	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	507
382	gezelter	115	end subroutine do_lj_pair
383	gezelter	507
384	chuckv	492	subroutine destroyLJTypes()
385	gezelter	572
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	492	haveMixingMap = .false.
400	gezelter	938
401	chuckv	492	end subroutine destroyLJTypes
402
403	gezelter	938	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	942	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	938	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	942	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	938	return
440			end subroutine getSoftFunc
441
442	gezelter	115	end module lj