rotrfmt.f90

Go to the documentation of this file.

 
! Copyright (C) 2015 J. K. Dewhurst, S. Sharma and E. K. U. Gross.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.
 
subroutine rotrfmt(rot,nr,nri,rfmt1,rfmt2)
use modmain
implicit none
! arguments
real(8), intent(in) :: rot(3,3)
integer, intent(in) :: nr,nri
real(8), intent(in) :: rfmt1(*)
real(8), intent(out) :: rfmt2(*)
! local variables
integer i
! inner part of muffin-tin
call rotrflm(rot,lmaxi,nri,lmmaxi,rfmt1,rfmt2)
! outer part of muffin-tin
i=lmmaxi*nri+1
call rotrflm(rot,lmaxo,nr-nri,lmmaxo,rfmt1(i),rfmt2(i))
return
 
contains
 
!BOP
! !ROUTINE: rotrflm
! !INTERFACE:
subroutine rotrflm(rot,lmax,n,ld,rflm1,rflm2)
! !INPUT/OUTPUT PARAMETERS:
!   rot   : rotation matrix (in,real(3,3))
!   lmax  : maximum angular momentum (in,integer)
!   n     : number of functions to rotate (in,integer)
!   ld    : leading dimension (in,integer)
!   rflm1 : coefficients of the real spherical harmonic expansion for each
!           function (in,real(ld,n))
!   rflm2 : coefficients of rotated functions (out,complex(ld,n))
! !DESCRIPTION:
!   Rotates a set of real functions
!   $$ f_i({\bf r})=\sum_{lm}f_{lm}^iR_{lm}(\hat{\bf r}) $$
!   for all $i$, given the coefficients $f_{lm}^i$ and a rotation matrix $R$.
!   This is done by first the computing the Euler angles $(\alpha,\beta,\gamma)$
!   of $R^{-1}$ (see routine {\tt roteuler}) and then applying the spherical
!   harmonic rotation matrix generated by the routine {\tt rlmrot}.
!
! !REVISION HISTORY:
!   Created December 2008 (JKD)
!EOP
!BOC
implicit none
! arguments
real(8), intent(in) :: rot(3,3)
integer, intent(in) :: lmax,n,ld
real(8), intent(in) :: rflm1(ld,*)
real(8), intent(out) :: rflm2(ld,*)
! local variables
integer l,lm,nm,p
real(8) det,ang(3),angi(3)
! automatic arrays
real(8) d(ld,ld)
if (n == 0) return
! find the determinant
det=rot(1,1)*(rot(2,2)*rot(3,3)-rot(3,2)*rot(2,3)) &
   +rot(2,1)*(rot(3,2)*rot(1,3)-rot(1,2)*rot(3,3)) &
   +rot(3,1)*(rot(1,2)*rot(2,3)-rot(2,2)*rot(1,3))
! calculate the Euler angles of the proper rotation
if (det > 0.d0) then
  p=1
  call roteuler(rot,ang)
else
  p=-1
  call roteuler(-rot(:,:),ang)
end if
! inverse rotation: the function is to be rotated, not the spherical harmonics
angi(1)=-ang(3)
angi(2)=-ang(2)
angi(3)=-ang(1)
! determine the rotation matrix for real spherical harmonics
call rlmrot(p,angi,lmax,ld,d)
! apply rotation matrix
do l=0,lmax
  nm=2*l+1
  lm=l**2+1
  call dgemm('N','N',nm,n,nm,1.d0,d(lm,lm),ld,rflm1(lm,1),ld,0.d0,rflm2(lm,1), &
   ld)
end do
end subroutine
!EOC
 
!BOP
! !ROUTINE: rlmrot
! !INTERFACE:
subroutine rlmrot(p,ang,lmax,ld,d)
! !INPUT/OUTPUT PARAMETERS:
!   p    : if p=-1 then the rotation matrix is improper (in,integer)
!   ang  : Euler angles; alpha, beta, gamma (in,real(3))
!   lmax : maximum angular momentum (in,integer)
!   ld   : leading dimension (in,integer)
!   d    : real spherical harmonic rotation matrix (out,real(ld,*))
! !DESCRIPTION:
!   Returns the rotation matrix in the basis of real spherical harmonics given
!   the three Euler angles, $(\alpha,\beta,\gamma)$, and the parity, $p$, of the
!   rotation. The matrix is determined using the formula of V. V. Nechaev,
!   [{\it J. Struct. Chem.} {\bf 35}, 115 (1994)], suitably modified for our
!   definition of the real spherical harmonics ($m_1>0$, $m_2>0$):
!   \begin{align*}
!    &\Delta^l_{00}=d^l_{00}, \\
!    &\Delta^l_{m_10}=\sqrt{2}\,(-1)^{m_1}d^l_{0m_1}\cos(m_1\alpha), \\
!    &\Delta^l_{0m_2}=\sqrt{2}\,(-1)^{m_2}d^l_{m_20}\cos(m_2\gamma), \\
!    &\Delta^l_{-m_10}=-\sqrt{2}\,d^l_{0m_1}\sin(m_1\alpha), \\
!    &\Delta^l_{0-m_2}=\sqrt{2}\,d^l_{m_20}\sin(m_2\gamma), \\
!    &\Delta^l_{m_1m_2}=(-1)^{m_1}(-1)^{m_2}\{\cos(m_1\alpha)\cos(m_2\gamma)
!     [d_A+d_B]-\sin(m_1\alpha)\sin(m_2\gamma)[d_A-d_B]\}, \\
!    &\Delta^l_{m_1-m_2}=(-1)^{m_1}\{\sin(m_1\alpha)\cos(m_2\gamma)
!     [d_A-d_B]+\cos(m_1\alpha)\sin(m_2\gamma)[d_A+d_B]\}, \\
!    &\Delta^l_{-m_1m_2}=-(-1)^{m_2}\{\sin(m_1\alpha)\cos(m_2\gamma)
!     [d_A+d_B]+\cos(m_1\alpha)\sin(m_2\gamma)[d_A-d_B]\}, \\
!    &\Delta^l_{-m_1-m_2}=\cos(m_1\alpha)\cos(m_2\gamma)
!     [d_A-d_B]-\sin(m_1\alpha)\sin(m_2\gamma)[d_A+d_B],
!   \end{align*}
!   where $d_A\equiv d^l_{-m_1-m_2}$, $d_B\equiv(-1)^{m_1}d^l_{m_1-m_2}$ and
!   $d$ is the rotation matrix about the $y$-axis for complex spherical
!   harmonics. See the routines {\tt genrlm}, {\tt roteuler} and {\tt ylmroty}.
!
! !REVISION HISTORY:
!   Created December 2008 (JKD)
!EOP
!BOC
implicit none
! arguments
integer, intent(in) :: p
real(8), intent(in) :: ang(3)
integer, intent(in) :: lmax,ld
real(8), intent(out) :: d(ld,*)
! local variables
integer l,m1,m2,lm0,lm1,lm2
real(8), parameter :: sqtwo=1.4142135623730950488d0
real(8) s1,s2,t1,t2,t3,t4,t5,t6,t7,t8
! automatic arrays
integer lmi(-lmax:lmax)
real(8) ca(lmax),sa(lmax),cg(lmax),sg(lmax),dy(ld,ld)
! generate the complex spherical harmonic rotation matrix about the y-axis
call ylmroty(ang(2),lmax,ld,dy)
do m1=1,lmax
  ca(m1)=cos(m1*ang(1)); sa(m1)=sin(m1*ang(1))
  cg(m1)=cos(m1*ang(3)); sg(m1)=sin(m1*ang(3))
end do
lm1=0
do l=0,lmax
  do m1=-l,l
    lm1=lm1+1
    lmi(m1)=lm1
  end do
  lm0=lmi(0)
  d(lm0,lm0)=dy(lm0,lm0)
  do m1=1,l
    if (mod(m1,2) == 0) then
      s1=1.d0
    else
      s1=-1.d0
    end if
    t1=sqtwo*dy(lm0,lmi(m1))
    t2=sqtwo*dy(lmi(m1),lm0)
    d(lmi(m1),lm0)=s1*t1*ca(m1)
    d(lm0,lmi(m1))=s1*t2*cg(m1)
    d(lmi(-m1),lm0)=-t1*sa(m1)
    d(lm0,lmi(-m1))=t2*sg(m1)
    do m2=1,l
      if (mod(m2,2) == 0) then
        s2=1.d0
      else
        s2=-1.d0
      end if
      t1=ca(m1)*cg(m2)
      t2=sa(m1)*sg(m2)
      t3=sa(m1)*cg(m2)
      t4=ca(m1)*sg(m2)
      t5=dy(lmi(-m1),lmi(-m2))
      t6=s1*dy(lmi(m1),lmi(-m2))
      t7=t5+t6
      t8=t5-t6
      d(lmi(m1),lmi(m2))=s1*s2*(t1*t7-t2*t8)
      d(lmi(m1),lmi(-m2))=s1*(t3*t8+t4*t7)
      d(lmi(-m1),lmi(m2))=-s2*(t3*t7+t4*t8)
      d(lmi(-m1),lmi(-m2))=t1*t8-t2*t7
    end do
  end do
end do
! apply inversion if required
if (p == -1) then
  do l=1,lmax,2
    lm1=l**2+1
    lm2=lm1+2*l
    d(lm1:lm2,lm1:lm2)=-d(lm1:lm2,lm1:lm2)
  end do
end if
end subroutine
!EOC
 
end subroutine