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.23 2006-04-21 19:32:07 chrisfen Exp $, $Date: 2006-04-21 19:32:07 $, $Name: not supported by cvs2svn $, $Revision: 1.23 $


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) :: 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 = dsqrt(s1 * s2)
          else
             MixingMap(i,j)%sigma = 0.5_dp * (s1 + s2)
          endif
          
          MixingMap(i,j)%epsilon = dsqrt(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

    !write(*,*) rij, ros, rcos, myPot, myDeriv, myPotC

    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:	2727
Committed:	Fri Apr 21 19:32:07 2006 UTC (18 years, 4 months ago) by chrisfen
File size:	12732 byte(s)
Log Message:	Removed sqrt splines and shape stuff (doesn't work, so why was it in there?). Changed some of the water samples to use shifted_force. Probably should set the defaults to use the damped version now that we sped it up.
#	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	chrisfen	2727	!! @version $Id: LJ.F90,v 1.23 2006-04-21 19:32:07 chrisfen Exp $, $Date: 2006-04-21 19:32:07 $, $Name: not supported by cvs2svn $, $Revision: 1.23 $
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			type, private :: LJtype
75			integer :: atid
76	gezelter	1628	real(kind=dp) :: sigma
77			real(kind=dp) :: epsilon
78	gezelter	2271	logical :: isSoftCore = .false.
79			end type LJtype
80	gezelter	2204
81	gezelter	2271	type, private :: LJList
82	chuckv	2319	integer :: Nljtypes = 0
83	gezelter	2271	integer :: currentLJtype = 0
84			type(LJtype), pointer :: LJtypes(:) => null()
85			integer, pointer :: atidToLJtype(:) => null()
86			end type LJList
87	gezelter	2204
88	gezelter	2271	type(LJList), save :: LJMap
89	gezelter	2204
90	gezelter	1628	type :: MixParameters
91			real(kind=DP) :: sigma
92			real(kind=DP) :: epsilon
93	gezelter	2717	real(kind=dp) :: sigmai
94	gezelter	2271	real(kind=dp) :: rCut
95			logical :: rCutWasSet = .false.
96			logical :: shiftedPot
97			logical :: isSoftCore = .false.
98	gezelter	1628	end type MixParameters
99	gezelter	2204
100	gezelter	1628	type(MixParameters), dimension(:,:), allocatable :: MixingMap
101	gezelter	2204
102	gezelter	2271	public :: newLJtype
103			public :: setLJDefaultCutoff
104			public :: getSigma
105			public :: getEpsilon
106	gezelter	1628	public :: do_lj_pair
107	gezelter	2271	public :: destroyLJtypes
108	gezelter	2204
109	gezelter	1608	contains
110
111	gezelter	2271	subroutine newLJtype(c_ident, sigma, epsilon, isSoftCore, status)
112	gezelter	1930	integer,intent(in) :: c_ident
113	gezelter	1628	real(kind=dp),intent(in) :: sigma
114			real(kind=dp),intent(in) :: epsilon
115	gezelter	2271	integer, intent(in) :: isSoftCore
116	chuckv	1624	integer,intent(out) :: status
117	gezelter	2271	integer :: nLJTypes, ntypes, myATID
118	gezelter	1930	integer, pointer :: MatchList(:) => null()
119	gezelter	2271	integer :: current
120	gezelter	1628
121	chuckv	1624	status = 0
122	gezelter	2271	! check to see if this is the first time into this routine...
123			if (.not.associated(LJMap%LJtypes)) then
124	gezelter	2204
125	gezelter	2271	call getMatchingElementList(atypes, "is_LennardJones", .true., &
126			nLJTypes, MatchList)
127
128			LJMap%nLJtypes = nLJTypes
129
130			allocate(LJMap%LJtypes(nLJTypes))
131
132			ntypes = getSize(atypes)
133
134			allocate(LJMap%atidToLJtype(ntypes))
135			end if
136
137			LJMap%currentLJtype = LJMap%currentLJtype + 1
138			current = LJMap%currentLJtype
139
140	gezelter	1930	myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
141	gezelter	2271	LJMap%atidToLJtype(myATID) = current
142			LJMap%LJtypes(current)%atid = myATID
143			LJMap%LJtypes(current)%sigma = sigma
144			LJMap%LJtypes(current)%epsilon = epsilon
145			if (isSoftCore .eq. 1) then
146			LJMap%LJtypes(current)%isSoftCore = .true.
147			else
148			LJMap%LJtypes(current)%isSoftCore = .false.
149			endif
150			end subroutine newLJtype
151	chuckv	1624
152	gezelter	2271	subroutine setLJDefaultCutoff(thisRcut, shiftedPot)
153			real(kind=dp), intent(in) :: thisRcut
154			logical, intent(in) :: shiftedPot
155			defaultCutoff = thisRcut
156			defaultShift = shiftedPot
157			haveDefaultCutoff = .true.
158	chrisfen	2727	!we only want to build LJ Mixing map if LJ is being used.
159	chuckv	2319	if(LJMap%nLJTypes /= 0) then
160			call createMixingMap()
161			end if
162	gezelter	2717
163	gezelter	2271	end subroutine setLJDefaultCutoff
164	gezelter	2204
165	gezelter	1633	function getSigma(atid) result (s)
166			integer, intent(in) :: atid
167	gezelter	2271	integer :: ljt1
168	gezelter	1633	real(kind=dp) :: s
169	gezelter	2204
170	gezelter	2271	if (LJMap%currentLJtype == 0) then
171			call handleError("LJ", "No members in LJMap")
172	gezelter	1633	return
173			end if
174	gezelter	2204
175	gezelter	2271	ljt1 = LJMap%atidToLJtype(atid)
176			s = LJMap%LJtypes(ljt1)%sigma
177
178	gezelter	1633	end function getSigma
179
180			function getEpsilon(atid) result (e)
181			integer, intent(in) :: atid
182	gezelter	2271	integer :: ljt1
183	gezelter	1633	real(kind=dp) :: e
184	gezelter	2204
185	gezelter	2271	if (LJMap%currentLJtype == 0) then
186			call handleError("LJ", "No members in LJMap")
187	gezelter	1633	return
188			end if
189	gezelter	2204
190	gezelter	2271	ljt1 = LJMap%atidToLJtype(atid)
191			e = LJMap%LJtypes(ljt1)%epsilon
192
193	gezelter	1633	end function getEpsilon
194
195	gezelter	2271	subroutine createMixingMap()
196			integer :: nLJtypes, i, j
197			real ( kind = dp ) :: s1, s2, e1, e2
198			real ( kind = dp ) :: rcut6, tp6, tp12
199	gezelter	2289	logical :: isSoftCore1, isSoftCore2, doShift
200	gezelter	1628
201	gezelter	2271	if (LJMap%currentLJtype == 0) then
202			call handleError("LJ", "No members in LJMap")
203	gezelter	2086	return
204	gezelter	2271	end if
205	gezelter	2204
206	gezelter	2271	nLJtypes = LJMap%nLJtypes
207	gezelter	2204
208	gezelter	2086	if (.not. allocated(MixingMap)) then
209	gezelter	2271	allocate(MixingMap(nLJtypes, nLJtypes))
210	gezelter	2086	endif
211
212	chuckv	2497	useGeometricDistanceMixing = usesGeometricDistanceMixing()
213	gezelter	2271	do i = 1, nLJtypes
214	gezelter	2204
215	gezelter	2271	s1 = LJMap%LJtypes(i)%sigma
216			e1 = LJMap%LJtypes(i)%epsilon
217			isSoftCore1 = LJMap%LJtypes(i)%isSoftCore
218	gezelter	2204
219	gezelter	2271	do j = i, nLJtypes
220
221			s2 = LJMap%LJtypes(j)%sigma
222			e2 = LJMap%LJtypes(j)%epsilon
223			isSoftCore2 = LJMap%LJtypes(j)%isSoftCore
224
225			MixingMap(i,j)%isSoftCore = isSoftCore1 .or. isSoftCore2
226	gezelter	2204
227	gezelter	2271	! only the distance parameter uses different mixing policies
228			if (useGeometricDistanceMixing) then
229			MixingMap(i,j)%sigma = dsqrt(s1 * s2)
230			else
231			MixingMap(i,j)%sigma = 0.5_dp * (s1 + s2)
232			endif
233
234			MixingMap(i,j)%epsilon = dsqrt(e1 * e2)
235	gezelter	2204
236	gezelter	2717	MixingMap(i,j)%sigmai = 1.0_DP / (MixingMap(i,j)%sigma)
237	gezelter	2204
238	gezelter	2461	if (haveDefaultCutoff) then
239			MixingMap(i,j)%shiftedPot = defaultShift
240	gezelter	2271	else
241	gezelter	2461	MixingMap(i,j)%shiftedPot = defaultShift
242			endif
243	gezelter	2204
244	gezelter	2289	if (i.ne.j) then
245			MixingMap(j,i)%sigma = MixingMap(i,j)%sigma
246			MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
247	gezelter	2717	MixingMap(j,i)%sigmai = MixingMap(i,j)%sigmai
248	gezelter	2289	MixingMap(j,i)%rCut = MixingMap(i,j)%rCut
249			MixingMap(j,i)%rCutWasSet = MixingMap(i,j)%rCutWasSet
250			MixingMap(j,i)%shiftedPot = MixingMap(i,j)%shiftedPot
251			MixingMap(j,i)%isSoftCore = MixingMap(i,j)%isSoftCore
252			endif
253
254	gezelter	2271	enddo
255			enddo
256
257	chrisfen	1754	haveMixingMap = .true.
258	gezelter	2271
259	gezelter	1628	end subroutine createMixingMap
260	chrisfen	2727
261	gezelter	2461	subroutine do_lj_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
262	gezelter	1608	pot, f, do_pot)
263	gezelter	2461
264	gezelter	1608	integer, intent(in) :: atom1, atom2
265	gezelter	2271	integer :: atid1, atid2, ljt1, ljt2
266	gezelter	2461	real( kind = dp ), intent(in) :: rij, r2, rcut
267	gezelter	1608	real( kind = dp ) :: pot, sw, vpair
268			real( kind = dp ), dimension(3,nLocal) :: f
269			real( kind = dp ), intent(in), dimension(3) :: d
270			real( kind = dp ), intent(inout), dimension(3) :: fpair
271			logical, intent(in) :: do_pot
272
273			! local Variables
274			real( kind = dp ) :: drdx, drdy, drdz
275			real( kind = dp ) :: fx, fy, fz
276	gezelter	2717	real( kind = dp ) :: myPot, myPotC, myDeriv, myDerivC, ros, rcos
277	gezelter	1608	real( kind = dp ) :: pot_temp, dudr
278	gezelter	2717	real( kind = dp ) :: sigmai
279	gezelter	1608	real( kind = dp ) :: epsilon
280	gezelter	2271	logical :: isSoftCore, shiftedPot
281	gezelter	1628	integer :: id1, id2, localError
282	gezelter	1608
283	gezelter	1628	if (.not.haveMixingMap) then
284	gezelter	2271	call createMixingMap()
285	gezelter	1628	endif
286
287	gezelter	1608	! Look up the correct parameters in the mixing matrix
288			#ifdef IS_MPI
289	gezelter	2271	atid1 = atid_Row(atom1)
290			atid2 = atid_Col(atom2)
291	gezelter	1608	#else
292	gezelter	2271	atid1 = atid(atom1)
293			atid2 = atid(atom2)
294	gezelter	1608	#endif
295
296	gezelter	2271	ljt1 = LJMap%atidToLJtype(atid1)
297			ljt2 = LJMap%atidToLJtype(atid2)
298
299	gezelter	2717	sigmai = MixingMap(ljt1,ljt2)%sigmai
300	gezelter	2271	epsilon = MixingMap(ljt1,ljt2)%epsilon
301			isSoftCore = MixingMap(ljt1,ljt2)%isSoftCore
302			shiftedPot = MixingMap(ljt1,ljt2)%shiftedPot
303
304	gezelter	2717	ros = rij * sigmai
305			myPotC = 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			call getLJfunc(rcos, myPotC, myDerivC)
323	gezelter	2204	endif
324	gezelter	2271
325	gezelter	1608	endif
326	gezelter	2204
327	gezelter	2717	!write(,) rij, ros, rcos, myPot, myDeriv, myPotC
328
329			pot_temp = epsilon * (myPot - myPotC)
330			vpair = vpair + pot_temp
331			dudr = sw * epsilon * myDeriv * sigmai
332
333	gezelter	1608	drdx = d(1) / rij
334			drdy = d(2) / rij
335			drdz = d(3) / rij
336	gezelter	2204
337	gezelter	1608	fx = dudr * drdx
338			fy = dudr * drdy
339			fz = dudr * drdz
340	gezelter	2204
341	gezelter	1608	#ifdef IS_MPI
342			if (do_pot) then
343	gezelter	2361	pot_Row(VDW_POT,atom1) = pot_Row(VDW_POT,atom1) + swpot_temp0.5
344			pot_Col(VDW_POT,atom2) = pot_Col(VDW_POT,atom2) + swpot_temp0.5
345	gezelter	1608	endif
346	gezelter	2204
347	gezelter	1608	f_Row(1,atom1) = f_Row(1,atom1) + fx
348			f_Row(2,atom1) = f_Row(2,atom1) + fy
349			f_Row(3,atom1) = f_Row(3,atom1) + fz
350	gezelter	2204
351	gezelter	1608	f_Col(1,atom2) = f_Col(1,atom2) - fx
352			f_Col(2,atom2) = f_Col(2,atom2) - fy
353			f_Col(3,atom2) = f_Col(3,atom2) - fz
354	gezelter	2204
355	gezelter	1608	#else
356			if (do_pot) pot = pot + sw*pot_temp
357
358			f(1,atom1) = f(1,atom1) + fx
359			f(2,atom1) = f(2,atom1) + fy
360			f(3,atom1) = f(3,atom1) + fz
361	gezelter	2204
362	gezelter	1608	f(1,atom2) = f(1,atom2) - fx
363			f(2,atom2) = f(2,atom2) - fy
364			f(3,atom2) = f(3,atom2) - fz
365			#endif
366	gezelter	2204
367	gezelter	1608	#ifdef IS_MPI
368			id1 = AtomRowToGlobal(atom1)
369			id2 = AtomColToGlobal(atom2)
370			#else
371			id1 = atom1
372			id2 = atom2
373			#endif
374
375			if (molMembershipList(id1) .ne. molMembershipList(id2)) then
376	gezelter	2204
377	gezelter	1608	fpair(1) = fpair(1) + fx
378			fpair(2) = fpair(2) + fy
379			fpair(3) = fpair(3) + fz
380
381			endif
382
383			return
384	gezelter	2204
385	gezelter	1608	end subroutine do_lj_pair
386	gezelter	2204
387	chuckv	2189	subroutine destroyLJTypes()
388	gezelter	2271
389			LJMap%nLJtypes = 0
390			LJMap%currentLJtype = 0
391
392			if (associated(LJMap%LJtypes)) then
393			deallocate(LJMap%LJtypes)
394			LJMap%LJtypes => null()
395			end if
396
397			if (associated(LJMap%atidToLJtype)) then
398			deallocate(LJMap%atidToLJtype)
399			LJMap%atidToLJtype => null()
400			end if
401
402	chuckv	2189	haveMixingMap = .false.
403	gezelter	2717
404	chuckv	2189	end subroutine destroyLJTypes
405
406	gezelter	2717	subroutine getLJfunc(r, myPot, myDeriv)
407
408			real(kind=dp), intent(in) :: r
409			real(kind=dp), intent(inout) :: myPot, myDeriv
410			real(kind=dp) :: ri, ri2, ri6, ri7, ri12, ri13
411			real(kind=dp) :: a, b, c, d, dx
412			integer :: j
413
414	chrisfen	2727	ri = 1.0_DP / r
415			ri2 = ri*ri
416			ri6 = ri2ri2ri2
417			ri7 = ri6*ri
418			ri12 = ri6*ri6
419			ri13 = ri12*ri
420
421			myPot = 4.0_DP * (ri12 - ri6)
422			myDeriv = 24.0_DP * (ri7 - 2.0_DP * ri13)
423
424	gezelter	2717	return
425			end subroutine getLJfunc
426
427			subroutine getSoftFunc(r, myPot, myDeriv)
428
429			real(kind=dp), intent(in) :: r
430			real(kind=dp), intent(inout) :: myPot, myDeriv
431			real(kind=dp) :: ri, ri2, ri6, ri7
432			real(kind=dp) :: a, b, c, d, dx
433			integer :: j
434
435	chrisfen	2727	ri = 1.0_DP / r
436			ri2 = ri*ri
437			ri6 = ri2ri2ri2
438			ri7 = ri6*ri
439			myPot = 4.0_DP * (ri6)
440			myDeriv = - 24.0_DP * ri7
441
442	gezelter	2717	return
443			end subroutine getSoftFunc
444
445	gezelter	1608	end module lj