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.6 2003-03-12 15:39:01 gezelter Exp $, $Date: 2003-03-12 15:39:01 $, $Name: not supported by cvs2svn $, $Revision: 1.6 $

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(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, 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, r2, pot, f)

    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

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