ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/LJ.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/LJ.F90 (file contents):
Revision 2355 by chuckv, Wed Oct 12 18:59:16 2005 UTC vs.
Revision 3129 by chrisfen, Fri Apr 20 18:15:48 2007 UTC

# Line 43 | Line 43
43   !! Calculates Long Range forces Lennard-Jones interactions.
44   !! @author Charles F. Vardeman II
45   !! @author Matthew Meineke
46 < !! @version $Id: LJ.F90,v 1.17 2005-10-12 18:59:16 chuckv Exp $, $Date: 2005-10-12 18:59:16 $, $Name: not supported by cvs2svn $, $Revision: 1.17 $
46 > !! @version $Id: LJ.F90,v 1.28 2007-04-20 18:15:47 chrisfen Exp $, $Date: 2007-04-20 18:15:47 $, $Name: not supported by cvs2svn $, $Revision: 1.28 $
47  
48  
49   module lj
50 +  use definitions
51    use atype_module
52    use vector_class
53    use simulation
54    use status
55 +  use fForceOptions
56   #ifdef IS_MPI
57    use mpiSimulation
58   #endif
# Line 61 | Line 63 | module lj
63   #define __FORTRAN90
64   #include "UseTheForce/DarkSide/fInteractionMap.h"
65  
64 –  integer, parameter :: DP = selected_real_kind(15)
65 –
66    logical, save :: useGeometricDistanceMixing = .false.
67    logical, save :: haveMixingMap = .false.
68  
69    real(kind=DP), save :: defaultCutoff = 0.0_DP
70 <  logical, save :: defaultShift = .false.
70 >  logical, save :: defaultShiftPot = .false.
71 >  logical, save :: defaultShiftFrc = .false.
72    logical, save :: haveDefaultCutoff = .false.
73  
73 –
74    type, private :: LJtype
75       integer       :: atid
76       real(kind=dp) :: sigma
# Line 90 | Line 90 | module lj
90    type :: MixParameters
91       real(kind=DP) :: sigma
92       real(kind=DP) :: epsilon
93 <     real(kind=dp) :: sigma6
93 >     real(kind=dp) :: sigmai
94       real(kind=dp) :: rCut
95 –     real(kind=dp) :: delta
95       logical       :: rCutWasSet = .false.
96       logical       :: shiftedPot
97       logical       :: isSoftCore = .false.
98 +     logical       :: shiftedFrc
99    end type MixParameters
100  
101    type(MixParameters), dimension(:,:), allocatable :: MixingMap
102  
103    public :: newLJtype
104    public :: setLJDefaultCutoff
105 –  public :: setLJUniformCutoff
106 –  public :: setLJCutoffByTypes
105    public :: getSigma
106    public :: getEpsilon
109 –  public :: useGeometricMixing
107    public :: do_lj_pair
108    public :: destroyLJtypes
109  
# Line 153 | Line 150 | contains
150      endif
151    end subroutine newLJtype
152  
153 <  subroutine setLJDefaultCutoff(thisRcut, shiftedPot)
153 >  subroutine setLJDefaultCutoff(thisRcut, shiftedPot, shiftedFrc)
154      real(kind=dp), intent(in) :: thisRcut
155      logical, intent(in) :: shiftedPot
156 +    logical, intent(in) :: shiftedFrc
157      defaultCutoff = thisRcut
158 <    defaultShift = shiftedPot
158 >    defaultShiftPot = shiftedPot
159 >    defaultShiftFrc = shiftedFrc
160      haveDefaultCutoff = .true.
161 +
162      !we only want to build LJ Mixing map if LJ is being used.
163      if(LJMap%nLJTypes /= 0) then
164         call createMixingMap()
165      end if
166 –  end subroutine setLJDefaultCutoff
166  
167 <  subroutine setLJUniformCutoff(thisRcut, shiftedPot)
169 <    real(kind=dp), intent(in) :: thisRcut
170 <    logical, intent(in) :: shiftedPot
171 <    integer :: nLJtypes, i, j
172 <
173 <    if (LJMap%currentLJtype == 0) then
174 <       call handleError("LJ", "No members in LJMap")
175 <       return
176 <    end if
177 <
178 <    nLJtypes = LJMap%nLJtypes
179 <    if (.not. allocated(MixingMap)) then
180 <       allocate(MixingMap(nLJtypes, nLJtypes))
181 <    endif
182 <
183 <    do i = 1, nLJtypes
184 <       do j = 1, nLJtypes
185 <          MixingMap(i,j)%rCut = thisRcut
186 <          MixingMap(i,j)%shiftedPot = shiftedPot
187 <          MixingMap(i,j)%rCutWasSet = .true.
188 <       enddo
189 <    enddo
190 <    call createMixingMap()
191 <  end subroutine setLJUniformCutoff
192 <
193 <  subroutine setLJCutoffByTypes(atid1, atid2, thisRcut, shiftedPot)
194 <    integer, intent(in) :: atid1, atid2
195 <    real(kind=dp), intent(in) :: thisRcut
196 <    logical, intent(in) :: shiftedPot
197 <    integer :: nLJtypes, ljt1, ljt2
198 <
199 <    if (LJMap%currentLJtype == 0) then
200 <       call handleError("LJ", "No members in LJMap")
201 <       return
202 <    end if
203 <
204 <    nLJtypes = LJMap%nLJtypes
205 <    if (.not. allocated(MixingMap)) then
206 <       allocate(MixingMap(nLJtypes, nLJtypes))
207 <    endif
208 <
209 <    ljt1 = LJMap%atidToLJtype(atid1)
210 <    ljt2 = LJMap%atidToLJtype(atid2)
167 >  end subroutine setLJDefaultCutoff
168  
212 –    MixingMap(ljt1,ljt2)%rCut = thisRcut
213 –    MixingMap(ljt1,ljt2)%shiftedPot = shiftedPot
214 –    MixingMap(ljt1,ljt2)%rCutWasSet = .true.
215 –    MixingMap(ljt2,ljt1)%rCut = thisRcut
216 –    MixingMap(ljt2,ljt1)%shiftedPot = shiftedPot
217 –    MixingMap(ljt2,ljt1)%rCutWasSet = .true.
218 –
219 –    call createMixingMap()
220 –  end subroutine setLJCutoffByTypes
221 –
169    function getSigma(atid) result (s)
170      integer, intent(in) :: atid
171      integer :: ljt1
# Line 249 | Line 196 | contains
196  
197    end function getEpsilon
198  
252 –  subroutine useGeometricMixing()
253 –    useGeometricDistanceMixing = .true.
254 –    haveMixingMap = .false.
255 –    return
256 –  end subroutine useGeometricMixing
257 –
199    subroutine createMixingMap()
200      integer :: nLJtypes, i, j
201      real ( kind = dp ) :: s1, s2, e1, e2
# Line 272 | Line 213 | contains
213         allocate(MixingMap(nLJtypes, nLJtypes))
214      endif
215  
216 +    useGeometricDistanceMixing = usesGeometricDistanceMixing()
217      do i = 1, nLJtypes
218  
219         s1 = LJMap%LJtypes(i)%sigma
# Line 288 | Line 230 | contains
230  
231            ! only the distance parameter uses different mixing policies
232            if (useGeometricDistanceMixing) then
233 <             MixingMap(i,j)%sigma = dsqrt(s1 * s2)
233 >             MixingMap(i,j)%sigma = sqrt(s1 * s2)
234            else
235               MixingMap(i,j)%sigma = 0.5_dp * (s1 + s2)
236            endif
237            
238 <          MixingMap(i,j)%epsilon = dsqrt(e1 * e2)
238 >          MixingMap(i,j)%epsilon = sqrt(e1 * e2)
239  
240 <          MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
240 >          MixingMap(i,j)%sigmai = 1.0_DP  / (MixingMap(i,j)%sigma)
241  
242 <          if (MixingMap(i,j)%rCutWasSet) then
243 <             rcut6 = (MixingMap(i,j)%rcut)**6
242 >          if (haveDefaultCutoff) then
243 >             MixingMap(i,j)%shiftedPot = defaultShiftPot
244 >             MixingMap(i,j)%shiftedFrc = defaultShiftFrc
245            else
246 <             if (haveDefaultCutoff) then
247 <                rcut6 = defaultCutoff**6
248 <                doShift = defaultShift
306 <             else
307 <                call handleError("LJ", "No specified or default cutoff value!")
308 <             endif
309 <          endif
310 <          
311 <          tp6    = MixingMap(i,j)%sigma6/rcut6
312 <          tp12    = tp6**2          
313 <          MixingMap(i,j)%delta =-4.0_DP*MixingMap(i,j)%epsilon*(tp12 - tp6)
314 <          MixingMap(i,j)%shiftedPot = doShift
246 >             MixingMap(i,j)%shiftedPot = defaultShiftPot
247 >             MixingMap(i,j)%shiftedFrc = defaultShiftFrc
248 >          endif          
249  
250            if (i.ne.j) then
251               MixingMap(j,i)%sigma      = MixingMap(i,j)%sigma
252               MixingMap(j,i)%epsilon    = MixingMap(i,j)%epsilon
253 <             MixingMap(j,i)%sigma6     = MixingMap(i,j)%sigma6
253 >             MixingMap(j,i)%sigmai     = MixingMap(i,j)%sigmai
254               MixingMap(j,i)%rCut       = MixingMap(i,j)%rCut
321 –             MixingMap(j,i)%delta      = MixingMap(i,j)%delta
255               MixingMap(j,i)%rCutWasSet = MixingMap(i,j)%rCutWasSet
256               MixingMap(j,i)%shiftedPot = MixingMap(i,j)%shiftedPot
257               MixingMap(j,i)%isSoftCore = MixingMap(i,j)%isSoftCore
258 +             MixingMap(j,i)%shiftedFrc = MixingMap(i,j)%shiftedFrc
259            endif
260  
261         enddo
262      enddo
263      
264      haveMixingMap = .true.
265 <    
265 >
266    end subroutine createMixingMap
267 <  
268 <  subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
267 >          
268 >  subroutine do_lj_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
269         pot, f, do_pot)
270 <
270 >    
271      integer, intent(in) ::  atom1, atom2
272      integer :: atid1, atid2, ljt1, ljt2
273 <    real( kind = dp ), intent(in) :: rij, r2
273 >    real( kind = dp ), intent(in) :: rij, r2, rcut
274      real( kind = dp ) :: pot, sw, vpair
275      real( kind = dp ), dimension(3,nLocal) :: f    
276      real( kind = dp ), intent(in), dimension(3) :: d
# Line 346 | Line 280 | contains
280      ! local Variables
281      real( kind = dp ) :: drdx, drdy, drdz
282      real( kind = dp ) :: fx, fy, fz
283 +    real( kind = dp ) :: myPot, myPotC, myDeriv, myDerivC, ros, rcos
284      real( kind = dp ) :: pot_temp, dudr
285 <    real( kind = dp ) :: sigma6
285 >    real( kind = dp ) :: sigmai
286      real( kind = dp ) :: epsilon
287 <    real( kind = dp ) :: r6
353 <    real( kind = dp ) :: t6
354 <    real( kind = dp ) :: t12
355 <    real( kind = dp ) :: delta
356 <    logical :: isSoftCore, shiftedPot
287 >    logical :: isSoftCore, shiftedPot, shiftedFrc
288      integer :: id1, id2, localError
289  
290      if (.not.haveMixingMap) then
# Line 372 | Line 303 | contains
303      ljt1 = LJMap%atidToLJtype(atid1)
304      ljt2 = LJMap%atidToLJtype(atid2)
305  
306 <    sigma6     = MixingMap(ljt1,ljt2)%sigma6
306 >    sigmai     = MixingMap(ljt1,ljt2)%sigmai
307      epsilon    = MixingMap(ljt1,ljt2)%epsilon
377 –    delta      = MixingMap(ljt1,ljt2)%delta
308      isSoftCore = MixingMap(ljt1,ljt2)%isSoftCore
309      shiftedPot = MixingMap(ljt1,ljt2)%shiftedPot
310 +    shiftedFrc = MixingMap(ljt1,ljt2)%shiftedFrc
311  
312 <    r6 = r2 * r2 * r2
312 >    ros = rij * sigmai
313  
383 –    t6  = sigma6/ r6
384 –    t12 = t6 * t6    
385 –
314      if (isSoftCore) then
315 <      
316 <       pot_temp = 4.0E0_DP * epsilon * t6
315 >
316 >       call getSoftFunc(ros, myPot, myDeriv)
317 >
318         if (shiftedPot) then
319 <          pot_temp = pot_temp + delta
319 >          rcos = rcut * sigmai
320 >          call getSoftFunc(rcos, myPotC, myDerivC)
321 >          myDerivC = 0.0_dp
322 >       elseif (shiftedFrc) then
323 >          rcos = rcut * sigmai
324 >          call getSoftFunc(rcos, myPotC, myDerivC)
325 >          myPotC = myPotC + myDerivC * (rij - rcut) * sigmai
326 >       else
327 >          myPotC = 0.0_dp
328 >          myDerivC = 0.0_dp
329         endif
330 <      
393 <       vpair = vpair + pot_temp
394 <      
395 <       dudr = -sw * 24.0E0_DP * epsilon * t6 / rij
396 <      
330 >              
331      else
332 <       pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
332 >
333 >       call getLJfunc(ros, myPot, myDeriv)
334 >
335         if (shiftedPot) then
336 <          pot_temp = pot_temp + delta
336 >          rcos = rcut * sigmai
337 >          call getLJfunc(rcos, myPotC, myDerivC)
338 >          myDerivC = 0.0_dp
339 >       elseif (shiftedFrc) then
340 >          rcos = rcut * sigmai
341 >          call getLJfunc(rcos, myPotC, myDerivC)
342 >          myPotC = myPotC + myDerivC * (rij - rcut) * sigmai
343 >       else
344 >          myPotC = 0.0_dp
345 >          myDerivC = 0.0_dp
346         endif
347        
403 –       vpair = vpair + pot_temp
404 –      
405 –       dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
348      endif
349  
350 +    pot_temp = epsilon * (myPot - myPotC)
351 +    vpair = vpair + pot_temp
352 +    dudr = sw * epsilon * (myDeriv - myDerivC) * sigmai
353 +
354      drdx = d(1) / rij
355      drdy = d(2) / rij
356      drdz = d(3) / rij
# Line 415 | Line 361 | contains
361  
362   #ifdef IS_MPI
363      if (do_pot) then
364 <       pot_Row(LJ_POT,atom1) = pot_Row(LJ_POT,atom1) + sw*pot_temp*0.5
365 <       pot_Col(LJ_POT,atom2) = pot_Col(LJ_POT,atom2) + sw*pot_temp*0.5
364 >       pot_Row(VDW_POT,atom1) = pot_Row(VDW_POT,atom1) + sw*pot_temp*0.5
365 >       pot_Col(VDW_POT,atom2) = pot_Col(VDW_POT,atom2) + sw*pot_temp*0.5
366      endif
367  
368      f_Row(1,atom1) = f_Row(1,atom1) + fx
# Line 475 | Line 421 | contains
421      end if
422      
423      haveMixingMap = .false.
424 +
425    end subroutine destroyLJTypes
426  
427 +  subroutine getLJfunc(r, myPot, myDeriv)
428 +
429 +    real(kind=dp), intent(in) :: r
430 +    real(kind=dp), intent(inout) :: myPot, myDeriv
431 +    real(kind=dp) :: ri, ri2, ri6, ri7, ri12, ri13
432 +    real(kind=dp) :: a, b, c, d, dx
433 +    integer :: j
434 +
435 +    ri = 1.0_DP / r
436 +    ri2 = ri*ri
437 +    ri6 = ri2*ri2*ri2
438 +    ri7 = ri6*ri
439 +    ri12 = ri6*ri6
440 +    ri13 = ri12*ri
441 +    
442 +    myPot = 4.0_DP * (ri12 - ri6)
443 +    myDeriv = 24.0_DP * (ri7 - 2.0_DP * ri13)
444 +    
445 +    return
446 +  end subroutine getLJfunc
447 +
448 +  subroutine getSoftFunc(r, myPot, myDeriv)
449 +    
450 +    real(kind=dp), intent(in) :: r
451 +    real(kind=dp), intent(inout) :: myPot, myDeriv
452 +    real(kind=dp) :: ri, ri2, ri6, ri7
453 +    real(kind=dp) :: a, b, c, d, dx
454 +    integer :: j
455 +    
456 +    ri = 1.0_DP / r    
457 +    ri2 = ri*ri
458 +    ri6 = ri2*ri2*ri2
459 +    ri7 = ri6*ri
460 +    myPot = 4.0_DP * (ri6)
461 +    myDeriv = - 24.0_DP * ri7
462 +    
463 +    return
464 +  end subroutine getSoftFunc
465 +
466   end module lj

Diff Legend

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