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 |
– |
|