| 234 |
|
e2 = GBMap%GBtypes(j)%eps |
| 235 |
|
er2 = GBMap%GBtypes(j)%eps_ratio |
| 236 |
|
dw2 = GBMap%GBtypes(j)%dw |
| 237 |
< |
|
| 238 |
< |
GBMixingMap(i,j)%sigma0 = sqrt(d1*d1 + d2*d2) |
| 237 |
> |
! write(*,*) 'd2 = ', d2,l2, d1,l2 |
| 238 |
> |
! GBMixingMap(i,j)%sigma0 = sqrt(d1*d1 + d2*d2) |
| 239 |
> |
GBMixingMap(i,j)%sigma0 = 0.5_dp*(d1 + d2) |
| 240 |
|
GBMixingMap(i,j)%xa2 = (l1*l1 - d1*d1)/(l1*l1 + d2*d2) |
| 241 |
|
GBMixingMap(i,j)%xai2 = (l2*l2 - d2*d2)/(l2*l2 + d1*d1) |
| 242 |
|
GBMixingMap(i,j)%x2 = (l1*l1 - d1*d1) * (l2*l2 - d2*d2) / & |
| 311 |
|
real (kind = dp), dimension(3) :: ul1, ul2, rxu1, rxu2, uxu, rhat |
| 312 |
|
|
| 313 |
|
real (kind = dp) :: sigma0, dw, eps0, x2, xa2, xai2, xp2, xpap2, xpapi2 |
| 314 |
< |
real (kind = dp) :: e1, e2, eps, sigma, s3, s03, au, bu, a, b, g, g2 |
| 314 |
> |
real (kind = dp) :: e1, e2, eps, sigma, s3, s03, au2, bu2, au, bu, a, b, g, g2 |
| 315 |
|
real (kind = dp) :: U, BigR, R3, R6, R7, R12, R13, H, Hp, fx, fy, fz |
| 316 |
|
real (kind = dp) :: dUdr, dUda, dUdb, dUdg, pref1, pref2 |
| 317 |
|
logical :: i_is_lj, j_is_lj |
| 381 |
|
else |
| 382 |
|
g = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3) |
| 383 |
|
endif |
| 383 |
– |
|
| 384 |
|
au = a / r |
| 385 |
|
bu = b / r |
| 386 |
+ |
|
| 387 |
+ |
au2 = au*au |
| 388 |
+ |
bu2 = bu*bu |
| 389 |
|
g2 = g*g |
| 390 |
|
|
| 391 |
< |
H = (xa2 * au + xai2 * bu - 2.0_dp*x2*au*bu*g) / (1.0_dp - x2*g2) |
| 392 |
< |
Hp = (xpap2*au + xpapi2*bu - 2.0_dp*xp2*au*bu*g) / (1.0_dp - xp2*g2) |
| 391 |
> |
|
| 392 |
> |
H = (xa2 * au2 + xai2 * bu2 - 2.0_dp*x2*au*bu*g) / (1.0_dp - x2*g2) |
| 393 |
> |
Hp = (xpap2*au2 + xpapi2*bu2 - 2.0_dp*xp2*au*bu*g) / (1.0_dp - xp2*g2) |
| 394 |
|
sigma = sigma0 / sqrt(1.0_dp - H) |
| 395 |
|
e1 = 1.0_dp / sqrt(1.0_dp - x2*g2) |
| 396 |
|
e2 = 1.0_dp - Hp |
| 410 |
|
|
| 411 |
|
pref1 = - 8.0_dp * eps * mu * (R12 - R6) / (e2 * r) |
| 412 |
|
|
| 413 |
+ |
|
| 414 |
|
pref2 = 8.0_dp * eps * s3 * (6.0_dp*R13 - 3.0_dp*R7) / (dw*r*s03) |
| 415 |
|
|
| 416 |
< |
dUdr = - (pref1 * Hp + pref2 * (sigma0*sigma0*r/s3 - H)) |
| 416 |
> |
dUdr = - (pref1 * Hp + pref2 * (sigma0*sigma0*r/s3 + H)) |
| 417 |
|
|
| 418 |
|
dUda = pref1 * (xpap2*au - xp2*bu*g) / (1.0_dp - xp2 * g2) & |
| 419 |
|
+ pref2 * (xa2 * au - x2 *bu*g) / (1.0_dp - x2 * g2) |
| 420 |
< |
|
| 420 |
> |
|
| 421 |
|
dUdb = pref1 * (xpapi2*bu - xp2*au*g) / (1.0_dp - xp2 * g2) & |
| 422 |
|
+ pref2 * (xai2 * bu - x2 *au*g) / (1.0_dp - x2 * g2) |
| 423 |
|
|
| 426 |
|
(1.0_dp - xp2 * g2) / e2 & |
| 427 |
|
+ 8.0_dp * eps * s3 * (3.0_dp * R7 - 6.0_dp * R13) * & |
| 428 |
|
(x2 * au * bu - H * x2 * g) / (1.0_dp - x2 * g2) / (dw * s03) |
| 429 |
+ |
|
| 430 |
|
|
| 431 |
|
rhat = d / r |
| 432 |
|
|
| 437 |
|
rxu1 = cross_product(d, ul1) |
| 438 |
|
rxu2 = cross_product(d, ul2) |
| 439 |
|
uxu = cross_product(ul1, ul2) |
| 440 |
< |
|
| 440 |
> |
|
| 441 |
> |
!!$ write(*,*) 'pref = ' , pref1, pref2 |
| 442 |
|
!!$ write(*,*) 'rxu1 = ' , rxu1(1), rxu1(2), rxu1(3) |
| 443 |
|
!!$ write(*,*) 'rxu2 = ' , rxu2(1), rxu2(2), rxu2(3) |
| 444 |
|
!!$ write(*,*) 'uxu = ' , uxu(1), uxu(2), uxu(3) |
| 445 |
< |
!!$ write(*,*) 'dUda = ', dUda, dudb, dudg |
| 446 |
< |
|
| 447 |
< |
|
| 445 |
> |
!!$ write(*,*) 'dUda = ', dUda, dudb, dudg, dudr |
| 446 |
> |
!!$ write(*,*) 'H = ', H,hp,sigma, e1, e2, BigR |
| 447 |
> |
!!$ write(*,*) 'chi = ', xa2, xai2, x2 |
| 448 |
> |
!!$ write(*,*) 'chip = ', xpap2, xpapi2, xp2 |
| 449 |
> |
!!$ write(*,*) 'eps = ', eps0, e1, e2, eps |
| 450 |
> |
!!$ write(*,*) 'U =', U, pref1, pref2 |
| 451 |
> |
!!$ write(*,*) 'f =', fx, fy, fz |
| 452 |
> |
!!$ write(*,*) 'au =', au, bu, g |
| 453 |
> |
!!$ |
| 454 |
> |
|
| 455 |
> |
|
| 456 |
|
#ifdef IS_MPI |
| 457 |
|
f_Row(1,atom1) = f_Row(1,atom1) + fx |
| 458 |
|
f_Row(2,atom1) = f_Row(2,atom1) + fy |
