orthoglo.F90 2.37 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
 MODULE m_orthoglo
      use m_juDFT
!*********************************************************************
! Each G-vector corresponds to a vector of C-coeff. These vectors must
! be linearly independent. This is checked by this soubroutine for an
! atom that doesn't have an inversion partner.
! Philipp Kurz 99/04
!*********************************************************************
      CONTAINS
      SUBROUTINE orthoglo(atoms, nkvec,lo,l,linindq,l_lo2, cwork, linind)
!
!*************** ABBREVIATIONS ***************************************
! cwork   : contains the vectors of C-coeff.
! l_lo2   : changes this routine to old 'orthoglo2': same as orthoglo, 
!           but for a pair of atoms that can be mapped onto eachother 
!           by inversion.
! CF Replaced (unstable) Gram-Schmidt by diagonalization.
!*********************************************************************
!
#include"cpp_double.h"
!
        USE m_types
      IMPLICIT NONE
      TYPE(t_atoms),INTENT(IN)   :: atoms
!     ..
!     .. Scalar Arguments ..
      INTEGER, INTENT (IN) :: l,lo,nkvec
      REAL,    INTENT (IN) :: linindq
      LOGICAL, INTENT (IN) :: l_lo2
      LOGICAL, INTENT (OUT) :: linind
!     ..
!     .. Array Arguments ..
      COMPLEX,INTENT (INOUT):: cwork(-2*atoms%llod:2*atoms%llod+1,2*(2*atoms%llod+1) ,atoms%nlod)
!     ..
!     .. Local Scalars ..
      INTEGER dim,low,i,j
!     ..
!     .. Local Arrays ..
      REAL eig(nkvec),rwork(3*nkvec)
#ifdef CPP_INVERSION
      REAL olap(nkvec,nkvec)
      EXTERNAL CPP_LAPACK_ssyev
#else
      COMPLEX olap(nkvec,nkvec),work(2*nkvec)
      EXTERNAL CPP_LAPACK_cheev
#endif
      
      IF (l_lo2) THEN
        dim = 2* (2*l+1)
        low = -2*l
      ELSE
        dim = 2*l+1
        low = -l
      ENDIF

      DO i = 1,nkvec
        DO j = 1,nkvec
          olap(i,j) = dot_product(cwork(low:low+dim-1,i,lo), cwork(low:low+dim-1,j,lo))
        ENDDO
      ENDDO
#ifdef CPP_INVERSION
      CALL CPP_LAPACK_ssyev('N','U',nkvec,olap,nkvec,eig, rwork,3*nkvec,i)
      IF(i/=0)  CALL juDFT_error("(S,D)SYEV failed.","orthoglo")
#else
      CALL CPP_LAPACK_cheev('N','U',nkvec,olap,nkvec,eig, work,2*nkvec,rwork,i)
      IF(i/=0)  CALL juDFT_error("(C,Z)HEEV failed.","orthoglo")
#endif
      IF(eig(1).LT.linindq) THEN
        linind = .false.
      ELSE
        linind = .true.
      ENDIF
      RETURN
      
      END SUBROUTINE orthoglo
      END MODULE m_orthoglo