LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-11-18 12:03:50 +01:00
Code Issues Releases Wiki Activity
85c8c2b58c
dft_tools/fortran/dmftproj/outputqmc.f
Olivier Parcollet 0e585ad9b4 First import. triqs 1.0 alpha1
2013-07-23 20:55:29 +02:00

1406 lines
59 KiB
Fortran
Raw Blame History

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 outqmc(elecn,qbbot)
C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C %% %%
C %% This subroutine creates the output files : %%
C %% - case.ctqmcout, with all the informations necessary for a %%
C %% CTQMC computation. %%
C %% - case.symqmc, describing all the symmetries of the system %%
c %% (necessary for CTQMC computation too). %%
C %% - case.parproj which gives the partial charge projectors for %%
C %% orbitals and for all atoms. %%
C %% - case.sympar which contains the symmetry matrices for all %%
C %% the included orbitals (for partial charge analysis) %%
C %% %%
C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C Definition of the variables :
C -----------------------------
USE almblm_data
USE common_data
USE file_names
USE prnt
USE reps
USE symm
USE projections
IMPLICIT NONE
C
INTEGER :: iorb, icrorb, irep, isrt
INTEGER :: l, m, is, i1, i2, isym
INTEGER :: ik, il, ib, ir, n
INTEGER :: ind1, ind2, iatom
INTEGER :: timeinvflag
REAL(KIND=8) :: qbbot, elecn, factor
COMPLEX(KIND=8) :: ephase
COMPLEX(KIND=8),DIMENSION(:,:), ALLOCATABLE :: hk
COMPLEX(KIND=8), DIMENSION(:,:), ALLOCATABLE :: spinrot
COMPLEX(KIND=8), DIMENSION(:,:), ALLOCATABLE :: densprint
COMPLEX(KIND=8), DIMENSION(:,:), ALLOCATABLE :: time_op
C
C
C =======================================
C Writing the file case.ctqmcout :
C =======================================
C
WRITE(buf,'(a)')'Writing the file case.ctqmcout...'
CALL printout(0)
C
C Definition of 1 electron-Volt
WRITE(ouctqmc,'(a)') '13.605698'
C
C ---------------------------------------
C General informations about the system :
C ---------------------------------------
C
C Number of k-points in the I-BZ
WRITE(ouctqmc,'(i6)') nk
C Definition of the spin-polarized flag ifSP
IF (ifSP) THEN
WRITE(ouctqmc,'(i6)') 1
ELSE
WRITE(ouctqmc,'(i6)') 0
ENDIF
C Definition of the Spin-orbit flag ifSO
IF (ifSO) THEN
WRITE(ouctqmc,'(i6)') 1
ELSE
WRITE(ouctqmc,'(i6)') 0
ENDIF
C The only possible combinations are :
C - (0,0) which stands for a paramagnetic computation without SO.
C - (1,0) which stands for computation without SO using spin-polarized input files.
C - (1,1) which stands for computation with SO using spin-polarized input files.
C
C Writing the total charge below the lower limit of the energy window (variable "qbbot")
WRITE(ouctqmc,*) qbbot
C Writing the total number of electrons (valence band+semicore) (variable "elecn").
C It is also the charge upto the Fermi level.
WRITE(ouctqmc,*) elecn
C
C ------------------------------------------
C Description of all the included orbitals :
C ------------------------------------------
C
C Definition of the number of included orbitals "norb"
WRITE(ouctqmc,'(i6)') norb
C Description of each orbital "iorb"
DO iorb=1,norb
IF (ifSO) THEN
WRITE(ouctqmc,'(4(i6,x))') orb(iorb)%atom, orb(iorb)%sort,
& orb(iorb)%l, 2*(2*orb(iorb)%l+1)
ELSE
WRITE(ouctqmc,'(4(i6,x))') orb(iorb)%atom, orb(iorb)%sort,
& orb(iorb)%l, 2*orb(iorb)%l+1
ENDIF
ENDDO
C an orbital "iorb" is described by :
C - the associated atom
C - the corresponding atomic sort
C - the considered orbital number l
C - the size of the corresponding matrices : 2*l+1 without SO ; 2*(2*l+1) with SO
C
C ----------------------------------------
C Description of the correlated orbitals :
C ----------------------------------------
C
C Definition of the number of correlated orbitals "ncrorb"
WRITE(ouctqmc,'(i6)') ncrorb
C Description of each correlated orbital "icrorb"
DO icrorb=1,ncrorb
l=crorb(icrorb)%l
isrt=crorb(icrorb)%sort
C
C The case l=0 is a particular case of "non-mixing" basis.
C --------------------------------------------------------
IF (l==0) THEN
C For the s-orbitals, the only irep possible is the matrix itself.
IF (ifSP.AND.ifSO) THEN
C If SO is taken into account, spinor rotation matrix is considered.
C The spin is not a good quantum number, so the whole representation is used.
WRITE(ouctqmc,'(6(i6,x))') crorb(icrorb)%atom, isrt, 0,
& 2, crorb(icrorb)%ifSOat, 1
ELSE
C Without SO, only the rotation matrix in orbital space is necessary.
WRITE(ouctqmc,'(6(i6,x))') crorb(icrorb)%atom, isrt, 0,
& 1, crorb(icrorb)%ifSOat, 1
ENDIF
C
C If the basis representation needs a complete spinor rotation approach (basis with "mixing" ).
C ---------------------------------------------------------------------------------------------
ELSEIF(reptrans(l,isrt)%ifmixing) THEN
C In this case, the SO is necessary considered, spinor rotation matrices are used.
IF (crorb(icrorb)%ifsplit) THEN
C If only an irep is correlated
DO irep=1,reptrans(l,isrt)%nreps
IF (crorb(icrorb)%correp(irep)) THEN
WRITE(ouctqmc,'(6(i6,x))') crorb(icrorb)%atom,
& isrt, l, reptrans(l,isrt)%dreps(irep),
& crorb(icrorb)%ifSOat, irep
ENDIF
ENDDO
ELSE
C If no particular irep is correlated
WRITE(ouctqmc,'(6(i6,x))') crorb(icrorb)%atom, isrt, l,
& 2*(2*l+1), crorb(icrorb)%ifSOat, 1
ENDIF
C
C If the basis representation can be reduce to the up/up block (basis without "mixing").
C --------------------------------------------------------------------------------------
ELSE
IF (ifSP.AND.ifSO) THEN
C If SO is taken into account, spinor rotation matrices are considered.
C The spin is not a good quantum number, so the whole representation is used.
WRITE(ouctqmc,'(6(i6,x))') crorb(icrorb)%atom, isrt, l,
& 2*(2*l+1), crorb(icrorb)%ifSOat, 1
ELSE
C Without SO, only the rotation matrix in orbital space is necessary.
IF (crorb(icrorb)%ifsplit) THEN
C If only an irep is correlated
DO irep=1,reptrans(l,isrt)%nreps
IF (crorb(icrorb)%correp(irep)) THEN
WRITE(ouctqmc,'(6(i6,x))') crorb(icrorb)%atom,
& isrt, l, reptrans(l,isrt)%dreps(irep),
& crorb(icrorb)%ifSOat, irep
ENDIF
ENDDO
ELSE
C If no particular irep is correlated
WRITE(ouctqmc,'(6(i6,x))') crorb(icrorb)%atom, isrt, l,
& (2*l+1), crorb(icrorb)%ifSOat, 1
END IF ! End of the ifsplit if-then-else
END IF ! End of the ifSO if-then-else
END IF ! End of the ifmixing if-then-else
END DO ! End of the icrorb loop
C an orbital "iorb" is described by :
C - the associated atom
C - the corresponding atomic sort
C - the considered orbital number l
C - the size of the "correlated" submatrix (can be the whole matrix)
C - the flag ifSOat which states that SO is considered for this orbital
C - the number of the irep
C
C ------------------------------------------------------------------------------------------------
C Description of the global to local coordinates transformation Rloc for each correlated orbital :
C ------------------------------------------------------------------------------------------------
C Description of each transformation Rloc
DO icrorb=1,ncrorb
l=crorb(icrorb)%l
isrt=crorb(icrorb)%sort
iatom=crorb(icrorb)%atom
C
C The case l=0 is a particular case of "non-mixing" basis.
C --------------------------------------------------------
IF(l==0) THEN
C For the s-orbitals, the only irep possible is the matrix itself.
IF (ifSP.AND.ifSO) THEN
C If SO is taken into account, spinor rotation matrix must be considered.
ALLOCATE(spinrot(1:2,1:2))
spinrot(:,:)=0.d0
C The spinor-rotation matrix is directly calculated from the Euler angles a,b and c.
IF (rotloc(iatom)%timeinv) THEN
factor=(rotloc(iatom)%a+rotloc(iatom)%g)/2.d0
spinrot(2,1)=EXP(CMPLX(0.d0,factor))*
& DCOS(rotloc(iatom)%b/2.d0)
spinrot(1,2)=-CONJG(spinrot(2,1))
C Up/dn and Dn/up terms
factor=-(rotloc(iatom)%a-rotloc(iatom)%g)/2.d0
spinrot(2,2)=-EXP(CMPLX(0.d0,factor))*
& DSIN(rotloc(iatom)%b/2.d0)
spinrot(1,1)=CONJG(spinrot(2,2))
ELSE
factor=(rotloc(iatom)%a+rotloc(iatom)%g)/2.d0
spinrot(1,1)=EXP(CMPLX(0.d0,factor))*
& DCOS(rotloc(iatom)%b/2.d0)
spinrot(2,2)=CONJG(spinrot(1,1))
C Up/dn and Dn/up terms
factor=-(rotloc(iatom)%a-rotloc(iatom)%g)/2.d0
spinrot(1,2)=EXP(CMPLX(0.d0,factor))*
& DSIN(rotloc(iatom)%b/2.d0)
spinrot(2,1)=-CONJG(spinrot(1,2))
ENDIF
C Printing the transformation informations
DO m=1,2
WRITE(ouctqmc,*) REAL(spinrot(m,1:2))
ENDDO
DO m=1,2
WRITE(ouctqmc,*) AIMAG(spinrot(m,1:2))
ENDDO
DEALLOCATE(spinrot)
C Without SO, only the rotation matrix in orbital space is necessary.
ELSE
C In this case, the Rloc matrix is merely identity.
WRITE(ouctqmc,*) 1.d0
WRITE(ouctqmc,*) 0.d0
ENDIF
C Printing the time inversion flag if the calculation is spin-polarized (ifSP=1)
IF (ifSP) THEN
timeinvflag=0
IF (rotloc(iatom)%timeinv) timeinvflag=1
WRITE(ouctqmc,'(i6)') timeinvflag
ENDIF
C
C If the basis representation needs a complete spinor rotation approach (basis with "mixing" ).
C ---------------------------------------------------------------------------------------------
ELSEIF (reptrans(l,isrt)%ifmixing) THEN
C In this case, the calculation is necessary spin-polarized with SO, spinor rotation matrices are used.
IF(crorb(icrorb)%ifsplit) THEN
C If only an irep is correlated
ind1=1
DO irep=1,reptrans(l,isrt)%nreps
IF(crorb(icrorb)%correp(irep)) THEN
ind2=ind1+reptrans(l,isrt)%dreps(irep)-1
DO m=ind1,ind2
WRITE(ouctqmc,*) REAL(rotloc(iatom)%
& rotrep(l)%mat(m,ind1:ind2))
ENDDO
DO m=ind1,ind2
WRITE(ouctqmc,*)AIMAG(rotloc(iatom)%
& rotrep(l)%mat(m,ind1:ind2))
ENDDO
ENDIF
ind1=ind1+reptrans(l,isrt)%dreps(irep)
ENDDO
ELSE
C If no particular irep is correlated
DO m=1,2*(2*l+1)
WRITE(ouctqmc,*) REAL(rotloc(iatom)%
& rotrep(l)%mat(m,1:2*(2*l+1)))
ENDDO
DO m=1,2*(2*l+1)
WRITE(ouctqmc,*) AIMAG(rotloc(iatom)%
& rotrep(l)%mat(m,1:2*(2*l+1)))
ENDDO
ENDIF
C Printing if the transformation included a time-reversal operation.
timeinvflag =0
IF (rotloc(iatom)%timeinv) timeinvflag=1
WRITE(ouctqmc,'(i6)') timeinvflag
C
C If the basis representation can be reduce to the up/up block (basis without "mixing").
C --------------------------------------------------------------------------------------
ELSE
IF (ifSP.AND.ifSO) THEN
C If SO is taken into account, spinor rotation matrices are considered.
C The spin is not a good quantum number, so the whole representation is used.
DO m=1,2*(2*l+1)
WRITE(ouctqmc,*) REAL(rotloc(iatom)%
& rotrep(l)%mat(m,1:2*(2*l+1)))
ENDDO
DO m=1,2*(2*l+1)
WRITE(ouctqmc,*) AIMAG(rotloc(iatom)%
& rotrep(l)%mat(m,1:2*(2*l+1)))
ENDDO
ELSE
C The calculation is either spin-polarized without SO or paramagnetic.
C The spin is a good quantum number and irep are possible.
IF(crorb(icrorb)%ifsplit) THEN
C If only an irep is correlated
ind1=-l
DO irep=1,reptrans(l,isrt)%nreps
IF(crorb(icrorb)%correp(irep)) THEN
ind2=ind1+reptrans(l,isrt)%dreps(irep)-1
DO m=ind1,ind2
WRITE(ouctqmc,*) REAL(rotloc(iatom)%
& rotrep(l)%mat(m,ind1:ind2))
ENDDO
DO m=ind1,ind2
WRITE(ouctqmc,*) AIMAG(rotloc(iatom)%
& rotrep(l)%mat(m,ind1:ind2))
ENDDO
ENDIF
ind1=ind1+reptrans(l,isrt)%dreps(irep)
ENDDO
ELSE
C If no particular irep is correlated
DO m=-l,l
WRITE(ouctqmc,*) REAL(rotloc(iatom)%
& rotrep(l)%mat(m,-l:l))
ENDDO
DO m=-l,l
WRITE(ouctqmc,*) AIMAG(rotloc(iatom)%
& rotrep(l)%mat(m,-l:l))
ENDDO
END IF ! End of the ifsplit if-then-else
END IF ! End of the ifSO if-then-else
C Printing if the transformation included a time-reversal operation.
IF (ifSP) THEN
timeinvflag =0
IF (rotloc(iatom)%timeinv) timeinvflag=1
WRITE(ouctqmc,'(i6)') timeinvflag
END IF
END IF ! End of the ifmixing if-then-else
END DO ! End of the icrorb loop
C for each correlated orbital icrorb, the transformation Rloc is described by :
C - the real part of the submatrix associated to "icrorb"
C - the imaginary part of the submatrix associated to "icrorb"
C - a flag which states if a time reversal operation is included in the transformation (if SP only )
C
C -------------------------------------------------------------------------------------------------------------
C Description of the transformation from complex harmonics to the basis associated to each correlated orbital :
C -------------------------------------------------------------------------------------------------------------
C
DO icrorb=1,ncrorb
l=crorb(icrorb)%l
isrt=crorb(icrorb)%sort
C The transformation is printed only for the first (representative) atom of each sort.
IF (crorb(icrorb)%first) THEN
C
C The case l=0 is a particular case of "non-mixing" basis.
C --------------------------------------------------------
IF (l==0) THEN
C For the s-orbitals, the only basis considered is the complex basis.
IF (ifSP.AND.ifSO) THEN
C If SO is taken into account, orbital+spin space is considered.
ALLOCATE(spinrot(1:2,1:2))
spinrot(:,:)=0.d0
spinrot(1,1)=1.d0
spinrot(2,2)=1.d0
C Printing the number of irep in the considered basis and the size of each irep
C in the considered basis (in this case, it's 1 and 2.)
WRITE(ouctqmc,*) 1, 2
C Printing the transformation matrix
DO m=1,2
WRITE(ouctqmc,*) REAL(spinrot(m,1:2))
ENDDO
DO m=1,2
WRITE(ouctqmc,*) AIMAG(spinrot(m,1:2))
ENDDO
C
DEALLOCATE(spinrot)
ELSE
C Without SO, only the matrix in orbital space is necessary.
C Printing the number of irep in the considered basis and the size of each irep
C in the considered basis (in this case, it's 1 and 1.)
WRITE(ouctqmc,'(2(i6,x))') 1, 1
C Printing the transformation matrix
WRITE(ouctqmc,*) 1.d0
WRITE(ouctqmc,*) 0.d0
ENDIF
C
C If the basis representation needs a complete spinor rotation approach (basis with "mixing" ).
C ---------------------------------------------------------------------------------------------
ELSEIF (reptrans(l,isrt)%ifmixing) THEN
C In this case, the SO is necessary considered, spinor rotation matrices are used.
IF (crorb(icrorb)%ifsplit) THEN
WRITE(ouctqmc,'(15(i6,x))') reptrans(l,isrt)%nreps,
& reptrans(l,isrt)%dreps(1:reptrans(l,isrt)%nreps)
ELSE
WRITE(ouctqmc,'(2(i6,x))') 1, 2*(2*l+1)
ENDIF
DO m=1,2*(2*l+1)
WRITE(ouctqmc,*) REAL(reptrans(l,isrt)%
& transmat(m,1:2*(2*l+1)) )
ENDDO
DO m=1,2*(2*l+1)
WRITE(ouctqmc,*) AIMAG(reptrans(l,isrt)%
& transmat(m,1:2*(2*l+1)) )
ENDDO
C
C If the basis representation can be reduce to the up/up block (basis without "mixing").
C --------------------------------------------------------------------------------------
ELSE
IF (ifSP.AND.ifSO) THEN
C If SO is taken into account, orbital+spin space is considered.
ALLOCATE(spinrot(1:2*(2*l+1),1:2*(2*l+1)))
spinrot(:,:)=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)
C Printing the number of irep in the considered basis and the size of each irep
C in the considered basis (in this case, it's 1 and 2*(2*l+1).)
WRITE(ouctqmc,*) 1, 2*(2*l+1)
C Printing the transformation matrix
DO m=1,2*(2*l+1)
WRITE(ouctqmc,*) REAL(spinrot(m,1:2*(2*l+1)) )
ENDDO
DO m=1,2*(2*l+1)
WRITE(ouctqmc,*) AIMAG(spinrot(m,1:2*(2*l+1)) )
ENDDO
C
DEALLOCATE(spinrot)
C Without SO, only the rotation matrix in orbital space is necessary.
ELSE
IF (crorb(icrorb)%ifsplit) THEN
WRITE(ouctqmc,'(8(i6,x))') reptrans(l,isrt)%nreps,
& reptrans(l,isrt)%dreps(1:reptrans(l,isrt)%nreps)
ELSE
WRITE(ouctqmc,'(2(i6,x))') 1, 2*l+1
ENDIF
DO m=-l,l
WRITE(ouctqmc,*) REAL(reptrans(l,isrt)%
& transmat(m,-l:l) )
ENDDO
DO m=-l,l
WRITE(ouctqmc,*) AIMAG(reptrans(l,isrt)%
& transmat(m,-l:l) )
END DO
END IF ! End of the ifSO if-then-else
END IF ! End of the ifmixing if-then-else
END IF ! End of the iffirst if-then-else
END DO ! End of the icrorb loop
C for each correlated orbital icrorb, the basis transformation is described by :
C - the number of irep in the new basis
C - the dimension of each irep in the new basis
C - the real part of the basis transformation
C - the imaginary part of the basis transformation
C
C -------------------------------------------------------------------------
C Description of the number of bands in the energy window at each k_point :
C -------------------------------------------------------------------------
C
DO is=1,ns
C If SO is considered, the number of up and dn bands are the same.
IF ((ifSP.AND.ifSO).and.(is.eq.2)) cycle
C Printing the number of included bands in the window for each k_point
DO ik=1,nk
WRITE(ouctqmc,'(i6)')
& ABS(kp(ik,is)%nb_top-kp(ik,is)%nb_bot+1)
ENDDO
ENDDO
C
C -----------------------------------------------------------
C Description of the projectors for the correlated orbitals :
C -----------------------------------------------------------
DO ik=1,nk
DO icrorb=1,ncrorb
l=crorb(icrorb)%l
isrt=crorb(icrorb)%sort
C
C The case l=0 is a particular case of "non-mixing" basis.
C --------------------------------------------------------
IF (l==0) THEN
C For the s-orbitals, the only irep possible is the matrix itself.
C if the calculation is spin-polarized, up and dn projectors are written the one above the other (orbital+spin-space of size 2)
DO is=1,ns
WRITE(ouctqmc,*)
& REAL(pr_crorb(icrorb,ik,is)%mat_rep(1,
& kp(ik,is)%nb_bot:kp(ik,is)%nb_top))
ENDDO
DO is=1,ns
WRITE(ouctqmc,*)
& AIMAG(pr_crorb(icrorb,ik,is)%mat_rep(1,
& kp(ik,is)%nb_bot:kp(ik,is)%nb_top))
ENDDO
C
C If the basis representation needs a complete spinor rotation approach (basis with "mixing" ).
C ---------------------------------------------------------------------------------------------
ELSEIF (reptrans(l,isrt)%ifmixing) THEN
C In this case, the calculation is necessary spin-polarized with SO, spinor rotation matrices are used.
IF(crorb(icrorb)%ifsplit) THEN
C If only 1 irep is correlated
ind1=1
DO irep=1,reptrans(l,isrt)%nreps
IF(crorb(icrorb)%correp(irep)) THEN
ind2=ind1+reptrans(l,isrt)%dreps(irep)-1
DO m=ind1,ind2
WRITE(ouctqmc,*)
& REAL(pr_crorb(icrorb,ik,1)%mat_rep(m,
& kp(ik,1)%nb_bot:kp(ik,1)%nb_top))
ENDDO
DO m=ind1,ind2
WRITE(ouctqmc,*)
& AIMAG(pr_crorb(icrorb,ik,1)%mat_rep(m,
& kp(ik,1)%nb_bot:kp(ik,1)%nb_top))
ENDDO
ENDIF
ind1=ind1+reptrans(l,isrt)%dreps(irep)
ENDDO
ELSE
C If no particular irep is correlated
DO m=1,2*(2*l+1)
WRITE(ouctqmc,*)
& REAL(pr_crorb(icrorb,ik,1)%mat_rep(m,
& kp(ik,1)%nb_bot:kp(ik,1)%nb_top))
ENDDO
DO m=1,2*(2*l+1)
WRITE(ouctqmc,*)
& AIMAG(pr_crorb(icrorb,ik,1)%mat_rep(m,
& kp(ik,1)%nb_bot:kp(ik,1)%nb_top))
ENDDO
ENDIF
C
C If the basis representation can be reduce to the up/up block (basis without "mixing").
C --------------------------------------------------------------------------------------
ELSE
IF ((.not.(ifSP.AND.ifSO)).AND.crorb(icrorb)%ifsplit) THEN
C If only 1 irep is correlated (case without SO)
ind1=-l
DO irep=1,reptrans(l,isrt)%nreps
IF(crorb(icrorb)%correp(irep)) THEN
ind2=ind1+reptrans(l,isrt)%dreps(irep)-1
DO is=1,ns
DO m=ind1,ind2
WRITE(ouctqmc,*)
& REAL(pr_crorb(icrorb,ik,is)%mat_rep(m,
& kp(ik,is)%nb_bot:kp(ik,is)%nb_top))
ENDDO
ENDDO
DO is=1,ns
DO m=ind1,ind2
WRITE(ouctqmc,*)
& AIMAG(pr_crorb(icrorb,ik,is)%mat_rep(m,
& kp(ik,is)%nb_bot:kp(ik,is)%nb_top))
ENDDO
ENDDO
ENDIF
ind1=ind1+reptrans(l,isrt)%dreps(irep)
ENDDO
ELSE
C If no particular irep is correlated (case with and without SO)
DO is=1,ns
DO m=-l,l
WRITE(ouctqmc,*)
& REAL(pr_crorb(icrorb,ik,is)%mat_rep(m,
& kp(ik,is)%nb_bot:kp(ik,is)%nb_top))
ENDDO
ENDDO
DO is=1,ns
DO m=-l,l
WRITE(ouctqmc,*)
& AIMAG(pr_crorb(icrorb,ik,is)%mat_rep(m,
& kp(ik,is)%nb_bot:kp(ik,is)%nb_top))
ENDDO
ENDDO
END IF ! End of the ifsplit if-then-else
END IF ! End of the ifmixing if-then-else
END DO ! End of the icrorb loop
END DO ! End of the ik loop
C for each k-point and each correlated orbital, the corresponding projector is described by :
C - the real part of the "correlated" submatrix
C - the imaginary part of the "correlated" submatrix
C
C ----------------------------------------------------------------------------
C Description of the weight of each k_point for the simple point integration :
C ----------------------------------------------------------------------------
DO ik=1,nk
WRITE(ouctqmc,*) kp(ik,1)%weight
C This is a geometrical factor independent of the spin value.
ENDDO
C
C -----------------------------------------------------
C Description of the Hamiltonian H(k) at each k_point :
C -----------------------------------------------------
DO is=1,ns
DO ik=1,nk
C If the calculation is spin-polarized with SO, the numbers for up and dn bands are the same.
IF ((ifSP.AND.ifSO).AND.(is.eq.2)) cycle
DO ib=kp(ik,is)%nb_bot,kp(ik,is)%nb_top
WRITE(ouctqmc,*) kp(ik,is)%eband(ib)
ENDDO
ENDDO
ENDDO
C for each spin value is and each k-point,
C - the energies of the band with spin is at point k
C
C
C ======================================
C Writing the file case.symqmc :
C ======================================
C
WRITE(buf,'(a)')'Writing the file case.symqmc...'
CALL printout(0)
C
C ----------------------------------------------------------
C Description of the general informations about the system :
C ----------------------------------------------------------
WRITE(ousymqmc,'(2(i6,x))') nsym, natom
C nysm is the total number of symmetries in the system
C natom is the total number of atom in the unit cell
C
C -------------------------------------------------------------------
C Description of the permutation matrix associated to each symmetry :
C -------------------------------------------------------------------
DO isym=1,nsym
WRITE(ousymqmc,'(100(i6,x))') srot(isym)%perm(1:natom)
ENDDO
C
C ------------------------------------------------------------
C Description of the time-reversal property for each symmetry :
C ------------------------------------------------------------
IF (ifSP) THEN
ALLOCATE(timeflag(nsym))
timeflag=0
DO isym=1,nsym
IF (srot(isym)%timeinv) timeflag(isym)= 1
ENDDO
WRITE(ousymqmc,'(100(i6,x))') timeflag(1:nsym)
DEALLOCATE(timeflag)
C When the calculation is spin-polarized (with SO), a flag which states
C if a time reversal operation is included in the transformation isym is written.
ENDIF
C
C -----------------------------------------------------------------------------------------
C Description of the representation matrices of each symmetry for each correlated orbital :
C -----------------------------------------------------------------------------------------
DO isym=1,nsym
DO icrorb=1,ncrorb
isrt=crorb(icrorb)%sort
l=crorb(icrorb)%l
C
C The case l=0 is a particular case of "non-mixing" basis.
C --------------------------------------------------------
IF(l==0) THEN
C For the s-orbitals, the only irep possible is the matrix itself.
IF (ifSP.AND.ifSO) THEN
C If SO is taken into account, spinor rotation matrix is considered.
ALLOCATE(spinrot(1:2,1:2))
spinrot(:,:)=0.d0
IF (srot(isym)%timeinv) THEN
C In this case, the Euler angle Beta is Pi. The spinor rotation matrix is block-diagonal,
C since the time-reversal operator is included in the definition of the transformation.
C
factor=srot(isym)%phase/2.d0
C We remind that the field phase is (g-a) if beta=Pi.
C Up/up and Dn/dn terms
spinrot(1,1)=EXP(CMPLX(0.d0,factor))
spinrot(2,2)=CONJG(spinrot(1,1))
C spinrot(1,1) = -exp(+i(alpha-gamma)/2) ; spinrot(2,2) = -exp(-i(alpha-gamma)/2)
C in good agreement with Wien conventions for the definition of this phase factor.
ELSE
C In this case, the Euler angle Beta is 0. The spinor rotation matrix is block-diagonal,
C No time-reversal treatment was applied to the transformation.
C
factor=srot(isym)%phase/2.d0
C We remind that the field phase is (a+g) if beta=0.
C Up/up and Dn/dn terms
spinrot(1,1)=EXP(CMPLX(0.d0,factor))
spinrot(2,2)=CONJG(spinrot(1,1))
C the field phase is 2pi-(alpha+gamma) in this case.
C spinrot(1,1) = -exp(-i(alpha+gamma)/2) ; spinrot(2,2) = -exp(i(alpha-gamma)/2)
C in good agreement with Wien conventions for the definition of this phase factor.
END IF
C Printing the transformation informations
DO m=1,2
WRITE(ousymqmc,*) REAL(spinrot(m,1:2))
ENDDO
DO m=1,2
WRITE(ousymqmc,*) AIMAG(spinrot(m,1:2))
ENDDO
C
DEALLOCATE(spinrot)
C Without SO, only the rotation matrix in orbital space is necessary.
ELSE
C In this case, the Rloc matrix is merely identity.
WRITE(ousymqmc,*) 1.d0
WRITE(ousymqmc,*) 0.d0
ENDIF
C
C If the basis representation needs a complete spinor rotation approach (basis with "mixing" ).
C ---------------------------------------------------------------------------------------------
ELSEIF (reptrans(l,isrt)%ifmixing) THEN
C In this case, the SO is necessary considered, spinor rotation matrices are used.
IF(crorb(icrorb)%ifsplit) THEN
C If only an irep is correlated
ind1=1
DO irep=1,reptrans(l,isrt)%nreps
IF(crorb(icrorb)%correp(irep)) THEN
ind2=ind1+reptrans(l,isrt)%dreps(irep)-1
DO m=ind1,ind2
WRITE(ousymqmc,*) REAL(srot(isym)%
& rotrep(l,isrt)%mat(m,ind1:ind2))
ENDDO
DO m=ind1,ind2
WRITE(ousymqmc,*)AIMAG(srot(isym)%
& rotrep(l,isrt)%mat(m,ind1:ind2))
ENDDO
ENDIF
ind1=ind1+reptrans(l,isrt)%dreps(irep)
ENDDO
ELSE
C If no particular irep is correlated
DO m=1,2*(2*l+1)
WRITE(ousymqmc,*) REAL(srot(isym)%
& rotrep(l,isrt)%mat(m,1:2*(2*l+1)))
ENDDO
DO m=1,2*(2*l+1)
WRITE(ousymqmc,*) AIMAG(srot(isym)%
& rotrep(l,isrt)%mat(m,1:2*(2*l+1)))
ENDDO
ENDIF
C
C If the basis representation can be reduce to the up/up block (basis without "mixing").
C --------------------------------------------------------------------------------------
ELSE
IF (ifSP.AND.ifSO) THEN
C If the calculation is spin-polarized with SO, spinor rotation matrices are considered.
C The spin is not a good quantum number, so the whole representation is used.
ALLOCATE(spinrot(1:2*(2*l+1),1:2*(2*l+1)))
spinrot(:,:)=0.d0
IF (srot(isym)%timeinv) THEN
C In this case, the Euler angle Beta is Pi. The spinor rotation matrix is block-diagonal,
C since the time-reversal operator is included in the definition of the transformation.
C
factor=srot(isym)%phase/2.d0
C We remind that the field phase is (g-a) in this case.
C Up/up block :
ephase=EXP(CMPLX(0.d0,factor))
C As a result, ephase = -exp(i(alpha-gamma)/2)
spinrot(1:2*l+1,1:2*l+1)=
= ephase*srot(isym)%rotrep(l,isrt)%mat(-l:l,-l:l)
C Dn/dn block :
ephase=CONJG(ephase)
C Now, ephase = -exp(-i(alpha-gamma)/2)
spinrot(2*l+2:2*(2*l+1),2*l+2:2*(2*l+1))=
= ephase*srot(isym)%rotrep(l,isrt)%mat(-l:l,-l:l)
ELSE
C In this case, the Euler angle Beta is 0. The spinor rotation matrix is block-diagonal,
C No time-reversal treatment was applied to the transformation.
C
factor=srot(isym)%phase/2.d0
C We remind that the field phase is (a+g) in this case.
C Up/up block :
ephase=EXP(CMPLX(0.d0,factor))
C As a result, ephase = -exp(-i(alpha+gamma)/2)
spinrot(1:2*l+1,1:2*l+1)=
= ephase*srot(isym)%rotrep(l,isrt)%mat(-l:l,-l:l)
C Dn/dn block :
ephase=CONJG(ephase)
C Now, ephase = -exp(i(alpha+gamma)/2)
spinrot(2*l+2:2*(2*l+1),2*l+2:2*(2*l+1))=
= ephase*srot(isym)%rotrep(l,isrt)%mat(-l:l,-l:l)
END IF
C Printing the transformation informations
DO m=1,2*(2*l+1)
WRITE(ousymqmc,*) REAL(spinrot(m,1:2*(2*l+1)) )
ENDDO
DO m=1,2*(2*l+1)
WRITE(ousymqmc,*) AIMAG(spinrot(m,1:2*(2*l+1)) )
ENDDO
C
DEALLOCATE(spinrot)
C In the other cases (spin-polarized without SO or paramagnetic), only the rotation matrix in orbital space is necessary.
ELSE
IF(crorb(icrorb)%ifsplit) THEN
C If only an irep is correlated
ind1=-l
DO irep=1,reptrans(l,isrt)%nreps
IF(crorb(icrorb)%correp(irep)) THEN
ind2=ind1+reptrans(l,isrt)%dreps(irep)-1
DO m=ind1,ind2
WRITE(ousymqmc,*) REAL(srot(isym)%
& rotrep(l,isrt)%mat(m,ind1:ind2))
ENDDO
DO m=ind1,ind2
WRITE(ousymqmc,*) AIMAG(srot(isym)%
& rotrep(l,isrt)%mat(m,ind1:ind2))
ENDDO
ENDIF
ind1=ind1+reptrans(l,isrt)%dreps(irep)
ENDDO
ELSE
C If no particular irep is correlated
DO m=-l,l
WRITE(ousymqmc,*) REAL(srot(isym)%
& rotrep(l,isrt)%mat(m,-l:l))
ENDDO
DO m=-l,l
WRITE(ousymqmc,*) AIMAG(srot(isym)%
& rotrep(l,isrt)%mat(m,-l:l))
ENDDO
END IF ! End of the ifsplit if-then-else
END IF ! End of the ifSO if-then-else
END IF ! End of the ifmixing if-then-else
END DO ! End of the icrorb loop
END DO ! End of the isym loop
C for each symmetry and each correlated orbital icrorb, is described :
C - the real part of the submatrix associated to "isym" for the "icrorb"
C - the imaginary part of the submatrix associated to "isym" for the "icrorb"
C
C -----------------------------------------------------------------------
C Description of the time reversal operator for each correlated orbital :
C -----------------------------------------------------------------------
C This description occurs only if the computation is paramagnetic. (ifSP=0)
IF (.not.ifSP) THEN
DO icrorb=1,ncrorb
l=crorb(icrorb)%l
isrt=crorb(icrorb)%sort
C
C The case l=0 is a particular case of "non-mixing" basis.
C --------------------------------------------------------
IF (l==0) THEN
C For the s-orbitals, the only basis considered is the complex basis.
WRITE(ousymqmc,*) 1.d0
WRITE(ousymqmc,*) 0.d0
C
C If the basis representation can be reduce to the up/up block (basis without "mixing").
C --------------------------------------------------------------------------------------
ELSE
C Calculation of the time-reversal operator
ALLOCATE(time_op(-l:l,-l:l))
time_op(:,:)=0.d0
DO m=-l,l
time_op(m,m)=1.d0
ENDDO
C time_op is Identity.
CALL timeinv_op(time_op,(2*l+1),l,isrt)
C time_op is now the time-reversal operator in the desired basis ({new_i})
C
IF(crorb(icrorb)%ifsplit) THEN
C If only an irep is correlated
ind1=-l
DO irep=1,reptrans(l,isrt)%nreps
IF(crorb(icrorb)%correp(irep)) THEN
ind2=ind1+reptrans(l,isrt)%dreps(irep)-1
DO m=ind1,ind2
WRITE(ousymqmc,*) REAL(time_op(m,ind1:ind2))
ENDDO
DO m=ind1,ind2
WRITE(ousymqmc,*) AIMAG(time_op(m,ind1:ind2))
ENDDO
ENDIF
ind1=ind1+reptrans(l,isrt)%dreps(irep)
ENDDO
ELSE
C If no particular irep is correlated
DO m=-l,l
WRITE(ousymqmc,*) REAL(time_op(m,-l:l))
ENDDO
DO m=-l,l
WRITE(ousymqmc,*) AIMAG(time_op(m,-l:l))
ENDDO
END IF ! End of the ifsplit if-then-else
DEALLOCATE(time_op)
END IF ! End of the l if-then-else
END DO ! End of the icrorb loop
END IF ! End of the ifsp if-then-else
C for each correlated orbital icrorb, the time-reversal operator is described by :
C - its real part
C - its imaginary part
C
C
C ===============================
C Writing the file case.parproj :
C ===============================
C
WRITE(buf,'(a)')'Writing the file case.parproj...'
CALL printout(0)
C
C ----------------------------------------------------------------
C Description of the size of the basis for each included orbital :
C ----------------------------------------------------------------
DO iorb=1,norb
WRITE(oupartial,'(3(i6))') norm_radf(iorb)%n
ENDDO
C There is not more than 1 LO for each orbital (hence n < 4 )
C
C The following descriptions are made for each included orbital.
DO iorb=1,norb
l=orb(iorb)%l
isrt=orb(iorb)%sort
C ------------------------------------
C Description of the Theta projector :
C ------------------------------------
C
C The case l=0 is a particular case of "non-mixing" basis.
C --------------------------------------------------------
IF (l==0) THEN
DO ik=1,nk
DO ir=1,norm_radf(iorb)%n
DO is=1,ns
WRITE(oupartial,*)
& REAL(pr_orb(iorb,ik,is)%matn_rep(1,
& kp(ik,is)%nb_bot:kp(ik,is)%nb_top,ir))
ENDDO
DO is=1,ns
WRITE(oupartial,*)
& AIMAG(pr_orb(iorb,ik,is)%matn_rep(1,
& kp(ik,is)%nb_bot:kp(ik,is)%nb_top,ir))
ENDDO
ENDDO ! End of the ir loop
ENDDO ! End of the ik loop
C
C If the basis representation needs a complete spinor rotation approach (basis with "mixing" ).
C ---------------------------------------------------------------------------------------------
ELSEIF (reptrans(l,isrt)%ifmixing) THEN
C In this case, the calculation is necessary spin-polarized with SO, spinor rotation matrices are used.
DO ik=1,nk
DO ir=1,norm_radf(iorb)%n
DO m=1,2*(2*l+1)
WRITE(oupartial,*)
& REAL(pr_orb(iorb,ik,1)%matn_rep(m,
& kp(ik,1)%nb_bot:kp(ik,1)%nb_top,ir))
ENDDO
DO m=1,2*(2*l+1)
WRITE(oupartial,*)
& AIMAG(pr_orb(iorb,ik,1)%matn_rep(m,
& kp(ik,1)%nb_bot:kp(ik,1)%nb_top,ir))
ENDDO
ENDDO ! End of the ir loop
ENDDO ! End of the ik loop
C
C If the basis representation can be reduce to the up/up block (basis without "mixing").
C --------------------------------------------------------------------------------------
ELSE
DO ik=1,nk
DO ir=1,norm_radf(iorb)%n
DO is=1,ns
DO m=-l,l
WRITE(oupartial,*)
& REAL(pr_orb(iorb,ik,is)%matn_rep(m,
& kp(ik,is)%nb_bot:kp(ik,is)%nb_top,ir))
ENDDO
ENDDO ! End of the is loop
DO is=1,ns
DO m=-l,l
WRITE(oupartial,*)
& AIMAG(pr_orb(iorb,ik,is)%matn_rep(m,
& kp(ik,is)%nb_bot:kp(ik,is)%nb_top,ir))
ENDDO
ENDDO ! End of the is loop
ENDDO ! End of the ir loop
ENDDO ! End of the ik loop
ENDIF ! End of the ifmixing if-then-else
C for each included orbital, for each k-point and each |phi_j> elmt,
C the corresponding Thetaprojector is described by :
C - the real part of the matrix
C - the imaginary part of the matrix
C
C -------------------------------------------------------------------------------
C Description of the density matrices below the lower limit e_bot of the window :
C -------------------------------------------------------------------------------
C
C The case l=0 is a particular case of "non-mixing" basis.
C --------------------------------------------------------
IF (l==0) THEN
C With SO, the density matrix is printed for the complete spin+orbital subspace.
C (in this case, this means a matrix of size 2)
IF (ifSP.AND.ifSO) THEN
ALLOCATE(densprint(2,2))
densprint(1,1)=densmat(1,iorb)%mat(1,1)
densprint(2,2)=densmat(2,iorb)%mat(1,1)
densprint(1,2)=densmat(3,iorb)%mat(1,1)
densprint(2,1)=densmat(4,iorb)%mat(1,1)
DO m=1,2
WRITE(oupartial,*) REAL(densprint(m,1:2))
ENDDO
DO m=1,2
WRITE(oupartial,*) AIMAG(densprint(m,1:2))
ENDDO
DEALLOCATE(densprint)
C In the other cases the density matrix is printed for each spin subspace independently.
C (in this case, this means two matrices of size 1)
ELSE
DO is=1,ns
WRITE(oupartial,*) REAL(densmat(is,iorb)%mat(1,1))
ENDDO
DO is=1,ns
WRITE(oupartial,*) AIMAG(densmat(is,iorb)%mat(1,1))
ENDDO
ENDIF ! End of the ifSO if-then-else
C
C If the basis representation needs a complete spinor rotation approach (basis with "mixing" ).
C ---------------------------------------------------------------------------------------------
ELSEIF (reptrans(l,isrt)%ifmixing) THEN
C In this case, the calculation is necessary spin-polarized with SO, spinor rotation matrices are used.
C As a result, the density matrix is printed for the complete spin+orbital subspace.
C (this means a matrix of size 2*(2*l+1))
DO m=1,2*(2*l+1)
WRITE(oupartial,*)
& REAL(densmat(1,iorb)%mat(m,1:2*(2*l+1)))
ENDDO
DO m=1,2*(2*l+1)
WRITE(oupartial,*)
& AIMAG(densmat(1,iorb)%mat(m,1:2*(2*l+1)))
ENDDO
C
C If the basis representation can be reduce to the up/up block (basis without "mixing").
C --------------------------------------------------------------------------------------
ELSE
C With SO, the density matrix is printed for the complete spin+orbital subspace.
C (this means a matrix of size 2*(2*l+1))
IF (ifSP.AND.ifSO) THEN
ALLOCATE(densprint(1:2*(2*l+1),1:2*(2*l+1)))
densprint(1:2*l+1,1:2*l+1)=
& densmat(1,iorb)%mat(-l:l,-l:l)
densprint(2*l+2:2*(2*l+1),2*l+2:2*(2*l+1))=
& densmat(2,iorb)%mat(-l:l,-l:l)
densprint(1:2*l+1,2*l+2:2*(2*l+1))=
& densmat(3,iorb)%mat(-l:l,-l:l)
densprint(2*l+2:2*(2*l+1),1:2*l+1)=
& densmat(4,iorb)%mat(-l:l,-l:l)
DO m=1,2*(2*l+1)
WRITE(oupartial,*) REAL(densprint(m,1:2*(2*l+1)))
ENDDO
DO m=1,2*(2*l+1)
WRITE(oupartial,*) AIMAG(densprint(m,1:2*(2*l+1)))
ENDDO
DEALLOCATE(densprint)
C In the other cases, the density matrix is printed for each spin subspace independently.
C (this means two matrices of size 2*l+1)
ELSE
DO is=1,ns
DO m=-l,l
WRITE(oupartial,*)
& REAL(densmat(is,iorb)%mat(m,-l:l))
ENDDO
ENDDO
DO is=1,ns
DO m=-l,l
WRITE(oupartial,*)
& AIMAG(densmat(is,iorb)%mat(m,-l:l))
ENDDO
ENDDO
ENDIF ! End of the ifSO if-then-else
ENDIF ! End of the ifmixing if-then-else
C for each included orbital, the corresponding density matrix is described by :
C - the real part of the matrix
C - the imaginary part of the matrix
C
C --------------------------------------------------------------------
C Description of the global to local coordinates transformation Rloc :
C --------------------------------------------------------------------
C Description of each transformation Rloc
iatom=orb(iorb)%atom
C
C The case l=0 is a particular case of "non-mixing" basis.
C --------------------------------------------------------
IF(l==0) THEN
C For the s-orbitals, the only irep possible is the matrix itself.
IF (ifSP.AND.ifSO) THEN
C If SO is taken into account, spinor rotation matrix must be considered.
ALLOCATE(spinrot(1:2,1:2))
spinrot(:,:)=0.d0
C The spinor-rotation matrix is directly calculated from the Euler angles a,b and c.
IF (rotloc(iatom)%timeinv) THEN
factor=(rotloc(iatom)%a+rotloc(iatom)%g)/2.d0
spinrot(2,1)=EXP(CMPLX(0.d0,factor))*
& DCOS(rotloc(iatom)%b/2.d0)
spinrot(1,2)=-CONJG(spinrot(2,1))
C Up/dn and Dn/up terms
factor=-(rotloc(iatom)%a-rotloc(iatom)%g)/2.d0
spinrot(2,2)=-EXP(CMPLX(0.d0,factor))*
& DSIN(rotloc(iatom)%b/2.d0)
spinrot(1,1)=CONJG(spinrot(2,2))
ELSE
factor=(rotloc(iatom)%a+rotloc(iatom)%g)/2.d0
spinrot(1,1)=EXP(CMPLX(0.d0,factor))*
& DCOS(rotloc(iatom)%b/2.d0)
spinrot(2,2)=CONJG(spinrot(1,1))
C Up/dn and Dn/up terms
factor=-(rotloc(iatom)%a-rotloc(iatom)%g)/2.d0
spinrot(1,2)=EXP(CMPLX(0.d0,factor))*
& DSIN(rotloc(iatom)%b/2.d0)
spinrot(2,1)=-CONJG(spinrot(1,2))
ENDIF
C Printing the transformation informations
DO m=1,2
WRITE(oupartial,*) REAL(spinrot(m,1:2))
ENDDO
DO m=1,2
WRITE(oupartial,*) AIMAG(spinrot(m,1:2))
ENDDO
DEALLOCATE(spinrot)
C Without SO, only the rotation matrix in orbital space is necessary.
ELSE
C In this case, the Rloc matrix is merely identity.
WRITE(oupartial,*) 1.d0
WRITE(oupartial,*) 0.d0
ENDIF
C Printing the time inversion flag if the calculation is spin-polarized (ifSP=1)
IF (ifSP) THEN
timeinvflag=0
IF (rotloc(iatom)%timeinv) timeinvflag=1
WRITE(oupartial,'(i6)') timeinvflag
ENDIF
C
C If the basis representation needs a complete spinor rotation approach (basis with "mixing" ).
C ---------------------------------------------------------------------------------------------
ELSEIF (reptrans(l,isrt)%ifmixing) THEN
C In this case, the calculation is necessary spin-polarized with SO, spinor rotation matrices are used.
DO m=1,2*(2*l+1)
WRITE(oupartial,*) REAL(rotloc(iatom)%
& rotrep(l)%mat(m,1:2*(2*l+1)))
ENDDO
DO m=1,2*(2*l+1)
WRITE(oupartial,*) AIMAG(rotloc(iatom)%
& rotrep(l)%mat(m,1:2*(2*l+1)))
ENDDO
C Printing if the transforamtion included a time-reversal operation.
timeinvflag =0
IF (rotloc(iatom)%timeinv) timeinvflag=1
WRITE(oupartial,'(i6)') timeinvflag
C
C If the basis representation can be reduce to the up/up block (basis without "mixing").
C --------------------------------------------------------------------------------------
ELSE
IF (ifSP.AND.ifSO) THEN
C If SO is taken into account, spinor rotation matrices are considered.
C The spin is not a good quantum number, so the whole representation is used.
DO m=1,2*(2*l+1)
WRITE(oupartial,*) REAL(rotloc(iatom)%
& rotrep(l)%mat(m,1:2*(2*l+1)))
ENDDO
DO m=1,2*(2*l+1)
WRITE(oupartial,*) AIMAG(rotloc(iatom)%
& rotrep(l)%mat(m,1:2*(2*l+1)))
ENDDO
ELSE
C The calculation is either spin-polarized without SO or paramagnetic.
C The spin is a good quantum number and irep are possible.
DO m=-l,l
WRITE(oupartial,*) REAL(rotloc(iatom)%
& rotrep(l)%mat(m,-l:l))
ENDDO
DO m=-l,l
WRITE(oupartial,*) AIMAG(rotloc(iatom)%
& rotrep(l)%mat(m,-l:l))
ENDDO
END IF ! End of the ifSO if-then-else
C Printing if the transformation included a time-reversal operation.
IF (ifSP) THEN
timeinvflag =0
IF (rotloc(iatom)%timeinv) timeinvflag=1
WRITE(oupartial,'(i6)') timeinvflag
END IF
END IF ! End of the ifmixing if-then-else
C for each included orbital iorb, the transformation Rloc is described by :
C - the real part of the submatrix associated to "iorb"
C - the imaginary part of the submatrix associated to "iorb"
C - a flag which states if a time reversal operation is included in the transformation (if SP only )
C
END DO ! End of the iorb loop
C
C
C ==============================
C Writing the file case.sympar :
C ==============================
C
WRITE(buf,'(a)')'Writing the file case.sympar...'
CALL printout(0)
C
C ----------------------------------------------------------
C Description of the general informations about the system :
C ----------------------------------------------------------
WRITE(ousympar,'(2(i6,x))') nsym, natom
C nysm is the total number of symmetries of the system
C natom is the total number of atom in the unit cell
C
C -------------------------------------------------------------------
C Description of the permutation matrix associated to each symmetry :
C -------------------------------------------------------------------
DO is=1,nsym
WRITE(ousympar,'(100(i6,x))') srot(is)%perm
ENDDO
C
C ------------------------------------------------------------
C Description of the time-reversal property for each symmetry :
C ------------------------------------------------------------
IF (ifSP) THEN
ALLOCATE(timeflag(nsym))
timeflag=0
DO isym=1,nsym
IF (srot(isym)%timeinv) timeflag(isym)= 1
ENDDO
WRITE(ousympar,'(100(i6,x))') timeflag(1:nsym)
DEALLOCATE(timeflag)
C When the calculation is spin-polarized (with SO), a flag which states
C if a time reversal operation is included in the transformation isym is written.
ENDIF
C
C -------------------------------------------------------------------------------------------
C Description of the representation matrices of each symmetry for all the included orbitals :
C -------------------------------------------------------------------------------------------
DO isym=1,nsym
DO iorb=1,norb
isrt=orb(iorb)%sort
l=orb(iorb)%l
C
C The case l=0 is a particular case of "non-mixing" basis.
C --------------------------------------------------------
IF(l==0) THEN
C For the s-orbitals, the only irep possible is the matrix itself.
IF (ifSP.AND.ifSO) THEN
C If SO is taken into account, spinor rotation matrix is considered.
ALLOCATE(spinrot(1:2,1:2))
spinrot(:,:)=0.d0
IF (srot(isym)%timeinv) THEN
C In this case, the Euler angle Beta is Pi. The spinor rotation matrix is block-diagonal,
C since the time-reversal operator is included in the definition of the transformation.
C
factor=srot(isym)%phase/2.d0
C We remind that the field phase is (g-a) if beta=Pi.
C Up/up and Dn/dn terms
spinrot(1,1)=EXP(CMPLX(0.d0,factor))
spinrot(2,2)=CONJG(spinrot(1,1))
C spinrot(1,1) = -exp(i(alpha-gamma)/2) ; spinrot(2,2) = -exp(-i(alpha-gamma)/2)
C in good agreement with Wien conventions for the definition of this phase factor.
ELSE
C In this case, the Euler angle Beta is 0. The spinor rotation matrix is block-diagonal,
C No time-reversal treatment was applied to the transformation.
C
factor=srot(isym)%phase/2.d0
C We remind that the field phase is (a+g) if beta=0.
C Up/up and Dn/dn terms
spinrot(1,1)=EXP(CMPLX(0.d0,factor))
spinrot(2,2)=CONJG(spinrot(1,1))
C spinrot(1,1) = -exp(-i(alpha+gamma)/2) ; spinrot(2,2) = -exp(i(alpha-gamma)/2)
C in good agreement with Wien conventions for the definition of this phase factor.
END IF
C Printing the transformation informations
DO m=1,2
WRITE(ousympar,*) REAL(spinrot(m,1:2))
ENDDO
DO m=1,2
WRITE(ousympar,*) AIMAG(spinrot(m,1:2))
ENDDO
C
DEALLOCATE(spinrot)
C Without SO, only the rotation matrix in orbital space is necessary.
ELSE
C In this case, the Rloc matrix is merely identity.
WRITE(ousympar,*) 1.d0
WRITE(ousympar,*) 0.d0
ENDIF
C
C If the basis representation needs a complete spinor rotation approach (basis with "mixing" ).
C ---------------------------------------------------------------------------------------------
ELSEIF (reptrans(l,isrt)%ifmixing) THEN
C In this case, the SO is necessary considered, spinor rotation matrices are used.
DO m=1,2*(2*l+1)
WRITE(ousympar,*) REAL(srot(isym)%
& rotrep(l,isrt)%mat(m,1:2*(2*l+1)))
ENDDO
DO m=1,2*(2*l+1)
WRITE(ousympar,*) AIMAG(srot(isym)%
& rotrep(l,isrt)%mat(m,1:2*(2*l+1)))
ENDDO
C
C If the basis representation can be reduce to the up/up block (basis without "mixing").
C --------------------------------------------------------------------------------------
ELSE
IF (ifSP.AND.ifSO) THEN
C If the calculation is spin-polarized with SO, spinor rotation matrices are considered.
C The spin is not a good quantum number, so the whole representation is used.
ALLOCATE(spinrot(1:2*(2*l+1),1:2*(2*l+1)))
spinrot(:,:)=0.d0
IF (srot(isym)%timeinv) THEN
C In this case, the Euler angle Beta is Pi. The spinor rotation matrix is block-diagonal,
C since the time-reversal operator is included in the definition of the transformation.
C
factor=srot(isym)%phase/2.d0
C We remind that the field phase is (g-a) in this case.
C Up/up block :
ephase=EXP(CMPLX(0.d0,factor))
C AS a result, ephase = -exp(i(alpha-gamma)/2)
spinrot(1:2*l+1,1:2*l+1)=
= ephase*srot(isym)%rotrep(l,isrt)%mat(-l:l,-l:l)
C Dn/dn block :
ephase=CONJG(ephase)
C Now, ephase = -exp(-i(alpha-gamma)/2)
spinrot(2*l+2:2*(2*l+1),2*l+2:2*(2*l+1))=
= ephase*srot(isym)%rotrep(l,isrt)%mat(-l:l,-l:l)
ELSE
C In this case, the Euler angle Beta is 0. The spinor rotation matrix is block-diagonal,
C No time-reversal treatment was applied to the transformation.
C
factor=srot(isym)%phase/2.d0
C We remind that the field phase is 2pi-(alpha+gamma) in this case.
C Up/up block :
ephase=EXP(CMPLX(0.d0,factor))
C As a result, ephase = -exp(-i(alpha+gamma)/2)
spinrot(1:2*l+1,1:2*l+1)=
= ephase*srot(isym)%rotrep(l,isrt)%mat(-l:l,-l:l)
C Dn/dn block :
ephase=CONJG(ephase)
C Now, ephase = -exp(i(alpha+gamma)/2)
spinrot(2*l+2:2*(2*l+1),2*l+2:2*(2*l+1))=
= ephase*srot(isym)%rotrep(l,isrt)%mat(-l:l,-l:l)
END IF
C Printing the transformation informations
DO m=1,2*(2*l+1)
WRITE(ousympar,*) REAL(spinrot(m,1:2*(2*l+1)) )
ENDDO
DO m=1,2*(2*l+1)
WRITE(ousympar,*) AIMAG(spinrot(m,1:2*(2*l+1)) )
ENDDO
C
DEALLOCATE(spinrot)
C In the other cases (spin-polarized without SO or paramagnetic), only the rotation matrix in orbital space is necessary.
ELSE
DO m=-l,l
WRITE(ousympar,*) REAL(srot(isym)%
& rotrep(l,isrt)%mat(m,-l:l))
ENDDO
DO m=-l,l
WRITE(ousympar,*) AIMAG(srot(isym)%
& rotrep(l,isrt)%mat(m,-l:l))
ENDDO
END IF ! End of the ifSO if-then-else
END IF ! End of the ifmixing if-then-else
END DO ! End of the iorb loop
END DO ! End of the isym loop
C for each symmetry and each included orbital iorb, is described :
C - the real part of the matrix associated to "isym" for the "iorb"
C - the imaginary part of the matrix associated to "isym" for the "iorb"
C
C ---------------------------------------------------------------------
C Description of the time reversal operator for each included orbital :
C ---------------------------------------------------------------------
C This description occurs only if the computation is paramagnetic (ifSP=0)
IF (.not.ifSP) THEN
DO iorb=1,norb
l=orb(iorb)%l
isrt=orb(iorb)%sort
C
C The case l=0 is a particular case of "non-mixing" basis.
C --------------------------------------------------------
IF (l==0) THEN
C For the s-orbitals, the only basis considered is the complex basis.
WRITE(ousympar,*) 1.d0
WRITE(ousympar,*) 0.d0
C
C If the basis representation can be reduce to the up/up block (basis without "mixing").
C --------------------------------------------------------------------------------------
ELSE
C Calculation of the time-reversal operator
ALLOCATE(time_op(-l:l,-l:l))
time_op(:,:)=0.d0
DO m=-l,l
time_op(m,m)=1.d0
ENDDO
C time_op is Identity.
CALL timeinv_op(time_op,(2*l+1),l,isrt)
C time_op is now the time-reversal operator in the desired basis ({new_i})
C
DO m=-l,l
WRITE(ousympar,*) REAL(time_op(m,-l:l))
ENDDO
DO m=-l,l
WRITE(ousympar,*) AIMAG(time_op(m,-l:l))
ENDDO
DEALLOCATE(time_op)
END IF ! End of the l if-then-else
END DO ! End of the icrorb loop
END IF ! End of the ifsp if-then-else
C for each included orbital iorb, the time-reversal operator is described by :
C - its real part
C - its imaginary part
C
RETURN
END
Reference in New Issue View Git Blame Copy Permalink