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 2716 by chuckv, Wed Dec 7 19:58:18 2005 UTC vs.
Revision 2717 by gezelter, Mon Apr 17 21:49:12 2006 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.20 2005-12-07 19:58:18 chuckv Exp $, $Date: 2005-12-07 19:58:18 $, $Name: not supported by cvs2svn $, $Revision: 1.20 $
46 > !! @version $Id: LJ.F90,v 1.21 2006-04-17 21:49:12 gezelter Exp $, $Date: 2006-04-17 21:49:12 $, $Name: not supported by cvs2svn $, $Revision: 1.21 $
47  
48  
49   module lj
# Line 52 | Line 52 | module lj
52    use simulation
53    use status
54    use fForceOptions
55 +  use interpolation
56   #ifdef IS_MPI
57    use mpiSimulation
58   #endif
# Line 66 | Line 67 | module lj
67  
68    logical, save :: useGeometricDistanceMixing = .false.
69    logical, save :: haveMixingMap = .false.
70 +  logical, save :: useSplines = .false.
71  
72    real(kind=DP), save :: defaultCutoff = 0.0_DP
73    logical, save :: defaultShift = .false.
74    logical, save :: haveDefaultCutoff = .false.
75  
74 –
76    type, private :: LJtype
77       integer       :: atid
78       real(kind=dp) :: sigma
# Line 91 | Line 92 | module lj
92    type :: MixParameters
93       real(kind=DP) :: sigma
94       real(kind=DP) :: epsilon
95 <     real(kind=dp) :: sigma6
95 >     real(kind=dp) :: sigmai
96       real(kind=dp) :: rCut
96 –     real(kind=dp) :: delta
97       logical       :: rCutWasSet = .false.
98       logical       :: shiftedPot
99       logical       :: isSoftCore = .false.
# Line 101 | Line 101 | module lj
101  
102    type(MixParameters), dimension(:,:), allocatable :: MixingMap
103  
104 +  type(cubicSpline), save :: vLJspline
105 +  type(cubicSpline), save :: vLJpspline
106 +  type(cubicSpline), save :: vSoftSpline
107 +  type(cubicSpline), save :: vSoftpSpline
108 +
109    public :: newLJtype
110    public :: setLJDefaultCutoff
111    public :: getSigma
112    public :: getEpsilon
113    public :: do_lj_pair
114    public :: destroyLJtypes
115 +  public :: setLJsplineRmax
116  
117   contains
118  
# Line 157 | Line 163 | contains
163      defaultCutoff = thisRcut
164      defaultShift = shiftedPot
165      haveDefaultCutoff = .true.
166 <    !we only want to build LJ Mixing map if LJ is being used.
166 >    !we only want to build LJ Mixing map and spline if LJ is being used.
167      if(LJMap%nLJTypes /= 0) then
168         call createMixingMap()
169 +       call setLJsplineRmax(defaultCutoff)
170      end if
171 +
172    end subroutine setLJDefaultCutoff
173  
174    function getSigma(atid) result (s)
# Line 234 | Line 242 | contains
242            
243            MixingMap(i,j)%epsilon = dsqrt(e1 * e2)
244  
245 <          MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
245 >          MixingMap(i,j)%sigmai = 1.0_DP  / (MixingMap(i,j)%sigma)
246  
247            if (haveDefaultCutoff) then
240 –             rcut6 = defaultCutoff**6
241 –             tp6    = MixingMap(i,j)%sigma6/rcut6
242 –             tp12    = tp6**2          
243 –             MixingMap(i,j)%delta =-4.0_DP*MixingMap(i,j)%epsilon*(tp12 - tp6)
248               MixingMap(i,j)%shiftedPot = defaultShift
249            else
246 –             MixingMap(i,j)%delta = 0.0_DP
250               MixingMap(i,j)%shiftedPot = defaultShift
251            endif          
252  
253            if (i.ne.j) then
254               MixingMap(j,i)%sigma      = MixingMap(i,j)%sigma
255               MixingMap(j,i)%epsilon    = MixingMap(i,j)%epsilon
256 <             MixingMap(j,i)%sigma6     = MixingMap(i,j)%sigma6
256 >             MixingMap(j,i)%sigmai     = MixingMap(i,j)%sigmai
257               MixingMap(j,i)%rCut       = MixingMap(i,j)%rCut
255 –             MixingMap(j,i)%delta      = MixingMap(i,j)%delta
258               MixingMap(j,i)%rCutWasSet = MixingMap(i,j)%rCutWasSet
259               MixingMap(j,i)%shiftedPot = MixingMap(i,j)%shiftedPot
260               MixingMap(j,i)%isSoftCore = MixingMap(i,j)%isSoftCore
# Line 265 | Line 267 | contains
267      
268    end subroutine createMixingMap
269    
270 +  subroutine setLJsplineRmax(largestRcut)
271 +    real( kind = dp ), intent(in) :: largestRcut
272 +    real( kind = dp ) :: s, bigS, smallS, rmax, rmin
273 +    integer :: np, i
274 +
275 +    if (LJMap%nLJtypes .ne. 0) then
276 +
277 +       !
278 +       ! find the largest and smallest values of sigma that we'll need
279 +       !
280 +       bigS = 0.0_DP
281 +       smallS = 1.0e9
282 +       do i = 1, LJMap%nLJtypes
283 +          s = LJMap%LJtypes(i)%sigma
284 +          if (s .gt. bigS) bigS = s
285 +          if (s .lt. smallS) smallS = s
286 +       enddo
287 +      
288 +       !
289 +       ! give ourselves a 20% margin just in case
290 +       !
291 +       rmax = 1.2 * largestRcut / smallS    
292 +       !
293 +       ! assume atoms will never get closer than 1 angstrom
294 +       !
295 +       rmin = 1 / bigS
296 +       !
297 +       ! assume 500 points is enough
298 +       !
299 +       np = 500
300 +      
301 +       write(*,*) 'calling setupSplines with rmin = ', rmin, ' rmax = ', rmax, &
302 +            ' np = ', np
303 +      
304 +       call setupSplines(rmin, rmax, np)
305 +
306 +    endif
307 +    return
308 +  end subroutine setLJsplineRmax
309 +        
310    subroutine do_lj_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
311         pot, f, do_pot)
312      
# Line 280 | Line 322 | contains
322      ! local Variables
323      real( kind = dp ) :: drdx, drdy, drdz
324      real( kind = dp ) :: fx, fy, fz
325 +    real( kind = dp ) :: myPot, myPotC, myDeriv, myDerivC, ros, rcos
326      real( kind = dp ) :: pot_temp, dudr
327 <    real( kind = dp ) :: sigma6
327 >    real( kind = dp ) :: sigmai
328      real( kind = dp ) :: epsilon
286 –    real( kind = dp ) :: r6, rc6
287 –    real( kind = dp ) :: t6, tp6
288 –    real( kind = dp ) :: t12, tp12
289 –    real( kind = dp ) :: delta
329      logical :: isSoftCore, shiftedPot
330      integer :: id1, id2, localError
331  
# Line 306 | Line 345 | contains
345      ljt1 = LJMap%atidToLJtype(atid1)
346      ljt2 = LJMap%atidToLJtype(atid2)
347  
348 <    sigma6     = MixingMap(ljt1,ljt2)%sigma6
348 >    sigmai     = MixingMap(ljt1,ljt2)%sigmai
349      epsilon    = MixingMap(ljt1,ljt2)%epsilon
311 –    delta      = MixingMap(ljt1,ljt2)%delta
350      isSoftCore = MixingMap(ljt1,ljt2)%isSoftCore
351      shiftedPot = MixingMap(ljt1,ljt2)%shiftedPot
352  
353 <    r6 = r2 * r2 * r2
353 >    ros = rij * sigmai
354 >    myPotC = 0.0_DP
355  
317 –    t6  = sigma6/ r6
318 –    t12 = t6 * t6    
319 –
356      if (isSoftCore) then
357 <      
358 <       pot_temp = 4.0E0_DP * epsilon * t6
357 >
358 >       call getSoftFunc(ros, myPot, myDeriv)
359 >
360         if (shiftedPot) then
361 <          rc6 = rcut**6
362 <          tp6 = sigma6 / rc6
326 <          delta =-4.0_DP*epsilon*(tp6)
327 <          pot_temp = pot_temp + delta
361 >          rcos = rcut * sigmai
362 >          call getSoftFunc(rcos, myPotC, myDerivC)
363         endif
364 <      
330 <       vpair = vpair + pot_temp
331 <      
332 <       dudr = -sw * 24.0E0_DP * epsilon * t6 / rij
333 <      
364 >              
365      else
366 <       pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
366 >
367 >       call getLJfunc(ros, myPot, myDeriv)
368 >
369         if (shiftedPot) then
370 <          rc6 = rcut**6
371 <          tp6 = sigma6 / rc6
339 <          tp12 = tp6*tp6
340 <          delta =-4.0_DP*epsilon*(tp12 - tp6)
341 <          pot_temp = pot_temp + delta
370 >          rcos = rcut * sigmai
371 >          call getLJfunc(rcos, myPotC, myDerivC)
372         endif
373        
344 –       vpair = vpair + pot_temp
345 –      
346 –       dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
374      endif
375  
376 +    !write(*,*) rij, ros, rcos, myPot, myDeriv, myPotC
377 +
378 +    pot_temp = epsilon * (myPot - myPotC)
379 +    vpair = vpair + pot_temp
380 +    dudr = sw * epsilon * myDeriv * sigmai
381 +
382      drdx = d(1) / rij
383      drdy = d(2) / rij
384      drdz = d(3) / rij
# Line 416 | Line 449 | contains
449      end if
450      
451      haveMixingMap = .false.
452 +
453 +    call deleteSpline(vLJspline)
454 +    call deleteSpline(vLJpspline)
455 +    call deleteSpline(vSoftSpline)
456 +    call deleteSpline(vSoftpSpline)
457 +
458    end subroutine destroyLJTypes
459 +
460 +  subroutine getLJfunc(r, myPot, myDeriv)
461  
462 +    real(kind=dp), intent(in) :: r
463 +    real(kind=dp), intent(inout) :: myPot, myDeriv
464 +    real(kind=dp) :: ri, ri2, ri6, ri7, ri12, ri13
465 +    real(kind=dp) :: a, b, c, d, dx
466 +    integer :: j
467 +
468 +    if (useSplines) then
469 +       j = MAX(1, MIN(vLJspline%np, idint((r-vLJspline%x(1)) * vLJspline%dx_i) + 1))
470 +      
471 +       dx = r - vLJspline%x(j)
472 +      
473 +       a = vLJspline%c(1,j)
474 +       b = vLJspline%c(2,j)
475 +       c = vLJspline%c(3,j)
476 +       d = vLJspline%c(4,j)
477 +      
478 +       myPot = c + dx * d
479 +       myPot = b + dx * myPot
480 +       myPot = a + dx * myPot
481 +
482 +       a = vLJpspline%c(1,j)
483 +       b = vLJpspline%c(2,j)
484 +       c = vLJpspline%c(3,j)
485 +       d = vLJpspline%c(4,j)
486 +      
487 +       myDeriv = c + dx * d
488 +       myDeriv = b + dx * myDeriv  
489 +       myDeriv = a + dx * myDeriv
490 +      
491 +    else
492 +       ri = 1.0_DP / r
493 +       ri2 = ri*ri
494 +       ri6 = ri2*ri2*ri2
495 +       ri7 = ri6*ri
496 +       ri12 = ri6*ri6
497 +       ri13 = ri12*ri
498 +      
499 +       myPot = 4.0_DP * (ri12 - ri6)
500 +       myDeriv = 24.0_DP * (ri7 - 2.0_DP * ri13)
501 +    endif
502 +
503 +    return
504 +  end subroutine getLJfunc
505 +
506 +  subroutine getSoftFunc(r, myPot, myDeriv)
507 +    
508 +    real(kind=dp), intent(in) :: r
509 +    real(kind=dp), intent(inout) :: myPot, myDeriv
510 +    real(kind=dp) :: ri, ri2, ri6, ri7
511 +    real(kind=dp) :: a, b, c, d, dx
512 +    integer :: j
513 +    
514 +    if (useSplines) then
515 +       j = MAX(1, MIN(vSoftSpline%np, idint((r-vSoftSpline%x(1)) * vSoftSpline%dx_i) + 1))
516 +      
517 +       dx = r - vSoftSpline%x(j)
518 +      
519 +       a = vSoftSpline%c(1,j)
520 +       b = vSoftSpline%c(2,j)
521 +       c = vSoftSpline%c(3,j)
522 +       d = vSoftSpline%c(4,j)
523 +      
524 +       myPot = c + dx * d
525 +       myPot = b + dx * myPot
526 +       myPot = a + dx * myPot
527 +
528 +       a = vSoftPspline%c(1,j)
529 +       b = vSoftPspline%c(2,j)
530 +       c = vSoftPspline%c(3,j)
531 +       d = vSoftPspline%c(4,j)
532 +      
533 +       myDeriv = c + dx * d
534 +       myDeriv = b + dx * myDeriv  
535 +       myDeriv = a + dx * myDeriv
536 +      
537 +    else
538 +       ri = 1.0_DP / r    
539 +       ri2 = ri*ri
540 +       ri6 = ri2*ri2*ri2
541 +       ri7 = ri6*ri
542 +       myPot = 4.0_DP * (ri6)
543 +       myDeriv = - 24.0_DP * ri7
544 +    endif
545 +
546 +    return
547 +  end subroutine getSoftFunc
548 +
549 +  subroutine setupSplines(rmin, rmax, np)
550 +    real( kind = dp ), intent(in) :: rmin, rmax
551 +    integer, intent(in) :: np
552 +    real( kind = dp ) :: rvals(np), vLJ(np), vLJp(np), vSoft(np), vSoftp(np)
553 +    real( kind = dp ) :: dr, r, ri, ri2, ri6, ri7, ri12, ri13
554 +    real( kind = dp ) :: vljpp1, vljppn, vsoftpp1, vsoftppn
555 +    integer :: i
556 +
557 +    dr = (rmax-rmin) / float(np-1)
558 +    
559 +    do i = 1, np
560 +       r = rmin + dble(i-1)*dr
561 +       ri = 1.0_DP / r
562 +       ri2 = ri*ri
563 +       ri6 = ri2*ri2*ri2
564 +       ri7 = ri6*ri
565 +       ri12 = ri6*ri6
566 +       ri13 = ri12*ri
567 +
568 +       rvals(i) = r
569 +       vLJ(i) = 4.0_DP * (ri12 - ri6)
570 +       vLJp(i) = 24.0_DP * (ri7 - 2.0_DP * ri13)
571 +
572 +       vSoft(i) = 4.0_DP * (ri6)
573 +       vSoftp(i) = - 24.0_DP * ri7      
574 +    enddo
575 +
576 +    vljpp1 = 624.0_DP / (rmin)**(14)  - 168.0_DP / (rmin)**(8)
577 +    vljppn = 624.0_DP / (rmax)**(14)  - 168.0_DP / (rmax)**(8)
578 +
579 +    vsoftpp1 = 168.0_DP / (rmin)**(8)
580 +    vsoftppn = 168.0_DP / (rmax)**(8)
581 +
582 +    call newSpline(vLJspline, rvals, vLJ, vLJp(1), vLJp(np), .true.)
583 +    call newSpline(vLJpspline, rvals, vLJp, vljpp1, vljppn, .true.)
584 +    call newSpline(vSoftSpline, rvals, vSoft, vSoftp(1), vSoftp(np), .true.)
585 +    call newSpline(vSoftpSpline, rvals, vSoftp, vsoftpp1, vsoftppn, .true.)
586 +
587 +    return
588 +  end subroutine setupSplines
589   end module lj

Diff Legend

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