UseTheForce/DarkSide/LJ.F90

!! Calculates Long Range forces Lennard-Jones interactions.
!! @author Charles F. Vardeman II
!! @author Matthew Meineke
!! @version $Id: LJ.F90,v 1.4 2004-10-22 20:22:47 gezelter Exp $, $Date: 2004-10-22 20:22:47 $, $Name: not supported by cvs2svn $, $Revision: 1.4 $

module lj
  use atype_module
  use switcheroo
  use vector_class
  use simulation
  use status
#ifdef IS_MPI
  use mpiSimulation
#endif
  use force_globals

  implicit none
  PRIVATE
  
  integer, parameter :: DP = selected_real_kind(15)
  
  type, private :: LjType
     integer :: ident
     real(kind=dp) :: sigma
     real(kind=dp) :: epsilon
  end type LjType
  
  type(LjType), dimension(:), allocatable :: ParameterMap
  
  logical, save :: haveMixingMap = .false.
  
  type :: MixParameters
     real(kind=DP) :: sigma
     real(kind=DP) :: epsilon
     real(kind=dp)  :: sigma6
     real(kind=dp)  :: tp6
     real(kind=dp)  :: tp12
     real(kind=dp)  :: delta 
  end type MixParameters
  
  type(MixParameters), dimension(:,:), allocatable :: MixingMap
  
  real(kind=DP), save :: LJ_rcut
  logical, save :: have_rcut = .false.
  logical, save :: LJ_do_shift = .false.
  logical, save :: useGeometricDistanceMixing = .false.
   
  !! Public methods and data
  
  public :: setCutoffLJ
  public :: useGeometricMixing
  public :: do_lj_pair
  public :: newLJtype  
  public :: getSigma
  public :: getEpsilon
  
contains

  subroutine newLJtype(ident, sigma, epsilon, status)
    integer,intent(in) :: ident
    real(kind=dp),intent(in) :: sigma
    real(kind=dp),intent(in) :: epsilon
    integer,intent(out) :: status
    integer :: nAtypes

    status = 0
    
    !! Be simple-minded and assume that we need a ParameterMap that
    !! is the same size as the total number of atom types

    if (.not.allocated(ParameterMap)) then
       
       nAtypes = getSize(atypes)
    
       if (nAtypes == 0) then
          status = -1
          return
       end if
       
       if (.not. allocated(ParameterMap)) then
          allocate(ParameterMap(nAtypes))
       endif
       
    end if

    if (ident .gt. size(ParameterMap)) then
       status = -1
       return
    endif
    
    ! set the values for ParameterMap for this atom type:

    ParameterMap(ident)%ident = ident
    ParameterMap(ident)%epsilon = epsilon
    ParameterMap(ident)%sigma = sigma
    
  end subroutine newLJtype

  function getSigma(atid) result (s)
    integer, intent(in) :: atid
    integer :: localError
    real(kind=dp) :: s
    
    if (.not.allocated(ParameterMap)) then
       call handleError("LJ", "no ParameterMap was present before first call of getSigma!")
       return
    end if
    
    s = ParameterMap(atid)%sigma
  end function getSigma

  function getEpsilon(atid) result (e)
    integer, intent(in) :: atid
    integer :: localError
    real(kind=dp) :: e
    
    if (.not.allocated(ParameterMap)) then
       call handleError("dipole-dipole", "no ParameterMap was present before first call of getEpsilon!")
       return
    end if
    
    e = ParameterMap(atid)%epsilon
  end function getEpsilon


  subroutine setCutoffLJ(rcut, do_shift, status)
    logical, intent(in):: do_shift
    integer :: status, myStatus
    real(kind=dp) :: rcut

#define __FORTRAN90
#include "UseTheForce/fSwitchingFunction.h"

    status = 0

    LJ_rcut = rcut
    LJ_do_shift = do_shift
    call set_switch(LJ_SWITCH, rcut, rcut)
    have_rcut = .true.
    
    return
  end subroutine setCutoffLJ

  subroutine useGeometricMixing() 
    useGeometricDistanceMixing = .true.
    haveMixingMap = .false.
    return
  end subroutine useGeometricMixing
  
  subroutine createMixingMap(status)
    integer :: nAtypes
    integer :: status
    integer :: i
    integer :: j
    real ( kind = dp ) :: Sigma_i, Sigma_j
    real ( kind = dp ) :: Epsilon_i, Epsilon_j
    real ( kind = dp ) :: rcut6

    status = 0
    
    nAtypes = size(ParameterMap)
    
    if (nAtypes == 0) then
       status = -1
       return
    end if

    if (.not.have_rcut) then
       status = -1
       return
    endif
    
    if (.not. allocated(MixingMap)) then
       allocate(MixingMap(nAtypes, nAtypes))
    endif
    
    rcut6 = LJ_rcut**6
    
    ! This loops through all atypes, even those that don't support LJ forces.
    do i = 1, nAtypes
       
       Epsilon_i = ParameterMap(i)%epsilon
       Sigma_i = ParameterMap(i)%sigma
       
       ! do self mixing rule
       MixingMap(i,i)%sigma   = Sigma_i          
       MixingMap(i,i)%sigma6  = Sigma_i ** 6          
       MixingMap(i,i)%tp6     = (MixingMap(i,i)%sigma6)/rcut6          
       MixingMap(i,i)%tp12    = (MixingMap(i,i)%tp6) ** 2
       MixingMap(i,i)%epsilon = Epsilon_i          
       MixingMap(i,i)%delta   = -4.0_DP * MixingMap(i,i)%epsilon * &
            (MixingMap(i,i)%tp12 - MixingMap(i,i)%tp6)
       
       do j = i + 1, nAtypes
          
          Epsilon_j = ParameterMap(j)%epsilon
          Sigma_j = ParameterMap(j)%sigma
          
          ! only the distance parameter uses different mixing policies
          if (useGeometricDistanceMixing) then
             ! only for OPLS as far as we can tell
             MixingMap(i,j)%sigma = dsqrt(Sigma_i * Sigma_j)
          else
             ! everyone else
             MixingMap(i,j)%sigma = 0.5_dp * (Sigma_i + Sigma_j)
          endif
          
          ! energy parameter is always geometric mean:
          MixingMap(i,j)%epsilon = dsqrt(Epsilon_i * Epsilon_j)
                    
          MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
          MixingMap(i,j)%tp6    = MixingMap(i,j)%sigma6/rcut6
          MixingMap(i,j)%tp12    = (MixingMap(i,j)%tp6) ** 2
          
          MixingMap(i,j)%delta = -4.0_DP * MixingMap(i,j)%epsilon * &
               (MixingMap(i,j)%tp12 - MixingMap(i,j)%tp6)
          
          MixingMap(j,i)%sigma   = MixingMap(i,j)%sigma
          MixingMap(j,i)%sigma6  = MixingMap(i,j)%sigma6
          MixingMap(j,i)%tp6     = MixingMap(i,j)%tp6
          MixingMap(j,i)%tp12    = MixingMap(i,j)%tp12
          MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
          MixingMap(j,i)%delta   = MixingMap(i,j)%delta
          
       end do
    end do
    
  end subroutine createMixingMap
        
  subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
       pot, f, do_pot)

    integer, intent(in) ::  atom1, atom2
    real( kind = dp ), intent(in) :: rij, r2
    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
    real( kind = dp ) :: t6
    real( kind = dp ) :: t12
    real( kind = dp ) :: delta
    integer :: id1, id2, localError

    if (.not.haveMixingMap) then
       localError = 0
       call createMixingMap(localError)
       if ( localError .ne. 0 ) then
          call handleError("LJ", "MixingMap creation failed!")
          return
       end if
    endif

    ! Look up the correct parameters in the mixing matrix
#ifdef IS_MPI
    sigma6   = MixingMap(atid_Row(atom1),atid_Col(atom2))%sigma6
    epsilon  = MixingMap(atid_Row(atom1),atid_Col(atom2))%epsilon
    delta    = MixingMap(atid_Row(atom1),atid_Col(atom2))%delta
#else
    sigma6   = MixingMap(atid(atom1),atid(atom2))%sigma6
    epsilon  = MixingMap(atid(atom1),atid(atom2))%epsilon
    delta    = MixingMap(atid(atom1),atid(atom2))%delta
#endif

    r6 = r2 * r2 * r2
    
    t6  = sigma6/ r6
    t12 = t6 * t6     
   
    pot_temp = 4.0E0_DP * epsilon * (t12 - t6) 
    if (LJ_do_shift) then
       pot_temp = pot_temp + delta
    endif

    vpair = vpair + pot_temp
       
    dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
       
    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(atom1) = pot_Row(atom1) + sw*pot_temp*0.5
       pot_Col(atom2) = pot_Col(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
  
  
  !! Calculates the mixing for sigma or epslon
    
end module lj

subroutine newLJtype(ident, sigma, epsilon, status)
  use lj, ONLY : module_newLJtype => newLJtype
  integer, parameter :: DP = selected_real_kind(15)
  integer,intent(inout) :: ident
  real(kind=dp),intent(inout) :: sigma
  real(kind=dp),intent(inout) :: epsilon
  integer,intent(inout) :: status
  
  call module_newLJtype(ident, sigma, epsilon, status)
  
end subroutine newLJtype

subroutine useGeometricMixing()
  use lj, ONLY: module_useGeometricMixing => useGeometricMixing
  
  call module_useGeometricMixing()
  return
end subroutine useGeometricMixing
Revision:	1633
Committed:	Fri Oct 22 20:22:48 2004 UTC (19 years, 10 months ago) by gezelter
File size:	9845 byte(s)
Log Message:	Added un-busticated fortran files and c/Fortran interfaces
#	User	Rev	Content
1	gezelter	1608	!! Calculates Long Range forces Lennard-Jones interactions.
2			!! @author Charles F. Vardeman II
3			!! @author Matthew Meineke
4	gezelter	1633	!! @version $Id: LJ.F90,v 1.4 2004-10-22 20:22:47 gezelter Exp $, $Date: 2004-10-22 20:22:47 $, $Name: not supported by cvs2svn $, $Revision: 1.4 $
5	gezelter	1608
6			module lj
7			use atype_module
8			use switcheroo
9			use vector_class
10			use simulation
11	gezelter	1628	use status
12	gezelter	1608	#ifdef IS_MPI
13			use mpiSimulation
14			#endif
15			use force_globals
16
17			implicit none
18			PRIVATE
19	chuckv	1624
20			integer, parameter :: DP = selected_real_kind(15)
21	gezelter	1608
22	gezelter	1628	type, private :: LjType
23			integer :: ident
24			real(kind=dp) :: sigma
25			real(kind=dp) :: epsilon
26			end type LjType
27	gezelter	1608
28	gezelter	1628	type(LjType), dimension(:), allocatable :: ParameterMap
29	gezelter	1608
30	gezelter	1628	logical, save :: haveMixingMap = .false.
31	gezelter	1608
32	gezelter	1628	type :: MixParameters
33			real(kind=DP) :: sigma
34			real(kind=DP) :: epsilon
35			real(kind=dp) :: sigma6
36			real(kind=dp) :: tp6
37			real(kind=dp) :: tp12
38			real(kind=dp) :: delta
39			end type MixParameters
40	chuckv	1624
41	gezelter	1628	type(MixParameters), dimension(:,:), allocatable :: MixingMap
42	chuckv	1624
43	gezelter	1628	real(kind=DP), save :: LJ_rcut
44			logical, save :: have_rcut = .false.
45			logical, save :: LJ_do_shift = .false.
46			logical, save :: useGeometricDistanceMixing = .false.
47
48			!! Public methods and data
49	chuckv	1624
50	gezelter	1628	public :: setCutoffLJ
51			public :: useGeometricMixing
52			public :: do_lj_pair
53			public :: newLJtype
54	gezelter	1633	public :: getSigma
55			public :: getEpsilon
56	chuckv	1624
57	gezelter	1608	contains
58
59	gezelter	1628	subroutine newLJtype(ident, sigma, epsilon, status)
60	chuckv	1624	integer,intent(in) :: ident
61	gezelter	1628	real(kind=dp),intent(in) :: sigma
62			real(kind=dp),intent(in) :: epsilon
63	chuckv	1624	integer,intent(out) :: status
64			integer :: nAtypes
65	gezelter	1628
66	chuckv	1624	status = 0
67
68	gezelter	1628	!! Be simple-minded and assume that we need a ParameterMap that
69			!! is the same size as the total number of atom types
70	chuckv	1624
71	gezelter	1628	if (.not.allocated(ParameterMap)) then
72
73			nAtypes = getSize(atypes)
74
75	chuckv	1624	if (nAtypes == 0) then
76	gezelter	1628	status = -1
77			return
78	chuckv	1624	end if
79	gezelter	1628
80			if (.not. allocated(ParameterMap)) then
81			allocate(ParameterMap(nAtypes))
82			endif
83
84	chuckv	1624	end if
85
86	gezelter	1628	if (ident .gt. size(ParameterMap)) then
87			status = -1
88			return
89			endif
90	chuckv	1624
91	gezelter	1628	! set the values for ParameterMap for this atom type:
92
93			ParameterMap(ident)%ident = ident
94			ParameterMap(ident)%epsilon = epsilon
95			ParameterMap(ident)%sigma = sigma
96	chuckv	1624
97			end subroutine newLJtype
98	gezelter	1633
99			function getSigma(atid) result (s)
100			integer, intent(in) :: atid
101			integer :: localError
102			real(kind=dp) :: s
103	gezelter	1608
104	gezelter	1633	if (.not.allocated(ParameterMap)) then
105			call handleError("LJ", "no ParameterMap was present before first call of getSigma!")
106			return
107			end if
108
109			s = ParameterMap(atid)%sigma
110			end function getSigma
111
112			function getEpsilon(atid) result (e)
113			integer, intent(in) :: atid
114			integer :: localError
115			real(kind=dp) :: e
116
117			if (.not.allocated(ParameterMap)) then
118			call handleError("dipole-dipole", "no ParameterMap was present before first call of getEpsilon!")
119			return
120			end if
121
122			e = ParameterMap(atid)%epsilon
123			end function getEpsilon
124
125
126	gezelter	1608	subroutine setCutoffLJ(rcut, do_shift, status)
127			logical, intent(in):: do_shift
128			integer :: status, myStatus
129			real(kind=dp) :: rcut
130
131			#define __FORTRAN90
132			#include "UseTheForce/fSwitchingFunction.h"
133
134			status = 0
135
136			LJ_rcut = rcut
137			LJ_do_shift = do_shift
138			call set_switch(LJ_SWITCH, rcut, rcut)
139	gezelter	1628	have_rcut = .true.
140	gezelter	1608
141			return
142			end subroutine setCutoffLJ
143	gezelter	1628
144			subroutine useGeometricMixing()
145			useGeometricDistanceMixing = .true.
146			haveMixingMap = .false.
147			return
148			end subroutine useGeometricMixing
149	gezelter	1608
150	gezelter	1628	subroutine createMixingMap(status)
151	gezelter	1608	integer :: nAtypes
152			integer :: status
153			integer :: i
154			integer :: j
155	gezelter	1628	real ( kind = dp ) :: Sigma_i, Sigma_j
156			real ( kind = dp ) :: Epsilon_i, Epsilon_j
157	gezelter	1608	real ( kind = dp ) :: rcut6
158	gezelter	1628
159	gezelter	1608	status = 0
160
161	gezelter	1628	nAtypes = size(ParameterMap)
162
163	gezelter	1608	if (nAtypes == 0) then
164			status = -1
165			return
166			end if
167	gezelter	1628
168			if (.not.have_rcut) then
169			status = -1
170			return
171	gezelter	1608	endif
172	gezelter	1628
173			if (.not. allocated(MixingMap)) then
174			allocate(MixingMap(nAtypes, nAtypes))
175			endif
176
177	gezelter	1608	rcut6 = LJ_rcut**6
178	gezelter	1628
179			! This loops through all atypes, even those that don't support LJ forces.
180	gezelter	1608	do i = 1, nAtypes
181	gezelter	1628
182			Epsilon_i = ParameterMap(i)%epsilon
183			Sigma_i = ParameterMap(i)%sigma
184
185			! do self mixing rule
186			MixingMap(i,i)%sigma = Sigma_i
187			MixingMap(i,i)%sigma6 = Sigma_i ** 6
188			MixingMap(i,i)%tp6 = (MixingMap(i,i)%sigma6)/rcut6
189			MixingMap(i,i)%tp12 = (MixingMap(i,i)%tp6) ** 2
190			MixingMap(i,i)%epsilon = Epsilon_i
191			MixingMap(i,i)%delta = -4.0_DP * MixingMap(i,i)%epsilon * &
192			(MixingMap(i,i)%tp12 - MixingMap(i,i)%tp6)
193
194			do j = i + 1, nAtypes
195	chuckv	1624
196	gezelter	1628	Epsilon_j = ParameterMap(j)%epsilon
197			Sigma_j = ParameterMap(j)%sigma
198	gezelter	1608
199	gezelter	1628	! only the distance parameter uses different mixing policies
200			if (useGeometricDistanceMixing) then
201			! only for OPLS as far as we can tell
202			MixingMap(i,j)%sigma = dsqrt(Sigma_i * Sigma_j)
203			else
204			! everyone else
205			MixingMap(i,j)%sigma = 0.5_dp * (Sigma_i + Sigma_j)
206			endif
207	gezelter	1608
208	gezelter	1628	! energy parameter is always geometric mean:
209			MixingMap(i,j)%epsilon = dsqrt(Epsilon_i * Epsilon_j)
210
211			MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
212			MixingMap(i,j)%tp6 = MixingMap(i,j)%sigma6/rcut6
213			MixingMap(i,j)%tp12 = (MixingMap(i,j)%tp6) ** 2
214	gezelter	1608
215	gezelter	1628	MixingMap(i,j)%delta = -4.0_DP * MixingMap(i,j)%epsilon * &
216			(MixingMap(i,j)%tp12 - MixingMap(i,j)%tp6)
217	gezelter	1608
218	gezelter	1628	MixingMap(j,i)%sigma = MixingMap(i,j)%sigma
219			MixingMap(j,i)%sigma6 = MixingMap(i,j)%sigma6
220			MixingMap(j,i)%tp6 = MixingMap(i,j)%tp6
221			MixingMap(j,i)%tp12 = MixingMap(i,j)%tp12
222			MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
223			MixingMap(j,i)%delta = MixingMap(i,j)%delta
224	gezelter	1608
225	gezelter	1628	end do
226	gezelter	1608	end do
227
228	gezelter	1628	end subroutine createMixingMap
229
230	gezelter	1608	subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
231			pot, f, do_pot)
232
233			integer, intent(in) :: atom1, atom2
234			real( kind = dp ), intent(in) :: rij, r2
235			real( kind = dp ) :: pot, sw, vpair
236			real( kind = dp ), dimension(3,nLocal) :: f
237			real( kind = dp ), intent(in), dimension(3) :: d
238			real( kind = dp ), intent(inout), dimension(3) :: fpair
239			logical, intent(in) :: do_pot
240
241			! local Variables
242			real( kind = dp ) :: drdx, drdy, drdz
243			real( kind = dp ) :: fx, fy, fz
244			real( kind = dp ) :: pot_temp, dudr
245			real( kind = dp ) :: sigma6
246			real( kind = dp ) :: epsilon
247			real( kind = dp ) :: r6
248			real( kind = dp ) :: t6
249			real( kind = dp ) :: t12
250			real( kind = dp ) :: delta
251	gezelter	1628	integer :: id1, id2, localError
252	gezelter	1608
253	gezelter	1628	if (.not.haveMixingMap) then
254			localError = 0
255			call createMixingMap(localError)
256			if ( localError .ne. 0 ) then
257			call handleError("LJ", "MixingMap creation failed!")
258			return
259			end if
260			endif
261
262	gezelter	1608	! Look up the correct parameters in the mixing matrix
263			#ifdef IS_MPI
264	gezelter	1628	sigma6 = MixingMap(atid_Row(atom1),atid_Col(atom2))%sigma6
265			epsilon = MixingMap(atid_Row(atom1),atid_Col(atom2))%epsilon
266			delta = MixingMap(atid_Row(atom1),atid_Col(atom2))%delta
267	gezelter	1608	#else
268	gezelter	1628	sigma6 = MixingMap(atid(atom1),atid(atom2))%sigma6
269			epsilon = MixingMap(atid(atom1),atid(atom2))%epsilon
270			delta = MixingMap(atid(atom1),atid(atom2))%delta
271	gezelter	1608	#endif
272
273			r6 = r2 * r2 * r2
274
275			t6 = sigma6/ r6
276			t12 = t6 * t6
277
278			pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
279			if (LJ_do_shift) then
280			pot_temp = pot_temp + delta
281			endif
282
283			vpair = vpair + pot_temp
284
285			dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
286
287			drdx = d(1) / rij
288			drdy = d(2) / rij
289			drdz = d(3) / rij
290
291			fx = dudr * drdx
292			fy = dudr * drdy
293			fz = dudr * drdz
294
295
296			#ifdef IS_MPI
297			if (do_pot) then
298			pot_Row(atom1) = pot_Row(atom1) + swpot_temp0.5
299			pot_Col(atom2) = pot_Col(atom2) + swpot_temp0.5
300			endif
301
302			f_Row(1,atom1) = f_Row(1,atom1) + fx
303			f_Row(2,atom1) = f_Row(2,atom1) + fy
304			f_Row(3,atom1) = f_Row(3,atom1) + fz
305
306			f_Col(1,atom2) = f_Col(1,atom2) - fx
307			f_Col(2,atom2) = f_Col(2,atom2) - fy
308			f_Col(3,atom2) = f_Col(3,atom2) - fz
309
310			#else
311			if (do_pot) pot = pot + sw*pot_temp
312
313			f(1,atom1) = f(1,atom1) + fx
314			f(2,atom1) = f(2,atom1) + fy
315			f(3,atom1) = f(3,atom1) + fz
316
317			f(1,atom2) = f(1,atom2) - fx
318			f(2,atom2) = f(2,atom2) - fy
319			f(3,atom2) = f(3,atom2) - fz
320			#endif
321
322			#ifdef IS_MPI
323			id1 = AtomRowToGlobal(atom1)
324			id2 = AtomColToGlobal(atom2)
325			#else
326			id1 = atom1
327			id2 = atom2
328			#endif
329
330			if (molMembershipList(id1) .ne. molMembershipList(id2)) then
331
332			fpair(1) = fpair(1) + fx
333			fpair(2) = fpair(2) + fy
334			fpair(3) = fpair(3) + fz
335
336			endif
337
338			return
339
340			end subroutine do_lj_pair
341
342
343			!! Calculates the mixing for sigma or epslon
344
345			end module lj
346	chuckv	1624
347	gezelter	1628	subroutine newLJtype(ident, sigma, epsilon, status)
348			use lj, ONLY : module_newLJtype => newLJtype
349			integer, parameter :: DP = selected_real_kind(15)
350			integer,intent(inout) :: ident
351			real(kind=dp),intent(inout) :: sigma
352			real(kind=dp),intent(inout) :: epsilon
353			integer,intent(inout) :: status
354
355			call module_newLJtype(ident, sigma, epsilon, status)
356
357			end subroutine newLJtype
358	chuckv	1624
359	gezelter	1628	subroutine useGeometricMixing()
360			use lj, ONLY: module_useGeometricMixing => useGeometricMixing
361
362			call module_useGeometricMixing()
363			return
364			end subroutine useGeometricMixing