TABLE OF CONTENTS
ABINIT/getspinrot [ Functions ]
NAME
getspinrot
FUNCTION
From the symmetry matrix symrel_conv expressed in the coordinate system rprimd, compute the components of the spinor rotation matrix
COPYRIGHT
Copyright (C) 2000-2018 ABINIT group (XG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
rprimd(3,3)=dimensional primitive translations for real space (bohr) symrel_conv(3,3)=symmetry operation in real space in terms of primitive translations rprimd
OUTPUT
spinrot(4)=components of the spinor rotation matrix : spinrot(1)=$\cos \phi / 2$ spinrot(2)=$\sin \phi / 2 \times u_x$ spinrot(3)=$\sin \phi / 2 \times u_y$ spinrot(4)=$\sin \phi / 2 \times u_z$ where $\phi$ is the angle of rotation, and $(u_x,u_y,u_z)$ is the normalized direction of the rotation axis
NOTES
Only the proper part of the symmetry operation is taken into account : pure rotations, while the inversion part is taken away, if present. The whole collection of symmetry matrices is call symrel(3,3,nsym) symrel1 contains just one of those matrices symrel1(3,3)
PARENTS
cg_rotate,m_crystal,wfconv
CHILDREN
mati3det,matr3inv
SOURCE
46 #if defined HAVE_CONFIG_H 47 #include "config.h" 48 #endif 49 50 #include "abi_common.h" 51 52 53 subroutine getspinrot(rprimd,spinrot,symrel_conv) 54 55 use defs_basis 56 use m_profiling_abi 57 use m_errors 58 59 !This section has been created automatically by the script Abilint (TD). 60 !Do not modify the following lines by hand. 61 #undef ABI_FUNC 62 #define ABI_FUNC 'getspinrot' 63 use interfaces_32_util 64 !End of the abilint section 65 66 implicit none 67 68 !Arguments ------------------------------------ 69 !arrays 70 integer,intent(in) :: symrel_conv(3,3) 71 real(dp),intent(in) :: rprimd(3,3) 72 real(dp),intent(out) :: spinrot(4) 73 74 !Local variables------------------------------- 75 !scalars 76 integer :: det,ii 77 real(dp) :: cos_phi,norminv,phi,scprod,sin_phi 78 !character(len=500) :: message 79 !arrays 80 integer :: identity(3,3),symrel1(3,3) 81 real(dp) :: axis(3),coord(3,3),coordinvt(3,3),matr1(3,3),matr2(3,3) 82 real(dp) :: rprimd_invt(3,3),vecta(3),vectb(3),vectc(3) 83 84 !************************************************************************** 85 86 symrel1(:,:)=symrel_conv(:,:) 87 88 !Compute determinant of the matrix 89 call mati3det(symrel1,det) 90 91 !Produce a rotation from an improper symmetry 92 if(det==-1)symrel1(:,:)=-symrel1(:,:) 93 94 !Test the possibility of the unit matrix 95 identity(:,:)=0 96 identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1 97 98 if( sum((symrel1(:,:)-identity(:,:))**2)/=0)then 99 100 ! Transform symmetry matrix in the system defined by rprimd 101 call matr3inv(rprimd,rprimd_invt) 102 do ii=1,3 103 coord(:,ii)=rprimd_invt(ii,:) 104 end do 105 call matr3inv(coord,coordinvt) 106 do ii=1,3 107 matr1(:,ii)=symrel1(:,1)*coord(1,ii)+& 108 & symrel1(:,2)*coord(2,ii)+& 109 & symrel1(:,3)*coord(3,ii) 110 end do 111 do ii=1,3 112 matr2(:,ii)=coordinvt(1,:)*matr1(1,ii)+& 113 & coordinvt(2,:)*matr1(2,ii)+& 114 & coordinvt(3,:)*matr1(3,ii) 115 end do 116 117 ! Find the eigenvector with unit eigenvalue of the 118 ! rotation matrix in cartesian coordinate, matr2 119 120 matr1(:,:)=matr2(:,:) 121 matr1(1,1)=matr1(1,1)-one 122 matr1(2,2)=matr1(2,2)-one 123 matr1(3,3)=matr1(3,3)-one 124 125 ! Compute the axis of rotation and the cos and sin of rotation angle 126 if(matr1(1,1)**2 + matr1(2,1)**2 + matr1(3,1)**2 < tol8 )then 127 ! The first direction is the axis 128 axis(1)=one ; axis(2)=zero ; axis(3)=zero 129 cos_phi=matr2(2,2) 130 sin_phi=matr2(3,2) 131 else if(matr1(1,2)**2 + matr1(2,2)**2 + matr1(3,2)**2 < tol8 )then 132 ! The second direction is the axis 133 axis(1)=zero ; axis(2)=one ; axis(3)=zero 134 cos_phi=matr2(3,3) 135 sin_phi=matr2(1,3) 136 else 137 ! In this case, try use the first and second vector to build the 138 ! rotation axis : compute their cross product 139 axis(1)=matr1(2,1)*matr1(3,2)-matr1(2,2)*matr1(3,1) 140 axis(2)=matr1(3,1)*matr1(1,2)-matr1(3,2)*matr1(1,1) 141 axis(3)=matr1(1,1)*matr1(2,2)-matr1(1,2)*matr1(2,1) 142 ! Then, try to normalize it 143 scprod=sum(axis(:)**2) 144 if(scprod<tol8)then 145 ! The first and second vectors were linearly dependent 146 ! Thus, use the first and third vectors 147 axis(1)=matr1(2,1)*matr1(3,3)-matr1(2,3)*matr1(3,1) 148 axis(2)=matr1(3,1)*matr1(1,3)-matr1(3,3)*matr1(1,1) 149 axis(3)=matr1(1,1)*matr1(2,3)-matr1(1,3)*matr1(2,1) 150 ! Normalize the vector 151 scprod=sum(axis(:)**2) 152 if(scprod<tol8)then 153 MSG_BUG('Cannot find the rotation axis.') 154 end if 155 end if 156 norminv=one/sqrt(scprod) 157 axis(:)=axis(:)*norminv 158 159 ! Project the axis vector out of the first unit vector, 160 ! and renormalize the projected vector 161 ! (the first vector cannot be the axis, as tested before) 162 vecta(1)=one-axis(1)**2 163 vecta(2)=-axis(1)*axis(2) 164 vecta(3)=-axis(1)*axis(3) 165 scprod=sum(vecta(:)**2) 166 norminv=one/sqrt(scprod) 167 vecta(:)=vecta(:)*norminv 168 ! Rotate the vector A, to get vector B 169 vectb(:)=matr2(:,1)*vecta(1)+matr2(:,2)*vecta(2)+matr2(:,3)*vecta(3) 170 ! Get dot product of vectors A and B, giving cos of the rotation angle 171 cos_phi=sum(vecta(:)*vectb(:)) 172 ! Compute the cross product of the axis and vector A 173 vectc(1)=axis(2)*vecta(3)-axis(3)*vecta(2) 174 vectc(2)=axis(3)*vecta(1)-axis(1)*vecta(3) 175 vectc(3)=axis(1)*vecta(2)-axis(2)*vecta(1) 176 ! Get dot product of vectors B and C, giving sin of the rotation angle 177 sin_phi=sum(vectb(:)*vectc(:)) 178 end if 179 180 ! Get the rotation angle, then the parameters of the spinor rotation 181 ! Here, treat possible inaccurate values of the cosine of phi 182 if(cos_phi> one-tol8 )cos_phi= one-tol8 183 if(cos_phi<-(one-tol8))cos_phi=-(one-tol8) 184 phi=acos(cos_phi) 185 if(sin_phi<zero)phi=-phi 186 ! Rectify the angle, such that its absolute values corresponds to 187 ! 180, 120, 90, 60, or 0 degrees 188 phi=(nint(six*phi/pi))/six*pi 189 ! Compute components of the spinor matrix 190 spinrot(1)=cos(half*phi) 191 spinrot(2)=axis(1)*sin(half*phi) 192 spinrot(3)=axis(2)*sin(half*phi) 193 spinrot(4)=axis(3)*sin(half*phi) 194 195 else 196 197 ! Here, the case of the unit matrix 198 axis(:)=zero 199 phi=zero 200 spinrot(1)=one 201 spinrot(2)=zero 202 spinrot(3)=zero 203 spinrot(4)=zero 204 205 end if ! the case of the identity matrix 206 207 !DEBUG 208 !write(std_out,*)' getspinrot :' 209 !write(std_out,*)' symrel_conv =',symrel_conv(:,:) 210 !write(std_out,*)' symrel =',symrel1(:,:) 211 !write(std_out,*)' rprimd =',rprimd(:,:) 212 !write(std_out,*)' matr2 =',matr2(:,:) 213 !write(std_out,*)' matr1 =',matr1(:,:) 214 !write(std_out,*)' phi (degree)=',phi*180._dp/pi 215 !write(std_out,'(a,3d16.6)' )' axis=',axis(:) 216 !write(std_out,*)' vecta=',vecta(:) 217 !stop 218 !ENDDEBUG 219 220 end subroutine getspinrot