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 2722 by gezelter, Thu Apr 20 18:24:24 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.22 2006-04-20 18:24:24 gezelter Exp $, $Date: 2006-04-20 18:24:24 $, $Name: not supported by cvs2svn $, $Revision: 1.22 $
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 +       call setLJsplineRmax(defaultCutoff)
170      end if
171 +
172    end subroutine setLJDefaultCutoff
173  
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
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)
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
357 >
358 >       call getSoftFunc(ros, myPot, myDeriv)
359 >
360         if (shiftedPot) then
361 <          pot_temp = pot_temp + delta
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 +
470 +       call lookupUniformSpline1d(vLJSpline, r, myPot, myDeriv)
471 +      
472 +    else
473 +       ri = 1.0_DP / r
474 +       ri2 = ri*ri
475 +       ri6 = ri2*ri2*ri2
476 +       ri7 = ri6*ri
477 +       ri12 = ri6*ri6
478 +       ri13 = ri12*ri
479 +      
480 +       myPot = 4.0_DP * (ri12 - ri6)
481 +       myDeriv = 24.0_DP * (ri7 - 2.0_DP * ri13)
482 +    endif
483 +
484 +    return
485 +  end subroutine getLJfunc
486 +
487 +  subroutine getSoftFunc(r, myPot, myDeriv)
488 +    
489 +    real(kind=dp), intent(in) :: r
490 +    real(kind=dp), intent(inout) :: myPot, myDeriv
491 +    real(kind=dp) :: ri, ri2, ri6, ri7
492 +    real(kind=dp) :: a, b, c, d, dx
493 +    integer :: j
494 +    
495 +    if (useSplines) then
496 +
497 +       call lookupUniformSpline1d(vSoftSpline, r, myPot, myDeriv)
498 +      
499 +    else
500 +       ri = 1.0_DP / r    
501 +       ri2 = ri*ri
502 +       ri6 = ri2*ri2*ri2
503 +       ri7 = ri6*ri
504 +       myPot = 4.0_DP * (ri6)
505 +       myDeriv = - 24.0_DP * ri7
506 +    endif
507 +
508 +    return
509 +  end subroutine getSoftFunc
510 +
511 +  subroutine setupSplines(rmin, rmax, np)
512 +    real( kind = dp ), intent(in) :: rmin, rmax
513 +    integer, intent(in) :: np
514 +    real( kind = dp ) :: rvals(np), vLJ(np), vLJp(np), vSoft(np), vSoftp(np)
515 +    real( kind = dp ) :: dr, r, ri, ri2, ri6, ri7, ri12, ri13
516 +    integer :: i
517 +
518 +    dr = (rmax-rmin) / float(np-1)
519 +    
520 +    do i = 1, np
521 +       r = rmin + dble(i-1)*dr
522 +       ri = 1.0_DP / r
523 +       ri2 = ri*ri
524 +       ri6 = ri2*ri2*ri2
525 +       ri7 = ri6*ri
526 +       ri12 = ri6*ri6
527 +       ri13 = ri12*ri
528 +
529 +       rvals(i) = r
530 +       vLJ(i) = 4.0_DP * (ri12 - ri6)
531 +       vLJp(i) = 24.0_DP * (ri7 - 2.0_DP * ri13)
532 +
533 +       vSoft(i) = 4.0_DP * (ri6)
534 +       vSoftp(i) = - 24.0_DP * ri7      
535 +    enddo
536 +
537 +    call newSpline(vLJspline, rvals, vLJ,  .true.)
538 +    call newSpline(vLJpspline, rvals, vLJp, .true.)
539 +    call newSpline(vSoftSpline, rvals, vSoft,  .true.)
540 +    call newSpline(vSoftpSpline, rvals, vSoftp, .true.)
541 +
542 +    return
543 +  end subroutine setupSplines
544   end module lj

Diff Legend

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