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 482 by chuckv, Tue Apr 8 22:38:43 2003 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.3 2003-04-08 22:38:43 chuckv Exp $, $Date: 2003-04-08 22:38:43 $, $Name: not supported by cvs2svn $, $Revision: 1.3 $
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
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 115 | Line 132 | contains
132          
133      if (.not. associated(ljMixed)) then
134         allocate(ljMixed(nAtypes, nAtypes))
135 <    else
119 <       status = -1
120 <       return
121 <    end if
135 >    endif
136  
137      rcut6 = LJ_rcut**6
138  
139 + ! This loops through all atypes, even those that don't support LJ forces.
140      do i = 1, nAtypes
141  
142         call getElementProperty(atypes, i, "lj_epsilon", myEpsilon_i)
# Line 129 | Line 144 | contains
144         ! do self mixing rule
145         ljMixed(i,i)%sigma   = mySigma_i
146        
147 <       ljMixed(i,i)%sigma6  = (ljMixed(i,i)%sigma) ** 6                                                                                                  
133 <       ljMixed(i,i)%tp6     = ljMixed(i,i)%sigma6/rcut6
147 >       ljMixed(i,i)%sigma6  = (ljMixed(i,i)%sigma) ** 6
148  
149 +       ljMixed(i,i)%tp6     = (ljMixed(i,i)%sigma6)/rcut6
150 +
151         ljMixed(i,i)%tp12    = (ljMixed(i,i)%tp6) ** 2
152  
153 +
154         ljMixed(i,i)%epsilon = myEpsilon_i
155  
156         ljMixed(i,i)%delta = -4.0_DP * ljMixed(i,i)%epsilon * &
# Line 176 | Line 193 | contains
193      
194    end subroutine createMixingList
195    
196 <  subroutine do_lj_pair(atom1, atom2, d, rij, r2, pot, f, do_pot, do_stress)
196 >  subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, pot, f, &
197 >       do_pot, do_stress)
198  
199      integer, intent(in) ::  atom1, atom2
200      real( kind = dp ), intent(in) :: rij, r2
201 <    real( kind = dp ) :: pot
202 <    real( kind = dp ), dimension(3,getNlocal()) :: f    
201 >    real( kind = dp ) :: pot, sw, vpair
202 >    real( kind = dp ), dimension(3,nLocal) :: f    
203      real( kind = dp ), intent(in), dimension(3) :: d
204      logical, intent(in) :: do_pot, do_stress
205  
# Line 195 | Line 213 | contains
213      real( kind = dp ) :: t6
214      real( kind = dp ) :: t12
215      real( kind = dp ) :: delta
216 +    integer :: id1, id2
217  
218 <
200 <    if (rij.lt.LJ_rcut)  then
201 <
202 <       ! 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) + delta      
235 <      
236 <       dudr = 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
237 <      
238 <       drdx = d(1) / rij
239 <       drdy = d(2) / rij
224 <       drdz = d(3) / rij
225 <      
226 <       fx = dudr * drdx
227 <       fy = dudr * drdy
228 <       fz = dudr * drdz
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 +    dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
242 +      
243 +    drdx = d(1) / rij
244 +    drdy = d(2) / rij
245 +    drdz = d(3) / rij
246 +      
247 +    fx = dudr * drdx
248 +    fy = dudr * drdy
249 +    fz = dudr * drdz
250 +    
251 +      
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 <
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 + pot_temp
267 >    if (do_pot) pot = pot + sw*pot_temp
268  
269 <       f(1,atom1) = f(1,atom1) + fx
270 <       f(2,atom1) = f(2,atom1) + fy
271 <       f(3,atom1) = f(3,atom1) + 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
279  
280 <       f(1,atom2) = f(1,atom2) - fx
281 <       f(2,atom2) = f(2,atom2) - fy
282 <       f(3,atom2) = f(3,atom2) - fz
280 > #ifdef IS_MPI
281 >       id1 = tagRow(atom1)
282 >       id2 = tagColumn(atom2)
283 > #else
284 >       id1 = atom1
285 >       id2 = atom2
286   #endif
287        
288 <       if (do_stress) then
289 <          tau_Temp(1) = tau_Temp(1) + fx * d(1)
290 <          tau_Temp(2) = tau_Temp(2) + fx * d(2)
291 <          tau_Temp(3) = tau_Temp(3) + fx * d(3)
292 <          tau_Temp(4) = tau_Temp(4) + fy * d(1)
293 <          tau_Temp(5) = tau_Temp(5) + fy * d(2)
294 <          tau_Temp(6) = tau_Temp(6) + fy * d(3)
295 <          tau_Temp(7) = tau_Temp(7) + fz * d(1)
296 <          tau_Temp(8) = tau_Temp(8) + fz * d(2)
297 <          tau_Temp(9) = tau_Temp(9) + fz * d(3)
298 <          virial_Temp = virial_Temp + (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
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
268 –      
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
# Line 286 | Line 326 | contains
326      
327      if (present(status)) status = 0
328      select case (LJ_Mixing_Policy)
329 <    case (LB_MIXING_RULE)
329 >    case (1)
330         select case (thisParam)
331         case ("sigma")
332            myMixParam = 0.5_dp * (param1 + param2)

Diff Legend

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