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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines