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 490 by gezelter, Fri Apr 11 15:16:59 2003 UTC vs.
Revision 1192 by gezelter, Mon May 24 21:03:30 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.5 2003-04-11 15:16:59 gezelter Exp $, $Date: 2003-04-11 15:16:59 $, $Name: not supported by cvs2svn $, $Revision: 1.5 $
5 > !! @version $Id: calc_LJ_FF.F90,v 1.20 2004-05-24 21:03:24 gezelter Exp $, $Date: 2004-05-24 21:03:24 $, $Name: not supported by cvs2svn $, $Revision: 1.20 $
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, fpair, &
197 >       pot, f, do_pot)
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
204 >    real( kind = dp ), intent(inout), dimension(3) :: fpair
205 >    logical, intent(in) :: do_pot
206  
207      ! local Variables
208      real( kind = dp ) :: drdx, drdy, drdz
# Line 197 | Line 216 | contains
216      real( kind = dp ) :: delta
217      integer :: id1, id2
218  
219 <
201 <    if (rij.lt.LJ_rcut)  then
202 <
203 <       ! Look up the correct parameters in the mixing matrix
219 >    ! Look up the correct parameters in the mixing matrix
220   #ifdef IS_MPI
221 <       sigma6   = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6
222 <       epsilon  = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon
223 <       delta    = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta
221 >    sigma6   = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6
222 >    epsilon  = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon
223 >    delta    = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta
224   #else
225 <       sigma6   = ljMixed(atid(atom1),atid(atom2))%sigma6
226 <       epsilon  = ljMixed(atid(atom1),atid(atom2))%epsilon
227 <       delta    = ljMixed(atid(atom1),atid(atom2))%delta
225 >    sigma6   = ljMixed(atid(atom1),atid(atom2))%sigma6
226 >    epsilon  = ljMixed(atid(atom1),atid(atom2))%epsilon
227 >    delta    = ljMixed(atid(atom1),atid(atom2))%delta
228   #endif
229 +
230 +    r6 = r2 * r2 * r2
231 +    
232 +    t6  = sigma6/ r6
233 +    t12 = t6 * t6    
234 +  
235 +    pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
236 +    if (LJ_do_shift) then
237 +       pot_temp = pot_temp + delta
238 +    endif
239 +
240 +    vpair = vpair + pot_temp
241        
242 <       r6 = r2 * r2 * r2
242 >    dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
243        
244 <       t6  = sigma6/ r6
245 <       t12 = t6 * t6    
246 <            
219 <       pot_temp = 4.0E0_DP * epsilon * (t12 - t6) + delta      
244 >    drdx = d(1) / rij
245 >    drdy = d(2) / rij
246 >    drdz = d(3) / rij
247        
248 <       dudr = 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
248 >    fx = dudr * drdx
249 >    fy = dudr * drdy
250 >    fz = dudr * drdz
251 >    
252        
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 –      
253   #ifdef IS_MPI
254 <       if (do_pot) then
255 <          pot_Row(atom1) = pot_Row(atom1) + pot_temp*0.5
256 <          pot_Col(atom2) = pot_Col(atom2) + pot_temp*0.5
257 <       endif
258 <
259 <       f_Row(1,atom1) = f_Row(1,atom1) + fx
260 <       f_Row(2,atom1) = f_Row(2,atom1) + fy
261 <       f_Row(3,atom1) = f_Row(3,atom1) + fz
262 <
263 <       f_Col(1,atom2) = f_Col(1,atom2) - fx
264 <       f_Col(2,atom2) = f_Col(2,atom2) - fy
265 <       f_Col(3,atom2) = f_Col(3,atom2) - fz      
266 <
254 >    if (do_pot) then
255 >       pot_Row(atom1) = pot_Row(atom1) + sw*pot_temp*0.5
256 >       pot_Col(atom2) = pot_Col(atom2) + sw*pot_temp*0.5
257 >    endif
258 >    
259 >    f_Row(1,atom1) = f_Row(1,atom1) + fx
260 >    f_Row(2,atom1) = f_Row(2,atom1) + fy
261 >    f_Row(3,atom1) = f_Row(3,atom1) + fz
262 >    
263 >    f_Col(1,atom2) = f_Col(1,atom2) - fx
264 >    f_Col(2,atom2) = f_Col(2,atom2) - fy
265 >    f_Col(3,atom2) = f_Col(3,atom2) - fz      
266 >    
267   #else
268 <       if (do_pot) pot = pot + pot_temp
268 >    if (do_pot) pot = pot + sw*pot_temp
269  
270 <       f(1,atom1) = f(1,atom1) + fx
271 <       f(2,atom1) = f(2,atom1) + fy
272 <       f(3,atom1) = f(3,atom1) + fz
273 <
274 <       f(1,atom2) = f(1,atom2) - fx
275 <       f(2,atom2) = f(2,atom2) - fy
276 <       f(3,atom2) = f(3,atom2) - fz
270 >    f(1,atom1) = f(1,atom1) + fx
271 >    f(2,atom1) = f(2,atom1) + fy
272 >    f(3,atom1) = f(3,atom1) + fz
273 >    
274 >    f(1,atom2) = f(1,atom2) - fx
275 >    f(2,atom2) = f(2,atom2) - fy
276 >    f(3,atom2) = f(3,atom2) - fz
277   #endif
278 <      
257 <       if (do_stress) then
258 <
278 >        
279   #ifdef IS_MPI
280 <          id1 = tagRow(atom1)
281 <          id2 = tagColumn(atom2)
280 >    id1 = tagRow(atom1)
281 >    id2 = tagColumn(atom2)
282   #else
283 <          id1 = atom1
284 <          id2 = atom2
283 >    id1 = atom1
284 >    id2 = atom2
285   #endif
286  
287 <          if (molMembershipList(id1) .ne. molMembershipList(id2)) then
268 <             tau_Temp(1) = tau_Temp(1) + fx * d(1)
269 <             tau_Temp(2) = tau_Temp(2) + fx * d(2)
270 <             tau_Temp(3) = tau_Temp(3) + fx * d(3)
271 <             tau_Temp(4) = tau_Temp(4) + fy * d(1)
272 <             tau_Temp(5) = tau_Temp(5) + fy * d(2)
273 <             tau_Temp(6) = tau_Temp(6) + fy * d(3)
274 <             tau_Temp(7) = tau_Temp(7) + fz * d(1)
275 <             tau_Temp(8) = tau_Temp(8) + fz * d(2)
276 <             tau_Temp(9) = tau_Temp(9) + fz * d(3)
277 <             virial_Temp = virial_Temp + &
278 <                  (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
279 <          else
280 <             write(0,*) 'skipping ', id1, id2, molMembershipList(id1), molMembershipList(id2)
281 <          endif
282 <
283 <       endif
287 >    if (molMembershipList(id1) .ne. molMembershipList(id2)) then
288        
289 +       fpair(1) = fpair(1) + fx
290 +       fpair(2) = fpair(2) + fy
291 +       fpair(3) = fpair(3) + fz
292 +
293      endif
294 +
295      return    
296 <      
296 >    
297    end subroutine do_lj_pair
298 <
299 <
298 >  
299 >  
300    !! Calculates the mixing for sigma or epslon
301 <
301 >  
302    function calcLJMix(thisParam,param1,param2,status) result(myMixParam)
303      character(len=*) :: thisParam
304      real(kind = dp)  :: param1
# Line 302 | Line 311 | contains
311      
312      if (present(status)) status = 0
313      select case (LJ_Mixing_Policy)
314 <    case (LB_MIXING_RULE)
314 >    case (1)
315         select case (thisParam)
316         case ("sigma")
317            myMixParam = 0.5_dp * (param1 + param2)

Diff Legend

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