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.29 2008-07-14 12:35:55 gezelter Exp $, $Date: 2008-07-14 12:35:55 $, $Name: not supported by cvs2svn $, $Revision: 1.29 $


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 :: 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, 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, shiftedFrc
    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
    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 = epsilon * (myPot - myPotC)
    vpair = vpair + pot_temp
    dudr = sw * 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

#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:	3432
Committed:	Mon Jul 14 12:35:58 2008 UTC (17 years ago) by gezelter
File size:	13582 byte(s)
Log Message:	Changes for implementing Amber force field: Added Inversions and worked on BaseAtomTypes so that they'd function with the fortran side.
#	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	3432	!! @version $Id: LJ.F90,v 1.29 2008-07-14 12:35:55 gezelter Exp $, $Date: 2008-07-14 12:35:55 $, $Name: not supported by cvs2svn $, $Revision: 1.29 $
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	chrisfen	3129	logical, save :: defaultShiftPot = .false.
71			logical, save :: defaultShiftFrc = .false.
72	gezelter	2271	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	chrisfen	3129	logical :: shiftedFrc
99	gezelter	1628	end type MixParameters
100	gezelter	2204
101	gezelter	1628	type(MixParameters), dimension(:,:), allocatable :: MixingMap
102	gezelter	2204
103	gezelter	2271	public :: newLJtype
104			public :: setLJDefaultCutoff
105			public :: getSigma
106			public :: getEpsilon
107	gezelter	1628	public :: do_lj_pair
108	gezelter	2271	public :: destroyLJtypes
109	gezelter	2204
110	gezelter	1608	contains
111
112	gezelter	2271	subroutine newLJtype(c_ident, sigma, epsilon, isSoftCore, status)
113	gezelter	1930	integer,intent(in) :: c_ident
114	gezelter	1628	real(kind=dp),intent(in) :: sigma
115			real(kind=dp),intent(in) :: epsilon
116	gezelter	2271	integer, intent(in) :: isSoftCore
117	chuckv	1624	integer,intent(out) :: status
118	gezelter	2271	integer :: nLJTypes, ntypes, myATID
119	gezelter	1930	integer, pointer :: MatchList(:) => null()
120	gezelter	2271	integer :: current
121	gezelter	1628
122	chuckv	1624	status = 0
123	gezelter	2271	! check to see if this is the first time into this routine...
124			if (.not.associated(LJMap%LJtypes)) then
125	gezelter	2204
126	gezelter	2271	call getMatchingElementList(atypes, "is_LennardJones", .true., &
127			nLJTypes, MatchList)
128
129			LJMap%nLJtypes = nLJTypes
130
131			allocate(LJMap%LJtypes(nLJTypes))
132
133			ntypes = getSize(atypes)
134
135			allocate(LJMap%atidToLJtype(ntypes))
136			end if
137
138			LJMap%currentLJtype = LJMap%currentLJtype + 1
139			current = LJMap%currentLJtype
140
141	gezelter	1930	myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
142	gezelter	3432
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	chrisfen	3129	subroutine setLJDefaultCutoff(thisRcut, shiftedPot, shiftedFrc)
155	gezelter	2271	real(kind=dp), intent(in) :: thisRcut
156			logical, intent(in) :: shiftedPot
157	chrisfen	3129	logical, intent(in) :: shiftedFrc
158	gezelter	2271	defaultCutoff = thisRcut
159	chrisfen	3129	defaultShiftPot = shiftedPot
160			defaultShiftFrc = shiftedFrc
161	gezelter	2271	haveDefaultCutoff = .true.
162	chrisfen	3129
163	chrisfen	2727	!we only want to build LJ Mixing map if LJ is being used.
164	chuckv	2319	if(LJMap%nLJTypes /= 0) then
165			call createMixingMap()
166			end if
167	gezelter	2717
168	gezelter	2271	end subroutine setLJDefaultCutoff
169	gezelter	2204
170	gezelter	1633	function getSigma(atid) result (s)
171			integer, intent(in) :: atid
172	gezelter	2271	integer :: ljt1
173	gezelter	1633	real(kind=dp) :: s
174	gezelter	2204
175	gezelter	2271	if (LJMap%currentLJtype == 0) then
176			call handleError("LJ", "No members in LJMap")
177	gezelter	1633	return
178			end if
179	gezelter	2204
180	gezelter	2271	ljt1 = LJMap%atidToLJtype(atid)
181			s = LJMap%LJtypes(ljt1)%sigma
182
183	gezelter	1633	end function getSigma
184
185			function getEpsilon(atid) result (e)
186			integer, intent(in) :: atid
187	gezelter	2271	integer :: ljt1
188	gezelter	1633	real(kind=dp) :: e
189	gezelter	2204
190	gezelter	2271	if (LJMap%currentLJtype == 0) then
191			call handleError("LJ", "No members in LJMap")
192	gezelter	1633	return
193			end if
194	gezelter	2204
195	gezelter	2271	ljt1 = LJMap%atidToLJtype(atid)
196			e = LJMap%LJtypes(ljt1)%epsilon
197
198	gezelter	1633	end function getEpsilon
199
200	gezelter	2271	subroutine createMixingMap()
201			integer :: nLJtypes, i, j
202			real ( kind = dp ) :: s1, s2, e1, e2
203			real ( kind = dp ) :: rcut6, tp6, tp12
204	gezelter	2289	logical :: isSoftCore1, isSoftCore2, doShift
205	gezelter	1628
206	gezelter	2271	if (LJMap%currentLJtype == 0) then
207			call handleError("LJ", "No members in LJMap")
208	gezelter	2086	return
209	gezelter	2271	end if
210	gezelter	2204
211	gezelter	2271	nLJtypes = LJMap%nLJtypes
212	gezelter	2204
213	gezelter	2086	if (.not. allocated(MixingMap)) then
214	gezelter	2271	allocate(MixingMap(nLJtypes, nLJtypes))
215	gezelter	2086	endif
216
217	chuckv	2497	useGeometricDistanceMixing = usesGeometricDistanceMixing()
218	gezelter	2271	do i = 1, nLJtypes
219	gezelter	2204
220	gezelter	2271	s1 = LJMap%LJtypes(i)%sigma
221			e1 = LJMap%LJtypes(i)%epsilon
222			isSoftCore1 = LJMap%LJtypes(i)%isSoftCore
223	gezelter	2204
224	gezelter	2271	do j = i, nLJtypes
225
226			s2 = LJMap%LJtypes(j)%sigma
227			e2 = LJMap%LJtypes(j)%epsilon
228			isSoftCore2 = LJMap%LJtypes(j)%isSoftCore
229
230			MixingMap(i,j)%isSoftCore = isSoftCore1 .or. isSoftCore2
231	gezelter	2204
232	gezelter	2271	! only the distance parameter uses different mixing policies
233			if (useGeometricDistanceMixing) then
234	gezelter	2756	MixingMap(i,j)%sigma = sqrt(s1 * s2)
235	gezelter	2271	else
236			MixingMap(i,j)%sigma = 0.5_dp * (s1 + s2)
237			endif
238
239	gezelter	2756	MixingMap(i,j)%epsilon = sqrt(e1 * e2)
240	gezelter	2204
241	gezelter	2717	MixingMap(i,j)%sigmai = 1.0_DP / (MixingMap(i,j)%sigma)
242	gezelter	2204
243	gezelter	2461	if (haveDefaultCutoff) then
244	chrisfen	3129	MixingMap(i,j)%shiftedPot = defaultShiftPot
245			MixingMap(i,j)%shiftedFrc = defaultShiftFrc
246	gezelter	2271	else
247	chrisfen	3129	MixingMap(i,j)%shiftedPot = defaultShiftPot
248			MixingMap(i,j)%shiftedFrc = defaultShiftFrc
249	gezelter	2461	endif
250	gezelter	2204
251	gezelter	2289	if (i.ne.j) then
252			MixingMap(j,i)%sigma = MixingMap(i,j)%sigma
253			MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
254	gezelter	2717	MixingMap(j,i)%sigmai = MixingMap(i,j)%sigmai
255	gezelter	2289	MixingMap(j,i)%rCut = MixingMap(i,j)%rCut
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	chrisfen	3129	MixingMap(j,i)%shiftedFrc = MixingMap(i,j)%shiftedFrc
260	gezelter	2289	endif
261
262	gezelter	2271	enddo
263			enddo
264
265	chrisfen	1754	haveMixingMap = .true.
266	chrisfen	3129
267	gezelter	1628	end subroutine createMixingMap
268	chrisfen	2727
269	gezelter	2461	subroutine do_lj_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
270	gezelter	1608	pot, f, do_pot)
271	gezelter	2461
272	gezelter	1608	integer, intent(in) :: atom1, atom2
273	gezelter	2271	integer :: atid1, atid2, ljt1, ljt2
274	gezelter	2461	real( kind = dp ), intent(in) :: rij, r2, rcut
275	gezelter	1608	real( kind = dp ) :: pot, sw, vpair
276			real( kind = dp ), dimension(3,nLocal) :: f
277			real( kind = dp ), intent(in), dimension(3) :: d
278			real( kind = dp ), intent(inout), dimension(3) :: fpair
279			logical, intent(in) :: do_pot
280
281			! local Variables
282			real( kind = dp ) :: drdx, drdy, drdz
283			real( kind = dp ) :: fx, fy, fz
284	gezelter	2717	real( kind = dp ) :: myPot, myPotC, myDeriv, myDerivC, ros, rcos
285	gezelter	1608	real( kind = dp ) :: pot_temp, dudr
286	gezelter	2717	real( kind = dp ) :: sigmai
287	gezelter	1608	real( kind = dp ) :: epsilon
288	chrisfen	3129	logical :: isSoftCore, shiftedPot, shiftedFrc
289	gezelter	1628	integer :: id1, id2, localError
290	gezelter	1608
291	gezelter	1628	if (.not.haveMixingMap) then
292	gezelter	2271	call createMixingMap()
293	gezelter	1628	endif
294
295	gezelter	1608	! Look up the correct parameters in the mixing matrix
296			#ifdef IS_MPI
297	gezelter	2271	atid1 = atid_Row(atom1)
298			atid2 = atid_Col(atom2)
299	gezelter	1608	#else
300	gezelter	2271	atid1 = atid(atom1)
301			atid2 = atid(atom2)
302	gezelter	1608	#endif
303
304	gezelter	2271	ljt1 = LJMap%atidToLJtype(atid1)
305			ljt2 = LJMap%atidToLJtype(atid2)
306
307	gezelter	2717	sigmai = MixingMap(ljt1,ljt2)%sigmai
308	gezelter	2271	epsilon = MixingMap(ljt1,ljt2)%epsilon
309			isSoftCore = MixingMap(ljt1,ljt2)%isSoftCore
310			shiftedPot = MixingMap(ljt1,ljt2)%shiftedPot
311	chrisfen	3129	shiftedFrc = MixingMap(ljt1,ljt2)%shiftedFrc
312	gezelter	2271
313	gezelter	2717	ros = rij * sigmai
314	gezelter	2204
315	gezelter	2717	if (isSoftCore) then
316	tim	2062
317	gezelter	2717	call getSoftFunc(ros, myPot, myDeriv)
318
319	gezelter	2271	if (shiftedPot) then
320	gezelter	2717	rcos = rcut * sigmai
321			call getSoftFunc(rcos, myPotC, myDerivC)
322	chrisfen	3129	myDerivC = 0.0_dp
323			elseif (shiftedFrc) then
324			rcos = rcut * sigmai
325			call getSoftFunc(rcos, myPotC, myDerivC)
326			myPotC = myPotC + myDerivC * (rij - rcut) * sigmai
327			else
328			myPotC = 0.0_dp
329			myDerivC = 0.0_dp
330	gezelter	2271	endif
331	gezelter	2717
332	tim	2062	else
333	gezelter	2717
334			call getLJfunc(ros, myPot, myDeriv)
335
336	gezelter	2271	if (shiftedPot) then
337	gezelter	2717	rcos = rcut * sigmai
338	gezelter	3126	call getLJfunc(rcos, myPotC, myDerivC)
339	chrisfen	3129	myDerivC = 0.0_dp
340			elseif (shiftedFrc) then
341			rcos = rcut * sigmai
342			call getLJfunc(rcos, myPotC, myDerivC)
343			myPotC = myPotC + myDerivC * (rij - rcut) * sigmai
344			else
345			myPotC = 0.0_dp
346			myDerivC = 0.0_dp
347	gezelter	2204	endif
348	gezelter	2271
349	gezelter	1608	endif
350	gezelter	2204
351	chrisfen	3129	pot_temp = epsilon * (myPot - myPotC)
352			vpair = vpair + pot_temp
353	gezelter	3127	dudr = sw * epsilon * (myDeriv - myDerivC) * sigmai
354	gezelter	2717
355	gezelter	1608	drdx = d(1) / rij
356			drdy = d(2) / rij
357			drdz = d(3) / rij
358	gezelter	2204
359	gezelter	1608	fx = dudr * drdx
360			fy = dudr * drdy
361			fz = dudr * drdz
362	gezelter	2204
363	gezelter	1608	#ifdef IS_MPI
364			if (do_pot) then
365	gezelter	2361	pot_Row(VDW_POT,atom1) = pot_Row(VDW_POT,atom1) + swpot_temp0.5
366			pot_Col(VDW_POT,atom2) = pot_Col(VDW_POT,atom2) + swpot_temp0.5
367	gezelter	1608	endif
368	gezelter	2204
369	gezelter	1608	f_Row(1,atom1) = f_Row(1,atom1) + fx
370			f_Row(2,atom1) = f_Row(2,atom1) + fy
371			f_Row(3,atom1) = f_Row(3,atom1) + fz
372	gezelter	2204
373	gezelter	1608	f_Col(1,atom2) = f_Col(1,atom2) - fx
374			f_Col(2,atom2) = f_Col(2,atom2) - fy
375			f_Col(3,atom2) = f_Col(3,atom2) - fz
376	gezelter	2204
377	gezelter	1608	#else
378			if (do_pot) pot = pot + sw*pot_temp
379
380			f(1,atom1) = f(1,atom1) + fx
381			f(2,atom1) = f(2,atom1) + fy
382			f(3,atom1) = f(3,atom1) + fz
383	gezelter	2204
384	gezelter	1608	f(1,atom2) = f(1,atom2) - fx
385			f(2,atom2) = f(2,atom2) - fy
386			f(3,atom2) = f(3,atom2) - fz
387			#endif
388	gezelter	2204
389	gezelter	1608	#ifdef IS_MPI
390			id1 = AtomRowToGlobal(atom1)
391			id2 = AtomColToGlobal(atom2)
392			#else
393			id1 = atom1
394			id2 = atom2
395			#endif
396
397			if (molMembershipList(id1) .ne. molMembershipList(id2)) then
398	gezelter	2204
399	gezelter	1608	fpair(1) = fpair(1) + fx
400			fpair(2) = fpair(2) + fy
401			fpair(3) = fpair(3) + fz
402
403			endif
404
405			return
406	gezelter	2204
407	gezelter	1608	end subroutine do_lj_pair
408	gezelter	2204
409	chuckv	2189	subroutine destroyLJTypes()
410	gezelter	2271
411			LJMap%nLJtypes = 0
412			LJMap%currentLJtype = 0
413
414			if (associated(LJMap%LJtypes)) then
415			deallocate(LJMap%LJtypes)
416			LJMap%LJtypes => null()
417			end if
418
419			if (associated(LJMap%atidToLJtype)) then
420			deallocate(LJMap%atidToLJtype)
421			LJMap%atidToLJtype => null()
422			end if
423
424	chuckv	2189	haveMixingMap = .false.
425	gezelter	2717
426	chuckv	2189	end subroutine destroyLJTypes
427
428	gezelter	2717	subroutine getLJfunc(r, myPot, myDeriv)
429
430			real(kind=dp), intent(in) :: r
431			real(kind=dp), intent(inout) :: myPot, myDeriv
432			real(kind=dp) :: ri, ri2, ri6, ri7, ri12, ri13
433			real(kind=dp) :: a, b, c, d, dx
434			integer :: j
435
436	chrisfen	2727	ri = 1.0_DP / r
437			ri2 = ri*ri
438			ri6 = ri2ri2ri2
439			ri7 = ri6*ri
440			ri12 = ri6*ri6
441			ri13 = ri12*ri
442
443			myPot = 4.0_DP * (ri12 - ri6)
444			myDeriv = 24.0_DP * (ri7 - 2.0_DP * ri13)
445
446	gezelter	2717	return
447			end subroutine getLJfunc
448
449			subroutine getSoftFunc(r, myPot, myDeriv)
450
451			real(kind=dp), intent(in) :: r
452			real(kind=dp), intent(inout) :: myPot, myDeriv
453			real(kind=dp) :: ri, ri2, ri6, ri7
454			real(kind=dp) :: a, b, c, d, dx
455			integer :: j
456
457	chrisfen	2727	ri = 1.0_DP / r
458			ri2 = ri*ri
459			ri6 = ri2ri2ri2
460			ri7 = ri6*ri
461			myPot = 4.0_DP * (ri6)
462			myDeriv = - 24.0_DP * ri7
463
464	gezelter	2717	return
465			end subroutine getSoftFunc
466
467	gezelter	1608	end module lj