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.4 2003-03-11 23:13:06 gezelter Exp $, $Date: 2003-03-11 23:13:06 $, $Name: not supported by cvs2svn $, $Revision: 1.4 $


module lj
  use simulation
  use definitions
  use forceGlobals
  use atype_typedefs
  use vector_class
#ifdef IS_MPI
  use mpiSimulation
#endif
  implicit none
  PRIVATE
  
  !! Logical has lj force field module been initialized?
  logical, save :: isljFFinit = .false.
  
  
  !! Public methods and data
  public :: init_ljFF
  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 Squared
     real ( kind = dp )  :: sigma2 = 0.0_dp
     !! Lennard-Jones Sigma to sixth
     real ( kind = dp )  :: sigma6 = 0.0_dp
  end type lj_mixed_params
  
  type (lj_mixed_params), dimension(:,:), pointer :: ljMixed
  
contains

  subroutine init_ljFF(ListHead,status)
    type (atype),pointer :: ListHead
    integer :: nTypes
    integer, intent(out) :: status
    
    integer :: myStatus
    
    status = 0
    call createMixingList(ListHead,myStatus)
    if (myStatus /= 0) then
       status = -1
       return
    end if
    
  end subroutine init_ljFF
  
  subroutine createMixingList(ListHead,status)
    type (atype), pointer :: ListHead
    integer :: listSize
    integer :: status
    integer :: i
    integer :: j
    type (atype), pointer :: tmpPtr_i => null()
    type (atype), pointer :: tmpPtr_j => null()
    status = 0
    
    listSize = get_this_ListLen(ListHead)
    if (listSize == 0) then
       status = -1
       return
    end if
    
    
    if (.not. associated(ljMixed)) then
       allocate(ljMixed(listSize,listSize))
    else
       status = -1
       return
    end if
    
    
    tmpPtr_i => ListHead
    tmpPtr_j => tmpPtr_i%next
    do while (associated(tmpPtr_i))
       i = i + 1
       ! do self mixing rule
       ljMixed(i,i)%sigma   = tmpPtr_i%lj_sigma
       
       ljMixed(i,i)%sigma2  = (ljMixed(i,i)%sigma) ** 2 
       
       ljMixed(i,i)%sigma6  = (ljMixed(i,i)%sigma) ** 6                                                                                                  
       ljMixed(i,i)%epsilon = tmpPtr_i%lj_epsilon
       
       j = i + 1
       do while (associated(tmpPtr_j))
          
          ljMixed(i,j)%sigma  =  &
               calcLJMix("sigma",tmpPtr_i%lj_sigma, &
               tmpPtr_j%lj_sigma)
          
          ljMixed(i,j)%sigma2  = &
               (ljMixed(i,j)%sigma)**2
          
          ljMixed(i,j)%sigma6 = &
               (ljMixed(i,j)%sigma)**6
          
          ljMixed(i,j)%epsilon = &
               calcLJMix("epsilon",tmpPtr_i%lj_epsilon, &
               tmpPtr_j%lj_epsilon)
          ljMixed(j,i)%sigma   = ljMixed(i,j)%sigma
          ljMixed(j,i)%sigma2  = ljMixed(i,j)%sigma2
          ljMixed(j,i)%sigma6  = ljMixed(i,j)%sigma6
          ljMixed(j,i)%epsilon = ljMixed(i,j)%epsilon
          
          
          tmpPtr_j => tmpPtr_j%next
          j = j + 1
       end do
       ! advance pointers 
       tmpPtr_i => tmpPtr_i%next
       if (associated(tmpPtr_i)) then
          tmpPtr_j => tmpPtr_i%next
       endif
       
    end do
    
  end subroutine createMixingList
  

  subroutine do_lj_pair(atom1, atom2, d, rij, pot, f)

    integer, intent(in) ::  atom1, atom2
    real( kind = dp ), intent(in) :: rij
    real( kind = dp ) :: pot
    real( kind = dp ), dimension(3, getNlocal()) :: f    
    real( kind = dp ), intent(in), dimension(3) :: d

    ! local Variables
    real( kind = dp ) :: drdx, drdy, drdz
    real( kind = dp ) :: fx, fy, fz
    real( kind = dp ) :: pot_temp, dudr
    real( kind = dp ) :: sigma
    real( kind = dp ) :: sigma2
    real( kind = dp ) :: sigma6
    real( kind = dp ) :: epsilon

    real( kind = dp ) :: rcut
    real( kind = dp ) :: rcut2
    real( kind = dp ) :: rcut6
    real( kind = dp ) :: r2
    real( kind = dp ) :: r6

    real( kind = dp ) :: t6
    real( kind = dp ) :: t12
    real( kind = dp ) :: tp6
    real( kind = dp ) :: tp12
    real( kind = dp ) :: delta

    ! Look up the correct parameters in the mixing matrix
    sigma    = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma
    sigma2   = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma2
    sigma6   = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma6
    epsilon  = ljMixed(atype1%atype_ident,atype2%atype_ident)%epsilon
   
    call getRcut(rcut,rcut2=rcut2,rcut6=rcut6)
    
    r2 = rij * rij
    r6 = r2 * r2 * r2

    t6  = sigma6/ r6
    t12 = t6 * t6
                                                                               
    tp6 = sigma6 / rcut6
    tp12 = tp6*tp6
                                                                               
    delta = -4.0_DP * epsilon * (tp12 - tp6)
    
    if (rij.le.rcut) then
       
       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
       pot_Row(atom1) = pot_Row(atom1) + pot_temp*0.5
       pot_Col(atom2) = pot_Col(atom2) + pot_temp*0.5

       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
       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 (doStress()) 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 based on character
  !! string for initialzition of mixing array.
  function calcLJMix(thisParam,param1,param2,status) result(myMixParam)
    character(len=*) :: thisParam
    real(kind = dp)  :: param1
    real(kind = dp)  :: param2
    real(kind = dp ) :: myMixParam
    character(len = getStringLen()) :: thisMixingRule
    integer, optional :: status
    
    !! get the mixing rules from the simulation
    thisMixingRule = returnMixingRules()
    myMixParam = 0.0_dp
    
    if (present(status)) status = 0
    select case (thisMixingRule)
    case ("standard")
       select case (thisParam)
       case ("sigma")
          myMixParam = 0.5_dp * (param1 + param2)
       case ("epsilon")
          myMixParam = sqrt(param1 * param2)
       case default
          status = -1
       end select
    case("LJglass")
    case default
       status = -1
    end select
  end function calcLJMix
  
end module lj
Revision:	317
Committed:	Tue Mar 11 23:13:06 2003 UTC (21 years, 3 months ago) by gezelter
File size:	7639 byte(s)
Log Message:	Bug fixes
#	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	317	!! @version $Id: calc_LJ_FF.F90,v 1.4 2003-03-11 23:13:06 gezelter Exp $, $Date: 2003-03-11 23:13:06 $, $Name: not supported by cvs2svn $, $Revision: 1.4 $
6	chuckv	298
7
8
9	gezelter	309	module lj
10	chuckv	298	use simulation
11			use definitions
12	gezelter	309	use forceGlobals
13			use atype_typedefs
14	gezelter	317	use vector_class
15	chuckv	298	#ifdef IS_MPI
16			use mpiSimulation
17			#endif
18			implicit none
19			PRIVATE
20	gezelter	309
21			!! Logical has lj force field module been initialized?
22	chuckv	298	logical, save :: isljFFinit = .false.
23	gezelter	309
24
25			!! Public methods and data
26	chuckv	298	public :: init_ljFF
27	gezelter	309	public :: do_lj_pair
28
29	chuckv	306	type :: lj_mixed_params
30	gezelter	309	!! Lennard-Jones epsilon
31	chuckv	298	real ( kind = dp ) :: epsilon = 0.0_dp
32	gezelter	309	!! Lennard-Jones Sigma
33	chuckv	298	real ( kind = dp ) :: sigma = 0.0_dp
34	gezelter	309	!! Lennard-Jones Sigma Squared
35	chuckv	298	real ( kind = dp ) :: sigma2 = 0.0_dp
36	gezelter	309	!! Lennard-Jones Sigma to sixth
37	chuckv	298	real ( kind = dp ) :: sigma6 = 0.0_dp
38	chuckv	306	end type lj_mixed_params
39	chuckv	298
40	gezelter	309	type (lj_mixed_params), dimension(:,:), pointer :: ljMixed
41
42	chuckv	298	contains
43
44			subroutine init_ljFF(ListHead,status)
45			type (atype),pointer :: ListHead
46			integer :: nTypes
47			integer, intent(out) :: status
48	gezelter	309
49	chuckv	298	integer :: myStatus
50	gezelter	309
51	chuckv	298	status = 0
52	chuckv	306	call createMixingList(ListHead,myStatus)
53	chuckv	298	if (myStatus /= 0) then
54			status = -1
55			return
56			end if
57	gezelter	309
58	chuckv	298	end subroutine init_ljFF
59	gezelter	309
60	chuckv	298	subroutine createMixingList(ListHead,status)
61			type (atype), pointer :: ListHead
62			integer :: listSize
63			integer :: status
64			integer :: i
65			integer :: j
66	chuckv	306	type (atype), pointer :: tmpPtr_i => null()
67			type (atype), pointer :: tmpPtr_j => null()
68	chuckv	298	status = 0
69	gezelter	309
70			listSize = get_this_ListLen(ListHead)
71	chuckv	298	if (listSize == 0) then
72			status = -1
73			return
74			end if
75	gezelter	309
76
77	chuckv	298	if (.not. associated(ljMixed)) then
78			allocate(ljMixed(listSize,listSize))
79			else
80			status = -1
81			return
82			end if
83
84	gezelter	309
85
86	chuckv	306	tmpPtr_i => ListHead
87			tmpPtr_j => tmpPtr_i%next
88			do while (associated(tmpPtr_i))
89			i = i + 1
90	gezelter	309	! do self mixing rule
91			ljMixed(i,i)%sigma = tmpPtr_i%lj_sigma
92
93	chuckv	306	ljMixed(i,i)%sigma2 = (ljMixed(i,i)%sigma) ** 2
94	gezelter	309
95	chuckv	306	ljMixed(i,i)%sigma6 = (ljMixed(i,i)%sigma) ** 6
96	gezelter	309	ljMixed(i,i)%epsilon = tmpPtr_i%lj_epsilon
97
98	chuckv	306	j = i + 1
99			do while (associated(tmpPtr_j))
100	chuckv	298
101	chuckv	306	ljMixed(i,j)%sigma = &
102	gezelter	309	calcLJMix("sigma",tmpPtr_i%lj_sigma, &
103			tmpPtr_j%lj_sigma)
104	chuckv	298
105	chuckv	306	ljMixed(i,j)%sigma2 = &
106			(ljMixed(i,j)%sigma)**2
107	chuckv	298
108	chuckv	306	ljMixed(i,j)%sigma6 = &
109			(ljMixed(i,j)%sigma)**6
110	chuckv	298
111	chuckv	306	ljMixed(i,j)%epsilon = &
112	gezelter	309	calcLJMix("epsilon",tmpPtr_i%lj_epsilon, &
113			tmpPtr_j%lj_epsilon)
114	chuckv	306	ljMixed(j,i)%sigma = ljMixed(i,j)%sigma
115			ljMixed(j,i)%sigma2 = ljMixed(i,j)%sigma2
116			ljMixed(j,i)%sigma6 = ljMixed(i,j)%sigma6
117			ljMixed(j,i)%epsilon = ljMixed(i,j)%epsilon
118	gezelter	309
119
120	chuckv	306	tmpPtr_j => tmpPtr_j%next
121			j = j + 1
122	chuckv	298	end do
123	gezelter	309	! advance pointers
124	chuckv	306	tmpPtr_i => tmpPtr_i%next
125			if (associated(tmpPtr_i)) then
126			tmpPtr_j => tmpPtr_i%next
127	chuckv	298	endif
128
129			end do
130	gezelter	309
131	chuckv	298	end subroutine createMixingList
132	gezelter	309
133	chuckv	298
134	gezelter	317	subroutine do_lj_pair(atom1, atom2, d, rij, pot, f)
135	chuckv	298
136	gezelter	309	integer, intent(in) :: atom1, atom2
137	gezelter	317	real( kind = dp ), intent(in) :: rij
138	gezelter	309	real( kind = dp ) :: pot
139			real( kind = dp ), dimension(3, getNlocal()) :: f
140	gezelter	317	real( kind = dp ), intent(in), dimension(3) :: d
141	chuckv	298
142	gezelter	309	! local Variables
143			real( kind = dp ) :: drdx, drdy, drdz
144			real( kind = dp ) :: fx, fy, fz
145			real( kind = dp ) :: pot_temp, dudr
146	chuckv	298	real( kind = dp ) :: sigma
147			real( kind = dp ) :: sigma2
148			real( kind = dp ) :: sigma6
149			real( kind = dp ) :: epsilon
150
151			real( kind = dp ) :: rcut
152			real( kind = dp ) :: rcut2
153			real( kind = dp ) :: rcut6
154			real( kind = dp ) :: r2
155			real( kind = dp ) :: r6
156
157			real( kind = dp ) :: t6
158			real( kind = dp ) :: t12
159			real( kind = dp ) :: tp6
160			real( kind = dp ) :: tp12
161			real( kind = dp ) :: delta
162
163	gezelter	309	! Look up the correct parameters in the mixing matrix
164	chuckv	298	sigma = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma
165			sigma2 = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma2
166			sigma6 = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma6
167			epsilon = ljMixed(atype1%atype_ident,atype2%atype_ident)%epsilon
168	gezelter	309
169			call getRcut(rcut,rcut2=rcut2,rcut6=rcut6)
170	chuckv	298
171	gezelter	309	r2 = rij * rij
172	chuckv	298	r6 = r2 * r2 * r2
173
174			t6 = sigma6/ r6
175			t12 = t6 * t6
176
177			tp6 = sigma6 / rcut6
178			tp12 = tp6*tp6
179
180	gezelter	309	delta = -4.0_DP * epsilon * (tp12 - tp6)
181
182			if (rij.le.rcut) then
183
184			pot_temp = 4.0E0_DP * epsilon * (t12 - t6) + delta
185	chuckv	298
186	gezelter	309	dudr = 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
187	chuckv	298
188	gezelter	317	drdx = -d(1) / rij
189			drdy = -d(2) / rij
190			drdz = -d(3) / rij
191	gezelter	309
192			fx = dudr * drdx
193			fy = dudr * drdy
194			fz = dudr * drdz
195
196			#ifdef IS_MPI
197			pot_Row(atom1) = pot_Row(atom1) + pot_temp*0.5
198			pot_Col(atom2) = pot_Col(atom2) + pot_temp*0.5
199	chuckv	298
200	gezelter	309	f_Row(1,atom1) = f_Row(1,atom1) + fx
201			f_Row(2,atom1) = f_Row(2,atom1) + fy
202			f_Row(3,atom1) = f_Row(3,atom1) + fz
203	chuckv	298
204	gezelter	309	f_Col(1,atom2) = f_Col(1,atom2) - fx
205			f_Col(2,atom2) = f_Col(2,atom2) - fy
206			f_Col(3,atom2) = f_Col(3,atom2) - fz
207	chuckv	298
208	gezelter	309	#else
209			pot = pot + pot_temp
210
211			f(1,atom1) = f(1,atom1) + fx
212			f(2,atom1) = f(2,atom1) + fy
213			f(3,atom1) = f(3,atom1) + fz
214
215			f(1,atom2) = f(1,atom2) - fx
216			f(2,atom2) = f(2,atom2) - fy
217			f(3,atom2) = f(3,atom2) - fz
218			#endif
219
220			if (doStress()) then
221	gezelter	317	tau_Temp(1) = tau_Temp(1) + fx * d(1)
222			tau_Temp(2) = tau_Temp(2) + fx * d(2)
223			tau_Temp(3) = tau_Temp(3) + fx * d(3)
224			tau_Temp(4) = tau_Temp(4) + fy * d(1)
225			tau_Temp(5) = tau_Temp(5) + fy * d(2)
226			tau_Temp(6) = tau_Temp(6) + fy * d(3)
227			tau_Temp(7) = tau_Temp(7) + fz * d(1)
228			tau_Temp(8) = tau_Temp(8) + fz * d(2)
229			tau_Temp(9) = tau_Temp(9) + fz * d(3)
230	gezelter	309	virial_Temp = virial_Temp + (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
231			endif
232
233			endif
234			return
235			end subroutine do_lj_pair
236
237
238			!! Calculates the mixing for sigma or epslon based on character
239			!! string for initialzition of mixing array.
240	chuckv	298	function calcLJMix(thisParam,param1,param2,status) result(myMixParam)
241			character(len=*) :: thisParam
242			real(kind = dp) :: param1
243			real(kind = dp) :: param2
244			real(kind = dp ) :: myMixParam
245			character(len = getStringLen()) :: thisMixingRule
246			integer, optional :: status
247	gezelter	309
248			!! get the mixing rules from the simulation
249	chuckv	298	thisMixingRule = returnMixingRules()
250			myMixParam = 0.0_dp
251	gezelter	309
252	chuckv	298	if (present(status)) status = 0
253			select case (thisMixingRule)
254	gezelter	309	case ("standard")
255			select case (thisParam)
256			case ("sigma")
257			myMixParam = 0.5_dp * (param1 + param2)
258			case ("epsilon")
259			myMixParam = sqrt(param1 * param2)
260			case default
261			status = -1
262			end select
263	chuckv	298	case("LJglass")
264			case default
265			status = -1
266			end select
267			end function calcLJMix
268	gezelter	309
269			end module lj