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.5 2004-11-18 15:57:16 chrisfen Exp $, $Date: 2004-11-18 15:57:16 $, $Name: not supported by cvs2svn $, $Revision: 1.5 $

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
    
    haveMixingMap = .true.

  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:	1754
Committed:	Thu Nov 18 15:57:16 2004 UTC (19 years, 7 months ago) by chrisfen
File size:	9873 byte(s)
Log Message:	Fixed a mixing list bug that slowed down the force loop
#	User	Rev	Content
1	gezelter	1608	!! Calculates Long Range forces Lennard-Jones interactions.
2			!! @author Charles F. Vardeman II
3			!! @author Matthew Meineke
4	chrisfen	1754	!! @version $Id: LJ.F90,v 1.5 2004-11-18 15:57:16 chrisfen Exp $, $Date: 2004-11-18 15:57:16 $, $Name: not supported by cvs2svn $, $Revision: 1.5 $
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	chrisfen	1754	haveMixingMap = .true.
229
230	gezelter	1628	end subroutine createMixingMap
231
232	gezelter	1608	subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
233			pot, f, do_pot)
234
235			integer, intent(in) :: atom1, atom2
236			real( kind = dp ), intent(in) :: rij, r2
237			real( kind = dp ) :: pot, sw, vpair
238			real( kind = dp ), dimension(3,nLocal) :: f
239			real( kind = dp ), intent(in), dimension(3) :: d
240			real( kind = dp ), intent(inout), dimension(3) :: fpair
241			logical, intent(in) :: do_pot
242
243			! local Variables
244			real( kind = dp ) :: drdx, drdy, drdz
245			real( kind = dp ) :: fx, fy, fz
246			real( kind = dp ) :: pot_temp, dudr
247			real( kind = dp ) :: sigma6
248			real( kind = dp ) :: epsilon
249			real( kind = dp ) :: r6
250			real( kind = dp ) :: t6
251			real( kind = dp ) :: t12
252			real( kind = dp ) :: delta
253	gezelter	1628	integer :: id1, id2, localError
254	gezelter	1608
255	gezelter	1628	if (.not.haveMixingMap) then
256			localError = 0
257			call createMixingMap(localError)
258			if ( localError .ne. 0 ) then
259			call handleError("LJ", "MixingMap creation failed!")
260			return
261			end if
262			endif
263
264	gezelter	1608	! Look up the correct parameters in the mixing matrix
265			#ifdef IS_MPI
266	gezelter	1628	sigma6 = MixingMap(atid_Row(atom1),atid_Col(atom2))%sigma6
267			epsilon = MixingMap(atid_Row(atom1),atid_Col(atom2))%epsilon
268			delta = MixingMap(atid_Row(atom1),atid_Col(atom2))%delta
269	gezelter	1608	#else
270	gezelter	1628	sigma6 = MixingMap(atid(atom1),atid(atom2))%sigma6
271			epsilon = MixingMap(atid(atom1),atid(atom2))%epsilon
272			delta = MixingMap(atid(atom1),atid(atom2))%delta
273	gezelter	1608	#endif
274
275			r6 = r2 * r2 * r2
276
277			t6 = sigma6/ r6
278			t12 = t6 * t6
279
280			pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
281			if (LJ_do_shift) then
282			pot_temp = pot_temp + delta
283			endif
284
285			vpair = vpair + pot_temp
286
287			dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
288
289			drdx = d(1) / rij
290			drdy = d(2) / rij
291			drdz = d(3) / rij
292
293			fx = dudr * drdx
294			fy = dudr * drdy
295			fz = dudr * drdz
296
297
298			#ifdef IS_MPI
299			if (do_pot) then
300			pot_Row(atom1) = pot_Row(atom1) + swpot_temp0.5
301			pot_Col(atom2) = pot_Col(atom2) + swpot_temp0.5
302			endif
303
304			f_Row(1,atom1) = f_Row(1,atom1) + fx
305			f_Row(2,atom1) = f_Row(2,atom1) + fy
306			f_Row(3,atom1) = f_Row(3,atom1) + fz
307
308			f_Col(1,atom2) = f_Col(1,atom2) - fx
309			f_Col(2,atom2) = f_Col(2,atom2) - fy
310			f_Col(3,atom2) = f_Col(3,atom2) - fz
311
312			#else
313			if (do_pot) pot = pot + sw*pot_temp
314
315			f(1,atom1) = f(1,atom1) + fx
316			f(2,atom1) = f(2,atom1) + fy
317			f(3,atom1) = f(3,atom1) + fz
318
319			f(1,atom2) = f(1,atom2) - fx
320			f(2,atom2) = f(2,atom2) - fy
321			f(3,atom2) = f(3,atom2) - fz
322			#endif
323
324			#ifdef IS_MPI
325			id1 = AtomRowToGlobal(atom1)
326			id2 = AtomColToGlobal(atom2)
327			#else
328			id1 = atom1
329			id2 = atom2
330			#endif
331
332			if (molMembershipList(id1) .ne. molMembershipList(id2)) then
333
334			fpair(1) = fpair(1) + fx
335			fpair(2) = fpair(2) + fy
336			fpair(3) = fpair(3) + fz
337
338			endif
339
340			return
341
342			end subroutine do_lj_pair
343
344
345			!! Calculates the mixing for sigma or epslon
346
347			end module lj
348	chuckv	1624
349	gezelter	1628	subroutine newLJtype(ident, sigma, epsilon, status)
350			use lj, ONLY : module_newLJtype => newLJtype
351			integer, parameter :: DP = selected_real_kind(15)
352			integer,intent(inout) :: ident
353			real(kind=dp),intent(inout) :: sigma
354			real(kind=dp),intent(inout) :: epsilon
355			integer,intent(inout) :: status
356
357			call module_newLJtype(ident, sigma, epsilon, status)
358
359			end subroutine newLJtype
360	chuckv	1624
361	gezelter	1628	subroutine useGeometricMixing()
362			use lj, ONLY: module_useGeometricMixing => useGeometricMixing
363
364			call module_useGeometricMixing()
365			return
366			end subroutine useGeometricMixing