| 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 |
| 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 |
| 234 |
|
e2 = GBMap%GBtypes(j)%eps |
| 235 |
|
er2 = GBMap%GBtypes(j)%eps_ratio |
| 236 |
|
dw2 = GBMap%GBtypes(j)%dw |
| 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) |
| 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) |
| 244 |
|
GBMixingMap(i,j)%x2 = (l1*l1 - d1*d1) * (l2*l2 - d2*d2) / & |
| 257 |
|
GBMixingMap(i,j)%xp2 = xp*xp |
| 258 |
|
GBMixingMap(i,j)%xpap2 = xp*ap2 |
| 259 |
|
GBMixingMap(i,j)%xpapi2 = xp/ap2 |
| 258 |
– |
|
| 259 |
– |
if (i.ne.j) then |
| 260 |
– |
GBMixingMap(j,i)%sigma0 = GBMixingMap(i,j)%sigma0 |
| 261 |
– |
GBMixingMap(j,i)%dw = GBMixingMap(i,j)%dw |
| 262 |
– |
GBMixingMap(j,i)%eps0 = GBMixingMap(i,j)%eps0 |
| 263 |
– |
GBMixingMap(j,i)%x2 = GBMixingMap(i,j)%x2 |
| 264 |
– |
GBMixingMap(j,i)%xa2 = GBMixingMap(i,j)%xa2 |
| 265 |
– |
GBMixingMap(j,i)%xai2 = GBMixingMap(i,j)%xai2 |
| 266 |
– |
GBMixingMap(j,i)%xp2 = GBMixingMap(i,j)%xp2 |
| 267 |
– |
GBMixingMap(j,i)%xpap2 = GBMixingMap(i,j)%xpap2 |
| 268 |
– |
GBMixingMap(j,i)%xpapi2 = GBMixingMap(i,j)%xpapi2 |
| 269 |
– |
endif |
| 260 |
|
enddo |
| 261 |
|
enddo |
| 262 |
|
haveMixingMap = .true. |
| 281 |
|
gbt1 = GBMap%atidToGBtype(atomID) |
| 282 |
|
l = GBMap%GBtypes(gbt1)%l |
| 283 |
|
d = GBMap%GBtypes(gbt1)%d |
| 294 |
– |
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, & |
| 374 |
|
else |
| 375 |
|
g = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3) |
| 376 |
|
endif |
| 377 |
+ |
|
| 378 |
|
au = a / r |
| 379 |
|
bu = b / r |
| 380 |
|
|
| 381 |
< |
au2 = au*au |
| 382 |
< |
bu2 = bu*bu |
| 383 |
< |
g2 = g*g |
| 381 |
> |
au2 = au * au |
| 382 |
> |
bu2 = bu * bu |
| 383 |
> |
g2 = g * g |
| 384 |
|
|
| 391 |
– |
|
| 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 |
| 404 |
|
|
| 405 |
|
pref1 = - 8.0_dp * eps * mu * (R12 - R6) / (e2 * r) |
| 406 |
|
|
| 413 |
– |
|
| 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)) |
