OOPSE/libmdtools/calc_LJ_FF.F90

!! Calculates Long Range forces Lennard-Jones interactions.
!! Corresponds to the force field defined in lj_FF.cpp 
!! @author Charles F. Vardeman II
!! @author Matthew Meineke
!! @version $Id: calc_LJ_FF.F90,v 1.18 2004-05-07 21:35:04 gezelter Exp $, $Date: 2004-05-07 21:35:04 $, $Name: not supported by cvs2svn $, $Revision: 1.18 $

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

  implicit none
  PRIVATE

#define __FORTRAN90
#include "fForceField.h"

  integer, save :: LJ_Mixing_Policy
  real(kind=DP), save :: LJ_rcut
  logical, save :: havePolicy = .false., haveCut = .false.

  
  !! Logical has lj force field module been initialized?
  
  logical, save :: LJ_FF_initialized = .false.
  
  !! Public methods and data
  public :: init_LJ_FF
  public :: setCutoffLJ
  public :: do_lj_pair
  
  type :: lj_mixed_params
     !! Lennard-Jones epsilon
     real ( kind = dp )  :: epsilon = 0.0_dp
     !! Lennard-Jones Sigma
     real ( kind = dp )  :: sigma = 0.0_dp
     !! Lennard-Jones Sigma to sixth
     real ( kind = dp )  :: sigma6 = 0.0_dp
     !! 
     real ( kind = dp )  :: tp6
     real ( kind = dp )  :: tp12

     real ( kind = dp )  :: delta  = 0.0_dp


  end type lj_mixed_params
  
  type (lj_mixed_params), dimension(:,:), pointer :: ljMixed
  
contains

  subroutine init_LJ_FF(mix_Policy, status)
    integer, intent(in) :: mix_Policy
    integer, intent(out) :: status
    integer :: myStatus
    
    if (mix_Policy == LB_MIXING_RULE) then
       LJ_Mixing_Policy = LB_MIXING_RULE
    else
       if (mix_Policy == EXPLICIT_MIXING_RULE) then
          LJ_Mixing_Policy = EXPLICIT_MIXING_RULE
       else
          write(*,*) 'Unknown Mixing Policy!'
          status = -1
          return
       endif
    endif

    havePolicy = .true.

    if (haveCut) then
       status = 0
       call createMixingList(myStatus)
       if (myStatus /= 0) then
          status = -1
          return
       end if
       
       LJ_FF_initialized = .true.
    end if
  
  end subroutine init_LJ_FF
  
  subroutine setCutoffLJ(rcut, status)
    integer :: status, myStatus
    real(kind=dp) :: rcut

#define __FORTRAN90
#include "fSwitchingFunction.h"

    status = 0

    LJ_rcut = rcut
    call set_switch(LJ_SWITCH, rcut, rcut)

!!$    ! ATTENTION! This is a hardwiring of rcut!
!!$    LJ_rcut = 9.0d0
    haveCut = .true.

    if (havePolicy) then
       status = 0
       call createMixingList(myStatus)
       if (myStatus /= 0) then
          status = -1
          return
       end if
       
       LJ_FF_initialized = .true.
    end if    
    
    return
  end subroutine setCutoffLJ
  
  subroutine createMixingList(status)
    integer :: nAtypes
    integer :: status
    integer :: i
    integer :: j
    real ( kind = dp ) :: mySigma_i,mySigma_j
    real ( kind = dp ) :: myEpsilon_i,myEpsilon_j
    real ( kind = dp ) :: rcut6
    status = 0
    
    nAtypes = getSize(atypes)
    if (nAtypes == 0) then
       status = -1
       return
    end if
        
    if (.not. associated(ljMixed)) then
       allocate(ljMixed(nAtypes, nAtypes))
    endif

    rcut6 = LJ_rcut**6

! This loops through all atypes, even those that don't support LJ forces.
    do i = 1, nAtypes

       call getElementProperty(atypes, i, "lj_epsilon", myEpsilon_i)
       call getElementProperty(atypes, i, "lj_sigma",   mySigma_i)
       ! do self mixing rule
       ljMixed(i,i)%sigma   = mySigma_i
       
       ljMixed(i,i)%sigma6  = (ljMixed(i,i)%sigma) ** 6

       ljMixed(i,i)%tp6     = (ljMixed(i,i)%sigma6)/rcut6

       ljMixed(i,i)%tp12    = (ljMixed(i,i)%tp6) ** 2


       ljMixed(i,i)%epsilon = myEpsilon_i

       ljMixed(i,i)%delta = -4.0_DP * ljMixed(i,i)%epsilon * &
            (ljMixed(i,i)%tp12 - ljMixed(i,i)%tp6)
      
       do j = i + 1, nAtypes
          call getElementProperty(atypes,j,"lj_epsilon",myEpsilon_j)
          call getElementProperty(atypes,j,"lj_sigma",  mySigma_j)
          
          ljMixed(i,j)%sigma  =  &
               calcLJMix("sigma",mySigma_i, &
               mySigma_j)
          
          ljMixed(i,j)%sigma6 = &
               (ljMixed(i,j)%sigma)**6
          
          
          ljMixed(i,j)%tp6     = ljMixed(i,j)%sigma6/rcut6
          
          ljMixed(i,j)%tp12    = (ljMixed(i,j)%tp6) ** 2
          
          
          ljMixed(i,j)%epsilon = &
               calcLJMix("epsilon",myEpsilon_i, &
               myEpsilon_j)
          
          ljMixed(i,j)%delta = -4.0_DP * ljMixed(i,j)%epsilon * &
               (ljMixed(i,j)%tp12 - ljMixed(i,j)%tp6)
          
          
          ljMixed(j,i)%sigma   = ljMixed(i,j)%sigma
          ljMixed(j,i)%sigma6  = ljMixed(i,j)%sigma6
          ljMixed(j,i)%tp6     = ljMixed(i,j)%tp6
          ljMixed(j,i)%tp12    = ljMixed(i,j)%tp12
          ljMixed(j,i)%epsilon = ljMixed(i,j)%epsilon
          ljMixed(j,i)%delta   = ljMixed(i,j)%delta
          
       end do
    end do
    
  end subroutine createMixingList
  
  subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, pot, f, &
       do_pot, do_stress)

    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
    logical, intent(in) :: do_pot, do_stress

    ! 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

    if (rij.lt.LJ_rcut)  then

       ! Look up the correct parameters in the mixing matrix
#ifdef IS_MPI
       sigma6   = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6
       epsilon  = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon
       !delta    = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta
#else
       sigma6   = ljMixed(atid(atom1),atid(atom2))%sigma6
       epsilon  = ljMixed(atid(atom1),atid(atom2))%epsilon
       !delta    = ljMixed(atid(atom1),atid(atom2))%delta
#endif
       
       r6 = r2 * r2 * r2
       
       t6  = sigma6/ r6
       t12 = t6 * t6     
             
       !pot_temp = 4.0E0_DP * epsilon * (t12 - t6) + delta       
       pot_temp = 4.0E0_DP * epsilon * (t12 - t6) * sw
       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) + pot_temp*0.5
          pot_Col(atom2) = pot_Col(atom2) + 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 + 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
       
       if (do_stress) then

#ifdef IS_MPI
          id1 = tagRow(atom1)
          id2 = tagColumn(atom2)
#else
          id1 = atom1
          id2 = atom2
#endif

          if (molMembershipList(id1) .ne. molMembershipList(id2)) then
 
             ! because the d vector is the rj - ri vector, and 
             ! because fx, fy, fz are the force on atom i, we need a
             ! negative sign here:

             tau_Temp(1) = tau_Temp(1) - d(1) * fx
             tau_Temp(2) = tau_Temp(2) - d(1) * fy
             tau_Temp(3) = tau_Temp(3) - d(1) * fz
             tau_Temp(4) = tau_Temp(4) - d(2) * fx
             tau_Temp(5) = tau_Temp(5) - d(2) * fy
             tau_Temp(6) = tau_Temp(6) - d(2) * fz
             tau_Temp(7) = tau_Temp(7) - d(3) * fx
             tau_Temp(8) = tau_Temp(8) - d(3) * fy
             tau_Temp(9) = tau_Temp(9) - d(3) * fz
             
             virial_Temp = virial_Temp + &
                  (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))

          endif
       endif
       
    endif
    return    
       
  end subroutine do_lj_pair


  !! Calculates the mixing for sigma or epslon

  function calcLJMix(thisParam,param1,param2,status) result(myMixParam)
    character(len=*) :: thisParam
    real(kind = dp)  :: param1
    real(kind = dp)  :: param2
    real(kind = dp ) :: myMixParam

    integer, optional :: status   

    myMixParam = 0.0_dp
    
    if (present(status)) status = 0
    select case (LJ_Mixing_Policy)
    case (1)
       select case (thisParam)
       case ("sigma")
          myMixParam = 0.5_dp * (param1 + param2)
       case ("epsilon")
          myMixParam = sqrt(param1 * param2)
       case default
          status = -1
       end select
    case default
       status = -1
    end select
  end function calcLJMix
  
end module lj
Revision:	1150
Committed:	Fri May 7 21:35:05 2004 UTC (20 years, 2 months ago) by gezelter
File size:	9313 byte(s)
Log Message:	Many changes to get group-based cutoffs to work
#	User	Rev	Content
1	mmeineke	377	!! Calculates Long Range forces Lennard-Jones interactions.
2			!! Corresponds to the force field defined in lj_FF.cpp
3			!! @author Charles F. Vardeman II
4			!! @author Matthew Meineke
5	gezelter	1150	!! @version $Id: calc_LJ_FF.F90,v 1.18 2004-05-07 21:35:04 gezelter Exp $, $Date: 2004-05-07 21:35:04 $, $Name: not supported by cvs2svn $, $Revision: 1.18 $
6	mmeineke	377
7			module lj
8			use definitions
9			use atype_module
10	gezelter	1150	use switcheroo
11	mmeineke	377	use vector_class
12	chuckv	460	use simulation
13	mmeineke	377	#ifdef IS_MPI
14			use mpiSimulation
15			#endif
16			use force_globals
17
18			implicit none
19			PRIVATE
20
21			#define __FORTRAN90
22			#include "fForceField.h"
23
24			integer, save :: LJ_Mixing_Policy
25	mmeineke	619	real(kind=DP), save :: LJ_rcut
26	mmeineke	626	logical, save :: havePolicy = .false., haveCut = .false.
27
28	mmeineke	377
29			!! Logical has lj force field module been initialized?
30
31			logical, save :: LJ_FF_initialized = .false.
32
33			!! Public methods and data
34			public :: init_LJ_FF
35	mmeineke	626	public :: setCutoffLJ
36	mmeineke	377	public :: do_lj_pair
37
38			type :: lj_mixed_params
39			!! Lennard-Jones epsilon
40			real ( kind = dp ) :: epsilon = 0.0_dp
41			!! Lennard-Jones Sigma
42			real ( kind = dp ) :: sigma = 0.0_dp
43			!! Lennard-Jones Sigma to sixth
44			real ( kind = dp ) :: sigma6 = 0.0_dp
45			!!
46			real ( kind = dp ) :: tp6
47			real ( kind = dp ) :: tp12
48
49			real ( kind = dp ) :: delta = 0.0_dp
50
51	mmeineke	626
52	mmeineke	377	end type lj_mixed_params
53
54			type (lj_mixed_params), dimension(:,:), pointer :: ljMixed
55
56			contains
57
58	mmeineke	626	subroutine init_LJ_FF(mix_Policy, status)
59	mmeineke	377	integer, intent(in) :: mix_Policy
60			integer, intent(out) :: status
61			integer :: myStatus
62
63	gezelter	834	if (mix_Policy == LB_MIXING_RULE) then
64			LJ_Mixing_Policy = LB_MIXING_RULE
65	mmeineke	377	else
66	gezelter	834	if (mix_Policy == EXPLICIT_MIXING_RULE) then
67			LJ_Mixing_Policy = EXPLICIT_MIXING_RULE
68	mmeineke	377	else
69			write(,) 'Unknown Mixing Policy!'
70			status = -1
71			return
72			endif
73			endif
74
75	mmeineke	626	havePolicy = .true.
76
77			if (haveCut) then
78			status = 0
79			call createMixingList(myStatus)
80			if (myStatus /= 0) then
81			status = -1
82			return
83			end if
84
85			LJ_FF_initialized = .true.
86	mmeineke	377	end if
87	mmeineke	626
88	mmeineke	377	end subroutine init_LJ_FF
89
90	mmeineke	626	subroutine setCutoffLJ(rcut, status)
91	mmeineke	377	integer :: status, myStatus
92			real(kind=dp) :: rcut
93	gezelter	1150
94			#define __FORTRAN90
95			#include "fSwitchingFunction.h"
96
97	mmeineke	377	status = 0
98
99	chrisfen	1041	LJ_rcut = rcut
100	gezelter	1150	call set_switch(LJ_SWITCH, rcut, rcut)
101
102	chrisfen	1041	!!$ ! ATTENTION! This is a hardwiring of rcut!
103			!!$ LJ_rcut = 9.0d0
104	mmeineke	626	haveCut = .true.
105
106			if (havePolicy) then
107			status = 0
108	mmeineke	619	call createMixingList(myStatus)
109			if (myStatus /= 0) then
110			status = -1
111			return
112			end if
113	mmeineke	626
114			LJ_FF_initialized = .true.
115			end if
116	mmeineke	619
117	mmeineke	377	return
118	mmeineke	626	end subroutine setCutoffLJ
119	mmeineke	377
120			subroutine createMixingList(status)
121			integer :: nAtypes
122			integer :: status
123			integer :: i
124			integer :: j
125			real ( kind = dp ) :: mySigma_i,mySigma_j
126			real ( kind = dp ) :: myEpsilon_i,myEpsilon_j
127			real ( kind = dp ) :: rcut6
128			status = 0
129
130			nAtypes = getSize(atypes)
131			if (nAtypes == 0) then
132			status = -1
133			return
134			end if
135
136			if (.not. associated(ljMixed)) then
137			allocate(ljMixed(nAtypes, nAtypes))
138	chuckv	623	endif
139	mmeineke	377
140			rcut6 = LJ_rcut**6
141
142	chuckv	900	! This loops through all atypes, even those that don't support LJ forces.
143	mmeineke	377	do i = 1, nAtypes
144
145			call getElementProperty(atypes, i, "lj_epsilon", myEpsilon_i)
146			call getElementProperty(atypes, i, "lj_sigma", mySigma_i)
147			! do self mixing rule
148			ljMixed(i,i)%sigma = mySigma_i
149
150	mmeineke	790	ljMixed(i,i)%sigma6 = (ljMixed(i,i)%sigma) ** 6
151	mmeineke	377
152	mmeineke	790	ljMixed(i,i)%tp6 = (ljMixed(i,i)%sigma6)/rcut6
153
154	mmeineke	377	ljMixed(i,i)%tp12 = (ljMixed(i,i)%tp6) ** 2
155
156	mmeineke	790
157	mmeineke	377	ljMixed(i,i)%epsilon = myEpsilon_i
158
159			ljMixed(i,i)%delta = -4.0_DP * ljMixed(i,i)%epsilon * &
160			(ljMixed(i,i)%tp12 - ljMixed(i,i)%tp6)
161
162			do j = i + 1, nAtypes
163			call getElementProperty(atypes,j,"lj_epsilon",myEpsilon_j)
164			call getElementProperty(atypes,j,"lj_sigma", mySigma_j)
165
166			ljMixed(i,j)%sigma = &
167			calcLJMix("sigma",mySigma_i, &
168			mySigma_j)
169
170			ljMixed(i,j)%sigma6 = &
171			(ljMixed(i,j)%sigma)**6
172
173
174			ljMixed(i,j)%tp6 = ljMixed(i,j)%sigma6/rcut6
175
176			ljMixed(i,j)%tp12 = (ljMixed(i,j)%tp6) ** 2
177
178
179			ljMixed(i,j)%epsilon = &
180			calcLJMix("epsilon",myEpsilon_i, &
181			myEpsilon_j)
182
183			ljMixed(i,j)%delta = -4.0_DP * ljMixed(i,j)%epsilon * &
184			(ljMixed(i,j)%tp12 - ljMixed(i,j)%tp6)
185
186
187			ljMixed(j,i)%sigma = ljMixed(i,j)%sigma
188			ljMixed(j,i)%sigma6 = ljMixed(i,j)%sigma6
189			ljMixed(j,i)%tp6 = ljMixed(i,j)%tp6
190			ljMixed(j,i)%tp12 = ljMixed(i,j)%tp12
191			ljMixed(j,i)%epsilon = ljMixed(i,j)%epsilon
192			ljMixed(j,i)%delta = ljMixed(i,j)%delta
193
194			end do
195			end do
196
197			end subroutine createMixingList
198
199	gezelter	1150	subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, pot, f, &
200			do_pot, do_stress)
201	mmeineke	377
202			integer, intent(in) :: atom1, atom2
203			real( kind = dp ), intent(in) :: rij, r2
204	gezelter	1150	real( kind = dp ) :: pot, sw, vpair
205	chuckv	898	real( kind = dp ), dimension(3,nLocal) :: f
206	mmeineke	377	real( kind = dp ), intent(in), dimension(3) :: d
207			logical, intent(in) :: do_pot, do_stress
208
209			! local Variables
210			real( kind = dp ) :: drdx, drdy, drdz
211			real( kind = dp ) :: fx, fy, fz
212			real( kind = dp ) :: pot_temp, dudr
213			real( kind = dp ) :: sigma6
214			real( kind = dp ) :: epsilon
215			real( kind = dp ) :: r6
216			real( kind = dp ) :: t6
217			real( kind = dp ) :: t12
218			real( kind = dp ) :: delta
219	gezelter	490	integer :: id1, id2
220	mmeineke	377
221			if (rij.lt.LJ_rcut) then
222
223			! Look up the correct parameters in the mixing matrix
224			#ifdef IS_MPI
225			sigma6 = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6
226			epsilon = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon
227	gezelter	1150	!delta = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta
228	mmeineke	377	#else
229			sigma6 = ljMixed(atid(atom1),atid(atom2))%sigma6
230			epsilon = ljMixed(atid(atom1),atid(atom2))%epsilon
231	gezelter	1150	!delta = ljMixed(atid(atom1),atid(atom2))%delta
232	mmeineke	377	#endif
233
234			r6 = r2 * r2 * r2
235
236			t6 = sigma6/ r6
237			t12 = t6 * t6
238
239	gezelter	1150	!pot_temp = 4.0E0_DP * epsilon * (t12 - t6) + delta
240			pot_temp = 4.0E0_DP * epsilon * (t12 - t6) * sw
241			vpair = vpair + pot_temp
242	mmeineke	377
243	gezelter	1150	dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
244	mmeineke	377
245	chuckv	482	drdx = d(1) / rij
246			drdy = d(2) / rij
247			drdz = d(3) / rij
248	mmeineke	377
249			fx = dudr * drdx
250			fy = dudr * drdy
251			fz = dudr * drdz
252
253	gezelter	1150
254	mmeineke	377	#ifdef IS_MPI
255			if (do_pot) then
256			pot_Row(atom1) = pot_Row(atom1) + pot_temp*0.5
257			pot_Col(atom2) = pot_Col(atom2) + pot_temp*0.5
258			endif
259
260			f_Row(1,atom1) = f_Row(1,atom1) + fx
261			f_Row(2,atom1) = f_Row(2,atom1) + fy
262			f_Row(3,atom1) = f_Row(3,atom1) + fz
263
264			f_Col(1,atom2) = f_Col(1,atom2) - fx
265			f_Col(2,atom2) = f_Col(2,atom2) - fy
266			f_Col(3,atom2) = f_Col(3,atom2) - fz
267
268			#else
269			if (do_pot) pot = pot + pot_temp
270
271			f(1,atom1) = f(1,atom1) + fx
272			f(2,atom1) = f(2,atom1) + fy
273			f(3,atom1) = f(3,atom1) + fz
274
275			f(1,atom2) = f(1,atom2) - fx
276			f(2,atom2) = f(2,atom2) - fy
277			f(3,atom2) = f(3,atom2) - fz
278			#endif
279
280			if (do_stress) then
281	gezelter	483
282	gezelter	490	#ifdef IS_MPI
283			id1 = tagRow(atom1)
284			id2 = tagColumn(atom2)
285			#else
286			id1 = atom1
287			id2 = atom2
288			#endif
289
290			if (molMembershipList(id1) .ne. molMembershipList(id2)) then
291	gezelter	611
292			! because the d vector is the rj - ri vector, and
293			! because fx, fy, fz are the force on atom i, we need a
294			! negative sign here:
295	mmeineke	491
296	gezelter	611	tau_Temp(1) = tau_Temp(1) - d(1) * fx
297			tau_Temp(2) = tau_Temp(2) - d(1) * fy
298			tau_Temp(3) = tau_Temp(3) - d(1) * fz
299			tau_Temp(4) = tau_Temp(4) - d(2) * fx
300			tau_Temp(5) = tau_Temp(5) - d(2) * fy
301			tau_Temp(6) = tau_Temp(6) - d(2) * fz
302			tau_Temp(7) = tau_Temp(7) - d(3) * fx
303			tau_Temp(8) = tau_Temp(8) - d(3) * fy
304			tau_Temp(9) = tau_Temp(9) - d(3) * fz
305
306	gezelter	483	virial_Temp = virial_Temp + &
307			(tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
308	mmeineke	491
309	gezelter	483	endif
310	mmeineke	377	endif
311
312			endif
313			return
314
315			end subroutine do_lj_pair
316
317
318			!! Calculates the mixing for sigma or epslon
319
320			function calcLJMix(thisParam,param1,param2,status) result(myMixParam)
321			character(len=*) :: thisParam
322			real(kind = dp) :: param1
323			real(kind = dp) :: param2
324			real(kind = dp ) :: myMixParam
325
326			integer, optional :: status
327
328			myMixParam = 0.0_dp
329
330			if (present(status)) status = 0
331			select case (LJ_Mixing_Policy)
332	gezelter	829	case (1)
333	mmeineke	377	select case (thisParam)
334			case ("sigma")
335			myMixParam = 0.5_dp * (param1 + param2)
336			case ("epsilon")
337			myMixParam = sqrt(param1 * param2)
338			case default
339			status = -1
340			end select
341			case default
342			status = -1
343			end select
344			end function calcLJMix
345
346			end module lj