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 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.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.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 61 | Line 63 | module lj
63   #define __FORTRAN90
64   #include "UseTheForce/DarkSide/fInteractionMap.h"
65  
64 –  integer, parameter :: DP = selected_real_kind(15)
65 –
66    logical, save :: useGeometricDistanceMixing = .false.
67    logical, save :: haveMixingMap = .false.
68  
# Line 70 | Line 70 | module lj
70    logical, save :: defaultShift = .false.
71    logical, save :: haveDefaultCutoff = .false.
72  
73 –
73    type, private :: LJtype
74       integer       :: atid
75       real(kind=dp) :: sigma
# Line 90 | 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
95 –     real(kind=dp) :: delta
94       logical       :: rCutWasSet = .false.
95       logical       :: shiftedPot
96       logical       :: isSoftCore = .false.
# Line 102 | Line 100 | module lj
100  
101    public :: newLJtype
102    public :: setLJDefaultCutoff
105 –  public :: setLJUniformCutoff
106 –  public :: setLJCutoffByTypes
103    public :: getSigma
104    public :: getEpsilon
109 –  public :: useGeometricMixing
105    public :: do_lj_pair
106    public :: destroyLJtypes
107  
# Line 163 | Line 158 | contains
158      if(LJMap%nLJTypes /= 0) then
159         call createMixingMap()
160      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
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)
161  
162 <    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.
162 >  end subroutine setLJDefaultCutoff
163  
219 –    call createMixingMap()
220 –  end subroutine setLJCutoffByTypes
221 –
164    function getSigma(atid) result (s)
165      integer, intent(in) :: atid
166      integer :: ljt1
# Line 249 | Line 191 | contains
191  
192    end function getEpsilon
193  
252 –  subroutine useGeometricMixing()
253 –    useGeometricDistanceMixing = .true.
254 –    haveMixingMap = .false.
255 –    return
256 –  end subroutine useGeometricMixing
257 –
194    subroutine createMixingMap()
195      integer :: nLJtypes, i, j
196      real ( kind = dp ) :: s1, s2, e1, e2
# Line 272 | 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 288 | 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
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
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
321 –             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 330 | 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 346 | 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
352 –    real( kind = dp ) :: r6
353 –    real( kind = dp ) :: t6
354 –    real( kind = dp ) :: t12
355 –    real( kind = dp ) :: delta
279      logical :: isSoftCore, shiftedPot
280      integer :: id1, id2, localError
281  
# Line 372 | 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
377 –    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  
383 –    t6  = sigma6/ r6
384 –    t12 = t6 * t6    
385 –
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 <      
393 <       vpair = vpair + pot_temp
394 <      
395 <       dudr = -sw * 24.0E0_DP * epsilon * t6 / rij
396 <      
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        
403 –       vpair = vpair + pot_temp
404 –      
405 –       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 415 | Line 337 | contains
337  
338   #ifdef IS_MPI
339      if (do_pot) then
340 <       pot_Row(LJ_POT,atom1) = pot_Row(LJ_POT,atom1) + sw*pot_temp*0.5
341 <       pot_Col(LJ_POT,atom2) = pot_Col(LJ_POT,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 475 | 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