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.20 2005-12-07 19:58:18 chuckv Exp $, $Date: 2005-12-07 19:58:18 $, $Name: not supported by cvs2svn $, $Revision: 1.20 $


module lj
  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"

  integer, parameter :: DP = selected_real_kind(15)

  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) :: sigma6
     real(kind=dp) :: rCut
     real(kind=dp) :: delta
     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 = dsqrt(s1 * s2)
          else
             MixingMap(i,j)%sigma = 0.5_dp * (s1 + s2)
          endif
          
          MixingMap(i,j)%epsilon = dsqrt(e1 * e2)

          MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6

          if (haveDefaultCutoff) then
             rcut6 = defaultCutoff**6
             tp6    = MixingMap(i,j)%sigma6/rcut6
             tp12    = tp6**2          
             MixingMap(i,j)%delta =-4.0_DP*MixingMap(i,j)%epsilon*(tp12 - tp6)
             MixingMap(i,j)%shiftedPot = defaultShift
          else
             MixingMap(i,j)%delta = 0.0_DP
             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)%sigma6     = MixingMap(i,j)%sigma6
             MixingMap(j,i)%rCut       = MixingMap(i,j)%rCut
             MixingMap(j,i)%delta      = MixingMap(i,j)%delta
             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 ) :: pot_temp, dudr
    real( kind = dp ) :: sigma6
    real( kind = dp ) :: epsilon
    real( kind = dp ) :: r6, rc6
    real( kind = dp ) :: t6, tp6
    real( kind = dp ) :: t12, tp12
    real( kind = dp ) :: delta
    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)

    sigma6     = MixingMap(ljt1,ljt2)%sigma6
    epsilon    = MixingMap(ljt1,ljt2)%epsilon
    delta      = MixingMap(ljt1,ljt2)%delta
    isSoftCore = MixingMap(ljt1,ljt2)%isSoftCore
    shiftedPot = MixingMap(ljt1,ljt2)%shiftedPot

    r6 = r2 * r2 * r2

    t6  = sigma6/ r6
    t12 = t6 * t6     

    if (isSoftCore) then
       
       pot_temp = 4.0E0_DP * epsilon * t6
       if (shiftedPot) then
          rc6 = rcut**6
          tp6 = sigma6 / rc6
          delta =-4.0_DP*epsilon*(tp6)
          pot_temp = pot_temp + delta
       endif
       
       vpair = vpair + pot_temp
       
       dudr = -sw * 24.0E0_DP * epsilon * t6 / rij
       
    else
       pot_temp = 4.0E0_DP * epsilon * (t12 - t6) 
       if (shiftedPot) then
          rc6 = rcut**6
          tp6 = sigma6 / rc6
          tp12 = tp6*tp6
          delta =-4.0_DP*epsilon*(tp12 - tp6)
          pot_temp = pot_temp + delta
       endif
       
       vpair = vpair + pot_temp
       
       dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
    endif

    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

end module lj
Revision:	2497
Committed:	Wed Dec 7 19:58:18 2005 UTC (18 years, 9 months ago) by chuckv
File size:	12461 byte(s)
Log Message:	Removed geometric distance mixing from individual modules and now use force options to detemine what the deal is.
#	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	chuckv	2497	!! @version $Id: LJ.F90,v 1.20 2005-12-07 19:58:18 chuckv Exp $, $Date: 2005-12-07 19:58:18 $, $Name: not supported by cvs2svn $, $Revision: 1.20 $
47	gezelter	1608
48	gezelter	1930
49	gezelter	1608	module lj
50			use atype_module
51			use vector_class
52			use simulation
53	gezelter	1628	use status
54	chuckv	2497	use fForceOptions
55	gezelter	1608	#ifdef IS_MPI
56			use mpiSimulation
57			#endif
58			use force_globals
59
60			implicit none
61			PRIVATE
62	chuckv	2355	#define __FORTRAN90
63			#include "UseTheForce/DarkSide/fInteractionMap.h"
64	gezelter	2204
65	chuckv	1624	integer, parameter :: DP = selected_real_kind(15)
66	gezelter	2204
67	gezelter	2271	logical, save :: useGeometricDistanceMixing = .false.
68			logical, save :: haveMixingMap = .false.
69
70			real(kind=DP), save :: defaultCutoff = 0.0_DP
71			logical, save :: defaultShift = .false.
72			logical, save :: haveDefaultCutoff = .false.
73
74
75			type, private :: LJtype
76			integer :: atid
77	gezelter	1628	real(kind=dp) :: sigma
78			real(kind=dp) :: epsilon
79	gezelter	2271	logical :: isSoftCore = .false.
80			end type LJtype
81	gezelter	2204
82	gezelter	2271	type, private :: LJList
83	chuckv	2319	integer :: Nljtypes = 0
84	gezelter	2271	integer :: currentLJtype = 0
85			type(LJtype), pointer :: LJtypes(:) => null()
86			integer, pointer :: atidToLJtype(:) => null()
87			end type LJList
88	gezelter	2204
89	gezelter	2271	type(LJList), save :: LJMap
90	gezelter	2204
91	gezelter	1628	type :: MixParameters
92			real(kind=DP) :: sigma
93			real(kind=DP) :: epsilon
94	gezelter	2271	real(kind=dp) :: sigma6
95			real(kind=dp) :: rCut
96			real(kind=dp) :: delta
97			logical :: rCutWasSet = .false.
98			logical :: shiftedPot
99			logical :: isSoftCore = .false.
100	gezelter	1628	end type MixParameters
101	gezelter	2204
102	gezelter	1628	type(MixParameters), dimension(:,:), allocatable :: MixingMap
103	gezelter	2204
104	gezelter	2271	public :: newLJtype
105			public :: setLJDefaultCutoff
106			public :: getSigma
107			public :: getEpsilon
108	gezelter	1628	public :: do_lj_pair
109	gezelter	2271	public :: destroyLJtypes
110	gezelter	2204
111	gezelter	1608	contains
112
113	gezelter	2271	subroutine newLJtype(c_ident, sigma, epsilon, isSoftCore, status)
114	gezelter	1930	integer,intent(in) :: c_ident
115	gezelter	1628	real(kind=dp),intent(in) :: sigma
116			real(kind=dp),intent(in) :: epsilon
117	gezelter	2271	integer, intent(in) :: isSoftCore
118	chuckv	1624	integer,intent(out) :: status
119	gezelter	2271	integer :: nLJTypes, ntypes, myATID
120	gezelter	1930	integer, pointer :: MatchList(:) => null()
121	gezelter	2271	integer :: current
122	gezelter	1628
123	chuckv	1624	status = 0
124	gezelter	2271	! check to see if this is the first time into this routine...
125			if (.not.associated(LJMap%LJtypes)) then
126	gezelter	2204
127	gezelter	2271	call getMatchingElementList(atypes, "is_LennardJones", .true., &
128			nLJTypes, MatchList)
129
130			LJMap%nLJtypes = nLJTypes
131
132			allocate(LJMap%LJtypes(nLJTypes))
133
134			ntypes = getSize(atypes)
135
136			allocate(LJMap%atidToLJtype(ntypes))
137			end if
138
139			LJMap%currentLJtype = LJMap%currentLJtype + 1
140			current = LJMap%currentLJtype
141
142	gezelter	1930	myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
143	gezelter	2271	LJMap%atidToLJtype(myATID) = current
144			LJMap%LJtypes(current)%atid = myATID
145			LJMap%LJtypes(current)%sigma = sigma
146			LJMap%LJtypes(current)%epsilon = epsilon
147			if (isSoftCore .eq. 1) then
148			LJMap%LJtypes(current)%isSoftCore = .true.
149			else
150			LJMap%LJtypes(current)%isSoftCore = .false.
151			endif
152			end subroutine newLJtype
153	chuckv	1624
154	gezelter	2271	subroutine setLJDefaultCutoff(thisRcut, shiftedPot)
155			real(kind=dp), intent(in) :: thisRcut
156			logical, intent(in) :: shiftedPot
157			defaultCutoff = thisRcut
158			defaultShift = shiftedPot
159			haveDefaultCutoff = .true.
160	chuckv	2319	!we only want to build LJ Mixing map if LJ is being used.
161			if(LJMap%nLJTypes /= 0) then
162			call createMixingMap()
163			end if
164	gezelter	2271	end subroutine setLJDefaultCutoff
165	gezelter	2204
166	gezelter	1633	function getSigma(atid) result (s)
167			integer, intent(in) :: atid
168	gezelter	2271	integer :: ljt1
169	gezelter	1633	real(kind=dp) :: s
170	gezelter	2204
171	gezelter	2271	if (LJMap%currentLJtype == 0) then
172			call handleError("LJ", "No members in LJMap")
173	gezelter	1633	return
174			end if
175	gezelter	2204
176	gezelter	2271	ljt1 = LJMap%atidToLJtype(atid)
177			s = LJMap%LJtypes(ljt1)%sigma
178
179	gezelter	1633	end function getSigma
180
181			function getEpsilon(atid) result (e)
182			integer, intent(in) :: atid
183	gezelter	2271	integer :: ljt1
184	gezelter	1633	real(kind=dp) :: e
185	gezelter	2204
186	gezelter	2271	if (LJMap%currentLJtype == 0) then
187			call handleError("LJ", "No members in LJMap")
188	gezelter	1633	return
189			end if
190	gezelter	2204
191	gezelter	2271	ljt1 = LJMap%atidToLJtype(atid)
192			e = LJMap%LJtypes(ljt1)%epsilon
193
194	gezelter	1633	end function getEpsilon
195
196	gezelter	2271	subroutine createMixingMap()
197			integer :: nLJtypes, i, j
198			real ( kind = dp ) :: s1, s2, e1, e2
199			real ( kind = dp ) :: rcut6, tp6, tp12
200	gezelter	2289	logical :: isSoftCore1, isSoftCore2, doShift
201	gezelter	1628
202	gezelter	2271	if (LJMap%currentLJtype == 0) then
203			call handleError("LJ", "No members in LJMap")
204	gezelter	2086	return
205	gezelter	2271	end if
206	gezelter	2204
207	gezelter	2271	nLJtypes = LJMap%nLJtypes
208	gezelter	2204
209	gezelter	2086	if (.not. allocated(MixingMap)) then
210	gezelter	2271	allocate(MixingMap(nLJtypes, nLJtypes))
211	gezelter	2086	endif
212
213	chuckv	2497	useGeometricDistanceMixing = usesGeometricDistanceMixing()
214	gezelter	2271	do i = 1, nLJtypes
215	gezelter	2204
216	gezelter	2271	s1 = LJMap%LJtypes(i)%sigma
217			e1 = LJMap%LJtypes(i)%epsilon
218			isSoftCore1 = LJMap%LJtypes(i)%isSoftCore
219	gezelter	2204
220	gezelter	2271	do j = i, nLJtypes
221
222			s2 = LJMap%LJtypes(j)%sigma
223			e2 = LJMap%LJtypes(j)%epsilon
224			isSoftCore2 = LJMap%LJtypes(j)%isSoftCore
225
226			MixingMap(i,j)%isSoftCore = isSoftCore1 .or. isSoftCore2
227	gezelter	2204
228	gezelter	2271	! only the distance parameter uses different mixing policies
229			if (useGeometricDistanceMixing) then
230			MixingMap(i,j)%sigma = dsqrt(s1 * s2)
231			else
232			MixingMap(i,j)%sigma = 0.5_dp * (s1 + s2)
233			endif
234
235			MixingMap(i,j)%epsilon = dsqrt(e1 * e2)
236	gezelter	2204
237	gezelter	2271	MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
238	gezelter	2204
239	gezelter	2461	if (haveDefaultCutoff) then
240			rcut6 = defaultCutoff**6
241			tp6 = MixingMap(i,j)%sigma6/rcut6
242			tp12 = tp6**2
243			MixingMap(i,j)%delta =-4.0_DPMixingMap(i,j)%epsilon(tp12 - tp6)
244			MixingMap(i,j)%shiftedPot = defaultShift
245	gezelter	2271	else
246	gezelter	2461	MixingMap(i,j)%delta = 0.0_DP
247			MixingMap(i,j)%shiftedPot = defaultShift
248			endif
249	gezelter	2204
250	gezelter	2289	if (i.ne.j) then
251			MixingMap(j,i)%sigma = MixingMap(i,j)%sigma
252			MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
253			MixingMap(j,i)%sigma6 = MixingMap(i,j)%sigma6
254			MixingMap(j,i)%rCut = MixingMap(i,j)%rCut
255			MixingMap(j,i)%delta = MixingMap(i,j)%delta
256			MixingMap(j,i)%rCutWasSet = MixingMap(i,j)%rCutWasSet
257			MixingMap(j,i)%shiftedPot = MixingMap(i,j)%shiftedPot
258			MixingMap(j,i)%isSoftCore = MixingMap(i,j)%isSoftCore
259			endif
260
261	gezelter	2271	enddo
262			enddo
263
264	chrisfen	1754	haveMixingMap = .true.
265	gezelter	2271
266	gezelter	1628	end subroutine createMixingMap
267	gezelter	2271
268	gezelter	2461	subroutine do_lj_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
269	gezelter	1608	pot, f, do_pot)
270	gezelter	2461
271	gezelter	1608	integer, intent(in) :: atom1, atom2
272	gezelter	2271	integer :: atid1, atid2, ljt1, ljt2
273	gezelter	2461	real( kind = dp ), intent(in) :: rij, r2, rcut
274	gezelter	1608	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			real( kind = dp ) :: pot_temp, dudr
284			real( kind = dp ) :: sigma6
285			real( kind = dp ) :: epsilon
286	gezelter	2461	real( kind = dp ) :: r6, rc6
287			real( kind = dp ) :: t6, tp6
288			real( kind = dp ) :: t12, tp12
289	gezelter	1608	real( kind = dp ) :: delta
290	gezelter	2271	logical :: isSoftCore, shiftedPot
291	gezelter	1628	integer :: id1, id2, localError
292	gezelter	1608
293	gezelter	1628	if (.not.haveMixingMap) then
294	gezelter	2271	call createMixingMap()
295	gezelter	1628	endif
296
297	gezelter	1608	! Look up the correct parameters in the mixing matrix
298			#ifdef IS_MPI
299	gezelter	2271	atid1 = atid_Row(atom1)
300			atid2 = atid_Col(atom2)
301	gezelter	1608	#else
302	gezelter	2271	atid1 = atid(atom1)
303			atid2 = atid(atom2)
304	gezelter	1608	#endif
305
306	gezelter	2271	ljt1 = LJMap%atidToLJtype(atid1)
307			ljt2 = LJMap%atidToLJtype(atid2)
308
309			sigma6 = MixingMap(ljt1,ljt2)%sigma6
310			epsilon = MixingMap(ljt1,ljt2)%epsilon
311			delta = MixingMap(ljt1,ljt2)%delta
312			isSoftCore = MixingMap(ljt1,ljt2)%isSoftCore
313			shiftedPot = MixingMap(ljt1,ljt2)%shiftedPot
314
315	gezelter	1608	r6 = r2 * r2 * r2
316	gezelter	2204
317	gezelter	1608	t6 = sigma6/ r6
318			t12 = t6 * t6
319	tim	2062
320	gezelter	2271	if (isSoftCore) then
321
322			pot_temp = 4.0E0_DP * epsilon * t6
323			if (shiftedPot) then
324	gezelter	2461	rc6 = rcut**6
325			tp6 = sigma6 / rc6
326			delta =-4.0_DPepsilon(tp6)
327	gezelter	2271	pot_temp = pot_temp + delta
328			endif
329
330			vpair = vpair + pot_temp
331
332			dudr = -sw * 24.0E0_DP * epsilon * t6 / rij
333
334	tim	2062	else
335	gezelter	2204	pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
336	gezelter	2271	if (shiftedPot) then
337	gezelter	2461	rc6 = rcut**6
338			tp6 = sigma6 / rc6
339			tp12 = tp6*tp6
340			delta =-4.0_DPepsilon(tp12 - tp6)
341	gezelter	2204	pot_temp = pot_temp + delta
342			endif
343	gezelter	2271
344	gezelter	2204	vpair = vpair + pot_temp
345	gezelter	2271
346	gezelter	2204	dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
347	gezelter	1608	endif
348	gezelter	2204
349	gezelter	1608	drdx = d(1) / rij
350			drdy = d(2) / rij
351			drdz = d(3) / rij
352	gezelter	2204
353	gezelter	1608	fx = dudr * drdx
354			fy = dudr * drdy
355			fz = dudr * drdz
356	gezelter	2204
357	gezelter	1608	#ifdef IS_MPI
358			if (do_pot) then
359	gezelter	2361	pot_Row(VDW_POT,atom1) = pot_Row(VDW_POT,atom1) + swpot_temp0.5
360			pot_Col(VDW_POT,atom2) = pot_Col(VDW_POT,atom2) + swpot_temp0.5
361	gezelter	1608	endif
362	gezelter	2204
363	gezelter	1608	f_Row(1,atom1) = f_Row(1,atom1) + fx
364			f_Row(2,atom1) = f_Row(2,atom1) + fy
365			f_Row(3,atom1) = f_Row(3,atom1) + fz
366	gezelter	2204
367	gezelter	1608	f_Col(1,atom2) = f_Col(1,atom2) - fx
368			f_Col(2,atom2) = f_Col(2,atom2) - fy
369			f_Col(3,atom2) = f_Col(3,atom2) - fz
370	gezelter	2204
371	gezelter	1608	#else
372			if (do_pot) pot = pot + sw*pot_temp
373
374			f(1,atom1) = f(1,atom1) + fx
375			f(2,atom1) = f(2,atom1) + fy
376			f(3,atom1) = f(3,atom1) + fz
377	gezelter	2204
378	gezelter	1608	f(1,atom2) = f(1,atom2) - fx
379			f(2,atom2) = f(2,atom2) - fy
380			f(3,atom2) = f(3,atom2) - fz
381			#endif
382	gezelter	2204
383	gezelter	1608	#ifdef IS_MPI
384			id1 = AtomRowToGlobal(atom1)
385			id2 = AtomColToGlobal(atom2)
386			#else
387			id1 = atom1
388			id2 = atom2
389			#endif
390
391			if (molMembershipList(id1) .ne. molMembershipList(id2)) then
392	gezelter	2204
393	gezelter	1608	fpair(1) = fpair(1) + fx
394			fpair(2) = fpair(2) + fy
395			fpair(3) = fpair(3) + fz
396
397			endif
398
399			return
400	gezelter	2204
401	gezelter	1608	end subroutine do_lj_pair
402	gezelter	2204
403	chuckv	2189	subroutine destroyLJTypes()
404	gezelter	2271
405			LJMap%nLJtypes = 0
406			LJMap%currentLJtype = 0
407
408			if (associated(LJMap%LJtypes)) then
409			deallocate(LJMap%LJtypes)
410			LJMap%LJtypes => null()
411			end if
412
413			if (associated(LJMap%atidToLJtype)) then
414			deallocate(LJMap%atidToLJtype)
415			LJMap%atidToLJtype => null()
416			end if
417
418	chuckv	2189	haveMixingMap = .false.
419			end subroutine destroyLJTypes
420
421	gezelter	1608	end module lj