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 2317 by chrisfen, Wed Sep 21 17:20:14 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.15 2005-09-21 17:20:14 chrisfen Exp $, $Date: 2005-09-21 17:20:14 $, $Name: not supported by cvs2svn $, $Revision: 1.15 $
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 58 | Line 60 | module lj
60  
61    implicit none
62    PRIVATE
63 + #define __FORTRAN90
64 + #include "UseTheForce/DarkSide/fInteractionMap.h"
65  
66    integer, parameter :: DP = selected_real_kind(15)
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  
71 –
76    type, private :: LJtype
77       integer       :: atid
78       real(kind=dp) :: sigma
# Line 77 | Line 81 | module lj
81    end type LJtype
82  
83    type, private :: LJList
84 <     integer               :: nLJtypes = 0
84 >     integer               :: Nljtypes = 0
85       integer               :: currentLJtype = 0
86       type(LJtype), pointer :: LJtypes(:)      => null()
87       integer, pointer      :: atidToLJtype(:) => null()
# Line 88 | 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
93 –     real(kind=dp) :: delta
97       logical       :: rCutWasSet = .false.
98       logical       :: shiftedPot
99       logical       :: isSoftCore = .false.
# Line 98 | 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
103 –  public :: setLJUniformCutoff
104 –  public :: setLJCutoffByTypes
111    public :: getSigma
112    public :: getEpsilon
107 –  public :: useGeometricMixing
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 <    call createMixingMap()
167 <  end subroutine setLJDefaultCutoff
168 <
169 <  subroutine setLJUniformCutoff(thisRcut, shiftedPot)
164 <    real(kind=dp), intent(in) :: thisRcut
165 <    logical, intent(in) :: shiftedPot
166 <    integer :: nLJtypes, i, j
167 <
168 <    if (LJMap%currentLJtype == 0) then
169 <       call handleError("LJ", "No members in LJMap")
170 <       return
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 <    nLJtypes = LJMap%nLJtypes
174 <    if (.not. allocated(MixingMap)) then
175 <       allocate(MixingMap(nLJtypes, nLJtypes))
176 <    endif
172 >  end subroutine setLJDefaultCutoff
173  
178 –    do i = 1, nLJtypes
179 –       do j = 1, nLJtypes
180 –          MixingMap(i,j)%rCut = thisRcut
181 –          MixingMap(i,j)%shiftedPot = shiftedPot
182 –          MixingMap(i,j)%rCutWasSet = .true.
183 –       enddo
184 –    enddo
185 –    call createMixingMap()
186 –  end subroutine setLJUniformCutoff
187 –
188 –  subroutine setLJCutoffByTypes(atid1, atid2, thisRcut, shiftedPot)
189 –    integer, intent(in) :: atid1, atid2
190 –    real(kind=dp), intent(in) :: thisRcut
191 –    logical, intent(in) :: shiftedPot
192 –    integer :: nLJtypes, ljt1, ljt2
193 –
194 –    if (LJMap%currentLJtype == 0) then
195 –       call handleError("LJ", "No members in LJMap")
196 –       return
197 –    end if
198 –
199 –    nLJtypes = LJMap%nLJtypes
200 –    if (.not. allocated(MixingMap)) then
201 –       allocate(MixingMap(nLJtypes, nLJtypes))
202 –    endif
203 –
204 –    ljt1 = LJMap%atidToLJtype(atid1)
205 –    ljt2 = LJMap%atidToLJtype(atid2)
206 –
207 –    MixingMap(ljt1,ljt2)%rCut = thisRcut
208 –    MixingMap(ljt1,ljt2)%shiftedPot = shiftedPot
209 –    MixingMap(ljt1,ljt2)%rCutWasSet = .true.
210 –    MixingMap(ljt2,ljt1)%rCut = thisRcut
211 –    MixingMap(ljt2,ljt1)%shiftedPot = shiftedPot
212 –    MixingMap(ljt2,ljt1)%rCutWasSet = .true.
213 –
214 –    call createMixingMap()
215 –  end subroutine setLJCutoffByTypes
216 –
174    function getSigma(atid) result (s)
175      integer, intent(in) :: atid
176      integer :: ljt1
# Line 243 | Line 200 | contains
200      e = LJMap%LJtypes(ljt1)%epsilon
201  
202    end function getEpsilon
246 –
247 –  subroutine useGeometricMixing()
248 –    useGeometricDistanceMixing = .true.
249 –    haveMixingMap = .false.
250 –    return
251 –  end subroutine useGeometricMixing
203  
204    subroutine createMixingMap()
205      integer :: nLJtypes, i, j
# Line 267 | 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 290 | 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
300 <                doShift = defaultShift
301 <             else
302 <                call handleError("LJ", "No specified or default cutoff value!")
303 <             endif
304 <          endif
305 <          
306 <          tp6    = MixingMap(i,j)%sigma6/rcut6
307 <          tp12    = tp6**2          
308 <          MixingMap(i,j)%delta =-4.0_DP*MixingMap(i,j)%epsilon*(tp12 - tp6)
309 <          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
316 –             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 326 | 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 341 | 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
347 –    real( kind = dp ) :: r6
348 –    real( kind = dp ) :: t6
349 –    real( kind = dp ) :: t12
350 –    real( kind = dp ) :: delta
329      logical :: isSoftCore, shiftedPot
330      integer :: id1, id2, localError
331  
# Line 367 | 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
372 –    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  
378 –    t6  = sigma6/ r6
379 –    t12 = t6 * t6    
380 –
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 <      
388 <       vpair = vpair + pot_temp
389 <      
390 <       dudr = -sw * 24.0E0_DP * epsilon * t6 / rij
391 <      
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        
398 –       vpair = vpair + pot_temp
399 –      
400 –       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 410 | Line 389 | contains
389  
390   #ifdef IS_MPI
391      if (do_pot) then
392 <       pot_Row(atom1) = pot_Row(atom1) + sw*pot_temp*0.5
393 <       pot_Col(atom2) = pot_Col(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 470 | 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