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 1159 by gezelter, Fri May 7 21:35:05 2004 UTC vs.
Revision 1160 by gezelter, Tue May 11 21:31:15 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.18 2004-05-07 21:35:04 gezelter Exp $, $Date: 2004-05-07 21:35:04 $, $Name: not supported by cvs2svn $, $Revision: 1.18 $
5 > !! @version $Id: calc_LJ_FF.F90,v 1.19 2004-05-11 21:31:14 gezelter Exp $, $Date: 2004-05-11 21:31:14 $, $Name: not supported by cvs2svn $, $Revision: 1.19 $
6  
7   module lj
8    use definitions
# Line 23 | Line 23 | module lj
23  
24    integer, save :: LJ_Mixing_Policy
25    real(kind=DP), save :: LJ_rcut
26 <  logical, save :: havePolicy = .false., haveCut = .false.
27 <
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 45 | Line 46 | module lj
46       !!
47       real ( kind = dp )  :: tp6
48       real ( kind = dp )  :: tp12
48 –
49       real ( kind = dp )  :: delta  = 0.0_dp
50 –
51 –
50    end type lj_mixed_params
51    
52    type (lj_mixed_params), dimension(:,:), pointer :: ljMixed
# Line 87 | Line 85 | contains
85    
86    end subroutine init_LJ_FF
87    
88 <  subroutine setCutoffLJ(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  
# Line 97 | Line 96 | contains
96      status = 0
97  
98      LJ_rcut = rcut
99 +    LJ_do_shift = do_shift
100      call set_switch(LJ_SWITCH, rcut, rcut)
101 –
102 – !!$    ! ATTENTION! This is a hardwiring of rcut!
103 – !!$    LJ_rcut = 9.0d0
101      haveCut = .true.
102  
103      if (havePolicy) then
# Line 218 | Line 215 | contains
215      real( kind = dp ) :: delta
216      integer :: id1, id2
217  
218 <    if (rij.lt.LJ_rcut)  then
222 <
223 <       ! Look up the correct parameters in the mixing matrix
218 >    ! Look up the correct parameters in the mixing matrix
219   #ifdef IS_MPI
220 <       sigma6   = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6
221 <       epsilon  = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon
222 <       !delta    = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta
220 >    sigma6   = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6
221 >    epsilon  = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon
222 >    delta    = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta
223   #else
224 <       sigma6   = ljMixed(atid(atom1),atid(atom2))%sigma6
225 <       epsilon  = ljMixed(atid(atom1),atid(atom2))%epsilon
226 <       !delta    = ljMixed(atid(atom1),atid(atom2))%delta
224 >    sigma6   = ljMixed(atid(atom1),atid(atom2))%sigma6
225 >    epsilon  = ljMixed(atid(atom1),atid(atom2))%epsilon
226 >    delta    = ljMixed(atid(atom1),atid(atom2))%delta
227   #endif
228 +
229 +    r6 = r2 * r2 * r2
230 +    
231 +    t6  = sigma6/ r6
232 +    t12 = t6 * t6    
233 +  
234 +    pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
235 +    if (LJ_do_shift) then
236 +       pot_temp = pot_temp + delta
237 +    endif
238 +
239 +    vpair = vpair + pot_temp
240        
241 <       r6 = r2 * r2 * r2
241 >    dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
242        
243 <       t6  = sigma6/ r6
244 <       t12 = t6 * t6    
245 <            
239 <       !pot_temp = 4.0E0_DP * epsilon * (t12 - t6) + delta      
240 <       pot_temp = 4.0E0_DP * epsilon * (t12 - t6) * sw
241 <       vpair = vpair + pot_temp
243 >    drdx = d(1) / rij
244 >    drdy = d(2) / rij
245 >    drdz = d(3) / rij
246        
247 <       dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
247 >    fx = dudr * drdx
248 >    fy = dudr * drdy
249 >    fz = dudr * drdz
250 >    
251        
245 –       drdx = d(1) / rij
246 –       drdy = d(2) / rij
247 –       drdz = d(3) / rij
248 –      
249 –       fx = dudr * drdx
250 –       fy = dudr * drdy
251 –       fz = dudr * drdz
252 –      
253 –      
252   #ifdef IS_MPI
253 <       if (do_pot) then
254 <          pot_Row(atom1) = pot_Row(atom1) + pot_temp*0.5
255 <          pot_Col(atom2) = pot_Col(atom2) + pot_temp*0.5
256 <       endif
257 <
258 <       f_Row(1,atom1) = f_Row(1,atom1) + fx
259 <       f_Row(2,atom1) = f_Row(2,atom1) + fy
260 <       f_Row(3,atom1) = f_Row(3,atom1) + fz
261 <
262 <       f_Col(1,atom2) = f_Col(1,atom2) - fx
263 <       f_Col(2,atom2) = f_Col(2,atom2) - fy
264 <       f_Col(3,atom2) = f_Col(3,atom2) - fz      
265 <
266 < #else
267 <       if (do_pot) pot = pot + pot_temp
270 <
271 <       f(1,atom1) = f(1,atom1) + fx
272 <       f(2,atom1) = f(2,atom1) + fy
273 <       f(3,atom1) = f(3,atom1) + fz
253 >    if (do_pot) then
254 >       pot_Row(atom1) = pot_Row(atom1) + sw*pot_temp*0.5
255 >       pot_Col(atom2) = pot_Col(atom2) + sw*pot_temp*0.5
256 >    endif
257 >    
258 >    f_Row(1,atom1) = f_Row(1,atom1) + fx
259 >    f_Row(2,atom1) = f_Row(2,atom1) + fy
260 >    f_Row(3,atom1) = f_Row(3,atom1) + fz
261 >    
262 >    f_Col(1,atom2) = f_Col(1,atom2) - fx
263 >    f_Col(2,atom2) = f_Col(2,atom2) - fy
264 >    f_Col(3,atom2) = f_Col(3,atom2) - fz      
265 >    
266 > #else
267 >    if (do_pot) pot = pot + sw*pot_temp
268  
269 <       f(1,atom2) = f(1,atom2) - fx
270 <       f(2,atom2) = f(2,atom2) - fy
271 <       f(3,atom2) = f(3,atom2) - fz
269 >    f(1,atom1) = f(1,atom1) + fx
270 >    f(2,atom1) = f(2,atom1) + fy
271 >    f(3,atom1) = f(3,atom1) + fz
272 >    
273 >    f(1,atom2) = f(1,atom2) - fx
274 >    f(2,atom2) = f(2,atom2) - fy
275 >    f(3,atom2) = f(3,atom2) - fz
276   #endif
277        
278 <       if (do_stress) then
278 >    if (do_stress) then
279  
280   #ifdef IS_MPI
281 <          id1 = tagRow(atom1)
282 <          id2 = tagColumn(atom2)
281 >       id1 = tagRow(atom1)
282 >       id2 = tagColumn(atom2)
283   #else
284 <          id1 = atom1
285 <          id2 = atom2
284 >       id1 = atom1
285 >       id2 = atom2
286   #endif
289 –
290 –          if (molMembershipList(id1) .ne. molMembershipList(id2)) then
291 –
292 –             ! because the d vector is the rj - ri vector, and
293 –             ! because fx, fy, fz are the force on atom i, we need a
294 –             ! negative sign here:
295 –
296 –             tau_Temp(1) = tau_Temp(1) - d(1) * fx
297 –             tau_Temp(2) = tau_Temp(2) - d(1) * fy
298 –             tau_Temp(3) = tau_Temp(3) - d(1) * fz
299 –             tau_Temp(4) = tau_Temp(4) - d(2) * fx
300 –             tau_Temp(5) = tau_Temp(5) - d(2) * fy
301 –             tau_Temp(6) = tau_Temp(6) - d(2) * fz
302 –             tau_Temp(7) = tau_Temp(7) - d(3) * fx
303 –             tau_Temp(8) = tau_Temp(8) - d(3) * fy
304 –             tau_Temp(9) = tau_Temp(9) - d(3) * fz
305 –            
306 –             virial_Temp = virial_Temp + &
307 –                  (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
308 –
309 –          endif
310 –       endif
287        
288 +       if (molMembershipList(id1) .ne. molMembershipList(id2)) then
289 +          
290 +          ! because the d vector is the rj - ri vector, and
291 +          ! because fx, fy, fz are the force on atom i, we need a
292 +          ! negative sign here:
293 +          
294 +          tau_Temp(1) = tau_Temp(1) - d(1) * fx
295 +          tau_Temp(2) = tau_Temp(2) - d(1) * fy
296 +          tau_Temp(3) = tau_Temp(3) - d(1) * fz
297 +          tau_Temp(4) = tau_Temp(4) - d(2) * fx
298 +          tau_Temp(5) = tau_Temp(5) - d(2) * fy
299 +          tau_Temp(6) = tau_Temp(6) - d(2) * fz
300 +          tau_Temp(7) = tau_Temp(7) - d(3) * fx
301 +          tau_Temp(8) = tau_Temp(8) - d(3) * fy
302 +          tau_Temp(9) = tau_Temp(9) - d(3) * fz
303 +          
304 +          virial_Temp = virial_Temp + &
305 +               (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
306 +          
307 +       endif
308      endif
309 +        
310      return    
311 <      
311 >    
312    end subroutine do_lj_pair
313 <
314 <
313 >  
314 >  
315    !! Calculates the mixing for sigma or epslon
316 <
316 >  
317    function calcLJMix(thisParam,param1,param2,status) result(myMixParam)
318      character(len=*) :: thisParam
319      real(kind = dp)  :: param1

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines