mirror of
https://github.com/triqs/dft_tools
synced 2024-12-22 12:23:41 +01:00
226 lines
8.5 KiB
Fortran
226 lines
8.5 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 orthogonal_h(s1,ndim,inv)
|
|
C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
C %% %%
|
|
C %% This subroutine computes : %%
|
|
C %% - if inv = .FALSE. the square root of the Hermitian matrix s1 %%
|
|
C %% - if inv = .TRUE. the inverse of the square root of the %%
|
|
C %% Hermitian matrix s1 %%
|
|
C %% %%
|
|
C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
C Definiton of the variables :
|
|
C ----------------------------
|
|
USE prnt
|
|
IMPLICIT NONE
|
|
INTEGER :: ndim, INFO, lm, lm1
|
|
COMPLEX(KIND=8), DIMENSION(ndim) :: WORK
|
|
COMPLEX(KIND=8), DIMENSION(ndim,ndim) :: s1
|
|
INTEGER, DIMENSION(ndim,ndim) :: IPIV
|
|
LOGICAL :: inv
|
|
C
|
|
C Calculation of S1^(1/2) or S1^(-1/2):
|
|
C -------------------------------------
|
|
CALL sqrtm(s1,ndim,inv)
|
|
C The resulting matrix is stored in s1.
|
|
RETURN
|
|
END
|
|
|
|
SUBROUTINE orthogonal_r(s2,ndim,inv)
|
|
C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
C %% %%
|
|
C %% This subroutine computes : %%
|
|
C %% - if inv = .FALSE. the square root of s1 %%
|
|
C %% - if inv = .TRUE. the inverse of the square root of s2 %%
|
|
C %% where s2 is a real symmetric matrix. %%
|
|
C %% %%
|
|
C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
C Definiton of the variables :
|
|
C ----------------------------
|
|
USE prnt
|
|
IMPLICIT NONE
|
|
INTEGER :: ndim, INFO, lm, lm1
|
|
COMPLEX(KIND=8), DIMENSION(ndim) :: WORK
|
|
COMPLEX(KIND=8), DIMENSION(ndim,ndim) :: s1
|
|
REAL(KIND=8), DIMENSION(ndim,ndim) :: s2
|
|
INTEGER, DIMENSION(ndim,ndim) :: IPIV
|
|
LOGICAL :: inv
|
|
C
|
|
C Calculation of S2^(1/2) or S2^(-1/2):
|
|
C -------------------------------------
|
|
s1=s2
|
|
CALL sqrtm(s1,ndim,inv)
|
|
s2=REAL(s1)
|
|
C The resulting matrix is stored in s2.
|
|
RETURN
|
|
END
|
|
|
|
|
|
SUBROUTINE sqrtm(cmat,m,inv)
|
|
C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
C %% %%
|
|
C %% This subroutine calculates the square root of a positively %%
|
|
C %% defined Hermitian matrix A=cmat using the decomposition %%
|
|
C %% A=Z*D*Z^H %%
|
|
C %% where D is a diagonal matrix of eigenvalues of A, %%
|
|
C %% Z is matrix of orthonormal eigenvectors of A, %%
|
|
C %% Z^H is its Hermitian conjugate. %%
|
|
C %% Then A^(1/2)=Z*D^(1/2)*Z^H. %%
|
|
C %% Correction: the matrix A is allowed to be negatively defined. %%
|
|
C %% %%
|
|
C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
C Definiton of the variables :
|
|
C ----------------------------
|
|
IMPLICIT NONE
|
|
INTEGER :: m
|
|
COMPLEX(KIND=8), DIMENSION(m,m):: cmat, D, D1
|
|
LOGICAL :: inv
|
|
C Calculation of Z*D^(1/2):
|
|
C -------------------------
|
|
CALL sqrt_eigenvec(cmat,D1,m,inv)
|
|
WRITE(95,*) cmat
|
|
WRITE(95,*) ' '
|
|
WRITE(95,*) D1
|
|
WRITE(95,*) ' '
|
|
C Calculation of A^(1/2)=Z*D^(1/2)*Z^H:
|
|
C -------------------------------------
|
|
D=CONJG(cmat)
|
|
call ZGEMM('N','T',m,m,m,DCMPLX(1.D0,0.D0),D1,
|
|
& m,D,m,DCMPLX(0.D0,0.D0),cmat,m)
|
|
C The resulting matrix is stored in cmat.
|
|
RETURN
|
|
END
|
|
|
|
|
|
SUBROUTINE sqrt_eigenvec(cmat,D1,m,inv)
|
|
C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
C %% %%
|
|
C %% This subroutine computes : %%
|
|
C %% - if inv = .FALSE. Z*D^(1/2) %%
|
|
C %% - if inv = .TRUE. Z*D^(-1/2) %%
|
|
C %% where Z is a matrix of orthonormal eigenvectors of cmat and %%
|
|
C %% D is the diagonal matrix of cmat's eigenvalues. %%
|
|
C %% %%
|
|
C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
C Definiton of the variables :
|
|
C ----------------------------
|
|
USE prnt
|
|
IMPLICIT NONE
|
|
LOGICAL :: inv, ifwrite
|
|
INTEGER :: m, INFO, i, j
|
|
INTEGER, PARAMETER :: nwork=40
|
|
C
|
|
COMPLEX(KIND=8), allocatable, DIMENSION(:) :: WORK
|
|
COMPLEX(KIND=8), DIMENSION(m,m) :: cmat, D1
|
|
REAL(KIND=8), DIMENSION(m) :: W
|
|
COMPLEX(KIND=8), DIMENSION(m) :: W_comp
|
|
REAL(KIND=8), allocatable, DIMENSION(:) :: RWORK
|
|
C
|
|
C Finding the eigenvalues and the eigenvectors of cmat :
|
|
C ------------------------------------------------------
|
|
ALLOCATE(rwork(3*m-2))
|
|
ALLOCATE(work(2*m-1))
|
|
CALL ZHEEV('V', 'U', m, cmat, m, W, WORK,2*m-1,RWORK,INFO)
|
|
IF (info.ne.0) THEN
|
|
WRITE(buf,'(a)')
|
|
& 'The subroutine zheev ends with info = ',info
|
|
CALL printout(0)
|
|
WRITE(buf,'(a)')'In sqrt_eigenvec, a pbm occurs in zheev.'
|
|
CALL printout(0)
|
|
WRITE(buf,'(a)')'END OF THE PRGM'
|
|
CALL printout(0)
|
|
STOP
|
|
ENDIF
|
|
C W contains the eigenvalues of cmat.
|
|
W_comp=CMPLX(W,0d0)
|
|
C
|
|
C Checking of the validity of the computation :
|
|
C ---------------------------------------------
|
|
ifwrite=.FALSE.
|
|
DO j=1,m
|
|
C The warning is written only once in the file case.outdmftpr
|
|
IF (ifwrite) EXIT
|
|
C Checking if the eigenvalues are not negative.
|
|
IF (W(j).lt.0.d0) THEN
|
|
WRITE(buf,'(a,i2,a,a)')
|
|
& 'WARNING : An eigenvalue (',j,') of the ',
|
|
& 'overlap matrix is negative.'
|
|
CALL printout(0)
|
|
WRITE(buf,'(a,a)')' The result ',
|
|
& 'of the calculation may thus be wrong.'
|
|
CALL printout(1)
|
|
ifwrite=.TRUE.
|
|
ENDIF
|
|
IF (ABS(W(j)).lt.1.d-12) THEN
|
|
WRITE(buf,'(a,i2,a,a)')
|
|
& 'WARNING : An eigenvalue (',j,') of the ',
|
|
& 'overlap matrix is almost zero.'
|
|
CALL printout(0)
|
|
WRITE(buf,'(a,a)')' The result ',
|
|
& 'of the calculation may thus be wrong.'
|
|
CALL printout(1)
|
|
ifwrite=.TRUE.
|
|
ENDIF
|
|
ENDDO
|
|
C
|
|
C Calculation of Z*D^(1/2) :
|
|
C --------------------------
|
|
C The result is stored in D1.
|
|
IF(.NOT.inv) THEN
|
|
DO i=1,m
|
|
DO j=1,m
|
|
D1(i,j)=cmat(i,j)*SQRT(W_comp(j))
|
|
ENDDO
|
|
ENDDO
|
|
ELSE
|
|
C Calculation of Z*D^(-1/2) :
|
|
C ---------------------------
|
|
C The result is stored in D1.
|
|
DO i=1,m
|
|
DO j=1,m
|
|
IF (ABS(W(j))==0.d0) THEN
|
|
WRITE(buf,'(a,i2,a)')
|
|
& 'An eigenvalue (',j,') of the ',
|
|
& 'overlap matrix has the value 0.'
|
|
CALL printout(0)
|
|
WRITE(buf,'(a)')
|
|
& 'The calculation can not be performed further.'
|
|
CALL printout(0)
|
|
CALL printout(0)
|
|
WRITE(buf,'(a)')'END OF THE PRGM'
|
|
CALL printout(0)
|
|
STOP
|
|
ENDIF
|
|
D1(i,j)=cmat(i,j)/SQRT(W_comp(j))
|
|
ENDDO
|
|
ENDDO
|
|
ENDIF
|
|
C The resulting matrix is stored in D1 and cmat is now Z.
|
|
RETURN
|
|
END
|
|
|