C
C                   SUB-SPACE ROTATION
C
C        FOR VAX ONLY   -->  USE SUBROTC.FOR FOR CRAY
C
C======================================================================
C     FORM A SUBSPACE HAMILTONIAN MATRIX USING THE ORTHOGONAL FILLED
C     BANDS AS BASIS FUNCTIONS. THE ELEMENTS OF THIS MATRIX ARE SIMPLY
C     THE LAGRANGE MULTIPLIERS FOR ORTHONORMALITY. THE EIGENVECTORS OF 
C     THIS SMALL MATRIX WILL GIVE THE LINEAR COMBINATIONS OF THE TRIAL
C     VECTORS WHICH ARE EIGENVECTORS IN THEIR SUBSPACE. 
C
C     THIS SIMILARITY TRANSFORMATION CHANGES NEITHER THE ENERGY NOR 
C     THE DENSITY OF THE SYSTEM. THE PROCESS NOT ONLY SOLVES ANY 
C     PROBLEMS ASSOCIATED WITH BAND ORDERING, BUT ALSO GUARANTEES THAT
C     ANY SUBSEQUENT CHANGE TO ONE BANDWHICH LOWERS ITS ENERGY WILL BE
C     ORTHOGONAL TO ALL OF THE OTHER BANDS.
C     (TETER, ET AL. PHYS. REV. B40, p12261)        X. WENG 27-FEB-90
C======================================================================
      SUBROUTINE SUBROT(NBANDS,NKPTS,NPLWV,MPLWV, NRPLWV,NINDPW,
     &      NPLWKP,CV, CPTWFP,CPTWFL,VOLC,CELEN,
     &      VNL, NGPTAR, DATAKE,CPTOWR,CPTNWR, CWORK, CGRA,
     &      NKP, HR, HI, AUX, FV1, FV2, FV3, CH0, NSPEC, NIONS, 
     &      NIONSP, PSCALE, 
     &      DNLKG, CPHSGR,VGNL, CELFRC, CWRK20,
     &      CWRK21, CWRK22, CWRK23, IVPTYP, IVPTYN,IPRINT,NLPOT,
     &      NGX,NGY,NGZ,NRGRPT,NRLPPI,CPHGRD,NRLNL,NIONST,MXRLNL,
     &      IRLNL,PRLSCA,VRLGRD,NADGRD,MXRLSH,CESAVE)
      IMPLICIT COMPLEX (C)
      DIMENSION NINDPW(*)
      DIMENSION CV(*)
      DIMENSION CPTWFP(*)
      DIMENSION CPTWFL(*)
      DIMENSION CELEN(*)
      DIMENSION NGPTAR(*)
      DIMENSION DATAKE(*)
      DIMENSION NIONSP(*)
C=======================================================================
C
C                      DIMENSION STATEMENTS
C
C CPTACC(MPLWV) = THE PRODUCT OF THE TOTAL POTENTIAL AND THE REAL SPACE
C          WAVEFUNCTION, THIS IS THEN FOURIER TRANSFORMED TO RECIPROCAL
C          SPACE
C
C CWORK(MPLWV) = A WORK ARRAY USED IN THE FOURIER TRANSFORM
C=======================================================================
      DIMENSION CPTOWR(*)
      DIMENSION CPTNWR(*)
      DIMENSION CWORK(*)
      DIMENSION CGRA(*)
C=======================================================================
C DIMENSION FOR SUBROUTINE ROTATION (WENG, 27-FEB-90)
C=======================================================================
      DIMENSION HR(NBANDS,NBANDS),HI(NBANDS,NBANDS),AUX(NBANDS),
     &     FV1(NBANDS),FV2(NBANDS),FV3(NBANDS),CH0(NBANDS,NBANDS)
C=======================================================================
C DIMENSION STATEMENTS FOR THE NON-LOCAL CALCULATION
C=======================================================================
      DIMENSION PSCALE(0:2,NSPEC)
      DIMENSION DNLKG(NRPLWV,0:3,NKPTS)
      DIMENSION IVPTYN(NSPEC)
      DIMENSION CPHSGR(NRPLWV,NIONS,NSPEC)
      DIMENSION VGNL(NRPLWV,0:2,NSPEC)
      DIMENSION CELFRC(NRPLWV)
      DIMENSION CWRK20(NIONS),CWRK21(3,NIONS)
      DIMENSION CWRK22(NIONS),CWRK23(3,3,NIONS)
C=====================================================================
C REAL SPACE NON-LOCAL                                                
C=====================================================================
      DIMENSION NRLPPI(NIONST),CPHGRD(NRGRPT,NIONST)                  
      DIMENSION NRLNL(NSPEC),IRLNL(MXRLNL,NSPEC),PRLSCA(MXRLNL,NSPEC) 
      DIMENSION VRLGRD(NRGRPT,MXRLSH,NIONST),NADGRD(NRGRPT,NIONST)    
      DIMENSION VNL(NBANDS,NKPTS)                                     
      DIMENSION CESAVE(NIONST,36)
C
      RINPLW=1.0/NPLWV
C======================================================================
C
C BECAUSE H(IB,JB) MUST BE HERMITIAN, WE ONLY NEED TO CALCULATE 
C H(IB,JB) FOR JB.GE.IB
C
C======================================================================
C 1. CALCULATE H*PSI
C======================================================================
      DO 10 IB=1,NBANDS
      NINDW=NRPLWV*(IB-1)
      NINDKE=0
      DO 15 M=1,MPLWV
      CPTOWR(M)=(0.0,0.0)
      CPTNWR(M)=(0.0,0.0)
      CWORK(M)=(0.0,0.0)
  15  CONTINUE
      DO 20  M=1,NPLWKP
      CPTOWR(NINDPW(M))=CPTWFP(M+NINDW)
  20  CONTINUE
C=======================================================================
C TRANSFORM THE WAVEFUNCTION INTO REAL SPACE
C=======================================================================
      CALL FFT3D(CPTOWR,CWORK,NGPTAR,1)
      IF(IVPTYP.EQ.1 .AND. NLPOT.EQ.1) THEN
C=======================================================================
C CALL THE ROUTINE THAT CALCULATES THE PRODUCT OF THE NON-LOCAL
C POTENTIAL AND THE WAVEFUNCTION IN REAL SPACE
C=======================================================================
      CALL VSINL(VOLC,NGX,NGY,NGZ,MPLWV,NRPLWV,
     &   NRGRPT,NIONSP,NRLPPI,NSPEC,CWORK,CPTOWR,CPTNWR,   
     &   CPHGRD,NRLNL,NIONST,MXRLNL,IRLNL,PRLSCA,VRLGRD,NADGRD,MXRLSH,
     &   CESAVE)                      
      ENDIF                                                             
C=======================================================================
C CALCULATE (WAVEFUNCTION*POTENTIAL) AND FFT IT INTO RECIPROCAL SPACE
C=======================================================================
      DO 25 NNN=1,NPLWV
      CPTNWR(NNN)=(CPTNWR(NNN)+CV(NNN)*CPTOWR(NNN))*RINPLW
  25  CONTINUE
      CALL FFT3D(CPTNWR,CWORK,NGPTAR,-1)
C=======================================================================
C ADDED THE KINETIC ENERGY TERM IN RECIPROCAL SPACE
C=======================================================================
      DO 30 M=1,NPLWKP
      CGRA(M)=CPTNWR(NINDPW(M))+CPTWFP(M+NINDW)*DATAKE(1+(7*(M-1)))
  30  CONTINUE
C=======================================================================
C 2. PSI(DAGGER) H PSI
C=======================================================================
      DO 40 JB=IB,NBANDS
      NINKI=NRPLWV*(JB-1)
      CHD=(0.0,0.0)
      DO 45 KL=1,NPLWKP
      CHD = CHD + CONJG(CPTWFP(KL+NINKI))*CGRA(KL)
  45  CONTINUE
      HR(IB,JB)=REAL(CHD)
      HI(IB,JB)=AIMAG(CHD)
  40  CONTINUE      
  10  CONTINUE 
C=========================================================================
C 3. THE CONTRIBUTION TO THE HAMILTONIAN FROM THE NON-LOCAL PARTS
C    OF PSEUDOPOTENTAIL
C=========================================================================
      IF(IVPTYP.EQ.0 .OR. NLPOT.EQ.1) GOTO 70
      DO 50 IB=1,NBANDS
      CALL VNLWAV(NRPLWV,NIONS,NSPEC,NIONSP,NPLWKP,DNLKG(1,0,NKP),VOLC, 
     &            PSCALE, CPTWFP,CWRK20,CWRK21,CWRK22,CWRK23,
     &            CPHSGR, VGNL, CELFRC, CPTWFL,CGRA,IB,IVPTYN)
      DO 55 JB=IB,NBANDS
      NNINDW=(JB-1)*NRPLWV
      CVNL=(0.0,0.0)
      DO 60 M=1,NPLWKP
  60  CVNL=CVNL+ CONJG(CPTWFP(M+NNINDW))*CELFRC(M)
      HR(IB,JB) = HR(IB,JB) + REAL(CVNL)
      HI(IB,JB) = HI(IB,JB) + AIMAG(CVNL)
  55  CONTINUE
  50  CONTINUE
C=======================================================================
C 4. NOW OBTAIN THE LOWER PART OF H(IB,IB)
C=======================================================================
  70  CONTINUE
      DO 71 IB=1,NBANDS
      DO 71 JB=1,IB
      HR(IB,JB)= HR(JB,IB)
      HI(IB,JB)=-HI(JB,IB)
  71  CONTINUE
C=======================================================================
C WRITE OUT THE HAMILTONIAN
C=======================================================================
      IF(IPRINT.LT.2) GOTO 74
      WRITE(*,*)' '
      WRITE(*,*)' ------- THE HAMILTONIAN BEFORE DIAGONALISATION ------'
      WRITE(*,*)'  THE UPPER TRIANGLE IS THE  R E A L  PART'
      WRITE(*,*)'  THE LOWER TRIANGLE IS THE IMAGINARY PART'
      DO 72 I=1,NBANDS
      IF(I.EQ.1) THEN
      WRITE(*,73)                    (HR(I,J), J=I,NBANDS)
      GOTO 72
      END IF
      WRITE(*,73) (HI(I,J),J=1,I-1), (HR(I,J), J=I,NBANDS)
  72  CONTINUE
  73  FORMAT(1X,10F11.5)
  74  CONTINUE
C=======================================================================
C DIAGONALISE H(I,J)     
C THE WAVE FUNCTIONS ARE ROTATED, SO THEY ARE EIGEN STATES
C=======================================================================
      CALL ROT(NBANDS,NRPLWV,CPTWFP,NINDPW,CWORK, NPLWV,NPLWKP,NGPTAR,
     &         MPLWV,HR,HI,AUX, FV1,FV2,FV3,NKP,NKPTS, CH0)
C=======================================================================
C UPDATE CELEN = THE KINETIC ENERGY OF EACH BAND
C=======================================================================
      DO 80 NN=1,NBANDS
      CSUM=(0.0,0.0)
      NINDW=NRPLWV*(NN-1)
      DO 85 M=1,NPLWKP
      CSUM=CSUM
     &    +DATAKE(1+7*(M-1))*CPTWFP(M+NINDW)*CONJG(CPTWFP(M+NINDW))
  85  CONTINUE
      CELEN(NN+(NKP-1)*NBANDS)=CSUM     
      IF(IPRINT.GE.2) WRITE(*,*)' AFTER ROT, CELEN=', CSUM
  80  CONTINUE
C=======================================================================
C UPDATE VNL
C=======================================================================
      IF(IVPTYP.EQ.0) GOTO 99
      DO 90 NB=1,NBANDS
      NINDW=(NB-1)*NRPLWV
      IF (NLPOT.EQ.0) THEN
      CALL VNLWAV(NRPLWV,NIONS,NSPEC,NIONSP,NPLWKP,DNLKG(1,0,NKP),
     &     VOLC, PSCALE, CPTWFP,CWRK20,CWRK21,CWRK22,CWRK23,
     &     CPHSGR, VGNL, CELFRC, CPTWFL,CGRA,NB,IVPTYN)
      CVNL=(0.0,0.0)
      DO 95 M=1,NPLWKP
  95  CVNL=CVNL+ CONJG(CPTWFP(M+NINDW))*CELFRC(M)
      VNL(NB,NKP)=REAL(CVNL)
      ELSE
      DO 91 M=1,MPLWV
      CPTOWR(M)=(0.0,0.0)                                             
      CPTNWR(M)=(0.0,0.0)
      CWORK(M)=(0.0,0.0)                                              
  91  CONTINUE
      DO 92  M=1,NPLWKP                                                 
      CPTOWR(NINDPW(M))=CPTWFP(M+NINDW)
  92  CONTINUE                                                          
C=======================================================================
C TRANSFORM THE WAVEFUNCTION INTO REAL SPACE
C=======================================================================
      CALL FFT3D(CPTOWR,CWORK,NGPTAR,1)                                 
C=======================================================================
C CALL ROUTINE TO CALCULATE NON-LOCAL ENERGY IN REAL SPACE
C=======================================================================
      CALL ENRLNL(VOLC,NGX,NGY,NGZ,VNL,NB,NKP,MPLWV,NRPLWV,             
     &   NRGRPT,NIONSP,NRLPPI,NBANDS,NKPTS,NSPEC,CWORK,CPTOWR,
     &   CPHGRD,NRLNL,NIONST,MXRLNL,IRLNL,PRLSCA,VRLGRD,NADGRD,MXRLSH,  
     &   CESAVE)                      
      END IF
  90  CONTINUE
C
  99  CONTINUE
      IF(IPRINT.GE.1)WRITE(*,*)' RETURN FROM SUBROT (SUBSPACE ROTATION)'
      RETURN
      END
C
C
C=======================================================================
      SUBROUTINE ROT(NBANDS,NRPLWV,CPTWFP,NINDPW,CWORK, NPLWV,NPLWKP,
     &  NGPTAR,MPLWV,HR,HI,AUX,FV1,FV2,FV3,NKP,NKPTS, CH0)
C=======================================================================
C     This routine performs sub-space diagonalisation of the occupied
C     states of the hamiltonian, by unitary matrix.
C
C     27-Feb-90 X. Weng (based on Qteish's version)
C     28-Feb-90 eigenvectors HR and HI normalised.
C     28-Feb-90 (HR,-HI)*CPTWFP part corrected.
C     28-Feb-90 Tested on Si 8-atom unit cell. 
C     04-MAR-90 Further tested
C
      IMPLICIT COMPLEX (C)
      DIMENSION CPTWFP(NRPLWV*NBANDS)
      DIMENSION NGPTAR(3)
      DIMENSION NINDPW(NRPLWV,NKPTS)
      DIMENSION HR(NBANDS,NBANDS),HI(NBANDS,NBANDS),AUX(NBANDS),
     &          FV1(NBANDS),FV2(NBANDS),FV3(NBANDS)
      DIMENSION CWORK(NRPLWV)
      DIMENSION CH0(NBANDS,NBANDS)
C=======================================================================
C     ITEST=1, FOR TEST 
C          =0, REGULAR USE
C=======================================================================
      ITEST=0
      IF(ITEST.EQ.1) THEN
        WRITE(*,*)'>>> THE INPUT H MATRIX'
        DO 5 J=1,NBANDS
        DO 6 I=1,NBANDS
        CH0(I,J)=CMPLX(HR(I,J), HI(I,J))
6       CONTINUE
        WRITE(*,102) (CH0(I,J), I=1,NBANDS)
5       CONTINUE
102     FORMAT(1X, 10F8.2 )
      END IF
C=======================================================================
C CHNAG SOLVES THE EIGENVALUE PROBLEM USING NAG ROUTINES.
C=======================================================================
      CALL CHNAG(NBANDS,NBANDS,HR,HI,AUX,FV1,FV2,FV3,IERR)
      IF(ITEST.EQ.1)WRITE(*,*)'IN ROT: EIGENVALUES FROM DIAGONALISATION'
      IF(ITEST.EQ.1)WRITE(*,100) (AUX(I), I=1, NBANDS)
100   FORMAT(2X,5F15.6)  
      AUXSUM=0.0
      DO 10 J=1, NBANDS
      AUXSUM=AUXSUM+AUX(J)
10    CONTINUE
      IF(ITEST.EQ.1) WRITE(*,*)'   SUM OF EIGENVALUES =:', AUXSUM
C=======================================================================
C     FIND NORM OF BANDS, and the Unitary matrix CUNIT
C
C                         /  VECTOR 1........\
C                         |  ................|
C       CUNIT =(HR, HI) = |  VECTOR J .......|      CUNIT(DAGGAR)*CUNIT=I
C                         |  ................|
C                         \  VECTOR NBANDS ../
C
C       THE EIGENVECTORS SHOULD BE NORMALISED BY NAG ROUTINE
C
C=======================================================================
      IF(ITEST.NE.1) GOTO 200
        DO 40 J = 1,NBANDS
        ANORM=0.0
        DO 50 I = 1,NBANDS
        ANORM = ANORM + HR(I, J)*HR(I, J) + HI(I, J)*HI(I, J)         
50      CONTINUE
        WRITE(*,150) J, AUX(J), ANORM
150     FORMAT(/,' BAND=', I2,'   EIGENVALUE=',F13.6,'  NORM=',F13.6,/)
        WRITE(*,102) (HR(I,J), HI(I,J), I=1,NBANDS)
40      CONTINUE 
C=======================================================================
C     BELOW IS FOR TESTING THE SUBROUTINE
C
C     TEST ON THE SIMILARITY TRANSFORM  U (DAGGAR) H0 U = H1
C=======================================================================
      WRITE(*,*)' BELOW IS THE DIAGNOLISED H1 = U(DAGGAR) H U'
      DO 95 J=1,NBANDS
C
C     FOR JTH ROW OF H1: CWORK(I), I=1,NBANDS
C
      DO 96 I=1,NBANDS
      CWORK(I)=(0.0,0.0)
      DO 97 K=1,NBANDS
      DO 98 L=1,NBANDS
      CWORK(I) = CWORK(I)
     1 + CMPLX(HR(K,I),-HI(K,I)) * CH0(K,L) *CMPLX(HR(L,J), HI(L,J))
98    CONTINUE
97    CONTINUE
96    CONTINUE
      WRITE(*,102) (CWORK(I), I=1,NBANDS)
95    CONTINUE
C
 200  CONTINUE 
C======================================================================
C
C     FOR ITH VECTOR, JTH COMPONENT
C
C                    ______
C                    \        +                        I=1,NPLWKP
C     NEW PSI(J,I) =  >      U (J,K) *OLD PSI(K,I)     J=1,NBANDS
C                    /_____
C                   
C                  K=1,NBANDS
C  
C     WE FIND NEW PSI(I,J), J FIRST 
C
C======================================================================
      DO 60 I=1,NPLWKP
      DO 70 J=1,NBANDS
      CWORK(J)=CMPLX(0.0,0.0)
      DO 70 K=1,NBANDS
      IOFK=NRPLWV*(K-1)+I
      CWORK(J) = CWORK(J) + CMPLX(HR(K,J),-HI(K,J)) * CPTWFP(IOFK)
70    CONTINUE  
      DO 80 K=1,NBANDS
      IOFK=NRPLWV*(K-1)+I
      CPTWFP(IOFK) = CWORK(K)
80    CONTINUE     
60    CONTINUE
C
C     THE NORM OF WAVEFUNCTION BEFORE AND AFTER, SHOULD BE THE SAME
C
      IF(ITEST.EQ.1) THEN
         DO 310 J=1,NBANDS
         CNORM=CMPLX(0.0,0.0)
         NINJ=NRPLWV*(J-1)
         DO 312 I=1,NPLWKP
312      CNORM=CNORM+CPTWFP(NINJ+I)*CONJG(CPTWFP(NINJ+I))
310      WRITE(*,*)' NORM FOR NEW PSI', J, CNORM
         END IF
C     
      RETURN
      END
C=====================================================================
      SUBROUTINE CHNAG(NM,N,A,B,E,FV1,FV2,FV3,IERR)  
C
      DIMENSION A(NM,NM),B(NM,NM),FV1(NM),FV2(NM),FV3(NM),E(NM)
      EXTERNAL X02AJE, F01BCE, F02AYE
C
C     F01BCE: reduction by similarity transfrmaions
C             complex Hermitian matrix to real tridigonal form
C
C     F02AYE: complex Hermitian matrix, all eigenvalues and vectors
C             tridigonal form by F01BCE, QL algorithm
C
C     The eigenvector corresponding to the eigenvalue E(J) is 
C     CMPLX( HR(I,J), HI(I,J)), J=1,N
C
C     TO BE LINKED WITH NAGSP/LIB
C
C=====================================================================
C     REDUCE TO TRIDIAGONAL FORM
C     TOL NO LONGER NEEDED
C=====================================================================
      TOL = 0.0
      CALL F01BCE(N,TOL,A,NM,B,NM,E,FV1,FV2,FV3)
C=====================================================================
C     EIGENVALUES AND EIGENVECTORS OF ORIGINAL MATRIX
C     EPS=MACHINE PRECISION
C=====================================================================
      EPS = X02AJE()
      CALL F02AYE(N,EPS,E,FV1,A,NM,B,NM,IERR)
      RETURN
      END
