mirror of
https://github.com/triqs/dft_tools
synced 2024-12-21 20:03:41 +01:00
369 lines
17 KiB
Fortran
369 lines
17 KiB
Fortran
|
|
c ******************************************************************************
|
|
c
|
|
c TRIQS: a Toolbox for Research in Interacting Quantum Systems
|
|
c
|
|
c Copyright (C) 2011 by L. Pourovskii, V. Vildosola, C. Martins, M. Aichhorn
|
|
c
|
|
c TRIQS is free software: you can redistribute it and/or modify it under the
|
|
c terms of the GNU General Public License as published by the Free Software
|
|
c Foundation, either version 3 of the License, or (at your option) any later
|
|
c version.
|
|
c
|
|
c TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
|
|
c WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
|
c FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
|
|
c details.
|
|
c
|
|
c You should have received a copy of the GNU General Public License along with
|
|
c TRIQS. If not, see <http://www.gnu.org/licenses/>.
|
|
c
|
|
c *****************************************************************************/
|
|
|
|
SUBROUTINE set_rotloc
|
|
C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
C %% %%
|
|
C %% This subroutine sets up the Global->local coordinates %%
|
|
C %% rotational matrices for each atom of the system. %%
|
|
C %% These matrices will be used to create the projectors. %%
|
|
C %% (They are the SR matrices defined in the tutorial file.) %%
|
|
C %% %%
|
|
C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
C Definiton of the variables :
|
|
C ----------------------------
|
|
USE common_data
|
|
USE reps
|
|
USE symm
|
|
USE prnt
|
|
IMPLICIT NONE
|
|
COMPLEX(KIND=8), DIMENSION(:,:), ALLOCATABLE :: tmp_rot, spinrot
|
|
REAL(KIND=8) :: alpha, beta, gama, factor
|
|
INTEGER :: iatom, jatom, imu, isrt
|
|
INTEGER :: is, is1, isym, l, lm
|
|
INTEGER :: ind1, ind2, inof1, inof2
|
|
COMPLEX(KIND=8) :: ephase
|
|
C
|
|
C ====================================================
|
|
C Multiplication by an S matrix for equivalent sites :
|
|
C ====================================================
|
|
C Up to now, rotloc is the rotloc matrix (from Global to local coordinates rotation : (rotloc)_ij = <x_global_i | x_local_j >)
|
|
C The matrix S to go from the representative atom of the sort to another one must be introduced. That's what is done here-after.
|
|
DO isrt=1,nsort
|
|
iatom=SUM(nmult(0:isrt-1))+1
|
|
DO imu=1,nmult(isrt)
|
|
jatom=iatom+imu-1
|
|
DO isym=1,nsym
|
|
C If the symmetry operation isym transforms the representative atom iatom in the jatom,
|
|
C the matrix rotloc is multiplied by the corresponding srot matrix, for each orbital number l.
|
|
C if R[isym](iatom) = jatom, rotloc is multiplied by R[isym] and Rloc is finally R[isym] X rotloc = <x_global|x_sym><x_sym|x_local>
|
|
IF(srot(isym)%perm(iatom)==jatom) THEN
|
|
WRITE(17,*) ' For jatom = ',jatom, ', isym =', isym
|
|
rotloc(jatom)%srotnum=isym
|
|
C Calculation of krotm and iprop.
|
|
rotloc(jatom)%krotm(1:3,1:3)=
|
|
= MATMUL(srot(isym)%krotm(1:3,1:3),
|
|
& rotloc(jatom)%krotm(1:3,1:3))
|
|
rotloc(jatom)%iprop=rotloc(jatom)%iprop*
|
|
* srot(isym)%iprop
|
|
C Evaluation of the Euler angles of the final operation Rloc
|
|
CALL euler(TRANSPOSE(rotloc(jatom)%krotm(1:3,1:3)),
|
|
& alpha,beta,gama)
|
|
C According to Wien convention, euler takes in argument the transpose
|
|
C of the matrix rotloc(jatom)%krotm to give a,b anc c of rotloc(jatom).
|
|
rotloc(jatom)%a=alpha
|
|
rotloc(jatom)%b=beta
|
|
rotloc(jatom)%g=gama
|
|
C
|
|
C =============================================================================================================
|
|
C Calculation of the rotational matrices and evaluation of the fields timeinv and phase for the Rloc matrices :
|
|
C =============================================================================================================
|
|
IF(ifSP.AND.ifSO) THEN
|
|
C No time reversal operation is applied to rotloc (alone). If a time reversal operation must be applied,
|
|
C it comes from the symmetry operation R[isym]. That is why the field timeinv is the same as the one from srot.
|
|
rotloc(jatom)%timeinv=srot(isym)%timeinv
|
|
rotloc(jatom)%phase=0.d0
|
|
DO l=1,lmax
|
|
ALLOCATE(tmp_rot(1:2*(2*l+1),1:2*(2*l+1)))
|
|
tmp_rot=0.d0
|
|
C Whatever the value of beta (0 or Pi), the spinor rotation matrix of isym is block-diagonal.
|
|
C because the time-reversal operation have been applied if necessary.
|
|
factor=srot(isym)%phase/2.d0
|
|
ephase=EXP(CMPLX(0.d0,factor))
|
|
C We remind that the field phase is (g-a) if beta=Pi. As a result, ephase = exp(+i(g-a)/2) = -exp(+i(alpha-gamma)/2)
|
|
C We remind that the field phase is (a+g) if beta=0. As a result, ephase = exp(+i(a+g)/2)=-exp(-i(alpha+gamma)/2)
|
|
C in good agreement with Wien conventions for the definition of this phase factor.
|
|
C Up/up block :
|
|
tmp_rot(1:2*l+1,1:2*l+1)=ephase*
|
|
& srot(isym)%rotl(-l:l,-l:l,l)
|
|
C Dn/dn block :
|
|
ephase=CONJG(ephase)
|
|
C now, ephase = exp(+i(a-g)/2) = -exp(-i(alpha-gamma)/2) if beta=Pi
|
|
C now, ephase = exp(-i(a+g)/2) = -exp(+i(alpha+gamma)/2) if beta=0
|
|
tmp_rot(2*l+2:2*(2*l+1),2*l+2:2*(2*l+1))=
|
|
& ephase*srot(isym)%rotl(-l:l,-l:l,l)
|
|
IF (rotloc(jatom)%timeinv) THEN
|
|
C In this case, the time reversal operator was applied to srot.
|
|
rotloc(jatom)%rotl(1:2*(2*l+1),1:2*(2*l+1),l)=
|
|
& MATMUL(tmp_rot(1:2*(2*l+1),1:2*(2*l+1)),CONJG(
|
|
& rotloc(jatom)%rotl(1:2*(2*l+1),1:2*(2*l+1),l)))
|
|
C rotloc(jatom)%rotl now contains D(Rloc) = D(R[isym])*transpose[D(rotloc)].
|
|
ELSE
|
|
C In this case, no time reversal operator was applied to srot.
|
|
rotloc(jatom)%rotl(1:2*(2*l+1),1:2*(2*l+1),l)=
|
|
& MATMUL(tmp_rot(1:2*(2*l+1),1:2*(2*l+1)),
|
|
& rotloc(jatom)%rotl(1:2*(2*l+1),1:2*(2*l+1),l))
|
|
C rotloc(jatom)%rotl now contains D(Rloc) = D(R[isym])*D(rotloc).
|
|
ENDIF
|
|
DEALLOCATE(tmp_rot)
|
|
ENDDO
|
|
ELSE
|
|
C Calculation of the rotational matrices associated to Rloc
|
|
ALLOCATE(tmp_rot(1:2*lmax+1,1:2*lmax+1))
|
|
DO l=1,lmax
|
|
C Use of the subroutine dmat to compute the rotational matrix
|
|
C associated to the Rloc operation in a (2*l+1) space :
|
|
tmp_rot=0.d0
|
|
CALL dmat(l,rotloc(jatom)%a,rotloc(jatom)%b,
|
|
& rotloc(jatom)%g,
|
|
& REAL(rotloc(jatom)%iprop,KIND=8),tmp_rot,2*lmax+1)
|
|
rotloc(jatom)%rotl(-l:l,-l:l,l)=
|
|
= tmp_rot(1:2*l+1,1:2*l+1)
|
|
C rotloc(jatom)%rotl = table of the rotational matrices of the symmetry operation
|
|
C for the different l orbital (from 1 to lmax), in the usual complex basis : dmat = D(R[isym])_l
|
|
C rotloc(jatom)%rotl = D(Rloc[jatom])_{lm}
|
|
ENDDO
|
|
DEALLOCATE(tmp_rot)
|
|
ENDIF ! End of the "ifSO-ifSP" if-then-else
|
|
C
|
|
EXIT
|
|
C Only one symmetry operation is necessary to be applied to R to get the complete rotloc matrix.
|
|
C This EXIT enables to leave the loop as soon as a symmetry operation which transforms the representative atom in jatom is found.
|
|
ENDIF ! End of the "perm" if-then-else
|
|
ENDDO ! End of the isym loop
|
|
C
|
|
C
|
|
C ===========================================================
|
|
C Computation of the rotational matrices in each sort basis :
|
|
C ===========================================================
|
|
ALLOCATE(rotloc(jatom)%rotrep(lmax))
|
|
C
|
|
C Initialization of the rotloc(jatom)%rotrep field = D(Rloc)_{new_i}
|
|
C This field is a table of size lmax which contains the rotloc matrices
|
|
C in the representation basis associated to each included orbital of the jatom.
|
|
DO l=1,lmax
|
|
ALLOCATE(rotloc(jatom)%rotrep(l)%mat(1,1))
|
|
rotloc(jatom)%rotrep(l)%mat(1,1)=0.d0
|
|
ENDDO
|
|
C
|
|
C Computation of the elements 'mat' in rotloc(jatom)%rotrep(l)
|
|
DO l=1,lmax
|
|
C The considered orbital is not included, hence no computation
|
|
IF (lsort(l,isrt)==0) cycle
|
|
C The considered orbital is included
|
|
IF (ifSP.AND.ifSO) THEN
|
|
C In this case, the basis representation needs a complete spinor rotation approach (matrices of size 2*(2*l+1) )
|
|
C --------------------------------------------------------------------------------------------------------------
|
|
DEALLOCATE(rotloc(jatom)%rotrep(l)%mat)
|
|
ALLOCATE(rotloc(jatom)%rotrep(l)%mat
|
|
& (1:2*(2*l+1),1:2*(2*l+1)))
|
|
ALLOCATE(tmp_rot(1:2*(2*l+1),1:2*(2*l+1)))
|
|
C Computation of rotloc(jatom)%rotrep(l)%mat
|
|
IF (reptrans(l,isrt)%ifmixing) THEN
|
|
C In this case, the basis representation requires a complete spinor rotation approach too.
|
|
IF(rotloc(jatom)%timeinv) THEN
|
|
tmp_rot(1:2*(2*l+1),1:2*(2*l+1))=MATMUL(
|
|
& reptrans(l,isrt)%transmat(1:2*(2*l+1),1:2*(2*l+1)),
|
|
& rotloc(jatom)%rotl(1:2*(2*l+1),1:2*(2*l+1),l))
|
|
rotloc(jatom)%rotrep(l)%mat(1:2*(2*l+1),1:2*(2*l+1))=
|
|
= MATMUL(tmp_rot(1:2*(2*l+1),1:2*(2*l+1)),
|
|
& TRANSPOSE(reptrans(l,isrt)%transmat
|
|
& (1:2*(2*l+1),1:2*(2*l+1))))
|
|
C Since the operation is antilinear, the field rotloc(jatom)%rotrep(l)%mat = (reptrans)*spinrot(l)*conjugate(inverse(reptrans))
|
|
C rotloc(jatom)%rotrep(l)%mat = D(Rloc)_{new_i} = <new_i|lm> D(Rloc)_{lm} [<lm|new_i>]^*
|
|
C which is exactly the expression of the spinor rotation matrix in the new basis.
|
|
ELSE
|
|
tmp_rot(1:2*(2*l+1),1:2*(2*l+1))=MATMUL(
|
|
& reptrans(l,isrt)%transmat(1:2*(2*l+1),1:2*(2*l+1)),
|
|
& rotloc(jatom)%rotl(1:2*(2*l+1),1:2*(2*l+1),l))
|
|
rotloc(jatom)%rotrep(l)%mat(1:2*(2*l+1),1:2*(2*l+1))=
|
|
= MATMUL(tmp_rot(1:2*(2*l+1),1:2*(2*l+1)),
|
|
& TRANSPOSE(CONJG(reptrans(l,isrt)%transmat
|
|
& (1:2*(2*l+1),1:2*(2*l+1)))))
|
|
C Since the operation is linear, the field rotloc(jatom)%rotrep(l)%mat = (reptrans)*spinrot(l)*inverse(reptrans)
|
|
C rotloc(jatom)%rotrep(l)%mat = D(Rloc)_{new_i} = <new_i|lm> D(Rloc)_{lm} <lm|new_i>
|
|
C which is exactly the expression of the spinor rotation matrix in the new basis.
|
|
ENDIF
|
|
ELSE
|
|
C In this case, the basis representation is reduced to the up/up block and must be extended.
|
|
ALLOCATE(spinrot(1:2*(2*l+1),1:2*(2*l+1)))
|
|
spinrot(1:2*(2*l+1),1:2*(2*l+1))=0.d0
|
|
spinrot(1:2*l+1,1:2*l+1)=
|
|
& reptrans(l,isrt)%transmat(-l:l,-l:l)
|
|
spinrot(2*l+2:2*(2*l+1),2*l+2:2*(2*l+1))=
|
|
& reptrans(l,isrt)%transmat(-l:l,-l:l)
|
|
IF(rotloc(jatom)%timeinv) THEN
|
|
tmp_rot(1:2*(2*l+1),1:2*(2*l+1))=MATMUL(
|
|
& spinrot(1:2*(2*l+1),1:2*(2*l+1)),
|
|
& rotloc(jatom)%rotl(1:2*(2*l+1),1:2*(2*l+1),l))
|
|
rotloc(jatom)%rotrep(l)%mat(1:2*(2*l+1),1:2*(2*l+1))=
|
|
= MATMUL(tmp_rot(1:2*(2*l+1),1:2*(2*l+1)),
|
|
& TRANSPOSE(spinrot(1:2*(2*l+1),1:2*(2*l+1))))
|
|
C Since the operation is antilinear, the field rotloc(jatom)%rotrep(l)%mat = (reptrans)*spinrot(l)*conjugate(inverse(reptrans))
|
|
C rotloc(jatom)%rotrep(l)%mat = D(Rloc)_{new_i} = <new_i|lm> D(Rloc)_{lm} [<lm|new_i>]^*
|
|
C which is exactly the expression of the spinor rotation matrix in the new basis.
|
|
ELSE
|
|
tmp_rot(1:2*(2*l+1),1:2*(2*l+1))=MATMUL(
|
|
& spinrot(1:2*(2*l+1),1:2*(2*l+1)),
|
|
& rotloc(jatom)%rotl(1:2*(2*l+1),1:2*(2*l+1),l))
|
|
rotloc(jatom)%rotrep(l)%mat(1:2*(2*l+1),1:2*(2*l+1))=
|
|
= MATMUL(tmp_rot(1:2*(2*l+1),1:2*(2*l+1)),
|
|
& TRANSPOSE(CONJG(spinrot(1:2*(2*l+1),1:2*(2*l+1)))))
|
|
C Since the operation is linear, the field rotloc(jatom)%rotrep(l)%mat = (reptrans)*spinrot(l)*inverse(reptrans)
|
|
C rotloc(jatom)%rotrep(l)%mat = D(Rloc)_{new_i} = <new_i|lm> D(Rloc)_{lm} <lm|new_i>
|
|
C which is exactly the expression of the spinor rotation matrix in the new basis.
|
|
ENDIF
|
|
DEALLOCATE(spinrot)
|
|
ENDIF ! End of the if mixing if-then-else
|
|
DEALLOCATE(tmp_rot)
|
|
C
|
|
ELSE
|
|
C If the basis representation can be reduce to the up/up block (matrices of size (2*l+1) only)
|
|
C --------------------------------------------------------------------------------------------
|
|
DEALLOCATE(rotloc(jatom)%rotrep(l)%mat)
|
|
ALLOCATE(rotloc(jatom)%rotrep(l)%mat(-l:l,-l:l))
|
|
ALLOCATE(tmp_rot(-l:l,-l:l))
|
|
C Computation of rotloc(jatom)%rotrep(l)%mat
|
|
tmp_rot(-l:l,-l:l)=MATMUL(
|
|
& reptrans(l,isrt)%transmat(-l:l,-l:l),
|
|
& rotloc(jatom)%rotl(-l:l,-l:l,l))
|
|
rotloc(jatom)%rotrep(l)%mat(-l:l,-l:l)=
|
|
= MATMUL(tmp_rot(-l:l,-l:l),
|
|
& TRANSPOSE(CONJG(reptrans(l,isrt)%transmat(-l:l,-l:l))))
|
|
C the field rotloc(jatom)%rotrep(l)%mat = (reptrans)*rotl*inverse(reptrans)
|
|
C rotloc(jatom)%rotrep(l)%mat = D(Rloc)_{new_i} = <new_i|lm> D(Rloc)_{lm} <lm|new_i>
|
|
C which is exactly the expression of the rotation matrix for the up/up block in the new basis.
|
|
DEALLOCATE(tmp_rot)
|
|
ENDIF
|
|
ENDDO ! End of the l loop
|
|
ENDDO ! End of the jatom loop
|
|
ENDDO ! End of the isrt loop
|
|
C
|
|
RETURN
|
|
END
|
|
|
|
|
|
SUBROUTINE euler(Rot,a,b,c)
|
|
C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
C %% %%
|
|
C %% This subroutine calculates the Euler angles a, b and c of Rot. %%
|
|
C %% The result are stored in a,b,c. (same as in SRC_lapwdm/euler.f) %%
|
|
C %% %%
|
|
C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
C
|
|
IMPLICIT NONE
|
|
REAL(KIND=8) :: a,aa,b,bb,c,cc,zero,pi,y_norm,dot
|
|
REAL(KIND=8), DIMENSION(3,3) :: Rot, Rot_temp
|
|
REAL(KIND=8), DIMENSION(3) :: z,zz,y,yy,yyy,pom,x,xx
|
|
INTEGER :: i,j
|
|
C Definition of the constants
|
|
zero=0d0
|
|
pi=ACOS(-1d0)
|
|
C Definition of Rot_temp=Id
|
|
DO i=1,3
|
|
DO j=1,3
|
|
Rot_temp(i,j)=0
|
|
IF (i.EQ.j) Rot_temp(i,i)=1
|
|
ENDDO
|
|
ENDDO
|
|
C Initialization of y=e_y, z=e_z, yyy and zz
|
|
DO j=1,3
|
|
y(j)=Rot_temp(j,2)
|
|
yyy(j)=Rot(j,2)
|
|
z(j)=Rot_temp(j,3)
|
|
zz(j)=Rot(j,3)
|
|
ENDDO
|
|
C Calculation of yy
|
|
CALL vecprod(z,zz,yy)
|
|
y_norm=DSQRT(dot(yy,yy))
|
|
IF (y_norm.lt.1d-10) THEN
|
|
C If yy=0, this implies that b is zero or pi
|
|
IF (ABS(dot(y,yyy)).gt.1d0) THEN
|
|
aa=dot(y,yyy)/ABS(dot(y,yyy))
|
|
a=ACOS(aa)
|
|
ELSE
|
|
a=ACOS(dot(y,yyy))
|
|
ENDIF
|
|
C
|
|
IF (dot(z,zz).gt.zero) THEN
|
|
c=zero
|
|
b=zero
|
|
IF (yyy(1).gt.zero) a=2*pi-a
|
|
ELSE
|
|
c=a
|
|
a=zero
|
|
b=pi
|
|
IF (yyy(1).lt.zero) c=2*pi-c
|
|
ENDIF
|
|
ELSE
|
|
C If yy is not 0, then b belongs to ]0,pi[
|
|
DO j=1,3
|
|
yy(j)=yy(j)/y_norm
|
|
ENDDO
|
|
C
|
|
aa=dot(y,yy)
|
|
bb=dot(z,zz)
|
|
cc=dot(yy,yyy)
|
|
IF (ABS(aa).gt.1d0) aa=aa/ABS(aa)
|
|
IF (ABS(bb).gt.1d0) bb=bb/ABS(bb)
|
|
IF (ABS(cc).gt.1d0) cc=cc/ABS(cc)
|
|
b=ACOS(bb)
|
|
a=ACOS(aa)
|
|
c=ACOS(cc)
|
|
IF (yy(1).gt.zero) a=2*pi-a
|
|
CALL vecprod(yy,yyy,pom)
|
|
IF (dot(pom,zz).lt.zero) c=2*pi-c
|
|
ENDIF
|
|
C
|
|
END
|
|
|
|
|
|
SUBROUTINE vecprod(a,b,c)
|
|
C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
C %% %%
|
|
C %% This subroutine calculates the vector product of a and b. %%
|
|
C %% The result is stored in c. (same as in SRC_lapwdm/euler.f) %%
|
|
C %% %%
|
|
C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
C
|
|
IMPLICIT NONE
|
|
REAL(KIND=8), DIMENSION(3) :: a,b,c
|
|
C
|
|
c(1)=a(2)*b(3)-a(3)*b(2)
|
|
c(2)=a(3)*b(1)-a(1)*b(3)
|
|
c(3)=a(1)*b(2)-a(2)*b(1)
|
|
C
|
|
END
|
|
|
|
REAL(KIND=8) FUNCTION dot(a,b)
|
|
C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
C %% %%
|
|
C %% This function calculates the scalar product of a and b. %%
|
|
C %% The result is stored in dot. (same as in SRC_lapwdm/euler.f) %%
|
|
C %% %%
|
|
C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
C
|
|
IMPLICIT NONE
|
|
REAL(KIND=8) :: a,b
|
|
INTEGER :: i
|
|
dimension a(3),b(3)
|
|
dot=0
|
|
DO i=1,3
|
|
dot=dot+a(i)*b(i)
|
|
ENDDO
|
|
C
|
|
END
|
|
|
|
|
|
|