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 306 by chuckv, Mon Mar 10 19:26:45 2003 UTC vs.
Revision 328 by gezelter, Wed Mar 12 20:00:58 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.2 2003-03-10 19:26:45 chuckv Exp $, $Date: 2003-03-10 19:26:45 $, $Name: not supported by cvs2svn $, $Revision: 1.2 $
5 > !! @version $Id: calc_LJ_FF.F90,v 1.8 2003-03-12 20:00:58 gezelter Exp $, $Date: 2003-03-12 20:00:58 $, $Name: not supported by cvs2svn $, $Revision: 1.8 $
6  
7 <
8 <
9 < module lj_ff
7 > module lj
8    use simulation
9    use definitions
10 <  use generic_lists
11 <
10 >  use forceGlobals
11 >  use atype_module
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 <
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 <     !! Mass of Particle
30 <     real ( kind = dp )  :: mass = 0.0_dp
31 < !! Lennard-Jones epslon
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 +     !!
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(:,:) :: ljMixed
47 <
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 createMixingList(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 :: i = 0
85 <    integer :: j = 0
70 <    type (atype), pointer :: tmpPtr_i => null()
71 <    type (atype), pointer :: tmpPtr_j => 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 <    
104 <
105 <    tmpPtr_i => ListHead
106 <    tmpPtr_j => tmpPtr_i%next
107 <    do while (associated(tmpPtr_i))
93 <       i = i + 1
94 < ! do self mixing rule
95 <       ljMixed(i,i)%sigma   = tmpPtr_i%sigma
96 <                                                                                                  
97 <       ljMixed(i,i)%sigma2  = (ljMixed(i,i)%sigma) ** 2
98 <                                                                                                  
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)%epsilon = tmpPtr_i%epsilon
109 >       ljMixed(i,i)%tp6     = ljMixed(i,i)%sigma6/rcut6
110  
111 <       j = i + 1
112 <       do while (associated(tmpPtr_j))
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            
122            ljMixed(i,j)%sigma  =  &
123 <               calcLJMix("sigma",tmpPtr_i%sigma, &
124 <               tmpPtr_j%sigma)
123 >               calcLJMix("sigma",mySigma_i, &
124 >               mySigma_j)
125            
109          ljMixed(i,j)%sigma2  = &
110               (ljMixed(i,j)%sigma)**2
111          
126            ljMixed(i,j)%sigma6 = &
127                 (ljMixed(i,j)%sigma)**6
128            
129 +          
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",tmpPtr_i%epsilon, &
137 <               tmpPtr_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
119          ljMixed(j,i)%sigma2  = ljMixed(i,j)%sigma2
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 <
149 <
124 <          tmpPtr_j => tmpPtr_j%next
125 <          j = j + 1
148 >          ljMixed(j,i)%delta   = ljMixed(i,j)%delta
149 >          
150         end do
127 ! advance pointers
128       tmpPtr_i => tmpPtr_i%next
129       if (associated(tmpPtr_i)) then
130          tmpPtr_j => tmpPtr_i%next
131       endif
132      
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 <
159 <
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 < !! Calculates the potential between two lj particles based on two lj_atype pointers, optionally returns second
166 < !! derivatives.
167 <  subroutine getLjForce(r,pot,dudr,atype1,atype2,d2,status)
168 < ! arguments
145 < !! Length of vector between particles
146 <    real( kind = dp ), intent(in)  :: r
147 < !! Potential Energy
148 <    real( kind = dp ), intent(out) :: pot
149 < !! Derivatve wrt postion
150 <    real( kind = dp ), intent(out) :: dudr
151 < !! Second Derivative, optional, used mainly for normal mode calculations.
152 <    real( kind = dp ), intent(out), optional :: d2
153 <    
154 <    type (atype), pointer :: atype1
155 <    type (atype), pointer :: atype2
156 <
157 <    integer, intent(out), optional :: status
158 <
159 < ! local Variables
160 <    real( kind = dp ) :: sigma
161 <    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
164
171      real( kind = dp ) :: rcut
166    real( kind = dp ) :: rcut2
167    real( kind = dp ) :: rcut6
168    real( kind = dp ) :: r2
172      real( kind = dp ) :: r6
170
173      real( kind = dp ) :: t6
174      real( kind = dp ) :: t12
173    real( kind = dp ) :: tp6
174    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      
184    if (present(d2))     doSec = .true.
185    if (present(status)) status = 0
186
187 ! Look up the correct parameters in the mixing matrix
188    sigma    = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma
189    sigma2   = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma2
190    sigma6   = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma6
191    epsilon  = ljMixed(atype1%atype_ident,atype2%atype_ident)%epsilon
192
193
194    
195
196    call getRcut(rcut,rcut2=rcut2,rcut6=rcut6,status=errorStat)
197    
198    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
218 <       if(doSec) d2 = 0.0E0_DP
219 <    endif
220 <                                                                              
221 <  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
234    character(len = getStringLen()) :: thisMixingRule
235    integer, optional :: status
260  
261 < !! get the mixing rules from the simulation
238 <    thisMixingRule = returnMixingRules()
239 <    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
252 <    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 <
259 <
260 < end module lj_ff
280 >  
281 > end module lj

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines