ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/calc_LJ_FF.F90
(Generate patch)

Comparing trunk/OOPSE_old/src/mdtools/libmdCode/calc_LJ_FF.F90 (file contents):
Revision 298 by chuckv, Fri Mar 7 18:26:30 2003 UTC vs.
Revision 329 by gezelter, Wed Mar 12 22:27:59 2003 UTC

# Line 2 | Line 2
2   !! Corresponds to the force field defined in lj_FF.cpp
3   !! @author Charles F. Vardeman II
4   !! @author Matthew Meineke
5 < !! @version $Id: calc_LJ_FF.F90,v 1.1 2003-03-07 18:26:30 chuckv Exp $, $Date: 2003-03-07 18:26:30 $, $Name: not supported by cvs2svn $, $Revision: 1.1 $
5 > !! @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  
7 <
8 <
9 < module lj_ff
7 > module lj
8    use simulation
9    use definitions
10 <  use generic_atypes
11 <
10 >  use atype_module
11 >  use vector_class
12   #ifdef IS_MPI
13    use mpiSimulation
14   #endif
15    implicit none
16    PRIVATE
17  
18 < !! Logical has lj force field module been initialized?
19 <  logical, save :: isljFFinit = .false.
20 <
21 <
22 < !! Public methods and data
23 <  public :: getLjPot
24 <  public :: init_ljFF
25 <
26 <  type :: lj_mixed
27 <     !! Mass of Particle
28 <     real ( kind = dp )  :: mass = 0.0_dp
29 < !! Lennard-Jones epslon
18 >  integer, save :: LJ_Mixing_Policy
19 >  integer, parameter :: LB_Mixing_Rule = 1
20 >  integer, parameter :: Explicit_Mixing_Rule = 2
21 >  
22 >  !! Logical has lj force field module been initialized?
23 >  
24 >  logical, save :: LJ_FF_initialized = .false.
25 >  
26 >  !! Public methods and data
27 >  public :: init_LJ_FF
28 >  public :: do_lj_pair
29 >  
30 >  type :: lj_mixed_params
31 >     !! Lennard-Jones epsilon
32       real ( kind = dp )  :: epsilon = 0.0_dp
33 < !! Lennard-Jones Sigma
33 >     !! Lennard-Jones Sigma
34       real ( kind = dp )  :: sigma = 0.0_dp
35 < !! Lennard-Jones Sigma Squared
36 <     real ( kind = dp )  :: sigma2 = 0.0_dp
37 < !! Lennard-Jones Sigma to sixth
35 >     !! Lennard-Jones Sigma to sixth
36       real ( kind = dp )  :: sigma6 = 0.0_dp
37 <  end type lj_mixed
38 <  
37 >     !!
38 >     real ( kind = dp )  :: tp6
39 >     real ( kind = dp )  :: tp12
40  
41 +     real ( kind = dp )  :: delta  = 0.0_dp
42  
43 +  end type lj_mixed_params
44 +  
45 +  type (lj_mixed_params), dimension(:,:), pointer :: ljMixed
46 +  
47   contains
48  
49 <  subroutine init_ljFF(ListHead,status)
50 <    type (atype),pointer :: ListHead
47 <    integer :: nTypes
49 >  subroutine init_LJ_FF(mix_Policy, status)
50 >    character(len=100), intent(in) :: mix_Policy
51      integer, intent(out) :: status
49
52      integer :: myStatus
53 <
53 >    
54 >    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      status = 0
67 <    call createMixinList(ListHead,myStatus)
67 >    call createMixingList(myStatus)
68      if (myStatus /= 0) then
69         status = -1
70         return
71      end if
72 +    
73 +    LJ_FF_initialized = .true.
74  
75 <  end subroutine init_ljFF
76 <
77 <  subroutine createMixingList(ListHead,status)
78 <    type (atype), pointer :: ListHead
63 <    integer :: listSize
75 >  end subroutine init_LJ_FF
76 >  
77 >  subroutine createMixingList(status)
78 >    integer :: nAtypes
79      integer :: status
80      integer :: i
81      integer :: j
82 <
83 <    integer :: outerCounter = 0
84 <    integer :: innerCounter = 0
70 <    type (atype), pointer :: tmpPtrOuter => null()
71 <    type (atype), pointer :: tmpPtrInner => null()
82 >    real ( kind = dp ) :: mySigma_i,mySigma_j
83 >    real ( kind = dp ) :: myEpsilon_i,myEpsilon_j
84 >    real ( kind = dp ) :: rcut, rcut6
85      status = 0
86 <
87 <    listSize = getListLen(ListHead)
88 <    if (listSize == 0) then
86 >    
87 >    nAtypes = getSize(atypes)
88 >    if (nAtypes == 0) then
89         status = -1
90         return
91      end if
92 <  
80 <
92 >        
93      if (.not. associated(ljMixed)) then
94 <       allocate(ljMixed(listSize,listSize))
94 >       allocate(ljMixed(nAtypes, nAtypes))
95      else
96         status = -1
97         return
98      end if
99 +    call getRcut(rcut, rc6=rcut6)
100 +    do i = 1, nAtypes
101  
102 <    
102 >       call getElementProperty(atypes,i,"lj_epsilon",myEpsilon_i)
103 >       call getElementProperty(atypes,i,"lj_sigma",  mySigma_i)
104 >       ! do self mixing rule
105 >       ljMixed(i,i)%sigma   = mySigma_i
106 >      
107 >       ljMixed(i,i)%sigma6  = (ljMixed(i,i)%sigma) ** 6                                                                                                  
108 >       ljMixed(i,i)%tp6     = ljMixed(i,i)%sigma6/rcut6
109  
110 <    tmpPtrOuter => ljListHead
91 <    tmpPtrInner => tmpPtrOuter%next
92 <    do while (associated(tmpPtrOuter))
93 <       outerCounter = outerCounter + 1
94 < ! do self mixing rule
95 <       ljMixed(outerCounter,outerCounter)%sigma  = tmpPtrOuter%sigma
96 <                                                                                                  
97 <       ljMixed(outerCounter,outerCounter)%sigma2  = ljMixed(outerCounter,outerCounter)%sigma &
98 <            * ljMixed(outerCounter,outerCounter)%sigma
99 <                                                                                                  
100 <       ljMixed(outerCounter,outerCounter)%sigma6 = ljMixed(outerCounter,outerCounter)%sigma2 &
101 <            * ljMixed(outerCounter,outerCounter)%sigma2 &
102 <            * ljMixed(outerCounter,outerCounter)%sigma2
103 <                                                                                                  
104 <       ljMixed(outerCounter,outerCounter)%epsilon = tmpPtrOuter%epsilon
110 >       ljMixed(i,i)%tp12    = (ljMixed(i,i)%tp6) ** 2
111  
112 <       innerCounter = outerCounter + 1
113 <       do while (associated(tmpPtrInner))
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            
121 <          ljMixed(outerCounter,innerCounter)%sigma  =  &
122 <               calcLJMix("sigma",tmpPtrOuter%sigma, &
123 <               tmpPtrInner%sigma)
121 >          ljMixed(i,j)%sigma  =  &
122 >               calcLJMix("sigma",mySigma_i, &
123 >               mySigma_j)
124            
125 <          ljMixed(outerCounter,innerCounter)%sigma2  = &
126 <               (ljMixed(outerCounter,innerCounter)%sigma)**2
125 >          ljMixed(i,j)%sigma6 = &
126 >               (ljMixed(i,j)%sigma)**6
127            
116          ljMixed(outerCounter,innerCounter)%sigma6 = &
117               (ljMixed(outerCounter,innerCounter)%sigma)**6
128            
129 <          ljMixed(outerCounter,innerCounter)%epsilon = &
130 <               calcLJMix("epsilon",tmpPtrOuter%epsilon, &
131 <               tmpPtrInner%epsilon)
132 <          ljMixed(innerCounter,outerCounter)%sigma = ljMixed(outerCounter,innerCounter)%sigma
133 <          ljMixed(innerCounter,outerCounter)%sigma2 = ljMixed(outerCounter,innerCounter)%sigma2
134 <          ljMixed(innerCounter,outerCounter)%sigma6 = ljMixed(outerCounter,innerCounter)%sigma6
135 <          ljMixed(innerCounter,outerCounter)%epsilon = ljMixed(outerCounter,innerCounter)%epsilon
136 <
137 <
138 <          tmpPtrInner => tmpPtrInner%next
139 <          innerCounter = innerCounter + 1
129 >          ljMixed(i,j)%tp6     = ljMixed(i,j)%sigma6/rcut6
130 >          
131 >          ljMixed(i,j)%tp12    = (ljMixed(i,j)%tp6) ** 2
132 >          
133 >          
134 >          ljMixed(i,j)%epsilon = &
135 >               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 >          ljMixed(j,i)%sigma   = ljMixed(i,j)%sigma
143 >          ljMixed(j,i)%sigma6  = ljMixed(i,j)%sigma6
144 >          ljMixed(j,i)%tp6     = ljMixed(i,j)%tp6
145 >          ljMixed(j,i)%tp12    = ljMixed(i,j)%tp12
146 >          ljMixed(j,i)%epsilon = ljMixed(i,j)%epsilon
147 >          ljMixed(j,i)%delta   = ljMixed(i,j)%delta
148 >          
149         end do
131 ! advance pointers
132       tmpPtrOuter => tmpPtrOuter%next
133       if (associated(tmpPtrOuter)) then
134          tmpPtrInner => tmpPtrOuter%next
135       endif
136      
150      end do
151 <
151 >    
152    end subroutine createMixingList
153 +  
154  
155 +  subroutine do_lj_pair(atom1, atom2, d, rij, r2, pot, f, do_pot, do_stress)
156  
157 +    integer, intent(in) ::  atom1, atom2
158 +    real( kind = dp ), intent(in) :: rij, r2
159 +    real( kind = dp ) :: pot
160 +    real( kind = dp ), dimension(3, getNlocal()) :: f    
161 +    real( kind = dp ), intent(in), dimension(3) :: d
162 +    logical, intent(in) :: do_pot, do_stress
163  
164 <
165 <
166 < !! Calculates the potential between two lj particles based on two lj_atype pointers, optionally returns second
167 < !! derivatives.
147 <  subroutine getLjForce(r,pot,dudr,atype1,atype2,d2,status)
148 < ! arguments
149 < !! Length of vector between particles
150 <    real( kind = dp ), intent(in)  :: r
151 < !! Potential Energy
152 <    real( kind = dp ), intent(out) :: pot
153 < !! Derivatve wrt postion
154 <    real( kind = dp ), intent(out) :: dudr
155 < !! Second Derivative, optional, used mainly for normal mode calculations.
156 <    real( kind = dp ), intent(out), optional :: d2
157 <    
158 <    type (lj_atype), pointer :: atype1
159 <    type (lj_atype), pointer :: atype2
160 <
161 <    integer, intent(out), optional :: status
162 <
163 < ! local Variables
164 <    real( kind = dp ) :: sigma
165 <    real( kind = dp ) :: sigma2
164 >    ! local Variables
165 >    real( kind = dp ) :: drdx, drdy, drdz
166 >    real( kind = dp ) :: fx, fy, fz
167 >    real( kind = dp ) :: pot_temp, dudr
168      real( kind = dp ) :: sigma6
169      real( kind = dp ) :: epsilon
168
170      real( kind = dp ) :: rcut
170    real( kind = dp ) :: rcut2
171    real( kind = dp ) :: rcut6
172    real( kind = dp ) :: r2
171      real( kind = dp ) :: r6
174
172      real( kind = dp ) :: t6
173      real( kind = dp ) :: t12
177    real( kind = dp ) :: tp6
178    real( kind = dp ) :: tp12
174      real( kind = dp ) :: delta
175  
176 <    logical :: doSec = .false.
177 <
178 <    integer :: errorStat
179 <
180 < !! Optional Argument Checking
181 < ! Check to see if we need to do second derivatives
176 >    ! Look up the correct parameters in the mixing matrix
177 > #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 >  
187 >    call getRcut(rcut)
188      
188    if (present(d2))     doSec = .true.
189    if (present(status)) status = 0
190
191 ! Look up the correct parameters in the mixing matrix
192    sigma    = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma
193    sigma2   = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma2
194    sigma6   = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma6
195    epsilon  = ljMixed(atype1%atype_ident,atype2%atype_ident)%epsilon
196
197
198    
199
200    call getRcut(rcut,rcut2=rcut2,rcut6=rcut6,status=errorStat)
201    
202    r2 = r * r
189      r6 = r2 * r2 * r2
190  
191      t6  = sigma6/ r6
192      t12 = t6 * t6
193 +      
194 +    if (rij.le.rcut) then
195 +      
196 +       pot_temp = 4.0E0_DP * epsilon * (t12 - t6) + delta      
197  
198 +       dudr = 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
199  
200 <                                                                              
201 <    tp6 = sigma6 / rcut6
202 <    tp12 = tp6*tp6
203 <                                                                              
204 <    delta = -4.0E0_DP*epsilon * (tp12 - tp6)
205 <                                                                              
206 <    if (r.le.rcut) then
207 <       pot = 4.0E0_DP * epsilon * (t12 - t6) + delta
208 <       dudr = 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / r
209 <       if(doSec)  d2 = 24.0E0_DP * epsilon * (26.0E0_DP*t12 - 7.0E0_DP*t6)/r/r
210 <    else
211 <       pot = 0.0E0_DP
212 <       dudr = 0.0E0_DP
222 <       if(doSec) d2 = 0.0E0_DP
223 <    endif
224 <                                                                              
225 <  return
200 >       drdx = -d(1) / rij
201 >       drdy = -d(2) / rij
202 >       drdz = -d(3) / rij
203 >      
204 >       fx = dudr * drdx
205 >       fy = dudr * drdy
206 >       fz = dudr * drdz
207 >    
208 > #ifdef IS_MPI
209 >       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  
214 +       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  
218 +       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  
222 <  end subroutine getLJForce
222 > #else
223 >       if (do_pot) 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 < !! Calculates the mixing for sigma or epslon based on character string for initialzition of mixing array.
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 (do_stress) then
235 >          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 >          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
253 >
254    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
238    character(len = getStringLen()) :: thisMixingRule
239    integer, optional :: status
259  
260 < !! get the mixing rules from the simulation
242 <    thisMixingRule = returnMixingRules()
243 <    myMixParam = 0.0_dp
260 >    integer, optional :: status  
261  
262 +    myMixParam = 0.0_dp
263 +    
264      if (present(status)) status = 0
265 <    select case (thisMixingRule)
266 <       case ("standard")
267 <          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
256 <    case("LJglass")
265 >    select case (LJ_Mixing_Policy)
266 >    case (LB_Mixing_Rule)
267 >       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      case default
276         status = -1
277      end select
278    end function calcLJMix
279 <
280 <
263 <
264 < end module lj_ff
279 >  
280 > end module lj

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines