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:	3127
Committed:	Mon Apr 9 18:24:00 2007 UTC (17 years, 4 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	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	3127	!! @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	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	3127	myDerivC = 0.0_DP
306	gezelter	2204
307	gezelter	2717	if (isSoftCore) then
308	tim	2062
309	gezelter	2717	call getSoftFunc(ros, myPot, myDeriv)
310
311	gezelter	2271	if (shiftedPot) then
312	gezelter	2717	rcos = rcut * sigmai
313			call getSoftFunc(rcos, myPotC, myDerivC)
314	gezelter	2271	endif
315	gezelter	2717
316	tim	2062	else
317	gezelter	2717
318			call getLJfunc(ros, myPot, myDeriv)
319
320	gezelter	2271	if (shiftedPot) then
321	gezelter	2717	rcos = rcut * sigmai
322	gezelter	3126	call getLJfunc(rcos, myPotC, myDerivC)
323	gezelter	2204	endif
324	gezelter	2271
325	gezelter	1608	endif
326	gezelter	2204
327	gezelter	3127	!! 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	2717	vpair = vpair + pot_temp
337
338	gezelter	1608	drdx = d(1) / rij
339			drdy = d(2) / rij
340			drdz = d(3) / rij
341	gezelter	2204
342	gezelter	1608	fx = dudr * drdx
343			fy = dudr * drdy
344			fz = dudr * drdz
345	gezelter	2204
346	gezelter	1608	#ifdef IS_MPI
347			if (do_pot) then
348	gezelter	2361	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	1608	endif
351	gezelter	2204
352	gezelter	1608	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	2204
356	gezelter	1608	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	2204
360	gezelter	1608	#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	2204
367	gezelter	1608	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	2204
372	gezelter	1608	#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	2204
382	gezelter	1608	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	2204
390	gezelter	1608	end subroutine do_lj_pair
391	gezelter	2204
392	chuckv	2189	subroutine destroyLJTypes()
393	gezelter	2271
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	2189	haveMixingMap = .false.
408	gezelter	2717
409	chuckv	2189	end subroutine destroyLJTypes
410
411	gezelter	2717	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	2727	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	2717	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	2727	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	2717	return
448			end subroutine getSoftFunc
449
450	gezelter	1608	end module lj