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 2319 by chuckv, Wed Sep 21 23:45:48 2005 UTC vs.
Revision 2756 by gezelter, Wed May 17 15:37:15 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.16 2005-09-21 23:45:48 chuckv Exp $, $Date: 2005-09-21 23:45:48 $, $Name: not supported by cvs2svn $, $Revision: 1.16 $
46 > !! @version $Id: LJ.F90,v 1.25 2006-05-17 15:37:14 gezelter Exp $, $Date: 2006-05-17 15:37:14 $, $Name: not supported by cvs2svn $, $Revision: 1.25 $
47  
48  
49   module lj
50 +  use definitions
51    use atype_module
52    use vector_class
53    use simulation
54    use status
55 +  use fForceOptions
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  
62 –  integer, parameter :: DP = selected_real_kind(15)
63 –
66    logical, save :: useGeometricDistanceMixing = .false.
67    logical, save :: haveMixingMap = .false.
68  
# Line 68 | Line 70 | module lj
70    logical, save :: defaultShift = .false.
71    logical, save :: haveDefaultCutoff = .false.
72  
71 –
73    type, private :: LJtype
74       integer       :: atid
75       real(kind=dp) :: sigma
# Line 88 | Line 89 | module lj
89    type :: MixParameters
90       real(kind=DP) :: sigma
91       real(kind=DP) :: epsilon
92 <     real(kind=dp) :: sigma6
92 >     real(kind=dp) :: sigmai
93       real(kind=dp) :: rCut
93 –     real(kind=dp) :: delta
94       logical       :: rCutWasSet = .false.
95       logical       :: shiftedPot
96       logical       :: isSoftCore = .false.
# Line 100 | Line 100 | module lj
100  
101    public :: newLJtype
102    public :: setLJDefaultCutoff
103 –  public :: setLJUniformCutoff
104 –  public :: setLJCutoffByTypes
103    public :: getSigma
104    public :: getEpsilon
107 –  public :: useGeometricMixing
105    public :: do_lj_pair
106    public :: destroyLJtypes
107  
# Line 161 | Line 158 | contains
158      if(LJMap%nLJTypes /= 0) then
159         call createMixingMap()
160      end if
164 –  end subroutine setLJDefaultCutoff
161  
162 <  subroutine setLJUniformCutoff(thisRcut, shiftedPot)
167 <    real(kind=dp), intent(in) :: thisRcut
168 <    logical, intent(in) :: shiftedPot
169 <    integer :: nLJtypes, i, j
170 <
171 <    if (LJMap%currentLJtype == 0) then
172 <       call handleError("LJ", "No members in LJMap")
173 <       return
174 <    end if
175 <
176 <    nLJtypes = LJMap%nLJtypes
177 <    if (.not. allocated(MixingMap)) then
178 <       allocate(MixingMap(nLJtypes, nLJtypes))
179 <    endif
180 <
181 <    do i = 1, nLJtypes
182 <       do j = 1, nLJtypes
183 <          MixingMap(i,j)%rCut = thisRcut
184 <          MixingMap(i,j)%shiftedPot = shiftedPot
185 <          MixingMap(i,j)%rCutWasSet = .true.
186 <       enddo
187 <    enddo
188 <    call createMixingMap()
189 <  end subroutine setLJUniformCutoff
190 <
191 <  subroutine setLJCutoffByTypes(atid1, atid2, thisRcut, shiftedPot)
192 <    integer, intent(in) :: atid1, atid2
193 <    real(kind=dp), intent(in) :: thisRcut
194 <    logical, intent(in) :: shiftedPot
195 <    integer :: nLJtypes, ljt1, ljt2
196 <
197 <    if (LJMap%currentLJtype == 0) then
198 <       call handleError("LJ", "No members in LJMap")
199 <       return
200 <    end if
201 <
202 <    nLJtypes = LJMap%nLJtypes
203 <    if (.not. allocated(MixingMap)) then
204 <       allocate(MixingMap(nLJtypes, nLJtypes))
205 <    endif
206 <
207 <    ljt1 = LJMap%atidToLJtype(atid1)
208 <    ljt2 = LJMap%atidToLJtype(atid2)
162 >  end subroutine setLJDefaultCutoff
163  
210 –    MixingMap(ljt1,ljt2)%rCut = thisRcut
211 –    MixingMap(ljt1,ljt2)%shiftedPot = shiftedPot
212 –    MixingMap(ljt1,ljt2)%rCutWasSet = .true.
213 –    MixingMap(ljt2,ljt1)%rCut = thisRcut
214 –    MixingMap(ljt2,ljt1)%shiftedPot = shiftedPot
215 –    MixingMap(ljt2,ljt1)%rCutWasSet = .true.
216 –
217 –    call createMixingMap()
218 –  end subroutine setLJCutoffByTypes
219 –
164    function getSigma(atid) result (s)
165      integer, intent(in) :: atid
166      integer :: ljt1
# Line 247 | Line 191 | contains
191  
192    end function getEpsilon
193  
250 –  subroutine useGeometricMixing()
251 –    useGeometricDistanceMixing = .true.
252 –    haveMixingMap = .false.
253 –    return
254 –  end subroutine useGeometricMixing
255 –
194    subroutine createMixingMap()
195      integer :: nLJtypes, i, j
196      real ( kind = dp ) :: s1, s2, e1, e2
# Line 270 | Line 208 | contains
208         allocate(MixingMap(nLJtypes, nLJtypes))
209      endif
210  
211 +    useGeometricDistanceMixing = usesGeometricDistanceMixing()
212      do i = 1, nLJtypes
213  
214         s1 = LJMap%LJtypes(i)%sigma
# Line 286 | Line 225 | contains
225  
226            ! only the distance parameter uses different mixing policies
227            if (useGeometricDistanceMixing) then
228 <             MixingMap(i,j)%sigma = dsqrt(s1 * s2)
228 >             MixingMap(i,j)%sigma = sqrt(s1 * s2)
229            else
230               MixingMap(i,j)%sigma = 0.5_dp * (s1 + s2)
231            endif
232            
233 <          MixingMap(i,j)%epsilon = dsqrt(e1 * e2)
233 >          MixingMap(i,j)%epsilon = sqrt(e1 * e2)
234  
235 <          MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
235 >          MixingMap(i,j)%sigmai = 1.0_DP  / (MixingMap(i,j)%sigma)
236  
237 <          if (MixingMap(i,j)%rCutWasSet) then
238 <             rcut6 = (MixingMap(i,j)%rcut)**6
237 >          if (haveDefaultCutoff) then
238 >             MixingMap(i,j)%shiftedPot = defaultShift
239            else
240 <             if (haveDefaultCutoff) then
241 <                rcut6 = defaultCutoff**6
303 <                doShift = defaultShift
304 <             else
305 <                call handleError("LJ", "No specified or default cutoff value!")
306 <             endif
307 <          endif
308 <          
309 <          tp6    = MixingMap(i,j)%sigma6/rcut6
310 <          tp12    = tp6**2          
311 <          MixingMap(i,j)%delta =-4.0_DP*MixingMap(i,j)%epsilon*(tp12 - tp6)
312 <          MixingMap(i,j)%shiftedPot = doShift
240 >             MixingMap(i,j)%shiftedPot = defaultShift
241 >          endif          
242  
243            if (i.ne.j) then
244               MixingMap(j,i)%sigma      = MixingMap(i,j)%sigma
245               MixingMap(j,i)%epsilon    = MixingMap(i,j)%epsilon
246 <             MixingMap(j,i)%sigma6     = MixingMap(i,j)%sigma6
246 >             MixingMap(j,i)%sigmai     = MixingMap(i,j)%sigmai
247               MixingMap(j,i)%rCut       = MixingMap(i,j)%rCut
319 –             MixingMap(j,i)%delta      = MixingMap(i,j)%delta
248               MixingMap(j,i)%rCutWasSet = MixingMap(i,j)%rCutWasSet
249               MixingMap(j,i)%shiftedPot = MixingMap(i,j)%shiftedPot
250               MixingMap(j,i)%isSoftCore = MixingMap(i,j)%isSoftCore
# Line 328 | Line 256 | contains
256      haveMixingMap = .true.
257      
258    end subroutine createMixingMap
259 <  
260 <  subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
259 >          
260 >  subroutine do_lj_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
261         pot, f, do_pot)
262 <
262 >    
263      integer, intent(in) ::  atom1, atom2
264      integer :: atid1, atid2, ljt1, ljt2
265 <    real( kind = dp ), intent(in) :: rij, r2
265 >    real( kind = dp ), intent(in) :: rij, r2, rcut
266      real( kind = dp ) :: pot, sw, vpair
267      real( kind = dp ), dimension(3,nLocal) :: f    
268      real( kind = dp ), intent(in), dimension(3) :: d
# Line 344 | Line 272 | contains
272      ! local Variables
273      real( kind = dp ) :: drdx, drdy, drdz
274      real( kind = dp ) :: fx, fy, fz
275 +    real( kind = dp ) :: myPot, myPotC, myDeriv, myDerivC, ros, rcos
276      real( kind = dp ) :: pot_temp, dudr
277 <    real( kind = dp ) :: sigma6
277 >    real( kind = dp ) :: sigmai
278      real( kind = dp ) :: epsilon
350 –    real( kind = dp ) :: r6
351 –    real( kind = dp ) :: t6
352 –    real( kind = dp ) :: t12
353 –    real( kind = dp ) :: delta
279      logical :: isSoftCore, shiftedPot
280      integer :: id1, id2, localError
281  
# Line 370 | Line 295 | contains
295      ljt1 = LJMap%atidToLJtype(atid1)
296      ljt2 = LJMap%atidToLJtype(atid2)
297  
298 <    sigma6     = MixingMap(ljt1,ljt2)%sigma6
298 >    sigmai     = MixingMap(ljt1,ljt2)%sigmai
299      epsilon    = MixingMap(ljt1,ljt2)%epsilon
375 –    delta      = MixingMap(ljt1,ljt2)%delta
300      isSoftCore = MixingMap(ljt1,ljt2)%isSoftCore
301      shiftedPot = MixingMap(ljt1,ljt2)%shiftedPot
302  
303 <    r6 = r2 * r2 * r2
303 >    ros = rij * sigmai
304 >    myPotC = 0.0_DP
305  
381 –    t6  = sigma6/ r6
382 –    t12 = t6 * t6    
383 –
306      if (isSoftCore) then
307 <      
308 <       pot_temp = 4.0E0_DP * epsilon * t6
307 >
308 >       call getSoftFunc(ros, myPot, myDeriv)
309 >
310         if (shiftedPot) then
311 <          pot_temp = pot_temp + delta
311 >          rcos = rcut * sigmai
312 >          call getSoftFunc(rcos, myPotC, myDerivC)
313         endif
314 <      
391 <       vpair = vpair + pot_temp
392 <      
393 <       dudr = -sw * 24.0E0_DP * epsilon * t6 / rij
394 <      
314 >              
315      else
316 <       pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
316 >
317 >       call getLJfunc(ros, myPot, myDeriv)
318 >
319         if (shiftedPot) then
320 <          pot_temp = pot_temp + delta
320 >          rcos = rcut * sigmai
321 >          call getLJfunc(rcos, myPotC, myDerivC)
322         endif
323        
401 –       vpair = vpair + pot_temp
402 –      
403 –       dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
324      endif
325  
326 +    pot_temp = epsilon * (myPot - myPotC)
327 +    vpair = vpair + pot_temp
328 +    dudr = sw * epsilon * myDeriv * sigmai
329 +
330      drdx = d(1) / rij
331      drdy = d(2) / rij
332      drdz = d(3) / rij
# Line 413 | Line 337 | contains
337  
338   #ifdef IS_MPI
339      if (do_pot) then
340 <       pot_Row(atom1) = pot_Row(atom1) + sw*pot_temp*0.5
341 <       pot_Col(atom2) = pot_Col(atom2) + sw*pot_temp*0.5
340 >       pot_Row(VDW_POT,atom1) = pot_Row(VDW_POT,atom1) + sw*pot_temp*0.5
341 >       pot_Col(VDW_POT,atom2) = pot_Col(VDW_POT,atom2) + sw*pot_temp*0.5
342      endif
343  
344      f_Row(1,atom1) = f_Row(1,atom1) + fx
# Line 473 | Line 397 | contains
397      end if
398      
399      haveMixingMap = .false.
400 +
401    end subroutine destroyLJTypes
402  
403 +  subroutine getLJfunc(r, myPot, myDeriv)
404 +
405 +    real(kind=dp), intent(in) :: r
406 +    real(kind=dp), intent(inout) :: myPot, myDeriv
407 +    real(kind=dp) :: ri, ri2, ri6, ri7, ri12, ri13
408 +    real(kind=dp) :: a, b, c, d, dx
409 +    integer :: j
410 +
411 +    ri = 1.0_DP / r
412 +    ri2 = ri*ri
413 +    ri6 = ri2*ri2*ri2
414 +    ri7 = ri6*ri
415 +    ri12 = ri6*ri6
416 +    ri13 = ri12*ri
417 +    
418 +    myPot = 4.0_DP * (ri12 - ri6)
419 +    myDeriv = 24.0_DP * (ri7 - 2.0_DP * ri13)
420 +    
421 +    return
422 +  end subroutine getLJfunc
423 +
424 +  subroutine getSoftFunc(r, myPot, myDeriv)
425 +    
426 +    real(kind=dp), intent(in) :: r
427 +    real(kind=dp), intent(inout) :: myPot, myDeriv
428 +    real(kind=dp) :: ri, ri2, ri6, ri7
429 +    real(kind=dp) :: a, b, c, d, dx
430 +    integer :: j
431 +    
432 +    ri = 1.0_DP / r    
433 +    ri2 = ri*ri
434 +    ri6 = ri2*ri2*ri2
435 +    ri7 = ri6*ri
436 +    myPot = 4.0_DP * (ri6)
437 +    myDeriv = - 24.0_DP * ri7
438 +    
439 +    return
440 +  end subroutine getSoftFunc
441 +
442   end module lj

Diff Legend

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