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 2733 by gezelter, Tue Apr 25 02:09:01 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.24 2006-04-25 02:09:01 gezelter Exp $, $Date: 2006-04-25 02:09:01 $, $Name: not supported by cvs2svn $, $Revision: 1.24 $
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   #ifdef IS_MPI
56    use mpiSimulation
57   #endif
# Line 58 | Line 59 | module lj
59  
60    implicit none
61    PRIVATE
62 + #define __FORTRAN90
63 + #include "UseTheForce/DarkSide/fInteractionMap.h"
64  
65    integer, parameter :: DP = selected_real_kind(15)
66  
# Line 68 | Line 71 | module lj
71    logical, save :: defaultShift = .false.
72    logical, save :: haveDefaultCutoff = .false.
73  
71 –
74    type, private :: LJtype
75       integer       :: atid
76       real(kind=dp) :: sigma
# Line 77 | Line 79 | module lj
79    end type LJtype
80  
81    type, private :: LJList
82 <     integer               :: nLJtypes = 0
82 >     integer               :: Nljtypes = 0
83       integer               :: currentLJtype = 0
84       type(LJtype), pointer :: LJtypes(:)      => null()
85       integer, pointer      :: atidToLJtype(:) => null()
# Line 88 | 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
93 –     real(kind=dp) :: delta
95       logical       :: rCutWasSet = .false.
96       logical       :: shiftedPot
97       logical       :: isSoftCore = .false.
# Line 100 | Line 101 | module lj
101  
102    public :: newLJtype
103    public :: setLJDefaultCutoff
103 –  public :: setLJUniformCutoff
104 –  public :: setLJCutoffByTypes
104    public :: getSigma
105    public :: getEpsilon
107 –  public :: useGeometricMixing
106    public :: do_lj_pair
107    public :: destroyLJtypes
108  
# Line 157 | Line 155 | contains
155      defaultCutoff = thisRcut
156      defaultShift = shiftedPot
157      haveDefaultCutoff = .true.
158 <    call createMixingMap()
159 <  end subroutine setLJDefaultCutoff
160 <
163 <  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
158 >    !we only want to build LJ Mixing map if LJ is being used.
159 >    if(LJMap%nLJTypes /= 0) then
160 >       call createMixingMap()
161      end if
162  
163 <    nLJtypes = LJMap%nLJtypes
174 <    if (.not. allocated(MixingMap)) then
175 <       allocate(MixingMap(nLJtypes, nLJtypes))
176 <    endif
177 <
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
163 >  end subroutine setLJDefaultCutoff
164  
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 –
165    function getSigma(atid) result (s)
166      integer, intent(in) :: atid
167      integer :: ljt1
# Line 244 | Line 192 | contains
192  
193    end function getEpsilon
194  
247 –  subroutine useGeometricMixing()
248 –    useGeometricDistanceMixing = .true.
249 –    haveMixingMap = .false.
250 –    return
251 –  end subroutine useGeometricMixing
252 –
195    subroutine createMixingMap()
196      integer :: nLJtypes, i, j
197      real ( kind = dp ) :: s1, s2, e1, e2
# Line 267 | Line 209 | contains
209         allocate(MixingMap(nLJtypes, nLJtypes))
210      endif
211  
212 +    useGeometricDistanceMixing = usesGeometricDistanceMixing()
213      do i = 1, nLJtypes
214  
215         s1 = LJMap%LJtypes(i)%sigma
# Line 290 | Line 233 | contains
233            
234            MixingMap(i,j)%epsilon = dsqrt(e1 * e2)
235  
236 <          MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
236 >          MixingMap(i,j)%sigmai = 1.0_DP  / (MixingMap(i,j)%sigma)
237  
238 <          if (MixingMap(i,j)%rCutWasSet) then
239 <             rcut6 = (MixingMap(i,j)%rcut)**6
238 >          if (haveDefaultCutoff) then
239 >             MixingMap(i,j)%shiftedPot = defaultShift
240            else
241 <             if (haveDefaultCutoff) then
242 <                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
241 >             MixingMap(i,j)%shiftedPot = defaultShift
242 >          endif          
243  
244            if (i.ne.j) then
245               MixingMap(j,i)%sigma      = MixingMap(i,j)%sigma
246               MixingMap(j,i)%epsilon    = MixingMap(i,j)%epsilon
247 <             MixingMap(j,i)%sigma6     = MixingMap(i,j)%sigma6
247 >             MixingMap(j,i)%sigmai     = MixingMap(i,j)%sigmai
248               MixingMap(j,i)%rCut       = MixingMap(i,j)%rCut
316 –             MixingMap(j,i)%delta      = MixingMap(i,j)%delta
249               MixingMap(j,i)%rCutWasSet = MixingMap(i,j)%rCutWasSet
250               MixingMap(j,i)%shiftedPot = MixingMap(i,j)%shiftedPot
251               MixingMap(j,i)%isSoftCore = MixingMap(i,j)%isSoftCore
# Line 325 | Line 257 | contains
257      haveMixingMap = .true.
258      
259    end subroutine createMixingMap
260 <  
261 <  subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
260 >          
261 >  subroutine do_lj_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
262         pot, f, do_pot)
263 <
263 >    
264      integer, intent(in) ::  atom1, atom2
265      integer :: atid1, atid2, ljt1, ljt2
266 <    real( kind = dp ), intent(in) :: rij, r2
266 >    real( kind = dp ), intent(in) :: rij, r2, rcut
267      real( kind = dp ) :: pot, sw, vpair
268      real( kind = dp ), dimension(3,nLocal) :: f    
269      real( kind = dp ), intent(in), dimension(3) :: d
# Line 341 | Line 273 | contains
273      ! local Variables
274      real( kind = dp ) :: drdx, drdy, drdz
275      real( kind = dp ) :: fx, fy, fz
276 +    real( kind = dp ) :: myPot, myPotC, myDeriv, myDerivC, ros, rcos
277      real( kind = dp ) :: pot_temp, dudr
278 <    real( kind = dp ) :: sigma6
278 >    real( kind = dp ) :: sigmai
279      real( kind = dp ) :: epsilon
347 –    real( kind = dp ) :: r6
348 –    real( kind = dp ) :: t6
349 –    real( kind = dp ) :: t12
350 –    real( kind = dp ) :: delta
280      logical :: isSoftCore, shiftedPot
281      integer :: id1, id2, localError
282  
# Line 367 | Line 296 | contains
296      ljt1 = LJMap%atidToLJtype(atid1)
297      ljt2 = LJMap%atidToLJtype(atid2)
298  
299 <    sigma6     = MixingMap(ljt1,ljt2)%sigma6
299 >    sigmai     = MixingMap(ljt1,ljt2)%sigmai
300      epsilon    = MixingMap(ljt1,ljt2)%epsilon
372 –    delta      = MixingMap(ljt1,ljt2)%delta
301      isSoftCore = MixingMap(ljt1,ljt2)%isSoftCore
302      shiftedPot = MixingMap(ljt1,ljt2)%shiftedPot
303  
304 <    r6 = r2 * r2 * r2
304 >    ros = rij * sigmai
305 >    myPotC = 0.0_DP
306  
378 –    t6  = sigma6/ r6
379 –    t12 = t6 * t6    
380 –
307      if (isSoftCore) then
308 <      
309 <       pot_temp = 4.0E0_DP * epsilon * t6
308 >
309 >       call getSoftFunc(ros, myPot, myDeriv)
310 >
311         if (shiftedPot) then
312 <          pot_temp = pot_temp + delta
312 >          rcos = rcut * sigmai
313 >          call getSoftFunc(rcos, myPotC, myDerivC)
314         endif
315 <      
388 <       vpair = vpair + pot_temp
389 <      
390 <       dudr = -sw * 24.0E0_DP * epsilon * t6 / rij
391 <      
315 >              
316      else
317 <       pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
317 >
318 >       call getLJfunc(ros, myPot, myDeriv)
319 >
320         if (shiftedPot) then
321 <          pot_temp = pot_temp + delta
321 >          rcos = rcut * sigmai
322 >          call getLJfunc(rcos, myPotC, myDerivC)
323         endif
324        
398 –       vpair = vpair + pot_temp
399 –      
400 –       dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
325      endif
326  
327 +    pot_temp = epsilon * (myPot - myPotC)
328 +    vpair = vpair + pot_temp
329 +    dudr = sw * epsilon * myDeriv * sigmai
330 +
331      drdx = d(1) / rij
332      drdy = d(2) / rij
333      drdz = d(3) / rij
# Line 410 | Line 338 | contains
338  
339   #ifdef IS_MPI
340      if (do_pot) then
341 <       pot_Row(atom1) = pot_Row(atom1) + sw*pot_temp*0.5
342 <       pot_Col(atom2) = pot_Col(atom2) + sw*pot_temp*0.5
341 >       pot_Row(VDW_POT,atom1) = pot_Row(VDW_POT,atom1) + sw*pot_temp*0.5
342 >       pot_Col(VDW_POT,atom2) = pot_Col(VDW_POT,atom2) + sw*pot_temp*0.5
343      endif
344  
345      f_Row(1,atom1) = f_Row(1,atom1) + fx
# Line 470 | Line 398 | contains
398      end if
399      
400      haveMixingMap = .false.
401 +
402    end subroutine destroyLJTypes
403  
404 +  subroutine getLJfunc(r, myPot, myDeriv)
405 +
406 +    real(kind=dp), intent(in) :: r
407 +    real(kind=dp), intent(inout) :: myPot, myDeriv
408 +    real(kind=dp) :: ri, ri2, ri6, ri7, ri12, ri13
409 +    real(kind=dp) :: a, b, c, d, dx
410 +    integer :: j
411 +
412 +    ri = 1.0_DP / r
413 +    ri2 = ri*ri
414 +    ri6 = ri2*ri2*ri2
415 +    ri7 = ri6*ri
416 +    ri12 = ri6*ri6
417 +    ri13 = ri12*ri
418 +    
419 +    myPot = 4.0_DP * (ri12 - ri6)
420 +    myDeriv = 24.0_DP * (ri7 - 2.0_DP * ri13)
421 +    
422 +    return
423 +  end subroutine getLJfunc
424 +
425 +  subroutine getSoftFunc(r, myPot, myDeriv)
426 +    
427 +    real(kind=dp), intent(in) :: r
428 +    real(kind=dp), intent(inout) :: myPot, myDeriv
429 +    real(kind=dp) :: ri, ri2, ri6, ri7
430 +    real(kind=dp) :: a, b, c, d, dx
431 +    integer :: j
432 +    
433 +    ri = 1.0_DP / r    
434 +    ri2 = ri*ri
435 +    ri6 = ri2*ri2*ri2
436 +    ri7 = ri6*ri
437 +    myPot = 4.0_DP * (ri6)
438 +    myDeriv = - 24.0_DP * ri7
439 +    
440 +    return
441 +  end subroutine getSoftFunc
442 +
443   end module lj

Diff Legend

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