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 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.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.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 51 | Line 51 | module lj
51    use vector_class
52    use simulation
53    use status
54 +  use fForceOptions
55 +  use interpolation
56   #ifdef IS_MPI
57    use mpiSimulation
58   #endif
# Line 65 | 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  
73 –
76    type, private :: LJtype
77       integer       :: atid
78       real(kind=dp) :: sigma
# Line 90 | 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
95 –     real(kind=dp) :: delta
97       logical       :: rCutWasSet = .false.
98       logical       :: shiftedPot
99       logical       :: isSoftCore = .false.
# Line 100 | 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
105 –  public :: setLJUniformCutoff
106 –  public :: setLJCutoffByTypes
111    public :: getSigma
112    public :: getEpsilon
109 –  public :: useGeometricMixing
113    public :: do_lj_pair
114    public :: destroyLJtypes
115 +  public :: setLJsplineRmax
116  
117   contains
118  
# Line 159 | 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 <    end if
166 <  end subroutine setLJDefaultCutoff
167 <
168 <  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
169 >       call setLJsplineRmax(defaultCutoff)
170      end if
171  
172 <    nLJtypes = LJMap%nLJtypes
205 <    if (.not. allocated(MixingMap)) then
206 <       allocate(MixingMap(nLJtypes, nLJtypes))
207 <    endif
172 >  end subroutine setLJDefaultCutoff
173  
209 –    ljt1 = LJMap%atidToLJtype(atid1)
210 –    ljt2 = LJMap%atidToLJtype(atid2)
211 –
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 –
174    function getSigma(atid) result (s)
175      integer, intent(in) :: atid
176      integer :: ljt1
# Line 249 | Line 201 | contains
201  
202    end function getEpsilon
203  
252 –  subroutine useGeometricMixing()
253 –    useGeometricDistanceMixing = .true.
254 –    haveMixingMap = .false.
255 –    return
256 –  end subroutine useGeometricMixing
257 –
204    subroutine createMixingMap()
205      integer :: nLJtypes, i, j
206      real ( kind = dp ) :: s1, s2, e1, e2
# Line 272 | Line 218 | contains
218         allocate(MixingMap(nLJtypes, nLJtypes))
219      endif
220  
221 +    useGeometricDistanceMixing = usesGeometricDistanceMixing()
222      do i = 1, nLJtypes
223  
224         s1 = LJMap%LJtypes(i)%sigma
# Line 295 | 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 (MixingMap(i,j)%rCutWasSet) then
248 <             rcut6 = (MixingMap(i,j)%rcut)**6
247 >          if (haveDefaultCutoff) then
248 >             MixingMap(i,j)%shiftedPot = defaultShift
249            else
250 <             if (haveDefaultCutoff) then
251 <                rcut6 = defaultCutoff**6
305 <                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
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
321 –             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 331 | Line 267 | contains
267      
268    end subroutine createMixingMap
269    
270 <  subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
271 <       pot, f, do_pot)
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 +    
313      integer, intent(in) ::  atom1, atom2
314      integer :: atid1, atid2, ljt1, ljt2
315 <    real( kind = dp ), intent(in) :: rij, r2
315 >    real( kind = dp ), intent(in) :: rij, r2, rcut
316      real( kind = dp ) :: pot, sw, vpair
317      real( kind = dp ), dimension(3,nLocal) :: f    
318      real( kind = dp ), intent(in), dimension(3) :: d
# Line 346 | 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
352 –    real( kind = dp ) :: r6
353 –    real( kind = dp ) :: t6
354 –    real( kind = dp ) :: t12
355 –    real( kind = dp ) :: delta
329      logical :: isSoftCore, shiftedPot
330      integer :: id1, id2, localError
331  
# Line 372 | 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
377 –    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  
383 –    t6  = sigma6/ r6
384 –    t12 = t6 * t6    
385 –
356      if (isSoftCore) then
357 <      
358 <       pot_temp = 4.0E0_DP * epsilon * t6
359 <       if (shiftedPot) then
360 <          pot_temp = pot_temp + delta
357 >
358 >       call getSoftFunc(ros, myPot, myDeriv)
359 >
360 >       if (shiftedPot) then
361 >          rcos = rcut * sigmai
362 >          call getSoftFunc(rcos, myPotC, myDerivC)
363         endif
364 <      
393 <       vpair = vpair + pot_temp
394 <      
395 <       dudr = -sw * 24.0E0_DP * epsilon * t6 / rij
396 <      
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 <          pot_temp = pot_temp + delta
370 >          rcos = rcut * sigmai
371 >          call getLJfunc(rcos, myPotC, myDerivC)
372         endif
373        
403 –       vpair = vpair + pot_temp
404 –      
405 –       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 415 | Line 389 | contains
389  
390   #ifdef IS_MPI
391      if (do_pot) then
392 <       pot_Row(LJ_POT,atom1) = pot_Row(LJ_POT,atom1) + sw*pot_temp*0.5
393 <       pot_Col(LJ_POT,atom2) = pot_Col(LJ_POT,atom2) + sw*pot_temp*0.5
392 >       pot_Row(VDW_POT,atom1) = pot_Row(VDW_POT,atom1) + sw*pot_temp*0.5
393 >       pot_Col(VDW_POT,atom2) = pot_Col(VDW_POT,atom2) + sw*pot_temp*0.5
394      endif
395  
396      f_Row(1,atom1) = f_Row(1,atom1) + fx
# Line 475 | 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