| 57 |
|
real(kind=dp), parameter :: pre11 = 332.0637778_dp |
| 58 |
|
real(kind=dp), parameter :: pre12 = 69.13291783_dp |
| 59 |
|
real(kind=dp), parameter :: pre22 = 14.39289874_dp |
| 60 |
+ |
real(kind=dp), parameter :: pre14 = 0.0_dp |
| 61 |
|
|
| 62 |
|
public :: newElectrostaticType |
| 63 |
|
public :: setCharge |
| 318 |
|
integer :: me1, me2, id1, id2 |
| 319 |
|
real (kind=dp) :: q_i, q_j, mu_i, mu_j, d_i, d_j |
| 320 |
|
real (kind=dp) :: ct_i, ct_j, ct_ij, a1 |
| 321 |
< |
real (kind=dp) :: riji, ri2, ri3, ri4 |
| 321 |
> |
real (kind=dp) :: riji, ri, ri2, ri3, ri4 |
| 322 |
|
real (kind=dp) :: pref, vterm, epot, dudr |
| 323 |
+ |
real (kind=dp) :: xhat, yhat, zhat |
| 324 |
|
real (kind=dp) :: dudx, dudy, dudz |
| 325 |
|
real (kind=dp) :: drdxj, drdyj, drdzj |
| 326 |
|
real (kind=dp) :: duduix, duduiy, duduiz, dudujx, dudujy, dudujz |
| 327 |
< |
|
| 327 |
> |
real (kind=dp) :: scale, sc2, bigR |
| 328 |
|
|
| 329 |
|
if (.not.allocated(ElectrostaticMap)) then |
| 330 |
|
call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!") |
| 343 |
|
|
| 344 |
|
riji = 1.0d0 / rij |
| 345 |
|
|
| 346 |
< |
!! these are also useful as the unit vector of \vec{r} |
| 347 |
< |
!! \hat{r} = \vec{r} / r = {(x_j-x_i) / r, (y_j-y_i)/r, (z_j-z_i)/r} |
| 346 |
> |
xhat = d(1) * riji |
| 347 |
> |
yhat = d(2) * riji |
| 348 |
> |
zhat = d(3) * riji |
| 349 |
|
|
| 350 |
< |
drdxj = d(1) * riji |
| 351 |
< |
drdyj = d(2) * riji |
| 352 |
< |
drdzj = d(3) * riji |
| 350 |
> |
drdxj = xhat |
| 351 |
> |
drdyj = yhat |
| 352 |
> |
drdzj = zhat |
| 353 |
|
|
| 354 |
|
!! logicals |
| 355 |
|
|
| 439 |
|
|
| 440 |
|
if (j_is_Dipole) then |
| 441 |
|
|
| 442 |
< |
ri2 = riji * riji |
| 443 |
< |
ri3 = ri2 * riji |
| 442 |
> |
if (j_is_SplitDipole) then |
| 443 |
> |
BigR = sqrt(r2 + 0.25_dp * d_j * d_j) |
| 444 |
> |
ri = 1.0_dp / BigR |
| 445 |
> |
scale = rij * ri |
| 446 |
> |
else |
| 447 |
> |
ri = riji |
| 448 |
> |
scale = 1.0_dp |
| 449 |
> |
endif |
| 450 |
|
|
| 451 |
+ |
ri2 = ri * ri |
| 452 |
+ |
ri3 = ri2 * ri |
| 453 |
+ |
sc2 = scale * scale |
| 454 |
+ |
|
| 455 |
|
pref = pre12 * q_i * mu_j |
| 456 |
< |
vterm = pref * ct_j * riji * riji |
| 456 |
> |
vterm = pref * ct_j * ri2 * scale |
| 457 |
|
vpair = vpair + vterm |
| 458 |
|
epot = epot + sw * vterm |
| 459 |
|
|
| 460 |
< |
dudx = dudx + pref * sw * ri3 * ( ul_j(1) + 3.0d0 * ct_j * drdxj) |
| 461 |
< |
dudy = dudy + pref * sw * ri3 * ( ul_j(2) + 3.0d0 * ct_j * drdyj) |
| 462 |
< |
dudz = dudz + pref * sw * ri3 * ( ul_j(3) + 3.0d0 * ct_j * drdzj) |
| 460 |
> |
!! this has a + sign in the () because the rij vector is |
| 461 |
> |
!! r_j - r_i and the charge-dipole potential takes the origin |
| 462 |
> |
!! as the point dipole, which is atom j in this case. |
| 463 |
|
|
| 464 |
< |
dudujx = dudujx - pref * sw * ri2 * drdxj |
| 465 |
< |
dudujy = dudujy - pref * sw * ri2 * drdyj |
| 466 |
< |
dudujz = dudujz - pref * sw * ri2 * drdzj |
| 464 |
> |
dudx = dudx + pref * sw * ri3 * ( ul_j(1) + 3.0d0*ct_j*xhat*sc2) |
| 465 |
> |
dudy = dudy + pref * sw * ri3 * ( ul_j(2) + 3.0d0*ct_j*yhat*sc2) |
| 466 |
> |
dudz = dudz + pref * sw * ri3 * ( ul_j(3) + 3.0d0*ct_j*zhat*sc2) |
| 467 |
> |
|
| 468 |
> |
dudujx = dudujx - pref * sw * ri2 * xhat * scale |
| 469 |
> |
dudujy = dudujy - pref * sw * ri2 * yhat * scale |
| 470 |
> |
dudujz = dudujz - pref * sw * ri2 * zhat * scale |
| 471 |
|
|
| 472 |
|
endif |
| 473 |
+ |
|
| 474 |
|
endif |
| 475 |
|
|
| 476 |
|
if (i_is_Dipole) then |
| 477 |
|
|
| 478 |
|
if (j_is_Charge) then |
| 479 |
|
|
| 480 |
< |
ri2 = riji * riji |
| 481 |
< |
ri3 = ri2 * riji |
| 480 |
> |
if (i_is_SplitDipole) then |
| 481 |
> |
BigR = sqrt(r2 + 0.25_dp * d_i * d_i) |
| 482 |
> |
ri = 1.0_dp / BigR |
| 483 |
> |
scale = rij * ri |
| 484 |
> |
else |
| 485 |
> |
ri = riji |
| 486 |
> |
scale = 1.0_dp |
| 487 |
> |
endif |
| 488 |
|
|
| 489 |
+ |
ri2 = ri * ri |
| 490 |
+ |
ri3 = ri2 * ri |
| 491 |
+ |
sc2 = scale * scale |
| 492 |
+ |
|
| 493 |
|
pref = pre12 * q_j * mu_i |
| 494 |
< |
vterm = pref * ct_i * riji * riji |
| 494 |
> |
vterm = pref * ct_i * ri2 * scale |
| 495 |
|
vpair = vpair + vterm |
| 496 |
|
epot = epot + sw * vterm |
| 497 |
|
|
| 498 |
< |
dudx = dudx + pref * sw * ri3 * ( ul_i(1) - 3.0d0 * ct_i * drdxj) |
| 499 |
< |
dudy = dudy + pref * sw * ri3 * ( ul_i(2) - 3.0d0 * ct_i * drdyj) |
| 500 |
< |
dudz = dudz + pref * sw * ri3 * ( ul_i(3) - 3.0d0 * ct_i * drdzj) |
| 498 |
> |
dudx = dudx + pref * sw * ri3 * ( ul_i(1) - 3.0d0 * ct_i * xhat*sc2) |
| 499 |
> |
dudy = dudy + pref * sw * ri3 * ( ul_i(2) - 3.0d0 * ct_i * yhat*sc2) |
| 500 |
> |
dudz = dudz + pref * sw * ri3 * ( ul_i(3) - 3.0d0 * ct_i * zhat*sc2) |
| 501 |
|
|
| 502 |
< |
duduix = duduix + pref * sw * ri2 * drdxj |
| 503 |
< |
duduiy = duduiy + pref * sw * ri2 * drdyj |
| 504 |
< |
duduiz = duduiz + pref * sw * ri2 * drdzj |
| 502 |
> |
duduix = duduix + pref * sw * ri2 * xhat * scale |
| 503 |
> |
duduiy = duduiy + pref * sw * ri2 * yhat * scale |
| 504 |
> |
duduiz = duduiz + pref * sw * ri2 * zhat * scale |
| 505 |
|
endif |
| 506 |
|
|
| 507 |
|
if (j_is_Dipole) then |
| 508 |
|
|
| 509 |
+ |
if (i_is_SplitDipole) then |
| 510 |
+ |
if (j_is_SplitDipole) then |
| 511 |
+ |
BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j) |
| 512 |
+ |
else |
| 513 |
+ |
BigR = sqrt(r2 + 0.25_dp * d_i * d_i) |
| 514 |
+ |
endif |
| 515 |
+ |
ri = 1.0_dp / BigR |
| 516 |
+ |
scale = rij * ri |
| 517 |
+ |
else |
| 518 |
+ |
if (j_is_SplitDipole) then |
| 519 |
+ |
BigR = sqrt(r2 + 0.25_dp * d_j * d_j) |
| 520 |
+ |
ri = 1.0_dp / BigR |
| 521 |
+ |
scale = rij * ri |
| 522 |
+ |
else |
| 523 |
+ |
ri = riji |
| 524 |
+ |
scale = 1.0_dp |
| 525 |
+ |
endif |
| 526 |
+ |
endif |
| 527 |
+ |
|
| 528 |
|
ct_ij = ul_i(1)*ul_j(1) + ul_i(2)*ul_j(2) + ul_i(3)*ul_j(3) |
| 529 |
< |
ri2 = riji * riji |
| 530 |
< |
ri3 = ri2 * riji |
| 529 |
> |
|
| 530 |
> |
ri2 = ri * ri |
| 531 |
> |
ri3 = ri2 * ri |
| 532 |
|
ri4 = ri2 * ri2 |
| 533 |
+ |
sc2 = scale * scale |
| 534 |
|
|
| 535 |
|
pref = pre22 * mu_i * mu_j |
| 536 |
< |
vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j) |
| 536 |
> |
vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2) |
| 537 |
|
vpair = vpair + vterm |
| 538 |
|
epot = epot + sw * vterm |
| 539 |
|
|
| 540 |
< |
a1 = 5.0d0 * ct_i * ct_j - ct_ij |
| 540 |
> |
a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij |
| 541 |
|
|
| 542 |
< |
dudx = dudx + pref*sw*3.0d0*ri4*(a1*drdxj-ct_i*ul_j(1)-ct_j*ul_i(1)) |
| 543 |
< |
dudy = dudy + pref*sw*3.0d0*ri4*(a1*drdyj-ct_i*ul_j(2)-ct_j*ul_i(2)) |
| 544 |
< |
dudz = dudz + pref*sw*3.0d0*ri4*(a1*drdzj-ct_i*ul_j(3)-ct_j*ul_i(3)) |
| 542 |
> |
dudx=dudx+pref*sw*3.0d0*ri4*scale*(a1*xhat-ct_i*ul_j(1)-ct_j*ul_i(1)) |
| 543 |
> |
dudy=dudy+pref*sw*3.0d0*ri4*scale*(a1*yhat-ct_i*ul_j(2)-ct_j*ul_i(2)) |
| 544 |
> |
dudz=dudz+pref*sw*3.0d0*ri4*scale*(a1*zhat-ct_i*ul_j(3)-ct_j*ul_i(3)) |
| 545 |
|
|
| 546 |
< |
duduix = duduix + pref*sw*ri3*(ul_j(1) - 3.0d0*ct_j*drdxj) |
| 547 |
< |
duduiy = duduiy + pref*sw*ri3*(ul_j(2) - 3.0d0*ct_j*drdyj) |
| 548 |
< |
duduiz = duduiz + pref*sw*ri3*(ul_j(3) - 3.0d0*ct_j*drdzj) |
| 546 |
> |
duduix = duduix + pref*sw*ri3*(ul_j(1) - 3.0d0*ct_j*xhat*sc2) |
| 547 |
> |
duduiy = duduiy + pref*sw*ri3*(ul_j(2) - 3.0d0*ct_j*yhat*sc2) |
| 548 |
> |
duduiz = duduiz + pref*sw*ri3*(ul_j(3) - 3.0d0*ct_j*zhat*sc2) |
| 549 |
|
|
| 550 |
< |
dudujx = dudujx + pref*sw*ri3*(ul_i(1) - 3.0d0*ct_i*drdxj) |
| 551 |
< |
dudujy = dudujy + pref*sw*ri3*(ul_i(2) - 3.0d0*ct_i*drdyj) |
| 552 |
< |
dudujz = dudujz + pref*sw*ri3*(ul_i(3) - 3.0d0*ct_i*drdzj) |
| 550 |
> |
dudujx = dudujx + pref*sw*ri3*(ul_i(1) - 3.0d0*ct_i*xhat*sc2) |
| 551 |
> |
dudujy = dudujy + pref*sw*ri3*(ul_i(2) - 3.0d0*ct_i*yhat*sc2) |
| 552 |
> |
dudujz = dudujz + pref*sw*ri3*(ul_i(3) - 3.0d0*ct_i*zhat*sc2) |
| 553 |
|
endif |
| 554 |
|
|
| 555 |
|
endif |
| 626 |
|
end subroutine doElectrostaticPair |
| 627 |
|
|
| 628 |
|
end module electrostatic_module |
| 580 |
– |
|
