| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=rp), | intent(in) | :: | quad(6) |
point quadrupole stored as (xx, xy, yy, xz, yz, zz) |
||
| real(kind=rp), | intent(in) | :: | dr(3) |
Distance vector |
||
| real(kind=rp), | intent(in) | :: | kernel(:) |
Array of coulomb kernel (either damped or undamped) |
||
| logical, | intent(in) | :: | do_V |
Flags to enable/disable calculation of different electrostatic properties |
||
| real(kind=rp), | intent(inout) | :: | V |
Electric potential |
||
| logical, | intent(in) | :: | do_E |
Flags to enable/disable calculation of different electrostatic properties |
||
| real(kind=rp), | intent(inout) | :: | E(3) |
Electric potential |
||
| logical, | intent(in) | :: | do_grdE |
Flags to enable/disable calculation of different electrostatic properties |
||
| real(kind=rp), | intent(inout) | :: | grdE(6) |
Electric potential |
||
| logical, | intent(in) | :: | do_HE |
Flags to enable/disable calculation of different electrostatic properties |
||
| real(kind=rp), | intent(inout) | :: | HE(10) |
Electric potential |
subroutine quad_elec_prop(quad, dr, kernel, &
do_V, V, do_E, E, do_grdE, grdE, do_HE, HE)
implicit none
real(rp), intent(in) :: quad(6)
!! point quadrupole stored as (xx, xy, yy, xz, yz, zz)
real(rp), intent(in) :: dr(3)
!! Distance vector
real(rp), intent(in) :: kernel(:)
!! Array of coulomb kernel (either damped or undamped)
logical, intent(in) :: do_V, do_E, do_grdE, do_HE
!! Flags to enable/disable calculation of different electrostatic
!! properties
real(rp), intent(inout) :: V, E(3), grdE(6), HE(10)
!! Electric potential
real(rp) :: quadxr(3), quadxr_dot_r
quadxr(1) = quad(1)*dr(1) + quad(2)*dr(2) + quad(4)*dr(3)
quadxr(2) = quad(2)*dr(1) + quad(3)*dr(2) + quad(5)*dr(3)
quadxr(3) = quad(4)*dr(1) + quad(5)*dr(2) + quad(6)*dr(3)
quadxr_dot_r = dot_product(quadxr, dr)
if(do_V) then
V = V + 3.0_rp * quadxr_dot_r * kernel(3)
end if
if(do_E) then
E = E + 15.0_rp * quadxr_dot_r * dr * kernel(4) &
- 6.0_rp * quadxr * kernel(3)
end if
if(do_grdE) then
! xx
grdE(1) = grdE(1) + 105.0 * kernel(5) * quadxr_dot_r * dr(1)*dr(1) + &
6.0 * kernel(3) * quad(1) - &
15.0*kernel(4)*(quadxr_dot_r+4.0*quadxr(1)*dr(1))
! xy
grdE(2) = grdE(2) + 105.0 * kernel(5) * quadxr_dot_r * dr(1)*dr(2) + &
6.0 * kernel(3) * quad(2) - &
30.0*kernel(4)*(quadxr(1)*dr(2)+quadxr(2)*dr(1))
! yy
grdE(3) = grdE(3) + 105.0 * kernel(5) * quadxr_dot_r * dr(2)*dr(2) + &
6.0 * kernel(3) * quad(3) - &
15.0*kernel(4)*(quadxr_dot_r+4.0*quadxr(2)*dr(2))
! xz
grdE(4) = grdE(4) + 105.0 * kernel(5) * quadxr_dot_r * dr(1)*dr(3) + &
6.0 * kernel(3) * quad(4) - &
30.0*kernel(4)*(quadxr(1)*dr(3)+quadxr(3)*dr(1))
! yz
grdE(5) = grdE(5) + 105.0 * kernel(5) * quadxr_dot_r * dr(2)*dr(3) + &
6.0 * kernel(3) * quad(5) - &
30.0*kernel(4)*(quadxr(2)*dr(3)+quadxr(3)*dr(2))
! zz
grdE(6) = grdE(6) + 105.0 * kernel(5) * quadxr_dot_r * dr(3)*dr(3) + &
6.0 * kernel(3) * quad(6) - &
15.0*kernel(4)*(quadxr_dot_r+4.0*quadxr(3)*dr(3))
end if
if(do_HE) then
! 3
HE(_xxx_) = HE(_xxx_) + 945*kernel(6)*dr(_x_)**3 * quadxr_dot_r &
- 315*kernel(5)*dr(_x_)*(2*quadxr(_x_)*dr(_x_) + quadxr_dot_r) &
+ 90*kernel(4)*(quad(_xx_)*dr(_x_) + quadxr(_x_))
HE(_yyy_) = HE(_yyy_) + 945*kernel(6)*dr(_y_)**3 * quadxr_dot_r &
- 315*kernel(5)*dr(_y_)*(2*quadxr(_y_)*dr(_y_) + quadxr_dot_r) &
+ 90*kernel(4)*(quad(_yy_)*dr(_y_) + quadxr(_y_))
HE(_zzz_) = HE(_zzz_) + 945*kernel(6)*dr(_z_)**3 * quadxr_dot_r &
- 315*kernel(5)*dr(_z_)*(2*quadxr(_z_)*dr(_z_) + quadxr_dot_r) &
+ 90*kernel(4)*(quad(_zz_)*dr(_z_) + quadxr(_z_))
! 2 + 1
HE(_xxy_) = HE(_xxy_) + 945*kernel(6)*dr(_x_)**2*dr(_y_)*quadxr_dot_r &
- 105*kernel(5)*(4*quadxr(_x_)*dr(_x_)*dr(_y_) + &
2*quadxr(_y_)*dr(_x_)*dr(_x_) + &
dr(_y_)*quadxr_dot_r) &
+ 30*kernel(4)*(quad(_xx_)*dr(_y_) + 2*quad(_xy_)*dr(_x_) + quadxr(_y_))
HE(_xxz_) = HE(_xxz_) + 945*kernel(6)*dr(_x_)**2*dr(_z_)*quadxr_dot_r &
- 105*kernel(5)*(4*quadxr(_x_)*dr(_x_)*dr(_z_) + &
2*quadxr(_z_)*dr(_x_)*dr(_x_) + &
dr(_z_)*quadxr_dot_r) &
+ 30*kernel(4)*(quad(_xx_)*dr(_z_) + 2*quad(_xz_)*dr(_x_) + quadxr(_z_))
HE(_yyx_) = HE(_yyx_) + 945*kernel(6)*dr(_y_)**2*dr(_x_)*quadxr_dot_r &
- 105*kernel(5)*(4*quadxr(_y_)*dr(_y_)*dr(_x_) + &
2*quadxr(_x_)*dr(_y_)*dr(_y_) + &
dr(_x_)*quadxr_dot_r) &
+ 30*kernel(4)*(quad(_yy_)*dr(_x_) + 2*quad(_xy_)*dr(_y_) + quadxr(_x_))
HE(_yyz_) = HE(_yyz_) + 945*kernel(6)*dr(_y_)**2*dr(_z_)*quadxr_dot_r &
- 105*kernel(5)*(4*quadxr(_y_)*dr(_y_)*dr(_z_) + &
2*quadxr(_z_)*dr(_y_)*dr(_y_) + &
dr(_z_)*quadxr_dot_r) &
+ 30*kernel(4)*(quad(_yy_)*dr(_z_) + 2*quad(_zy_)*dr(_y_) + quadxr(_z_))
HE(_zzx_) = HE(_zzx_) + 945*kernel(6)*dr(_z_)**2*dr(_x_)*quadxr_dot_r &
- 105*kernel(5)*(4*quadxr(_z_)*dr(_z_)*dr(_x_) + &
2*quadxr(_x_)*dr(_z_)*dr(_z_) + &
dr(_x_)*quadxr_dot_r) &
+ 30*kernel(4)*(quad(_zz_)*dr(_x_) + 2*quad(_zx_)*dr(_z_) + quadxr(_x_))
HE(_zzy_) = HE(_zzy_) + 945*kernel(6)*dr(_z_)**2*dr(_y_)*quadxr_dot_r &
- 105*kernel(5)*(4*quadxr(_z_)*dr(_z_)*dr(_y_) + &
2*quadxr(_y_)*dr(_z_)*dr(_z_) + &
dr(_y_)*quadxr_dot_r) &
+ 30*kernel(4)*(quad(_zz_)*dr(_y_) + 2*quad(_zy_)*dr(_z_) + quadxr(_y_))
! 1 + 1 + 1
HE(_xyz_) = HE(_xyz_) + 945*kernel(6)*dr(_x_)*dr(_y_)*dr(_z_)*quadxr_dot_r &
- 210*kernel(5)*(quadxr(_x_)*dr(_y_)*dr(_z_) + &
quadxr(_y_)*dr(_x_)*dr(_z_) + &
quadxr(_z_)*dr(_x_)*dr(_y_)) &
+ 30*kernel(4)*(quad(_xy_)*dr(_z_) + quad(_xz_)*dr(_y_) + quad(_yz_)*dr(_x_))
end if
end subroutine quad_elec_prop