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.7 2003-03-12 19:31:55 gezelter Exp $, $Date: 2003-03-12 19:31:55 $, $Name: not supported by cvs2svn $, $Revision: 1.7 $

module lj
  use simulation
  use definitions
  use forceGlobals
  use atype_typedefs
  use vector_class
#ifdef IS_MPI
  use mpiSimulation
#endif
  implicit none
  PRIVATE

  integer, save :: LJ_Mixing_Policy
  integer, parameter :: LB_Mixing_Rule = 1
  integer, parameter :: Explicit_Mixing_Rule = 2
  
  !! Logical has lj force field module been initialized?
  
  logical, save :: LJ_FF_initialized = .false.
  
  !! Public methods and data
  public :: init_LJ_FF
  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)
    character(len=100), intent(in) :: mix_Policy
    integer, intent(out) :: status
    integer :: myStatus
    
    if ((mix_Policy == "standard") .or.  (mix_Policy == "LB")) then
       LJ_Mixing_Policy = LB_Mixing_Rule
    else
       if (mix_Policy == "explicit") then
          LJ_Mixing_Policy = Explicit_Mixing_Rule
       else
          write(*,*) 'Unknown Mixing Policy!'
          status = -1
          return
       endif
    endif
        
    status = 0
    call createMixingList(myStatus)
    if (myStatus /= 0) then
       status = -1
       return
    end if
    
    LJ_FF_initialized = .true.

  end subroutine init_LJ_FF
  
  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(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
    call getRcut(rcut, rc6=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, 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 ) :: rcut
    real( kind = dp ) :: r6
    real( kind = dp ) :: t6
    real( kind = dp ) :: t12
    real( kind = dp ) :: delta

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