Computes the electric potential, field and field gradients of static multipoles at all sites (polarizable sites are a subset of static ones)
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(ommp_electrostatics_type), | intent(inout) | :: | eel |
Electrostatics data structure |
||
| logical, | intent(in) | :: | do_V |
Flags to enable/disable the calculation of different components |
||
| logical, | intent(in) | :: | do_E |
Flags to enable/disable the calculation of different components |
||
| logical, | intent(in) | :: | do_Egrd |
Flags to enable/disable the calculation of different components |
||
| logical, | intent(in) | :: | do_EHes |
Flags to enable/disable the calculation of different components |
subroutine elec_prop_M2M(eel, do_V, do_E, do_Egrd, do_EHes)
!! Computes the electric potential, field and field gradients of
!! static multipoles at all sites (polarizable sites are a
!! subset of static ones)
implicit none
type(ommp_electrostatics_type), intent(inout) :: eel
!! Electrostatics data structure
logical, intent(in) :: do_V, do_E, do_Egrd, do_EHes
!! Flags to enable/disable the calculation of different components
real(rp) :: kernel(6), dr(3), tmpV, tmpE(3), tmpEgr(6), tmpHE(10), scalf
integer(ip) :: i, j, idx, ikernel
logical :: to_do, to_scale
type(ommp_topology_type), pointer :: top
top => eel%top
if(do_EHes) then
ikernel = 3
elseif(do_Egrd) then
ikernel = 2
elseif(do_E) then
ikernel = 1
elseif(do_V) then
ikernel = 0
else
return
end if
if(eel%amoeba) ikernel = ikernel + 2
if(eel%amoeba) then
!$omp parallel do default(shared) schedule(dynamic) &
!$omp private(i,j,idx,to_do,to_scale,scalf,dr,kernel,tmpV,tmpE,tmpEgr,tmpHE)
do j=1, top%mm_atoms
! loop on sources
do i=1, top%mm_atoms
if(j == i) cycle
!loop on target
to_do = .true.
to_scale = .false.
scalf = 1.0
! Check if the element should be scaled
do idx=eel%list_S_S%ri(i), eel%list_S_S%ri(i+1)-1
if(eel%list_S_S%ci(idx) == j) then
to_scale = .true.
exit
end if
end do
!If it should set the correct variables
if(to_scale) then
to_do = eel%todo_S_S(idx)
scalf = eel%scalef_S_S(idx)
end if
if(to_do) then
dr = top%cmm(:,j) - top%cmm(:, i)
call coulomb_kernel(dr, ikernel, kernel)
if(do_V) tmpV = 0.0_rp
if(do_E) tmpE = 0.0_rp
if(do_Egrd) tmpEgr = 0.0_rp
if(do_EHes) tmpHE = 0.0_rp
call q_elec_prop(eel%q(1,i), dr, kernel, &
do_V, tmpV, &
do_E, tmpE, &
do_Egrd, tmpEgr, &
do_EHes, tmpHE)
call mu_elec_prop(eel%q(2:4,i), dr, kernel, &
do_V, tmpV, &
do_E, tmpE, &
do_Egrd, tmpEgr, &
do_EHes, tmpHE)
call quad_elec_prop(eel%q(5:10,i), dr, kernel, &
do_V, tmpV, &
do_E, tmpE, &
do_Egrd, tmpEgr, &
do_EHes, tmpHE)
if(to_scale) then
if(do_V) eel%V_M2M(j) = eel%V_M2M(j) + tmpV * scalf
if(do_E) eel%E_M2M(:,j) = eel%E_M2M(:,j) + tmpE * scalf
if(do_Egrd) eel%Egrd_M2M(:,j) = eel%Egrd_M2M(:,j) + tmpEgr * scalf
if(do_EHes) eel%EHes_M2M(:,j) = eel%EHes_M2M(:,j) + tmpHE * scalf
else
if(do_V) eel%V_M2M(j) = eel%V_M2M(j) + tmpV
if(do_E) eel%E_M2M(:,j) = eel%E_M2M(:,j) + tmpE
if(do_Egrd) eel%Egrd_M2M(:,j) = eel%Egrd_M2M(:,j) + tmpEgr
if(do_EHes) eel%EHes_M2M(:,j) = eel%EHes_M2M(:,j) + tmpHE
end if
end if
end do
end do
else
!$omp parallel do default(shared) schedule(dynamic) &
!$omp private(i,j,idx,to_do,to_scale,scalf,dr,kernel,tmpV,tmpE,tmpEgr,tmpHE)
do j=1, top%mm_atoms
! loop on sources
do i=1, top%mm_atoms
if(j == i) cycle
!loop on target
to_do = .true.
to_scale = .false.
scalf = 1.0
! Check if the element should be scaled
do idx=eel%list_S_S%ri(i), eel%list_S_S%ri(i+1)-1
if(eel%list_S_S%ci(idx) == j) then
to_scale = .true.
exit
end if
end do
!If it should set the correct variables
if(to_scale) then
to_do = eel%todo_S_S(idx)
scalf = eel%scalef_S_S(idx)
end if
if(to_do) then
dr = top%cmm(:,j) - top%cmm(:, i)
call coulomb_kernel(dr, ikernel, kernel)
if(do_V) tmpV = 0.0_rp
if(do_E) tmpE = 0.0_rp
if(do_Egrd) tmpEgr = 0.0_rp
if(do_EHes) tmpHE = 0.0_rp
call q_elec_prop(eel%q(1,i), dr, kernel, &
do_V, tmpV, &
do_E, tmpE, &
do_Egrd, tmpEgr, &
do_EHes, tmpHE)
if(to_scale) then
if(do_V) eel%V_M2M(j) = eel%V_M2M(j) + tmpV * scalf
if(do_E) eel%E_M2M(:,j) = eel%E_M2M(:,j) + tmpE * scalf
if(do_Egrd) eel%Egrd_M2M(:,j) = eel%Egrd_M2M(:,j) + tmpEgr * scalf
if(do_EHes) eel%EHes_M2M(:,j) = eel%EHes_M2M(:,j) + tmpHE * scalf
else
if(do_V) eel%V_M2M(j) = eel%V_M2M(j) + tmpV
if(do_E) eel%E_M2M(:,j) = eel%E_M2M(:,j) + tmpE
if(do_Egrd) eel%Egrd_M2M(:,j) = eel%Egrd_M2M(:,j) + tmpEgr
if(do_EHes) eel%EHes_M2M(:,j) = eel%EHes_M2M(:,j) + tmpHE
end if
end if
end do
end do
end if
end subroutine elec_prop_M2M