ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/calc_LJ_FF.F90
(Generate patch)

Comparing trunk/OOPSE/libmdtools/calc_LJ_FF.F90 (file contents):
Revision 611 by gezelter, Tue Jul 15 17:10:50 2003 UTC vs.
Revision 1217 by gezelter, Tue Jun 1 21:45:22 2004 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.7 2003-07-15 17:10:50 gezelter Exp $, $Date: 2003-07-15 17:10:50 $, $Name: not supported by cvs2svn $, $Revision: 1.7 $
5 > !! @version $Id: calc_LJ_FF.F90,v 1.22 2004-06-01 21:45:22 gezelter Exp $, $Date: 2004-06-01 21:45:22 $, $Name: not supported by cvs2svn $, $Revision: 1.22 $
6  
7   module lj
8    use definitions
9    use atype_module
10 +  use switcheroo
11    use vector_class
12    use simulation
13   #ifdef IS_MPI
# Line 21 | Line 22 | module lj
22   #include "fForceField.h"
23  
24    integer, save :: LJ_Mixing_Policy
25 <  integer, save :: LJ_rcut
25 >  real(kind=DP), save :: LJ_rcut
26 >  logical, save :: havePolicy = .false.
27 >  logical, save :: haveCut = .false.
28 >  logical, save :: LJ_do_shift = .false.
29    
30    !! Logical has lj force field module been initialized?
31    
# Line 29 | Line 33 | module lj
33    
34    !! Public methods and data
35    public :: init_LJ_FF
36 <  public :: LJ_new_rcut
36 >  public :: setCutoffLJ
37    public :: do_lj_pair
38    
39    type :: lj_mixed_params
# Line 42 | Line 46 | module lj
46       !!
47       real ( kind = dp )  :: tp6
48       real ( kind = dp )  :: tp12
45
49       real ( kind = dp )  :: delta  = 0.0_dp
47
50    end type lj_mixed_params
51    
52    type (lj_mixed_params), dimension(:,:), pointer :: ljMixed
53    
54   contains
55  
56 <  subroutine init_LJ_FF(mix_Policy, rcut, status)
56 >  subroutine init_LJ_FF(mix_Policy, status)
57      integer, intent(in) :: mix_Policy
58      integer, intent(out) :: status
59      integer :: myStatus
58    real(kind=dp) :: rcut
60      
61      if (mix_Policy == LB_MIXING_RULE) then
62         LJ_Mixing_Policy = LB_MIXING_RULE
# Line 68 | Line 69 | contains
69            return
70         endif
71      endif
71    
72    LJ_rcut = rcut
72  
73 <    status = 0
74 <    call createMixingList(myStatus)
75 <    if (myStatus /= 0) then
76 <       status = -1
77 <       return
73 >    havePolicy = .true.
74 >
75 >    if (haveCut) then
76 >       status = 0
77 >       call createMixingList(myStatus)
78 >       if (myStatus /= 0) then
79 >          status = -1
80 >          return
81 >       end if
82 >      
83 >       LJ_FF_initialized = .true.
84      end if
85 <    
81 <    LJ_FF_initialized = .true.
82 <    
85 >  
86    end subroutine init_LJ_FF
87    
88 <  subroutine LJ_new_rcut(rcut, status)
88 >  subroutine setCutoffLJ(rcut, do_shift, status)
89 >    logical, intent(in):: do_shift
90      integer :: status, myStatus
91      real(kind=dp) :: rcut
92  
93 <    LJ_rcut = rcut
93 > #define __FORTRAN90
94 > #include "fSwitchingFunction.h"
95 >
96      status = 0
91    call createMixingList(myStatus)
92    if (myStatus /= 0) then
93       status = -1
94       return
95    end if
97  
98 +    LJ_rcut = rcut
99 +    LJ_do_shift = do_shift
100 +    call set_switch(LJ_SWITCH, rcut, rcut)
101 +    haveCut = .true.
102 +
103 +    if (havePolicy) then
104 +       status = 0
105 +       call createMixingList(myStatus)
106 +       if (myStatus /= 0) then
107 +          status = -1
108 +          return
109 +       end if
110 +      
111 +       LJ_FF_initialized = .true.
112 +    end if    
113 +    
114      return
115 <  end subroutine LJ_new_rcut
115 >  end subroutine setCutoffLJ
116    
117    subroutine createMixingList(status)
118      integer :: nAtypes
# Line 105 | Line 122 | contains
122      real ( kind = dp ) :: mySigma_i,mySigma_j
123      real ( kind = dp ) :: myEpsilon_i,myEpsilon_j
124      real ( kind = dp ) :: rcut6
125 +    logical :: I_isLJ, J_isLJ
126      status = 0
127      
128      nAtypes = getSize(atypes)
# Line 115 | Line 133 | contains
133          
134      if (.not. associated(ljMixed)) then
135         allocate(ljMixed(nAtypes, nAtypes))
136 <    else
119 <       status = -1
120 <       return
121 <    end if
136 >    endif
137  
138      rcut6 = LJ_rcut**6
139  
140 + ! This loops through all atypes, even those that don't support LJ forces.
141      do i = 1, nAtypes
142  
143 <       call getElementProperty(atypes, i, "lj_epsilon", myEpsilon_i)
128 <       call getElementProperty(atypes, i, "lj_sigma",   mySigma_i)
129 <       ! do self mixing rule
130 <       ljMixed(i,i)%sigma   = mySigma_i
131 <      
132 <       ljMixed(i,i)%sigma6  = (ljMixed(i,i)%sigma) ** 6                                                                                                  
133 <       ljMixed(i,i)%tp6     = ljMixed(i,i)%sigma6/rcut6
143 >       call getElementProperty(atypes, i, "is_LJ", I_isLJ)
144  
145 <       ljMixed(i,i)%tp12    = (ljMixed(i,i)%tp6) ** 2
145 >       if (I_isLJ) then
146  
147 <       ljMixed(i,i)%epsilon = myEpsilon_i
148 <
149 <       ljMixed(i,i)%delta = -4.0_DP * ljMixed(i,i)%epsilon * &
150 <            (ljMixed(i,i)%tp12 - ljMixed(i,i)%tp6)
141 <      
142 <       do j = i + 1, nAtypes
143 <          call getElementProperty(atypes,j,"lj_epsilon",myEpsilon_j)
144 <          call getElementProperty(atypes,j,"lj_sigma",  mySigma_j)
147 >          call getElementProperty(atypes, i, "lj_epsilon", myEpsilon_i)
148 >          call getElementProperty(atypes, i, "lj_sigma",   mySigma_i)
149 >          ! do self mixing rule
150 >          ljMixed(i,i)%sigma   = mySigma_i
151            
152 <          ljMixed(i,j)%sigma  =  &
147 <               calcLJMix("sigma",mySigma_i, &
148 <               mySigma_j)
152 >          ljMixed(i,i)%sigma6  = (ljMixed(i,i)%sigma) ** 6
153            
154 <          ljMixed(i,j)%sigma6 = &
151 <               (ljMixed(i,j)%sigma)**6
154 >          ljMixed(i,i)%tp6     = (ljMixed(i,i)%sigma6)/rcut6
155            
156 +          ljMixed(i,i)%tp12    = (ljMixed(i,i)%tp6) ** 2
157            
154          ljMixed(i,j)%tp6     = ljMixed(i,j)%sigma6/rcut6
158            
159 <          ljMixed(i,j)%tp12    = (ljMixed(i,j)%tp6) ** 2
159 >          ljMixed(i,i)%epsilon = myEpsilon_i
160            
161 +          ljMixed(i,i)%delta = -4.0_DP * ljMixed(i,i)%epsilon * &
162 +            (ljMixed(i,i)%tp12 - ljMixed(i,i)%tp6)
163            
164 <          ljMixed(i,j)%epsilon = &
165 <               calcLJMix("epsilon",myEpsilon_i, &
166 <               myEpsilon_j)
167 <          
168 <          ljMixed(i,j)%delta = -4.0_DP * ljMixed(i,j)%epsilon * &
169 <               (ljMixed(i,j)%tp12 - ljMixed(i,j)%tp6)
170 <          
171 <          
172 <          ljMixed(j,i)%sigma   = ljMixed(i,j)%sigma
173 <          ljMixed(j,i)%sigma6  = ljMixed(i,j)%sigma6
174 <          ljMixed(j,i)%tp6     = ljMixed(i,j)%tp6
175 <          ljMixed(j,i)%tp12    = ljMixed(i,j)%tp12
176 <          ljMixed(j,i)%epsilon = ljMixed(i,j)%epsilon
177 <          ljMixed(j,i)%delta   = ljMixed(i,j)%delta
178 <          
179 <       end do
164 >          do j = i + 1, nAtypes
165 >
166 >             call getElementProperty(atypes, j, "is_LJ", J_isLJ)
167 >            
168 >             if (J_isLJ) then                
169 >            
170 >                call getElementProperty(atypes,j,"lj_epsilon",myEpsilon_j)
171 >                call getElementProperty(atypes,j,"lj_sigma",  mySigma_j)
172 >                
173 >                ljMixed(i,j)%sigma  =  &
174 >                     calcLJMix("sigma",mySigma_i, &
175 >                     mySigma_j)
176 >                
177 >                ljMixed(i,j)%sigma6 = &
178 >                     (ljMixed(i,j)%sigma)**6
179 >                
180 >                
181 >                ljMixed(i,j)%tp6     = ljMixed(i,j)%sigma6/rcut6
182 >                
183 >                ljMixed(i,j)%tp12    = (ljMixed(i,j)%tp6) ** 2
184 >                
185 >                
186 >                ljMixed(i,j)%epsilon = &
187 >                     calcLJMix("epsilon",myEpsilon_i, &
188 >                     myEpsilon_j)
189 >                
190 >                ljMixed(i,j)%delta = -4.0_DP * ljMixed(i,j)%epsilon * &
191 >                     (ljMixed(i,j)%tp12 - ljMixed(i,j)%tp6)
192 >                
193 >                
194 >                ljMixed(j,i)%sigma   = ljMixed(i,j)%sigma
195 >                ljMixed(j,i)%sigma6  = ljMixed(i,j)%sigma6
196 >                ljMixed(j,i)%tp6     = ljMixed(i,j)%tp6
197 >                ljMixed(j,i)%tp12    = ljMixed(i,j)%tp12
198 >                ljMixed(j,i)%epsilon = ljMixed(i,j)%epsilon
199 >                ljMixed(j,i)%delta   = ljMixed(i,j)%delta
200 >             endif                
201 >          end do
202 >       endif
203      end do
204      
205    end subroutine createMixingList
206    
207 <  subroutine do_lj_pair(atom1, atom2, d, rij, r2, pot, f, do_pot, do_stress)
207 >  subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
208 >       pot, f, do_pot)
209  
210      integer, intent(in) ::  atom1, atom2
211      real( kind = dp ), intent(in) :: rij, r2
212 <    real( kind = dp ) :: pot
213 <    real( kind = dp ), dimension(3,getNlocal()) :: f    
212 >    real( kind = dp ) :: pot, sw, vpair
213 >    real( kind = dp ), dimension(3,nLocal) :: f    
214      real( kind = dp ), intent(in), dimension(3) :: d
215 <    logical, intent(in) :: do_pot, do_stress
215 >    real( kind = dp ), intent(inout), dimension(3) :: fpair
216 >    logical, intent(in) :: do_pot
217  
218      ! local Variables
219      real( kind = dp ) :: drdx, drdy, drdz
# Line 197 | Line 227 | contains
227      real( kind = dp ) :: delta
228      integer :: id1, id2
229  
230 <
201 <    if (rij.lt.LJ_rcut)  then
202 <
203 <       ! Look up the correct parameters in the mixing matrix
230 >    ! Look up the correct parameters in the mixing matrix
231   #ifdef IS_MPI
232 <       sigma6   = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6
233 <       epsilon  = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon
234 <       delta    = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta
232 >    sigma6   = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6
233 >    epsilon  = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon
234 >    delta    = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta
235   #else
236 <       sigma6   = ljMixed(atid(atom1),atid(atom2))%sigma6
237 <       epsilon  = ljMixed(atid(atom1),atid(atom2))%epsilon
238 <       delta    = ljMixed(atid(atom1),atid(atom2))%delta
236 >    sigma6   = ljMixed(atid(atom1),atid(atom2))%sigma6
237 >    epsilon  = ljMixed(atid(atom1),atid(atom2))%epsilon
238 >    delta    = ljMixed(atid(atom1),atid(atom2))%delta
239   #endif
240 +
241 +    r6 = r2 * r2 * r2
242 +    
243 +    t6  = sigma6/ r6
244 +    t12 = t6 * t6    
245 +  
246 +    pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
247 +    if (LJ_do_shift) then
248 +       pot_temp = pot_temp + delta
249 +    endif
250 +
251 +    vpair = vpair + pot_temp
252        
253 <       r6 = r2 * r2 * r2
253 >    dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
254        
255 <       t6  = sigma6/ r6
256 <       t12 = t6 * t6    
257 <            
219 <       pot_temp = 4.0E0_DP * epsilon * (t12 - t6) + delta      
255 >    drdx = d(1) / rij
256 >    drdy = d(2) / rij
257 >    drdz = d(3) / rij
258        
259 <       dudr = 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
259 >    fx = dudr * drdx
260 >    fy = dudr * drdy
261 >    fz = dudr * drdz
262 >    
263        
223       drdx = d(1) / rij
224       drdy = d(2) / rij
225       drdz = d(3) / rij
226      
227       fx = dudr * drdx
228       fy = dudr * drdy
229       fz = dudr * drdz
230      
264   #ifdef IS_MPI
265 <       if (do_pot) then
266 <          pot_Row(atom1) = pot_Row(atom1) + pot_temp*0.5
267 <          pot_Col(atom2) = pot_Col(atom2) + pot_temp*0.5
268 <       endif
269 <
270 <       f_Row(1,atom1) = f_Row(1,atom1) + fx
271 <       f_Row(2,atom1) = f_Row(2,atom1) + fy
272 <       f_Row(3,atom1) = f_Row(3,atom1) + fz
273 <
274 <       f_Col(1,atom2) = f_Col(1,atom2) - fx
275 <       f_Col(2,atom2) = f_Col(2,atom2) - fy
276 <       f_Col(3,atom2) = f_Col(3,atom2) - fz      
277 <
265 >    if (do_pot) then
266 >       pot_Row(atom1) = pot_Row(atom1) + sw*pot_temp*0.5
267 >       pot_Col(atom2) = pot_Col(atom2) + sw*pot_temp*0.5
268 >    endif
269 >    
270 >    f_Row(1,atom1) = f_Row(1,atom1) + fx
271 >    f_Row(2,atom1) = f_Row(2,atom1) + fy
272 >    f_Row(3,atom1) = f_Row(3,atom1) + fz
273 >    
274 >    f_Col(1,atom2) = f_Col(1,atom2) - fx
275 >    f_Col(2,atom2) = f_Col(2,atom2) - fy
276 >    f_Col(3,atom2) = f_Col(3,atom2) - fz      
277 >    
278   #else
279 <       if (do_pot) pot = pot + pot_temp
279 >    if (do_pot) pot = pot + sw*pot_temp
280  
281 <       f(1,atom1) = f(1,atom1) + fx
282 <       f(2,atom1) = f(2,atom1) + fy
283 <       f(3,atom1) = f(3,atom1) + fz
284 <
285 <       f(1,atom2) = f(1,atom2) - fx
286 <       f(2,atom2) = f(2,atom2) - fy
287 <       f(3,atom2) = f(3,atom2) - fz
281 >    f(1,atom1) = f(1,atom1) + fx
282 >    f(2,atom1) = f(2,atom1) + fy
283 >    f(3,atom1) = f(3,atom1) + fz
284 >    
285 >    f(1,atom2) = f(1,atom2) - fx
286 >    f(2,atom2) = f(2,atom2) - fy
287 >    f(3,atom2) = f(3,atom2) - fz
288   #endif
289 <      
257 <       if (do_stress) then
258 <
289 >        
290   #ifdef IS_MPI
291 <          id1 = tagRow(atom1)
292 <          id2 = tagColumn(atom2)
291 >    id1 = AtomRowToGlobal(atom1)
292 >    id2 = AtomColToGlobal(atom2)
293   #else
294 <          id1 = atom1
295 <          id2 = atom2
294 >    id1 = atom1
295 >    id2 = atom2
296   #endif
297  
298 <          if (molMembershipList(id1) .ne. molMembershipList(id2)) then
268 <
269 <             ! because the d vector is the rj - ri vector, and
270 <             ! because fx, fy, fz are the force on atom i, we need a
271 <             ! negative sign here:
272 <
273 <             tau_Temp(1) = tau_Temp(1) - d(1) * fx
274 <             tau_Temp(2) = tau_Temp(2) - d(1) * fy
275 <             tau_Temp(3) = tau_Temp(3) - d(1) * fz
276 <             tau_Temp(4) = tau_Temp(4) - d(2) * fx
277 <             tau_Temp(5) = tau_Temp(5) - d(2) * fy
278 <             tau_Temp(6) = tau_Temp(6) - d(2) * fz
279 <             tau_Temp(7) = tau_Temp(7) - d(3) * fx
280 <             tau_Temp(8) = tau_Temp(8) - d(3) * fy
281 <             tau_Temp(9) = tau_Temp(9) - d(3) * fz
282 <            
283 <             virial_Temp = virial_Temp + &
284 <                  (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
285 <
286 <          endif
287 <       endif
298 >    if (molMembershipList(id1) .ne. molMembershipList(id2)) then
299        
300 +       fpair(1) = fpair(1) + fx
301 +       fpair(2) = fpair(2) + fy
302 +       fpair(3) = fpair(3) + fz
303 +
304      endif
305 +
306      return    
307 <      
307 >    
308    end subroutine do_lj_pair
309 <
310 <
309 >  
310 >  
311    !! Calculates the mixing for sigma or epslon
312 <
312 >  
313    function calcLJMix(thisParam,param1,param2,status) result(myMixParam)
314      character(len=*) :: thisParam
315      real(kind = dp)  :: param1
# Line 306 | Line 322 | contains
322      
323      if (present(status)) status = 0
324      select case (LJ_Mixing_Policy)
325 <    case (LB_MIXING_RULE)
325 >    case (1)
326         select case (thisParam)
327         case ("sigma")
328            myMixParam = 0.5_dp * (param1 + param2)

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines