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.3 2003-03-10 23:19:23 gezelter Exp $, $Date: 2003-03-10 23:19:23 $, $Name: not supported by cvs2svn $, $Revision: 1.3 $


module lj
  use simulation
  use definitions
  use forceGlobals
  use atype_typedefs
  use generic_atypes
#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, atype1, atype2, dx, dy, dz, rij, &
       pot, f)

    integer, intent(in) ::  atom1, atom2
    real( kind = dp ), intent(in) :: dx, dy, dz, rij
    real( kind = dp ) :: pot
    real( kind = dp ), dimension(3, getNlocal()) :: f    
    type (atype), pointer :: atype1
    type (atype), pointer :: atype2

    ! 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 = -dx / rij
       drdy = -dy / rij
       drdz = -dz / 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 * dx
          tau_Temp(2) = tau_Temp(2) + fx * dy
          tau_Temp(3) = tau_Temp(3) + fx * dz
          tau_Temp(4) = tau_Temp(4) + fy * dx
          tau_Temp(5) = tau_Temp(5) + fy * dy
          tau_Temp(6) = tau_Temp(6) + fy * dz
          tau_Temp(7) = tau_Temp(7) + fz * dx
          tau_Temp(8) = tau_Temp(8) + fz * dy
          tau_Temp(9) = tau_Temp(9) + fz * dz
          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:	309
Committed:	Mon Mar 10 23:19:23 2003 UTC (21 years, 3 months ago) by gezelter
File size:	7682 byte(s)
Log Message:	Massive rewrite underway. This way be dragons.
#	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	309	!! @version $Id: calc_LJ_FF.F90,v 1.3 2003-03-10 23:19:23 gezelter Exp $, $Date: 2003-03-10 23:19:23 $, $Name: not supported by cvs2svn $, $Revision: 1.3 $
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			use generic_atypes
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	309	subroutine do_lj_pair(atom1, atom2, atype1, atype2, dx, dy, dz, rij, &
135			pot, f)
136	chuckv	298
137	gezelter	309	integer, intent(in) :: atom1, atom2
138			real( kind = dp ), intent(in) :: dx, dy, dz, rij
139			real( kind = dp ) :: pot
140			real( kind = dp ), dimension(3, getNlocal()) :: f
141	chuckv	306	type (atype), pointer :: atype1
142			type (atype), pointer :: atype2
143	chuckv	298
144	gezelter	309	! local Variables
145			real( kind = dp ) :: drdx, drdy, drdz
146			real( kind = dp ) :: fx, fy, fz
147			real( kind = dp ) :: pot_temp, dudr
148	chuckv	298	real( kind = dp ) :: sigma
149			real( kind = dp ) :: sigma2
150			real( kind = dp ) :: sigma6
151			real( kind = dp ) :: epsilon
152
153			real( kind = dp ) :: rcut
154			real( kind = dp ) :: rcut2
155			real( kind = dp ) :: rcut6
156			real( kind = dp ) :: r2
157			real( kind = dp ) :: r6
158
159			real( kind = dp ) :: t6
160			real( kind = dp ) :: t12
161			real( kind = dp ) :: tp6
162			real( kind = dp ) :: tp12
163			real( kind = dp ) :: delta
164
165	gezelter	309	! Look up the correct parameters in the mixing matrix
166	chuckv	298	sigma = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma
167			sigma2 = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma2
168			sigma6 = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma6
169			epsilon = ljMixed(atype1%atype_ident,atype2%atype_ident)%epsilon
170	gezelter	309
171			call getRcut(rcut,rcut2=rcut2,rcut6=rcut6)
172	chuckv	298
173	gezelter	309	r2 = rij * rij
174	chuckv	298	r6 = r2 * r2 * r2
175
176			t6 = sigma6/ r6
177			t12 = t6 * t6
178
179			tp6 = sigma6 / rcut6
180			tp12 = tp6*tp6
181
182	gezelter	309	delta = -4.0_DP * epsilon * (tp12 - tp6)
183
184			if (rij.le.rcut) then
185
186			pot_temp = 4.0E0_DP * epsilon * (t12 - t6) + delta
187	chuckv	298
188	gezelter	309	dudr = 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
189	chuckv	298
190	gezelter	309	drdx = -dx / rij
191			drdy = -dy / rij
192			drdz = -dz / rij
193
194			fx = dudr * drdx
195			fy = dudr * drdy
196			fz = dudr * drdz
197
198			#ifdef IS_MPI
199			pot_Row(atom1) = pot_Row(atom1) + pot_temp*0.5
200			pot_Col(atom2) = pot_Col(atom2) + pot_temp*0.5
201	chuckv	298
202	gezelter	309	f_Row(1,atom1) = f_Row(1,atom1) + fx
203			f_Row(2,atom1) = f_Row(2,atom1) + fy
204			f_Row(3,atom1) = f_Row(3,atom1) + fz
205	chuckv	298
206	gezelter	309	f_Col(1,atom2) = f_Col(1,atom2) - fx
207			f_Col(2,atom2) = f_Col(2,atom2) - fy
208			f_Col(3,atom2) = f_Col(3,atom2) - fz
209	chuckv	298
210	gezelter	309	#else
211			pot = pot + pot_temp
212
213			f(1,atom1) = f(1,atom1) + fx
214			f(2,atom1) = f(2,atom1) + fy
215			f(3,atom1) = f(3,atom1) + fz
216
217			f(1,atom2) = f(1,atom2) - fx
218			f(2,atom2) = f(2,atom2) - fy
219			f(3,atom2) = f(3,atom2) - fz
220			#endif
221
222			if (doStress()) then
223			tau_Temp(1) = tau_Temp(1) + fx * dx
224			tau_Temp(2) = tau_Temp(2) + fx * dy
225			tau_Temp(3) = tau_Temp(3) + fx * dz
226			tau_Temp(4) = tau_Temp(4) + fy * dx
227			tau_Temp(5) = tau_Temp(5) + fy * dy
228			tau_Temp(6) = tau_Temp(6) + fy * dz
229			tau_Temp(7) = tau_Temp(7) + fz * dx
230			tau_Temp(8) = tau_Temp(8) + fz * dy
231			tau_Temp(9) = tau_Temp(9) + fz * dz
232			virial_Temp = virial_Temp + (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
233			endif
234
235			endif
236			return
237			end subroutine do_lj_pair
238
239
240			!! Calculates the mixing for sigma or epslon based on character
241			!! string for initialzition of mixing array.
242	chuckv	298	function calcLJMix(thisParam,param1,param2,status) result(myMixParam)
243			character(len=*) :: thisParam
244			real(kind = dp) :: param1
245			real(kind = dp) :: param2
246			real(kind = dp ) :: myMixParam
247			character(len = getStringLen()) :: thisMixingRule
248			integer, optional :: status
249	gezelter	309
250			!! get the mixing rules from the simulation
251	chuckv	298	thisMixingRule = returnMixingRules()
252			myMixParam = 0.0_dp
253	gezelter	309
254	chuckv	298	if (present(status)) status = 0
255			select case (thisMixingRule)
256	gezelter	309	case ("standard")
257			select case (thisParam)
258			case ("sigma")
259			myMixParam = 0.5_dp * (param1 + param2)
260			case ("epsilon")
261			myMixParam = sqrt(param1 * param2)
262			case default
263			status = -1
264			end select
265	chuckv	298	case("LJglass")
266			case default
267			status = -1
268			end select
269			end function calcLJMix
270	gezelter	309
271			end module lj