76 |
|
public :: doElectrostaticPair |
77 |
|
public :: getCharge |
78 |
|
public :: getDipoleMoment |
79 |
+ |
public :: pre22 |
80 |
|
|
81 |
|
type :: Electrostatic |
82 |
|
integer :: c_ident |
339 |
|
real (kind=dp) :: pref, vterm, epot, dudr |
340 |
|
real (kind=dp) :: xhat, yhat, zhat |
341 |
|
real (kind=dp) :: dudx, dudy, dudz |
341 |
– |
real (kind=dp) :: drdxj, drdyj, drdzj |
342 |
|
real (kind=dp) :: scale, sc2, bigR |
343 |
|
|
344 |
|
if (.not.allocated(ElectrostaticMap)) then |
361 |
|
xhat = d(1) * riji |
362 |
|
yhat = d(2) * riji |
363 |
|
zhat = d(3) * riji |
364 |
– |
|
365 |
– |
drdxj = xhat |
366 |
– |
drdyj = yhat |
367 |
– |
drdzj = zhat |
364 |
|
|
365 |
|
!! logicals |
366 |
|
|
443 |
|
uz_j(2) = eFrame(6,atom2) |
444 |
|
uz_j(3) = eFrame(9,atom2) |
445 |
|
#endif |
446 |
< |
ct_j = uz_j(1)*drdxj + uz_j(2)*drdyj + uz_j(3)*drdzj |
446 |
> |
ct_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat |
447 |
|
|
448 |
|
if (j_is_SplitDipole) then |
449 |
|
d_j = ElectrostaticMap(me2)%split_dipole_distance |
503 |
|
|
504 |
|
dudr = - sw * vterm * riji |
505 |
|
|
506 |
< |
dudx = dudx + dudr * drdxj |
507 |
< |
dudy = dudy + dudr * drdyj |
508 |
< |
dudz = dudz + dudr * drdzj |
506 |
> |
dudx = dudx + dudr * xhat |
507 |
> |
dudy = dudy + dudr * yhat |
508 |
> |
dudz = dudz + dudr * zhat |
509 |
|
|
510 |
|
endif |
511 |
|
|
525 |
|
sc2 = scale * scale |
526 |
|
|
527 |
|
pref = pre12 * q_i * mu_j |
528 |
< |
vterm = pref * ct_j * ri2 * scale |
528 |
> |
vterm = - pref * ct_j * ri2 * scale |
529 |
|
vpair = vpair + vterm |
530 |
|
epot = epot + sw * vterm |
531 |
|
|
533 |
|
!! r_j - r_i and the charge-dipole potential takes the origin |
534 |
|
!! as the point dipole, which is atom j in this case. |
535 |
|
|
536 |
< |
dudx = dudx + pref * sw * ri3 * ( uz_j(1) + 3.0d0*ct_j*xhat*sc2) |
537 |
< |
dudy = dudy + pref * sw * ri3 * ( uz_j(2) + 3.0d0*ct_j*yhat*sc2) |
538 |
< |
dudz = dudz + pref * sw * ri3 * ( uz_j(3) + 3.0d0*ct_j*zhat*sc2) |
536 |
> |
dudx = dudx - pref * sw * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2) |
537 |
> |
dudy = dudy - pref * sw * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2) |
538 |
> |
dudz = dudz - pref * sw * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2) |
539 |
|
|
540 |
|
duduz_j(1) = duduz_j(1) - pref * sw * ri2 * xhat * scale |
541 |
|
duduz_j(2) = duduz_j(2) - pref * sw * ri2 * yhat * scale |
546 |
|
if (j_is_Quadrupole) then |
547 |
|
ri2 = riji * riji |
548 |
|
ri3 = ri2 * riji |
549 |
< |
ri4 = ri2 * ri4 |
549 |
> |
ri4 = ri2 * ri2 |
550 |
|
cx2 = cx_j * cx_j |
551 |
|
cy2 = cy_j * cy_j |
552 |
|
cz2 = cz_j * cz_j |
553 |
|
|
554 |
< |
pref = pre14 * q_i / 6.0_dp |
554 |
> |
|
555 |
> |
pref = pre14 * q_i / 3.0_dp |
556 |
|
vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + & |
557 |
|
qyy_j * (3.0_dp*cy2 - 1.0_dp) + & |
558 |
|
qzz_j * (3.0_dp*cz2 - 1.0_dp)) |
559 |
|
vpair = vpair + vterm |
560 |
|
epot = epot + sw * vterm |
561 |
|
|
562 |
< |
dudx = dudx - 5.0_dp*sw*vterm*riji*xhat - pref * sw * ri4 * ( & |
562 |
> |
dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + pref * sw * ri4 * ( & |
563 |
|
qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + & |
564 |
|
qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + & |
565 |
|
qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) ) |
566 |
< |
dudy = dudy - 5.0_dp*sw*vterm*riji*yhat - pref * sw * ri4 * ( & |
566 |
> |
dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + pref * sw * ri4 * ( & |
567 |
|
qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + & |
568 |
|
qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + & |
569 |
|
qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) ) |
570 |
< |
dudz = dudz - 5.0_dp*sw*vterm*riji*zhat - pref * sw * ri4 * ( & |
570 |
> |
dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + pref * sw * ri4 * ( & |
571 |
|
qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + & |
572 |
|
qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + & |
573 |
|
qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) ) |
673 |
|
|
674 |
|
ri2 = riji * riji |
675 |
|
ri3 = ri2 * riji |
676 |
< |
ri4 = ri2 * ri4 |
676 |
> |
ri4 = ri2 * ri2 |
677 |
|
cx2 = cx_i * cx_i |
678 |
|
cy2 = cy_i * cy_i |
679 |
|
cz2 = cz_i * cz_i |
680 |
|
|
681 |
< |
pref = pre14 * q_i / 6.0_dp |
681 |
> |
pref = pre14 * q_j / 3.0_dp |
682 |
|
vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + & |
683 |
|
qyy_i * (3.0_dp*cy2 - 1.0_dp) + & |
684 |
|
qzz_i * (3.0_dp*cz2 - 1.0_dp)) |
685 |
|
vpair = vpair + vterm |
686 |
|
epot = epot + sw * vterm |
687 |
|
|
688 |
< |
dudx = dudx - 5.0_dp*sw*vterm*riji*xhat - pref * sw * ri4 * ( & |
688 |
> |
dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + pref * sw * ri4 * ( & |
689 |
|
qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + & |
690 |
|
qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + & |
691 |
|
qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) ) |
692 |
< |
dudy = dudy - 5.0_dp*sw*vterm*riji*yhat - pref * sw * ri4 * ( & |
692 |
> |
dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + pref * sw * ri4 * ( & |
693 |
|
qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + & |
694 |
|
qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + & |
695 |
|
qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) ) |
696 |
< |
dudz = dudz - 5.0_dp*sw*vterm*riji*zhat - pref * sw * ri4 * ( & |
696 |
> |
dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + pref * sw * ri4 * ( & |
697 |
|
qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + & |
698 |
|
qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + & |
699 |
|
qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) ) |