mdtools/libmdCode/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.12 2003-03-18 16:46:47 gezelter Exp $, $Date: 2003-03-18 16:46:47 $, $Name: not supported by cvs2svn $, $Revision: 1.12 $

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

  implicit none
  PRIVATE

#define __FORTRAN90
#include "fForceField.h"

  integer, save :: LJ_Mixing_Policy
  integer, save :: LJ_rcut
  
  !! Logical has lj force field module been initialized?
  
  logical, save :: LJ_FF_initialized = .false.
  
  !! Public methods and data
  public :: init_LJ_FF
  public :: LJ_new_rcut
  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, rcut, status)
    integer, intent(in) :: mix_Policy
    integer, intent(out) :: status
    integer :: myStatus
    real(kind=dp) :: rcut
    
    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
    
    LJ_rcut = rcut

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

    LJ_rcut = rcut
    status = 0
    call createMixingList(myStatus)
    if (myStatus /= 0) then
       status = -1
       return
    end if

    return
  end subroutine LJ_new_rcut
  
  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))
    else
       status = -1
       return
    end if

    rcut6 = LJ_rcut**6

    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, pot, f, do_pot, do_stress)

    integer, intent(in) ::  atom1, atom2
    real( kind = dp ), intent(in) :: rij, r2
    real( kind = dp ) :: pot
    real( kind = dp ), dimension(3, getNlocal()) :: 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


    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       
       
       dudr = 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
          tau_Temp(1) = tau_Temp(1) + fx * d(1)
          tau_Temp(2) = tau_Temp(2) + fx * d(2)
          tau_Temp(3) = tau_Temp(3) + fx * d(3)
          tau_Temp(4) = tau_Temp(4) + fy * d(1)
          tau_Temp(5) = tau_Temp(5) + fy * d(2)
          tau_Temp(6) = tau_Temp(6) + fy * d(3)
          tau_Temp(7) = tau_Temp(7) + fz * d(1)
          tau_Temp(8) = tau_Temp(8) + fz * d(2)
          tau_Temp(9) = tau_Temp(9) + fz * d(3)
          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 (LB_MIXING_RULE)
       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:	360
Committed:	Tue Mar 18 16:46:47 2003 UTC (21 years, 3 months ago) by gezelter
File size:	8344 byte(s)
Log Message:	brought force_globals back from the dead to fix circular references
#	User	Rev	Content
1	chuckv	298	!! 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	360	!! @version $Id: calc_LJ_FF.F90,v 1.12 2003-03-18 16:46:47 gezelter Exp $, $Date: 2003-03-18 16:46:47 $, $Name: not supported by cvs2svn $, $Revision: 1.12 $
6	chuckv	298
7	gezelter	309	module lj
8	chuckv	298	use definitions
9	gezelter	328	use atype_module
10	gezelter	317	use vector_class
11	chuckv	298	#ifdef IS_MPI
12			use mpiSimulation
13			#endif
14	gezelter	360	use force_globals
15
16	chuckv	298	implicit none
17			PRIVATE
18	gezelter	326
19	gezelter	358	#define __FORTRAN90
20			#include "fForceField.h"
21
22	gezelter	326	integer, save :: LJ_Mixing_Policy
23	gezelter	360	integer, save :: LJ_rcut
24	gezelter	309
25			!! Logical has lj force field module been initialized?
26
27	gezelter	326	logical, save :: LJ_FF_initialized = .false.
28
29	gezelter	309	!! Public methods and data
30	gezelter	326	public :: init_LJ_FF
31	gezelter	360	public :: LJ_new_rcut
32	gezelter	309	public :: do_lj_pair
33
34	chuckv	306	type :: lj_mixed_params
35	gezelter	309	!! Lennard-Jones epsilon
36	chuckv	298	real ( kind = dp ) :: epsilon = 0.0_dp
37	gezelter	309	!! Lennard-Jones Sigma
38	chuckv	298	real ( kind = dp ) :: sigma = 0.0_dp
39	gezelter	309	!! Lennard-Jones Sigma to sixth
40	chuckv	298	real ( kind = dp ) :: sigma6 = 0.0_dp
41	chuckv	320	!!
42			real ( kind = dp ) :: tp6
43			real ( kind = dp ) :: tp12
44
45			real ( kind = dp ) :: delta = 0.0_dp
46
47	chuckv	306	end type lj_mixed_params
48	chuckv	298
49	gezelter	309	type (lj_mixed_params), dimension(:,:), pointer :: ljMixed
50
51	chuckv	298	contains
52
53	gezelter	360	subroutine init_LJ_FF(mix_Policy, rcut, status)
54	gezelter	358	integer, intent(in) :: mix_Policy
55	chuckv	298	integer, intent(out) :: status
56			integer :: myStatus
57	gezelter	360	real(kind=dp) :: rcut
58	gezelter	309
59	gezelter	358	if (mix_Policy == LB_MIXING_RULE) then
60			LJ_Mixing_Policy = LB_MIXING_RULE
61	gezelter	326	else
62	gezelter	358	if (mix_Policy == EXPLICIT_MIXING_RULE) then
63			LJ_Mixing_Policy = EXPLICIT_MIXING_RULE
64	gezelter	326	else
65			write(,) 'Unknown Mixing Policy!'
66			status = -1
67			return
68			endif
69			endif
70	gezelter	356
71	gezelter	360	LJ_rcut = rcut
72
73	chuckv	298	status = 0
74	chuckv	320	call createMixingList(myStatus)
75	chuckv	298	if (myStatus /= 0) then
76			status = -1
77			return
78			end if
79	gezelter	309
80	gezelter	326	LJ_FF_initialized = .true.
81	gezelter	358
82	gezelter	326	end subroutine init_LJ_FF
83	gezelter	309
84	gezelter	360	subroutine LJ_new_rcut(rcut, status)
85			integer :: status, myStatus
86			real(kind=dp) :: rcut
87
88			LJ_rcut = rcut
89			status = 0
90			call createMixingList(myStatus)
91			if (myStatus /= 0) then
92			status = -1
93			return
94			end if
95
96			return
97			end subroutine LJ_new_rcut
98
99	chuckv	320	subroutine createMixingList(status)
100			integer :: nAtypes
101	chuckv	298	integer :: status
102			integer :: i
103			integer :: j
104	chuckv	320	real ( kind = dp ) :: mySigma_i,mySigma_j
105			real ( kind = dp ) :: myEpsilon_i,myEpsilon_j
106	gezelter	360	real ( kind = dp ) :: rcut6
107	chuckv	298	status = 0
108	gezelter	309
109	gezelter	323	nAtypes = getSize(atypes)
110			if (nAtypes == 0) then
111	chuckv	298	status = -1
112			return
113			end if
114	gezelter	323
115	chuckv	298	if (.not. associated(ljMixed)) then
116	gezelter	323	allocate(ljMixed(nAtypes, nAtypes))
117	chuckv	298	else
118			status = -1
119			return
120			end if
121	gezelter	360
122			rcut6 = LJ_rcut**6
123
124	chuckv	320	do i = 1, nAtypes
125
126	gezelter	360	call getElementProperty(atypes, i, "lj_epsilon", myEpsilon_i)
127			call getElementProperty(atypes, i, "lj_sigma", mySigma_i)
128	gezelter	309	! do self mixing rule
129	chuckv	320	ljMixed(i,i)%sigma = mySigma_i
130	gezelter	309
131	chuckv	306	ljMixed(i,i)%sigma6 = (ljMixed(i,i)%sigma) ** 6
132	chuckv	320	ljMixed(i,i)%tp6 = ljMixed(i,i)%sigma6/rcut6
133
134			ljMixed(i,i)%tp12 = (ljMixed(i,i)%tp6) ** 2
135
136			ljMixed(i,i)%epsilon = myEpsilon_i
137
138			ljMixed(i,i)%delta = -4.0_DP * ljMixed(i,i)%epsilon * &
139			(ljMixed(i,i)%tp12 - ljMixed(i,i)%tp6)
140
141			do j = i + 1, nAtypes
142			call getElementProperty(atypes,j,"lj_epsilon",myEpsilon_j)
143			call getElementProperty(atypes,j,"lj_sigma", mySigma_j)
144	chuckv	298
145	chuckv	306	ljMixed(i,j)%sigma = &
146	chuckv	320	calcLJMix("sigma",mySigma_i, &
147			mySigma_j)
148	chuckv	298
149	chuckv	306	ljMixed(i,j)%sigma6 = &
150			(ljMixed(i,j)%sigma)**6
151	chuckv	298
152	chuckv	320
153			ljMixed(i,j)%tp6 = ljMixed(i,j)%sigma6/rcut6
154
155			ljMixed(i,j)%tp12 = (ljMixed(i,j)%tp6) ** 2
156
157
158	chuckv	306	ljMixed(i,j)%epsilon = &
159	chuckv	320	calcLJMix("epsilon",myEpsilon_i, &
160			myEpsilon_j)
161
162			ljMixed(i,j)%delta = -4.0_DP * ljMixed(i,j)%epsilon * &
163			(ljMixed(i,j)%tp12 - ljMixed(i,j)%tp6)
164
165
166	chuckv	306	ljMixed(j,i)%sigma = ljMixed(i,j)%sigma
167			ljMixed(j,i)%sigma6 = ljMixed(i,j)%sigma6
168	chuckv	320	ljMixed(j,i)%tp6 = ljMixed(i,j)%tp6
169			ljMixed(j,i)%tp12 = ljMixed(i,j)%tp12
170	chuckv	306	ljMixed(j,i)%epsilon = ljMixed(i,j)%epsilon
171	chuckv	320	ljMixed(j,i)%delta = ljMixed(i,j)%delta
172	gezelter	309
173	chuckv	298	end do
174			end do
175	gezelter	309
176	chuckv	298	end subroutine createMixingList
177	gezelter	309
178	gezelter	326	subroutine do_lj_pair(atom1, atom2, d, rij, r2, pot, f, do_pot, do_stress)
179	chuckv	298
180	gezelter	309	integer, intent(in) :: atom1, atom2
181	gezelter	323	real( kind = dp ), intent(in) :: rij, r2
182	gezelter	309	real( kind = dp ) :: pot
183			real( kind = dp ), dimension(3, getNlocal()) :: f
184	gezelter	317	real( kind = dp ), intent(in), dimension(3) :: d
185	gezelter	326	logical, intent(in) :: do_pot, do_stress
186	chuckv	298
187	gezelter	309	! local Variables
188			real( kind = dp ) :: drdx, drdy, drdz
189			real( kind = dp ) :: fx, fy, fz
190			real( kind = dp ) :: pot_temp, dudr
191	chuckv	298	real( kind = dp ) :: sigma6
192			real( kind = dp ) :: epsilon
193			real( kind = dp ) :: r6
194			real( kind = dp ) :: t6
195			real( kind = dp ) :: t12
196			real( kind = dp ) :: delta
197
198	gezelter	360
199			if (rij.lt.LJ_rcut) then
200
201			! Look up the correct parameters in the mixing matrix
202	chuckv	320	#ifdef IS_MPI
203	gezelter	360	sigma6 = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6
204			epsilon = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon
205			delta = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta
206	chuckv	320	#else
207	gezelter	360	sigma6 = ljMixed(atid(atom1),atid(atom2))%sigma6
208			epsilon = ljMixed(atid(atom1),atid(atom2))%epsilon
209			delta = ljMixed(atid(atom1),atid(atom2))%delta
210	chuckv	320	#endif
211	gezelter	309
212	gezelter	360	r6 = r2 * r2 * r2
213
214			t6 = sigma6/ r6
215			t12 = t6 * t6
216
217	gezelter	309	pot_temp = 4.0E0_DP * epsilon * (t12 - t6) + delta
218	gezelter	360
219	gezelter	309	dudr = 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
220	gezelter	360
221	gezelter	317	drdx = -d(1) / rij
222			drdy = -d(2) / rij
223			drdz = -d(3) / rij
224	gezelter	309
225			fx = dudr * drdx
226			fy = dudr * drdy
227			fz = dudr * drdz
228	gezelter	360
229	gezelter	309	#ifdef IS_MPI
230	gezelter	326	if (do_pot) then
231			pot_Row(atom1) = pot_Row(atom1) + pot_temp*0.5
232			pot_Col(atom2) = pot_Col(atom2) + pot_temp*0.5
233			endif
234	chuckv	298
235	gezelter	309	f_Row(1,atom1) = f_Row(1,atom1) + fx
236			f_Row(2,atom1) = f_Row(2,atom1) + fy
237			f_Row(3,atom1) = f_Row(3,atom1) + fz
238	chuckv	298
239	gezelter	309	f_Col(1,atom2) = f_Col(1,atom2) - fx
240			f_Col(2,atom2) = f_Col(2,atom2) - fy
241			f_Col(3,atom2) = f_Col(3,atom2) - fz
242	chuckv	298
243	gezelter	309	#else
244	gezelter	326	if (do_pot) pot = pot + pot_temp
245	gezelter	309
246			f(1,atom1) = f(1,atom1) + fx
247			f(2,atom1) = f(2,atom1) + fy
248			f(3,atom1) = f(3,atom1) + fz
249
250			f(1,atom2) = f(1,atom2) - fx
251			f(2,atom2) = f(2,atom2) - fy
252			f(3,atom2) = f(3,atom2) - fz
253			#endif
254
255	gezelter	326	if (do_stress) then
256	gezelter	317	tau_Temp(1) = tau_Temp(1) + fx * d(1)
257			tau_Temp(2) = tau_Temp(2) + fx * d(2)
258			tau_Temp(3) = tau_Temp(3) + fx * d(3)
259			tau_Temp(4) = tau_Temp(4) + fy * d(1)
260			tau_Temp(5) = tau_Temp(5) + fy * d(2)
261			tau_Temp(6) = tau_Temp(6) + fy * d(3)
262			tau_Temp(7) = tau_Temp(7) + fz * d(1)
263			tau_Temp(8) = tau_Temp(8) + fz * d(2)
264			tau_Temp(9) = tau_Temp(9) + fz * d(3)
265	gezelter	309	virial_Temp = virial_Temp + (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
266			endif
267	gezelter	323
268	gezelter	309	endif
269			return
270			end subroutine do_lj_pair
271
272
273	gezelter	326	!! Calculates the mixing for sigma or epslon
274
275	chuckv	298	function calcLJMix(thisParam,param1,param2,status) result(myMixParam)
276			character(len=*) :: thisParam
277			real(kind = dp) :: param1
278			real(kind = dp) :: param2
279			real(kind = dp ) :: myMixParam
280	gezelter	326
281			integer, optional :: status
282
283	chuckv	298	myMixParam = 0.0_dp
284	gezelter	309
285	chuckv	298	if (present(status)) status = 0
286	gezelter	326	select case (LJ_Mixing_Policy)
287	gezelter	358	case (LB_MIXING_RULE)
288	gezelter	309	select case (thisParam)
289			case ("sigma")
290			myMixParam = 0.5_dp * (param1 + param2)
291			case ("epsilon")
292			myMixParam = sqrt(param1 * param2)
293			case default
294			status = -1
295			end select
296	chuckv	298	case default
297			status = -1
298			end select
299			end function calcLJMix
300	gezelter	309
301			end module lj