ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/gb.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/gb.F90 (file contents):
Revision 2787 by gezelter, Mon Jun 5 18:24:45 2006 UTC vs.
Revision 2958 by xsun, Mon Jul 24 14:51:09 2006 UTC

# Line 49 | Line 49 | module gayberne
49    use linearalgebra
50    use status
51    use lj
52 +  use fForceOptions
53   #ifdef IS_MPI
54    use mpiSimulation
55   #endif
# Line 186 | Line 187 | contains
187            GBMap%atidToGBtype(myATID) = current
188            GBMap%GBtypes(current)%atid      = myATID      
189            GBMap%GBtypes(current)%isLJ      = .true.          
190 <          GBMap%GBtypes(current)%d         = getSigma(myATID)
190 >          GBMap%GBtypes(current)%d         = getSigma(myATID) / sqrt(2.0_dp)
191            GBMap%GBtypes(current)%l         = GBMap%GBtypes(current)%d
192            GBMap%GBtypes(current)%eps       = getEpsilon(myATID)
193            GBMap%GBtypes(current)%eps_ratio = 1.0_dp
# Line 197 | Line 198 | contains
198      end do
199      
200      haveGBMap = .true.
201 +
202      
203    end subroutine complete_GB_FF
204  
# Line 225 | Line 227 | contains
227         er1 = GBMap%GBtypes(i)%eps_ratio
228         dw1 = GBMap%GBtypes(i)%dw
229  
230 <       do j = i, nGBtypes
230 >       do j = 1, nGBtypes
231  
232            d2 = GBMap%GBtypes(j)%d
233            l2 = GBMap%GBtypes(j)%l
# Line 233 | Line 235 | contains
235            er2 = GBMap%GBtypes(j)%eps_ratio
236            dw2 = GBMap%GBtypes(j)%dw
237  
238 + !  Cleaver paper uses sqrt of squares to get sigma0 for
239 + !  mixed interactions.
240 +            
241            GBMixingMap(i,j)%sigma0 = sqrt(d1*d1 + d2*d2)
242            GBMixingMap(i,j)%xa2 = (l1*l1 - d1*d1)/(l1*l1 + d2*d2)
243            GBMixingMap(i,j)%xai2 = (l2*l2 - d2*d2)/(l2*l2 + d1*d1)
# Line 252 | Line 257 | contains
257            GBMixingMap(i,j)%xp2 = xp*xp
258            GBMixingMap(i,j)%xpap2 = xp*ap2
259            GBMixingMap(i,j)%xpapi2 = xp/ap2
255 –          
256 –          if (i.ne.j) then
257 –             GBMixingMap(j,i)%sigma0 = GBMixingMap(i,j)%sigma0
258 –             GBMixingMap(j,i)%dw     = GBMixingMap(i,j)%dw    
259 –             GBMixingMap(j,i)%eps0   = GBMixingMap(i,j)%eps0  
260 –             GBMixingMap(j,i)%x2     = GBMixingMap(i,j)%x2    
261 –             GBMixingMap(j,i)%xa2    = GBMixingMap(i,j)%xa2  
262 –             GBMixingMap(j,i)%xai2   = GBMixingMap(i,j)%xai2  
263 –             GBMixingMap(j,i)%xp2    = GBMixingMap(i,j)%xp2  
264 –             GBMixingMap(j,i)%xpap2  = GBMixingMap(i,j)%xpap2
265 –             GBMixingMap(j,i)%xpapi2 = GBMixingMap(i,j)%xpapi2
266 –          endif
260         enddo
261      enddo
262      haveMixingMap = .true.
263 <    
263 >    mu = getGayBerneMu()
264 >    nu = getGayBerneNu()    
265    end subroutine createGBMixingMap
266    
267  
# Line 287 | Line 281 | contains
281      gbt1 = GBMap%atidToGBtype(atomID)
282      l = GBMap%GBtypes(gbt1)%l
283      d = GBMap%GBtypes(gbt1)%d  
290 –    cutValue = 2.5_dp*max(l,d)
284  
285 +    ! sigma is actually sqrt(2)*l  for prolate ellipsoids
286 +    
287 +    cutValue = 2.5_dp*sqrt(2.0_dp)*max(l,d)
288 +
289    end function getGayBerneCut
290  
291    subroutine do_gb_pair(atom1, atom2, d, r, r2, sw, vpair, fpair, &
# Line 307 | Line 304 | contains
304      real (kind = dp), dimension(3) :: ul1, ul2, rxu1, rxu2, uxu, rhat
305  
306      real (kind = dp) :: sigma0, dw, eps0, x2, xa2, xai2, xp2, xpap2, xpapi2
307 <    real (kind = dp) :: e1, e2, eps, sigma, s3, s03, au, bu, a, b, g, g2
307 >    real (kind = dp) :: e1, e2, eps, sigma, s3, s03, au2, bu2, au, bu, a, b, g, g2
308      real (kind = dp) :: U, BigR, R3, R6, R7, R12, R13, H, Hp, fx, fy, fz
309      real (kind = dp) :: dUdr, dUda, dUdb, dUdg, pref1, pref2
310      logical :: i_is_lj, j_is_lj
# Line 380 | Line 377 | contains
377  
378      au = a / r
379      bu = b / r
383 –    g2 = g*g
380  
381 <    H  = (xa2 * au + xai2 * bu - 2.0_dp*x2*au*bu*g)  / (1.0_dp - x2*g2)
382 <    Hp = (xpap2*au + xpapi2*bu - 2.0_dp*xp2*au*bu*g) / (1.0_dp - xp2*g2)
381 >    au2 = au * au
382 >    bu2 = bu * bu
383 >    g2 = g * g
384 >
385 >    H  = (xa2 * au2 + xai2 * bu2 - 2.0_dp*x2*au*bu*g)  / (1.0_dp - x2*g2)
386 >    Hp = (xpap2*au2 + xpapi2*bu2 - 2.0_dp*xp2*au*bu*g) / (1.0_dp - xp2*g2)
387 >
388      sigma = sigma0 / sqrt(1.0_dp - H)
389      e1 = 1.0_dp / sqrt(1.0_dp - x2*g2)
390      e2 = 1.0_dp - Hp
# Line 402 | Line 403 | contains
403      s03 = sigma0*sigma0*sigma0
404  
405      pref1 = - 8.0_dp * eps * mu * (R12 - R6) / (e2 * r)
406 +
407      pref2 = 8.0_dp * eps * s3 * (6.0_dp*R13 - 3.0_dp*R7) / (dw*r*s03)
408  
409 <    dUdr = - (pref1 * Hp + pref2 * (sigma0*sigma0*r/s3 - H))
409 >    dUdr = - (pref1 * Hp + pref2 * (sigma0*sigma0*r/s3 + H))
410      
411      dUda = pref1 * (xpap2*au - xp2*bu*g) / (1.0_dp - xp2 * g2) &
412           + pref2 * (xa2 * au - x2 *bu*g) / (1.0_dp - x2 * g2)
413 <    
413 >
414      dUdb = pref1 * (xpapi2*bu - xp2*au*g) / (1.0_dp - xp2 * g2) &
415           + pref2 * (xai2 * bu - x2 *au*g) / (1.0_dp - x2 * g2)
416  
417      dUdg = 4.0_dp * eps * nu * (R12 - R6) * x2 * g / (1.0_dp - x2*g2) &
418           + 8.0_dp * eps * mu * (R12 - R6) * (xp2*au*bu - Hp*xp2*g) / &
419           (1.0_dp - xp2 * g2) / e2 &
420 <         + 8.0_dp * eps * s3 * (3.0_dp * R7 - 6.0_dp * R13) * &
420 >         + 8.0_dp * eps * s3 * (3.0_dp * R7 - 6.0_dp * R13) * &  
421           (x2 * au * bu - H * x2 * g) / (1.0_dp - x2 * g2) / (dw * s03)
422 +
423              
424      rhat = d / r
425  
426 <    fx = -dUdr * rhat(1) - dUda * ul1(1) - dUdb * ul2(1)
427 <    fy = -dUdr * rhat(2) - dUda * ul1(2) - dUdb * ul2(2)
428 <    fx = -dUdr * rhat(3) - dUda * ul1(3) - dUdb * ul2(3)
426 >    fx = dUdr * rhat(1) + dUda * ul1(1) + dUdb * ul2(1)
427 >    fy = dUdr * rhat(2) + dUda * ul1(2) + dUdb * ul2(2)
428 >    fz = dUdr * rhat(3) + dUda * ul1(3) + dUdb * ul2(3)    
429  
430      rxu1 = cross_product(d, ul1)
431      rxu2 = cross_product(d, ul2)    
432      uxu = cross_product(ul1, ul2)
433 <          
433 >    
434 > !!$    write(*,*) 'pref = ' , pref1, pref2
435 > !!$    write(*,*) 'rxu1 = ' , rxu1(1), rxu1(2), rxu1(3)
436 > !!$    write(*,*) 'rxu2 = ' , rxu2(1), rxu2(2), rxu2(3)
437 > !!$    write(*,*) 'uxu = ' , uxu(1), uxu(2), uxu(3)
438 > !!$    write(*,*) 'dUda = ', dUda, dudb, dudg, dudr
439 > !!$    write(*,*) 'H = ', H,hp,sigma, e1, e2, BigR
440 > !!$    write(*,*) 'chi = ', xa2, xai2, x2
441 > !!$    write(*,*) 'chip = ', xpap2, xpapi2, xp2
442 > !!$    write(*,*) 'eps = ', eps0, e1, e2, eps
443 > !!$    write(*,*) 'U =', U, pref1, pref2
444 > !!$    write(*,*) 'f =', fx, fy, fz
445 > !!$    write(*,*) 'au =', au, bu, g
446 > !!$    
447 >        
448 >  
449   #ifdef IS_MPI
450      f_Row(1,atom1) = f_Row(1,atom1) + fx
451      f_Row(2,atom1) = f_Row(2,atom1) + fy

Diff Legend

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