subroutine rotation_geomgrad(eel, E, Egrd, grad)
use mod_memory, only: ip, rp
use mod_electrostatics, only: ommp_electrostatics_type
implicit none
type(ommp_electrostatics_type), intent(in) :: eel
real(rp), intent(in) :: E(3, eel%top%mm_atoms), Egrd(6, eel%top%mm_atoms)
real(rp), dimension(3, eel%top%mm_atoms), intent(inout) :: grad
integer(ip) :: j, jx, jy, jz, k, l, m, n
real(rp), dimension(3) :: dip
real(rp), dimension(3,3) :: r, rt, qua, rqua, ddip
real(rp), dimension(3,3,3) :: dri, driz, drix, driy, dqua, dtmp
logical :: frozen_j, frozen_jx, frozen_jy, frozen_jz
! TODO prepare_fixedelec
! loop over the mm sites and build the derivatives of the rotation
! matrices with respect to the positions of all the relevant atoms.
do j = 1, eel%top%mm_atoms
jz = eel%iz(j)
if(jz == 0) jz = j
jx = eel%ix(j)
if(jx == 0) jx = j
jy = eel%iy(j)
if(jy == 0) jy = j
if(eel%top%use_frozen) then
if(eel%top%frozen(j) .and. &
eel%top%frozen(jx) .and. &
eel%top%frozen(jy) .and. &
eel%top%frozen(jz)) cycle
end if
frozen_j = .false.
frozen_jx = .false.
frozen_jy = .false.
frozen_jz = .false.
if(eel%top%use_frozen) then
frozen_j = eel%top%frozen(j)
frozen_jx = eel%top%frozen(jx)
frozen_jy = eel%top%frozen(jy)
frozen_jz = eel%top%frozen(jz)
end if
call rotation_matrix(.true., &
eel%top%cmm(:,j), eel%top%cmm(:,jx), &
eel%top%cmm(:,jy), eel%top%cmm(:,jz), &
eel%mol_frame(j), &
r, dri, driz, drix, driy)
rt = transpose(r)
! get the multipoles. we also need the rotated quadrupoles:
dip(_x_) = eel%q0(1+_x_,j)
dip(_y_) = eel%q0(1+_y_,j)
dip(_z_) = eel%q0(1+_z_,j)
qua(_x_,_x_) = eel%q0(4+_xx_,j)
qua(_x_,_y_) = eel%q0(4+_xy_,j)
qua(_x_,_z_) = eel%q0(4+_xz_,j)
qua(_y_,_x_) = eel%q0(4+_yx_,j)
qua(_y_,_y_) = eel%q0(4+_yy_,j)
qua(_y_,_z_) = eel%q0(4+_yz_,j)
qua(_z_,_x_) = eel%q0(4+_zx_,j)
qua(_z_,_y_) = eel%q0(4+_zy_,j)
qua(_z_,_z_) = eel%q0(4+_zz_,j)
rqua(_x_,_x_) = eel%q(4+_xx_,j)
rqua(_x_,_y_) = eel%q(4+_xy_,j)
rqua(_x_,_z_) = eel%q(4+_xz_,j)
rqua(_y_,_x_) = eel%q(4+_yx_,j)
rqua(_y_,_y_) = eel%q(4+_yy_,j)
rqua(_y_,_z_) = eel%q(4+_yz_,j)
rqua(_z_,_x_) = eel%q(4+_zx_,j)
rqua(_z_,_y_) = eel%q(4+_zy_,j)
rqua(_z_,_z_) = eel%q(4+_zz_,j)
! contributions to the forces on the j-th atoms:
ddip = 0.0_rp
dqua = 0.0_rp
dtmp = 0.0_rp
! compute the differentiated multipoles:
do k = 1, 3
do l = 1, 3
ddip(:,k) = ddip(:,k) + dri(:,k,l)*dip(l)
end do
end do
do k = 1, 3
do l = 1, 3
do m = 1, 3
do n = 1, 3
!dtmp(:,l,m)*rt(m,k)
dqua(:,k,l) = dqua(:,k,l) + (dri(:,k,m)*qua(m,n) - &
rqua(k,m)*dri(:,m,n))*rt(n,l)
end do
end do
end do
end do
if(.not. frozen_j) then
! increment the forces for the dipoles...
grad(:,j) = grad(:,j) - ddip(:,_x_)*E(_x_,j) &
- ddip(:,_y_)*E(_y_,j) &
- ddip(:,_z_)*E(_z_,j)
! ... and for the quadrupoles:
grad(:,j) = grad(:,j) + dqua(:,_x_,_x_)*Egrd(_xx_,j) &
+ dqua(:,_y_,_y_)*Egrd(_yy_,j) &
+ dqua(:,_z_,_z_)*Egrd(_zz_,j) &
+ 2*(dqua(:,_x_,_y_)*Egrd(_xy_,j) &
+ dqua(:,_x_,_z_)*Egrd(_xz_,j) &
+ dqua(:,_y_,_z_)*Egrd(_yz_,j))
end if
! do jx
ddip = 0.0_rp
dqua = 0.0_rp
dtmp = 0.0_rp
do k = 1, 3
do l = 1, 3
ddip(:,k) = ddip(:,k) + drix(:,k,l)*dip(l)
end do
end do
do k = 1, 3
do l = 1, 3
do m = 1, 3
do n = 1, 3
!dtmp(:,l,m)*rt(m,k)
dqua(:,k,l) = dqua(:,k,l) + (drix(:,k,m)*qua(m,n) -&
rqua(k,m)*drix(:,m,n))*rt(n,l)
end do
end do
end do
end do
if(.not. frozen_jx) then
! increment the forces for the dipoles...
grad(:,jx) = grad(:,jx) - ddip(:,_x_)*E(_x_,j) &
- ddip(:,_y_)*E(_y_,j) &
- ddip(:,_z_)*E(_z_,j)
! ... and for the quadrupoles:
grad(:,jx) = grad(:,jx) + dqua(:,_x_,_x_)*Egrd(_xx_,j) &
+ dqua(:,_y_,_y_)*Egrd(_yy_,j) &
+ dqua(:,_z_,_z_)*Egrd(_zz_,j) &
+ 2*(dqua(:,_x_,_y_)*Egrd(_xy_,j) &
+ dqua(:,_x_,_z_)*Egrd(_xz_,j) &
+ dqua(:,_y_,_z_)*Egrd(_yz_,j))
end if
! do jy
ddip = 0.0_rp
dqua = 0.0_rp
dtmp = 0.0_rp
do k = 1, 3
do l = 1, 3
ddip(:,k) = ddip(:,k) + driy(:,k,l)*dip(l)
end do
end do
do k = 1, 3
do l = 1, 3
do m = 1, 3
do n = 1, 3
! dtmp(:,l,m)*rt(m,k)
dqua(:,k,l) = dqua(:,k,l) + (driy(:,k,m)*qua(m,n) -&
rqua(k,m)*driy(:,m,n))*rt(n,l)
end do
end do
end do
end do
if(.not. frozen_jy) then
! increment the forces for the dipoles...
grad(:,jy) = grad(:,jy) - ddip(:,_x_)*E(_x_,j) &
- ddip(:,_y_)*E(_y_,j) &
- ddip(:,_z_)*E(_z_,j)
! ... and for the quadrupoles:
grad(:,jy) = grad(:,jy) + dqua(:,_x_,_x_)*Egrd(_xx_,j) &
+ dqua(:,_y_,_y_)*Egrd(_yy_,j) &
+ dqua(:,_z_,_z_)*Egrd(_zz_,j) &
+ 2*(dqua(:,_x_,_y_)*Egrd(_xy_,j) &
+ dqua(:,_x_,_z_)*Egrd(_xz_,j) &
+ dqua(:,_y_,_z_)*Egrd(_yz_,j))
end if
! Do jz
ddip = 0.0_rp
dqua = 0.0_rp
dtmp = 0.0_rp
do k = 1, 3
do l = 1, 3
ddip(:,k) = ddip(:,k) + driz(:,k,l)*dip(l)
end do
end do
do k = 1, 3
do l = 1, 3
do m = 1, 3
do n = 1, 3
! tmp(:,l,m)*rt(m,k)
dqua(:,k,l) = dqua(:,k,l) + (driz(:,k,m)*qua(m,n) -&
rqua(k,m)*driz(:,m,n))*rt(n,l)
end do
end do
end do
end do
if(.not. frozen_jz) then
! increment the forces for the dipoles...
grad(:,jz) = grad(:,jz) - ddip(:,_x_)*E(_x_,j) &
- ddip(:,_y_)*E(_y_,j) &
- ddip(:,_z_)*E(_z_,j)
! ... and for the quadrupoles:
grad(:,jz) = grad(:,jz) + dqua(:,_x_,_x_)*Egrd(_xx_,j) &
+ dqua(:,_y_,_y_)*Egrd(_yy_,j) &
+ dqua(:,_z_,_z_)*Egrd(_zz_,j) &
+ 2*(dqua(:,_x_,_y_)*Egrd(_xy_,j) &
+ dqua(:,_x_,_z_)*Egrd(_xz_,j) &
+ dqua(:,_y_,_z_)*Egrd(_yz_,j))
end if
end do
end subroutine rotation_geomgrad