73 |
|
|
74 |
|
!! variables to handle different summation methods for long-range electrostatics: |
75 |
|
integer, save :: summationMethod = NONE |
76 |
+ |
logical, save :: summationMethodChecked = .false. |
77 |
|
real(kind=DP), save :: defaultCutoff = 0.0_DP |
78 |
|
logical, save :: haveDefaultCutoff = .false. |
79 |
|
real(kind=DP), save :: dampingAlpha = 0.0_DP |
83 |
|
real(kind=DP), save :: constERFC = 0.0_DP |
84 |
|
real(kind=DP), save :: constEXP = 0.0_DP |
85 |
|
logical, save :: haveDWAconstants = .false. |
86 |
+ |
real(kind=dp), save :: rcuti = 0.0_dp |
87 |
+ |
real(kind=dp), save :: rcuti2 = 0.0_dp |
88 |
+ |
real(kind=dp), save :: rcuti3 = 0.0_dp |
89 |
+ |
real(kind=dp), save :: rcuti4 = 0.0_dp |
90 |
|
|
91 |
|
|
92 |
|
public :: setElectrostaticSummationMethod |
360 |
|
end function getDipoleMoment |
361 |
|
|
362 |
|
subroutine checkSummationMethod() |
363 |
+ |
|
364 |
+ |
if (.not.haveDefaultCutoff) then |
365 |
+ |
call handleError("checkSummationMethod", "no Default Cutoff set!") |
366 |
+ |
endif |
367 |
|
|
368 |
+ |
rcuti = 1.0d0 / defaultCutoff |
369 |
+ |
rcuti2 = rcuti*rcuti |
370 |
+ |
rcuti3 = rcuti2*rcuti |
371 |
+ |
rcuti4 = rcuti2*rcuti2 |
372 |
+ |
|
373 |
|
if (summationMethod .eq. DAMPED_WOLF) then |
374 |
|
if (.not.haveDWAconstants) then |
375 |
|
|
377 |
|
call handleError("checkSummationMethod", "no Damping Alpha set!") |
378 |
|
endif |
379 |
|
|
380 |
< |
if (.not.have....) |
381 |
< |
constEXP = |
382 |
< |
constERFC = |
380 |
> |
if (.not.haveDefaultCutoff) then |
381 |
> |
call handleError("checkSummationMethod", "no Default Cutoff set!") |
382 |
> |
endif |
383 |
> |
|
384 |
> |
constEXP = exp(-dampingAlpha*dampingAlpha*defaultCutoff*defaultCutoff) |
385 |
> |
constERFC = erfc(dampingAlpha*defaultCutoff) |
386 |
|
|
387 |
|
haveDWAconstants = .true. |
388 |
+ |
endif |
389 |
+ |
endif |
390 |
+ |
|
391 |
+ |
if (summationMethod .eq. REACTION_FIELD) then |
392 |
+ |
if (.not.haveDielectric) then |
393 |
+ |
call handleError("checkSummationMethod", "no reaction field Dielectric set!") |
394 |
|
endif |
395 |
|
endif |
396 |
|
|
397 |
+ |
summationMethodChecked = .true. |
398 |
|
end subroutine checkSummationMethod |
399 |
|
|
400 |
|
|
401 |
|
|
402 |
|
subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, & |
403 |
< |
vpair, fpair, pot, eFrame, f, t, do_pot, corrMethod, rcuti) |
403 |
> |
vpair, fpair, pot, eFrame, f, t, do_pot) |
404 |
|
|
405 |
|
logical, intent(in) :: do_pot |
406 |
|
|
407 |
|
integer, intent(in) :: atom1, atom2 |
408 |
|
integer :: localError |
385 |
– |
integer, intent(in) :: corrMethod |
409 |
|
|
410 |
< |
real(kind=dp), intent(in) :: rij, r2, sw, rcuti |
410 |
> |
real(kind=dp), intent(in) :: rij, r2, sw |
411 |
|
real(kind=dp), intent(in), dimension(3) :: d |
412 |
|
real(kind=dp), intent(inout) :: vpair |
413 |
|
real(kind=dp), intent(inout), dimension(3) :: fpair |
438 |
|
real (kind=dp) :: xhat, yhat, zhat |
439 |
|
real (kind=dp) :: dudx, dudy, dudz |
440 |
|
real (kind=dp) :: scale, sc2, bigR, switcher, dswitcher |
418 |
– |
real (kind=dp) :: rcuti2, rcuti3, rcuti4 |
441 |
|
|
442 |
|
if (.not.allocated(ElectrostaticMap)) then |
443 |
|
call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!") |
464 |
|
xhat = d(1) * riji |
465 |
|
yhat = d(2) * riji |
466 |
|
zhat = d(3) * riji |
445 |
– |
|
446 |
– |
rcuti2 = rcuti*rcuti |
447 |
– |
rcuti3 = rcuti2*rcuti |
448 |
– |
rcuti4 = rcuti2*rcuti2 |
467 |
|
|
468 |
|
swi = 1.0d0 / sw |
469 |
|
|
610 |
|
|
611 |
|
if (j_is_Charge) then |
612 |
|
|
613 |
< |
if (corrMethod .eq. 1) then |
613 |
> |
if (summationMethod .eq. 1) then |
614 |
|
vterm = pre11 * q_i * q_j * (riji - rcuti) |
615 |
|
|
616 |
|
vpair = vpair + vterm |
642 |
|
|
643 |
|
pref = sw * pre12 * q_i * mu_j |
644 |
|
|
645 |
< |
if (corrMethod .eq. 1) then |
645 |
> |
if (summationMethod .eq. 1) then |
646 |
|
ri2 = riji * riji |
647 |
|
ri3 = ri2 * riji |
648 |
|
|
709 |
|
|
710 |
|
pref = sw * pre14 * q_i / 3.0_dp |
711 |
|
|
712 |
< |
if (corrMethod .eq. 1) then |
712 |
> |
if (summationMethod .eq. 1) then |
713 |
|
vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + & |
714 |
|
qyy_j * (3.0_dp*cy2 - 1.0_dp) + & |
715 |
|
qzz_j * (3.0_dp*cz2 - 1.0_dp) ) |
807 |
|
|
808 |
|
pref = sw * pre12 * q_j * mu_i |
809 |
|
|
810 |
< |
if (corrMethod .eq. 1) then |
810 |
> |
if (summationMethod .eq. 1) then |
811 |
|
ri2 = riji * riji |
812 |
|
ri3 = ri2 * riji |
813 |
|
|
862 |
|
|
863 |
|
pref = sw * pre22 * mu_i * mu_j |
864 |
|
|
865 |
< |
if (corrMethod .eq. 1) then |
865 |
> |
if (summationMethod .eq. 1) then |
866 |
|
ri2 = riji * riji |
867 |
|
ri3 = ri2 * riji |
868 |
|
ri4 = ri2 * ri2 |
965 |
|
|
966 |
|
pref = sw * pre14 * q_j / 3.0_dp |
967 |
|
|
968 |
< |
if (corrMethod .eq. 1) then |
968 |
> |
if (summationMethod .eq. 1) then |
969 |
|
vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + & |
970 |
|
qyy_i * (3.0_dp*cy2 - 1.0_dp) + & |
971 |
|
qzz_i * (3.0_dp*cz2 - 1.0_dp) ) |