This subroutine computes the damped coulomb kernel between two atoms. Note that this only makes sense between two MM atoms, as it is only used to compute the field that induces the point dipoles!
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
type(ommp_electrostatics_type), | intent(in) | :: | eel |
Electrostatics data structure |
||
integer(kind=ip), | intent(in) | :: | i |
Index of atoms (in MM atom list) for which the kernel should be computed |
||
integer(kind=ip), | intent(in) | :: | j |
Index of atoms (in MM atom list) for which the kernel should be computed |
||
integer(kind=ip), | intent(in) | :: | maxder |
Maximum derivative to be computed |
||
real(kind=rp), | intent(out), | dimension(maxder+1) | :: | res |
Results vector |
|
real(kind=rp), | intent(out), | dimension(3) | :: | dr |
Distance vector between i and j |
subroutine damped_coulomb_kernel(eel, i, j, maxder, res, dr)
!! This subroutine computes the damped coulomb kernel between two atoms.
!! Note that this only makes sense between two MM atoms, as it is only used
!! to compute the field that induces the point dipoles!
use mod_memory, only: ip, rp
use mod_constants, only: eps_rp
implicit none
type(ommp_electrostatics_type), intent(in) :: eel
!! Electrostatics data structure
integer(ip), intent(in) :: i, j
!! Index of atoms (in MM atom list) for which the kernel should be
!! computed
integer(ip), intent(in) :: maxder
!! Maximum derivative to be computed
real(rp), intent(out), dimension(maxder+1) :: res
!! Results vector
real(rp), intent(out), dimension(3) :: dr
!! Distance vector between i and j
real(rp) :: s, u, u3, u4, fexp, eexp
s = eel%thole(i) * eel%thole(j)
! Compute undamped kernels
dr = eel%top%cmm(:,j) - eel%top%cmm(:,i)
call coulomb_kernel(dr, maxder, res)
if(abs(s) < eps_rp) then
! either thole(i) or thole(j) are zero, so the damped kernel
! is just the same as the regular one. Job is done.
return
end if
u = 1.0_rp / (res(1) * s)
u3 = u**3
if(eel%amoeba .and. u3 * 0.39_rp < 50.0_rp) then
! 1.0 / (res(1) * s)**3 * 0.39_rp < 50_rp TODO Remove 0.39, this is a
! FF constants, ask someone why this condition is here.
! Here basically we multiply the standard interactions kernel for
! dumping coefficients. The equations implemented here correspond to
! eq. (5) of 10.1021/jp027815
fexp = -0.39_rp * u3
eexp = exp(fexp)
if(maxder >= 1) res(2) = res(2) * (1.0_rp - eexp)
if(maxder >= 2) res(3) = res(3) * (1.0_rp - (1.0_rp - fexp) * eexp)
if(maxder >= 3) res(4) = res(4) * &
(1.0_rp - (1.0_rp - fexp + 0.6_rp * fexp * fexp) * eexp)
if(maxder >= 4) res(5) = res(5) * &
(1.0_rp - (1.0_rp - fexp + 18.0_rp/35.0_rp * fexp * fexp - &
9.0_rp/35.0_rp * fexp * fexp * fexp) * eexp)
if(maxder >= 5) then
call fatal_error("Damped Coulomb kernel (AMOEBA) only supports up to the 5th derivative")
end if
else if(.not. eel%amoeba .and. res(1) > 1_rp/s) then
! TODO Again it is not clear to me why condition res(1) > 1_rp/s is here.
u4 = u3*u
if(maxder >= 1) res(2) = res(2) * (4.0_rp * u3 - 3.0_rp * u4)
if(maxder >= 2) res(3) = res(3) * u4
if(maxder >= 3) res(4) = res(4) * 0.2_rp * u4
if(maxder >= 4) then
call fatal_error("Damped Coulomb kernel (WANG) only supports up to the 4th derivative")
end if
end if
end subroutine damped_coulomb_kernel