this routine rotates the atomic multipoles from the molecular frame where they are defined as force field parameters to the lab frame. if required, it also computes the contribution to the forces that stems from the derivatives of the rotation matrices, sometimes referred to as "torques". for the latter task, it uses the field and field gradient from at the multipoles, which is passed in def. for consistency reasons, def is dimensioned (ld_cder,mm_atoms).
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
type(ommp_electrostatics_type), | intent(inout) | :: | eel |
subroutine rotate_multipoles(eel)
!! this routine rotates the atomic multipoles from the molecular frame
!! where they are defined as force field parameters to the lab frame.
!! if required, it also computes the contribution to the forces that
!! stems from the derivatives of the rotation matrices, sometimes
!! referred to as "torques".
!! for the latter task, it uses the field and field gradient from
!! at the multipoles, which is passed in def.
!! for consistency reasons, def is dimensioned (ld_cder,mm_atoms).
use mod_memory, only: ip, rp
use mod_electrostatics, only: ommp_electrostatics_type
implicit none
type(ommp_electrostatics_type), intent(inout) :: eel
integer(ip) :: j, jx, jy, jz
real(rp), dimension(3,3) :: r, rt, qua, rqua, tmp
real(rp), dimension(3,3,3) :: dri, driz, drix, driy
! loop over the mm sites and build the rotation matrices.
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
call rotation_matrix(.false., &
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)
! copy the monopole:
eel%q(1,j) = eel%q0(1,j)
! rotate the dipole
eel%q(2:4,j) = matmul(r,eel%q0(2:4,j))
! exctract, rotate and put back the quadrupole:
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)
tmp = matmul(r,qua)
rqua = matmul(tmp,rt)
eel%q(4+_xx_,j) = rqua(_x_,_x_)
eel%q(4+_yy_,j) = rqua(_y_,_y_)
eel%q(4+_zz_,j) = rqua(_z_,_z_)
eel%q(4+_xy_,j) = rqua(_x_,_y_)
eel%q(4+_xz_,j) = rqua(_x_,_z_)
eel%q(4+_yz_,j) = rqua(_y_,_z_)
end do
end subroutine rotate_multipoles