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.9 2003-03-12 22:27:59 gezelter Exp $, $Date: 2003-03-12 22:27:59 $, $Name: not supported by cvs2svn $, $Revision: 1.9 $

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