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.5 2003-03-12 14:55:29 chuckv Exp $, $Date: 2003-03-12 14:55:29 $, $Name: not supported by cvs2svn $, $Revision: 1.5 $


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 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_ljFF(status)
    integer, intent(out) :: status
    integer :: myStatus
    
    status = 0
    call createMixingList(myStatus)
    if (myStatus /= 0) then
       status = -1
       return
    end if
    
  end subroutine init_ljFF
  
  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 ) :: rcut, rcut6
    status = 0
    
    nAtypes = getSize(atype)
    if (listSize == 0) then
       status = -1
       return
    end if
    
    
    if (.not. associated(ljMixed)) then
       allocate(ljMixed(listSize,listSize))
    else
       status = -1
       return
    end if
    call getRcut(rcut, rcut6=rcut6)
    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, 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
#ifdef IS_MPI
    sigma    = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma
    sigma6   = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6
    tp6      = ljMixed(atid_Row(atom1),atid_col(atom2))%tp6
    tp12     = ljMixed(atid_Row(atom1),atid_Col(atom2))%tp12
    epsilon  = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon
    delta    = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta
#else
    sigma    = ljMixed(atid(atom1),atid(atom2))%sigma
    sigma6   = ljMixed(atid(atom1),atid(atom2))%sigma6
    tp6      = ljMixed(atid(atom1),atid(atom2))%tp6
    tp12     = ljMixed(atid(atom1),atid(atom2))%tp12
    epsilon  = ljMixed(atid(atom1),atid(atom2))%epsilon
    delta    = ljMixed(atid(atom1),atid(atom2))%delta
#endif

   
    call getRcut(rcut)
    
    r2 = rij * rij
    r6 = r2 * r2 * r2

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