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 |