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.22 2004-06-01 21:45:22 gezelter Exp $, $Date: 2004-06-01 21:45:22 $, $Name: not supported by cvs2svn $, $Revision: 1.22 $

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.
  logical, save :: haveCut = .false.
  logical, save :: LJ_do_shift = .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, do_shift, status)
    logical, intent(in):: do_shift
    integer :: status, myStatus
    real(kind=dp) :: rcut

#define __FORTRAN90
#include "fSwitchingFunction.h"

    status = 0

    LJ_rcut = rcut
    LJ_do_shift = do_shift
    call set_switch(LJ_SWITCH, rcut, rcut)
    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
    logical :: I_isLJ, J_isLJ
    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, "is_LJ", I_isLJ)

       if (I_isLJ) then

          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, "is_LJ", J_isLJ)
             
             if (J_isLJ) then                
             
                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
             endif                
          end do
       endif
    end do
    
  end subroutine createMixingList
  
  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

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