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.19 2004-05-11 21:31:14 gezelter Exp $, $Date: 2004-05-11 21:31:14 $, $Name: not supported by cvs2svn $, $Revision: 1.19 $

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
    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

    ! 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
       
    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
        
    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:	1160
Committed:	Tue May 11 21:31:15 2004 UTC (20 years, 1 month ago) by gezelter
File size:	9190 byte(s)
Log Message:	Fortran-side changes for group-based cutoffs
#	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	1160	!! @version $Id: calc_LJ_FF.F90,v 1.19 2004-05-11 21:31:14 gezelter Exp $, $Date: 2004-05-11 21:31:14 $, $Name: not supported by cvs2svn $, $Revision: 1.19 $
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			status = 0
126
127			nAtypes = getSize(atypes)
128			if (nAtypes == 0) then
129			status = -1
130			return
131			end if
132
133			if (.not. associated(ljMixed)) then
134			allocate(ljMixed(nAtypes, nAtypes))
135	chuckv	623	endif
136	mmeineke	377
137			rcut6 = LJ_rcut**6
138
139	chuckv	900	! This loops through all atypes, even those that don't support LJ forces.
140	mmeineke	377	do i = 1, nAtypes
141
142			call getElementProperty(atypes, i, "lj_epsilon", myEpsilon_i)
143			call getElementProperty(atypes, i, "lj_sigma", mySigma_i)
144			! do self mixing rule
145			ljMixed(i,i)%sigma = mySigma_i
146
147	mmeineke	790	ljMixed(i,i)%sigma6 = (ljMixed(i,i)%sigma) ** 6
148	mmeineke	377
149	mmeineke	790	ljMixed(i,i)%tp6 = (ljMixed(i,i)%sigma6)/rcut6
150
151	mmeineke	377	ljMixed(i,i)%tp12 = (ljMixed(i,i)%tp6) ** 2
152
153	mmeineke	790
154	mmeineke	377	ljMixed(i,i)%epsilon = myEpsilon_i
155
156			ljMixed(i,i)%delta = -4.0_DP * ljMixed(i,i)%epsilon * &
157			(ljMixed(i,i)%tp12 - ljMixed(i,i)%tp6)
158
159			do j = i + 1, nAtypes
160			call getElementProperty(atypes,j,"lj_epsilon",myEpsilon_j)
161			call getElementProperty(atypes,j,"lj_sigma", mySigma_j)
162
163			ljMixed(i,j)%sigma = &
164			calcLJMix("sigma",mySigma_i, &
165			mySigma_j)
166
167			ljMixed(i,j)%sigma6 = &
168			(ljMixed(i,j)%sigma)**6
169
170
171			ljMixed(i,j)%tp6 = ljMixed(i,j)%sigma6/rcut6
172
173			ljMixed(i,j)%tp12 = (ljMixed(i,j)%tp6) ** 2
174
175
176			ljMixed(i,j)%epsilon = &
177			calcLJMix("epsilon",myEpsilon_i, &
178			myEpsilon_j)
179
180			ljMixed(i,j)%delta = -4.0_DP * ljMixed(i,j)%epsilon * &
181			(ljMixed(i,j)%tp12 - ljMixed(i,j)%tp6)
182
183
184			ljMixed(j,i)%sigma = ljMixed(i,j)%sigma
185			ljMixed(j,i)%sigma6 = ljMixed(i,j)%sigma6
186			ljMixed(j,i)%tp6 = ljMixed(i,j)%tp6
187			ljMixed(j,i)%tp12 = ljMixed(i,j)%tp12
188			ljMixed(j,i)%epsilon = ljMixed(i,j)%epsilon
189			ljMixed(j,i)%delta = ljMixed(i,j)%delta
190
191			end do
192			end do
193
194			end subroutine createMixingList
195
196	gezelter	1150	subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, pot, f, &
197			do_pot, do_stress)
198	mmeineke	377
199			integer, intent(in) :: atom1, atom2
200			real( kind = dp ), intent(in) :: rij, r2
201	gezelter	1150	real( kind = dp ) :: pot, sw, vpair
202	chuckv	898	real( kind = dp ), dimension(3,nLocal) :: f
203	mmeineke	377	real( kind = dp ), intent(in), dimension(3) :: d
204			logical, intent(in) :: do_pot, do_stress
205
206			! local Variables
207			real( kind = dp ) :: drdx, drdy, drdz
208			real( kind = dp ) :: fx, fy, fz
209			real( kind = dp ) :: pot_temp, dudr
210			real( kind = dp ) :: sigma6
211			real( kind = dp ) :: epsilon
212			real( kind = dp ) :: r6
213			real( kind = dp ) :: t6
214			real( kind = dp ) :: t12
215			real( kind = dp ) :: delta
216	gezelter	490	integer :: id1, id2
217	mmeineke	377
218	gezelter	1160	! Look up the correct parameters in the mixing matrix
219	mmeineke	377	#ifdef IS_MPI
220	gezelter	1160	sigma6 = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6
221			epsilon = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon
222			delta = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta
223	mmeineke	377	#else
224	gezelter	1160	sigma6 = ljMixed(atid(atom1),atid(atom2))%sigma6
225			epsilon = ljMixed(atid(atom1),atid(atom2))%epsilon
226			delta = ljMixed(atid(atom1),atid(atom2))%delta
227	mmeineke	377	#endif
228	gezelter	1160
229			r6 = r2 * r2 * r2
230
231			t6 = sigma6/ r6
232			t12 = t6 * t6
233
234			pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
235			if (LJ_do_shift) then
236			pot_temp = pot_temp + delta
237			endif
238
239			vpair = vpair + pot_temp
240	mmeineke	377
241	gezelter	1160	dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
242	mmeineke	377
243	gezelter	1160	drdx = d(1) / rij
244			drdy = d(2) / rij
245			drdz = d(3) / rij
246	mmeineke	377
247	gezelter	1160	fx = dudr * drdx
248			fy = dudr * drdy
249			fz = dudr * drdz
250
251	mmeineke	377
252			#ifdef IS_MPI
253	gezelter	1160	if (do_pot) then
254			pot_Row(atom1) = pot_Row(atom1) + swpot_temp0.5
255			pot_Col(atom2) = pot_Col(atom2) + swpot_temp0.5
256			endif
257
258			f_Row(1,atom1) = f_Row(1,atom1) + fx
259			f_Row(2,atom1) = f_Row(2,atom1) + fy
260			f_Row(3,atom1) = f_Row(3,atom1) + fz
261
262			f_Col(1,atom2) = f_Col(1,atom2) - fx
263			f_Col(2,atom2) = f_Col(2,atom2) - fy
264			f_Col(3,atom2) = f_Col(3,atom2) - fz
265
266	mmeineke	377	#else
267	gezelter	1160	if (do_pot) pot = pot + sw*pot_temp
268	mmeineke	377
269	gezelter	1160	f(1,atom1) = f(1,atom1) + fx
270			f(2,atom1) = f(2,atom1) + fy
271			f(3,atom1) = f(3,atom1) + fz
272
273			f(1,atom2) = f(1,atom2) - fx
274			f(2,atom2) = f(2,atom2) - fy
275			f(3,atom2) = f(3,atom2) - fz
276	mmeineke	377	#endif
277
278	gezelter	1160	if (do_stress) then
279	gezelter	483
280	gezelter	490	#ifdef IS_MPI
281	gezelter	1160	id1 = tagRow(atom1)
282			id2 = tagColumn(atom2)
283	gezelter	490	#else
284	gezelter	1160	id1 = atom1
285			id2 = atom2
286	gezelter	490	#endif
287	gezelter	1160
288			if (molMembershipList(id1) .ne. molMembershipList(id2)) then
289
290			! because the d vector is the rj - ri vector, and
291			! because fx, fy, fz are the force on atom i, we need a
292			! negative sign here:
293
294			tau_Temp(1) = tau_Temp(1) - d(1) * fx
295			tau_Temp(2) = tau_Temp(2) - d(1) * fy
296			tau_Temp(3) = tau_Temp(3) - d(1) * fz
297			tau_Temp(4) = tau_Temp(4) - d(2) * fx
298			tau_Temp(5) = tau_Temp(5) - d(2) * fy
299			tau_Temp(6) = tau_Temp(6) - d(2) * fz
300			tau_Temp(7) = tau_Temp(7) - d(3) * fx
301			tau_Temp(8) = tau_Temp(8) - d(3) * fy
302			tau_Temp(9) = tau_Temp(9) - d(3) * fz
303
304			virial_Temp = virial_Temp + &
305			(tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
306
307	mmeineke	377	endif
308			endif
309	gezelter	1160
310	mmeineke	377	return
311	gezelter	1160
312	mmeineke	377	end subroutine do_lj_pair
313	gezelter	1160
314
315	mmeineke	377	!! Calculates the mixing for sigma or epslon
316	gezelter	1160
317	mmeineke	377	function calcLJMix(thisParam,param1,param2,status) result(myMixParam)
318			character(len=*) :: thisParam
319			real(kind = dp) :: param1
320			real(kind = dp) :: param2
321			real(kind = dp ) :: myMixParam
322
323			integer, optional :: status
324
325			myMixParam = 0.0_dp
326
327			if (present(status)) status = 0
328			select case (LJ_Mixing_Policy)
329	gezelter	829	case (1)
330	mmeineke	377	select case (thisParam)
331			case ("sigma")
332			myMixParam = 0.5_dp * (param1 + param2)
333			case ("epsilon")
334			myMixParam = sqrt(param1 * param2)
335			case default
336			status = -1
337			end select
338			case default
339			status = -1
340			end select
341			end function calcLJMix
342
343			end module lj