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.31 2009-10-23 18:41:08 gezelter Exp $, $Date: 2009-10-23 18:41:08 $, $Name: not supported by cvs2svn $, $Revision: 1.31 $


module lj
  use definitions
  use atype_module
  use vector_class
  use simulation
  use status
  use fForceOptions
  use force_globals

  implicit none
  PRIVATE
#define __FORTRAN90
#include "UseTheForce/DarkSide/fInteractionMap.h"

  logical, save :: useGeometricDistanceMixing = .false.
  logical, save :: haveMixingMap = .false.

  real(kind=DP), save :: defaultCutoff = 0.0_DP
  logical, save :: defaultShiftPot = .false.
  logical, save :: defaultShiftFrc = .false.
  logical, save :: haveDefaultCutoff = .false.

  type, private :: LJtype
     integer       :: atid
     real(kind=dp) :: sigma
     real(kind=dp) :: epsilon
     logical       :: isSoftCore = .false.
  end type LJtype

  type, private :: LJList
     integer               :: Nljtypes = 0
     integer               :: currentLJtype = 0
     type(LJtype), pointer :: LJtypes(:)      => null()
     integer, pointer      :: atidToLJtype(:) => null()
  end type LJList

  type(LJList), save :: LJMap

  type :: MixParameters
     real(kind=DP) :: sigma
     real(kind=DP) :: epsilon
     real(kind=dp) :: sigmai
     real(kind=dp) :: rCut
     logical       :: rCutWasSet = .false.
     logical       :: shiftedPot
     logical       :: isSoftCore = .false.
     logical       :: shiftedFrc
  end type MixParameters

  type(MixParameters), dimension(:,:), allocatable :: MixingMap

  public :: newLJtype
  public :: setLJDefaultCutoff
  public :: getSigma
  public :: getEpsilon
  public :: do_lj_pair
  public :: destroyLJtypes

contains

  subroutine newLJtype(c_ident, sigma, epsilon, isSoftCore, status)
    integer,intent(in) :: c_ident
    real(kind=dp),intent(in) :: sigma
    real(kind=dp),intent(in) :: epsilon
    integer, intent(in) :: isSoftCore
    integer,intent(out) :: status
    integer :: nLJTypes, ntypes, myATID
    integer, pointer :: MatchList(:) => null()
    integer :: current

    status = 0
    ! check to see if this is the first time into this routine...
    if (.not.associated(LJMap%LJtypes)) then

       call getMatchingElementList(atypes, "is_LennardJones", .true., &
            nLJTypes, MatchList)
       
       LJMap%nLJtypes =  nLJTypes

       allocate(LJMap%LJtypes(nLJTypes))

       ntypes = getSize(atypes)

       allocate(LJMap%atidToLJtype(ntypes))
    end if

    LJMap%currentLJtype = LJMap%currentLJtype + 1
    current = LJMap%currentLJtype

    myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)

    LJMap%atidToLJtype(myATID)        = current
    LJMap%LJtypes(current)%atid       = myATID
    LJMap%LJtypes(current)%sigma      = sigma
    LJMap%LJtypes(current)%epsilon    = epsilon 
    if (isSoftCore .eq. 1) then
       LJMap%LJtypes(current)%isSoftCore = .true.
    else
       LJMap%LJtypes(current)%isSoftCore = .false.
    endif
  end subroutine newLJtype

  subroutine setLJDefaultCutoff(thisRcut, shiftedPot, shiftedFrc)
    real(kind=dp), intent(in) :: thisRcut
    logical, intent(in) :: shiftedPot
    logical, intent(in) :: shiftedFrc
    defaultCutoff = thisRcut
    defaultShiftPot = shiftedPot
    defaultShiftFrc = shiftedFrc
    haveDefaultCutoff = .true.

    !we only want to build LJ Mixing map if LJ is being used.
    if(LJMap%nLJTypes /= 0) then
       call createMixingMap()
    end if

  end subroutine setLJDefaultCutoff

  function getSigma(atid) result (s)
    integer, intent(in) :: atid
    integer :: ljt1
    real(kind=dp) :: s

    if (LJMap%currentLJtype == 0) then
       call handleError("LJ", "No members in LJMap")
       return
    end if

    ljt1 = LJMap%atidToLJtype(atid)
    s = LJMap%LJtypes(ljt1)%sigma

  end function getSigma

  function getEpsilon(atid) result (e)
    integer, intent(in) :: atid
    integer :: ljt1
    real(kind=dp) :: e

    if (LJMap%currentLJtype == 0) then
       call handleError("LJ", "No members in LJMap")
       return
    end if

    ljt1 = LJMap%atidToLJtype(atid)
    e = LJMap%LJtypes(ljt1)%epsilon

  end function getEpsilon

  subroutine createMixingMap()
    integer :: nLJtypes, i, j
    real ( kind = dp ) :: s1, s2, e1, e2
    real ( kind = dp ) :: rcut6, tp6, tp12
    logical :: isSoftCore1, isSoftCore2, doShift

    if (LJMap%currentLJtype == 0) then
       call handleError("LJ", "No members in LJMap")
       return
    end if

    nLJtypes = LJMap%nLJtypes

    if (.not. allocated(MixingMap)) then
       allocate(MixingMap(nLJtypes, nLJtypes))
    endif

    useGeometricDistanceMixing = usesGeometricDistanceMixing()
    do i = 1, nLJtypes

       s1 = LJMap%LJtypes(i)%sigma
       e1 = LJMap%LJtypes(i)%epsilon
       isSoftCore1 = LJMap%LJtypes(i)%isSoftCore

       do j = i, nLJtypes
          
          s2 = LJMap%LJtypes(j)%sigma
          e2 = LJMap%LJtypes(j)%epsilon
          isSoftCore2 = LJMap%LJtypes(j)%isSoftCore
          
          MixingMap(i,j)%isSoftCore = isSoftCore1 .or. isSoftCore2

          ! only the distance parameter uses different mixing policies
          if (useGeometricDistanceMixing) then
             MixingMap(i,j)%sigma = sqrt(s1 * s2)
          else
             MixingMap(i,j)%sigma = 0.5_dp * (s1 + s2)
          endif
          
          MixingMap(i,j)%epsilon = sqrt(e1 * e2)

          MixingMap(i,j)%sigmai = 1.0_DP  / (MixingMap(i,j)%sigma)

          if (haveDefaultCutoff) then
             MixingMap(i,j)%shiftedPot = defaultShiftPot
             MixingMap(i,j)%shiftedFrc = defaultShiftFrc
          else
             MixingMap(i,j)%shiftedPot = defaultShiftPot
             MixingMap(i,j)%shiftedFrc = defaultShiftFrc
          endif          

          if (i.ne.j) then
             MixingMap(j,i)%sigma      = MixingMap(i,j)%sigma
             MixingMap(j,i)%epsilon    = MixingMap(i,j)%epsilon
             MixingMap(j,i)%sigmai     = MixingMap(i,j)%sigmai
             MixingMap(j,i)%rCut       = MixingMap(i,j)%rCut
             MixingMap(j,i)%rCutWasSet = MixingMap(i,j)%rCutWasSet
             MixingMap(j,i)%shiftedPot = MixingMap(i,j)%shiftedPot
             MixingMap(j,i)%isSoftCore = MixingMap(i,j)%isSoftCore
             MixingMap(j,i)%shiftedFrc = MixingMap(i,j)%shiftedFrc
          endif

       enddo
    enddo
    
    haveMixingMap = .true.

  end subroutine createMixingMap
          
  subroutine do_lj_pair(atom1, atom2, atid1, atid2, d, rij, r2, rcut, sw, vdwMult, &
       vpair, fpair, pot, f1, do_pot)
    
    integer, intent(in) ::  atom1, atom2, atid1, atid2
    integer :: ljt1, ljt2
    real( kind = dp ), intent(in) :: rij, r2, rcut, vdwMult
    real( kind = dp ) :: pot, sw, vpair
    real( kind = dp ), intent(inout), dimension(3) :: f1
    real( kind = dp ), intent(in), dimension(3) :: d
    real( kind = dp ), intent(inout), dimension(3) :: fpair
    logical, intent(in) :: do_pot

    ! local Variables
    real( kind = dp ) :: drdx, drdy, drdz
    real( kind = dp ) :: fx, fy, fz
    real( kind = dp ) :: myPot, myPotC, myDeriv, myDerivC, ros, rcos
    real( kind = dp ) :: pot_temp, dudr
    real( kind = dp ) :: sigmai
    real( kind = dp ) :: epsilon
    logical :: isSoftCore, shiftedPot, shiftedFrc
    integer :: id1, id2, localError

    if (.not.haveMixingMap) then
       call createMixingMap()
    endif

    ljt1 = LJMap%atidToLJtype(atid1)
    ljt2 = LJMap%atidToLJtype(atid2)

    sigmai     = MixingMap(ljt1,ljt2)%sigmai
    epsilon    = MixingMap(ljt1,ljt2)%epsilon
    isSoftCore = MixingMap(ljt1,ljt2)%isSoftCore
    shiftedPot = MixingMap(ljt1,ljt2)%shiftedPot
    shiftedFrc = MixingMap(ljt1,ljt2)%shiftedFrc

    ros = rij * sigmai

    if (isSoftCore) then

       call getSoftFunc(ros, myPot, myDeriv)

       if (shiftedPot) then
          rcos = rcut * sigmai
          call getSoftFunc(rcos, myPotC, myDerivC)
          myDerivC = 0.0_dp
       elseif (shiftedFrc) then
          rcos = rcut * sigmai
          call getSoftFunc(rcos, myPotC, myDerivC)
          myPotC = myPotC + myDerivC * (rij - rcut) * sigmai
       else
          myPotC = 0.0_dp
          myDerivC = 0.0_dp
       endif
              
    else

       call getLJfunc(ros, myPot, myDeriv)

       if (shiftedPot) then
          rcos = rcut * sigmai
          call getLJfunc(rcos, myPotC, myDerivC) 
          myDerivC = 0.0_dp
       elseif (shiftedFrc) then
          rcos = rcut * sigmai
          call getLJfunc(rcos, myPotC, myDerivC)
          myPotC = myPotC + myDerivC * (rij - rcut) * sigmai
       else
          myPotC = 0.0_dp
          myDerivC = 0.0_dp
       endif
       
    endif

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

    drdx = d(1) / rij
    drdy = d(2) / rij
    drdz = d(3) / rij

    fx = dudr * drdx
    fy = dudr * drdy
    fz = dudr * drdz
    
    pot = pot + sw*pot_temp

    f1(1) = fx
    f1(2) = fy
    f1(3) = fz

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