**********************************************************************
C IN VALUES
C
C CV(r)    ... THE REAL SPACE TOTAL KOHN-SHAM POTENTIAL
C CVD()   .... EMPTY
C CDUM    .... EMPTY
C CDIR    .... EMPTY
C CPTWFP  .... WAVE FUNCTIONS FOR CURRENT K-POINT
C CPTWFL  .... EMPTY
C CPTOWR,CPTNWR,CWORK,CGRA,CORGR,COGRPI,CDIRPI,PRECON : EMPTY
C
C CELEN(1)   . THE KINETIC ENERGY FOR 1ST BAND OF K-POINT 1
C EIGEN() ... EMPTY
C WTKPT(1) ... WEIGHT OF K-POINT NO 1
C OCC(N,1) ... WEIGHT OF Nth BAND FOR K-POINT NO 1
C
C OUT VALUES
C
C EIGEN(1)... EIGENVALUE FOR 1ST BAND OF K-POINT NKP
C=======================================================================
      SUBROUTINE CONGRABS(NBANDS,NKPTS,NPLWV,MPLWV,         
     &    NRPLWV,NINDPW,NPLWKP,WTKPT,CV,CPTWFP,             
     &    CPTWFL,DIRC,RECC,VOLC,CELEN,NGPTAR,DATAKE,        
     &    CPTOWR,CPTNWR,CWORK,CGRA,CORGR,CDIR,COGRPI,CDIRPI,
     &    PRECON,CVD,NPKPT,CDUM,
     &    HR, HI, AUX, FV1, FV2, FV3, CH0, NITMAX,          
     &    NSPEC, NIONS, NIONSP,PSCALE,                      
     &    DNLKG, CPHSGR,VGNL, CELFRC, CWRK20,               
     &    CWRK21,CWRK22,CWRK23,VNL,IVPTYP,IPRINT,ISBROT,    
     &    OCC,EIGVAL,NLPOT,ICLOCK,IVPTYN,NITER,N,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 CVD(*)
      DIMENSION WTKPT(*)
      DIMENSION NIONSP(*)
      DIMENSION OCC(NBANDS,NKPTS), EIGVAL(NBANDS,NKPTS)
      DIMENSION DIRC(3,3),RECC(3,3)
C======================================================================
C ARRAYS NEEDED FOR REAL SPACE PROJECTION OF THE NON-LOCAL POTENTIALS
C======================================================================
      DIMENSION  NRLNL(NSPEC),PRLSCA(MXRLNL,NSPEC)                 
      DIMENSION  IRLNL(MXRLNL,NSPEC)
      DIMENSION  VRLGRD(NRGRPT,MXRLSH,NIONST),CPHGRD(NRGRPT,NIONST)
      DIMENSION  NADGRD(NRGRPT,NIONST),NRLPPI(NIONST)              
      DIMENSION CESAVE(NIONST,20) 
C
C=======================================================================
      DIMENSION CPTOWR(*)
      DIMENSION CPTNWR(*)
      DIMENSION CWORK(*)
      DIMENSION PRECON(*)
      DIMENSION CGRA(*),CORGR(*),CDIR(*),COGRPI(*),CDIRPI(*)
      DIMENSION CDUM (*)
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),VGNL(NRPLWV,0:2,NSPEC)
      DIMENSION CELFRC(NRPLWV)
      DIMENSION CWRK20(NIONS),CWRK21(3,NIONS)
      DIMENSION CWRK22(NIONS),CWRK23(3,3,NIONS)
      DIMENSION VNL(NBANDS, NKPTS)
C=======================================================================
C     HAMMER'S
C=======================================================================
C
      RINPLW=1.0/NPLWV
C=======================================================================
C
C START CONJUGATE GRADIENTS LOOP
C
C=======================================================================
      DO 2001 NN = 1 , NBANDS
        NINDW = NRPLWV * ( NN - 1 )
        IF (NN.LE.1) GO TO 7026
C=======================================================================
C    ORTHOGONALISE PRESENT WAVEFUNCTION TO LOWER BAND VECTORS
C=======================================================================
        DO 7022 NB = 1 , NN - 1
          NINDD = NRPLWV * ( NB - 1 )
          COVERL = (0.0,0.0)
          DO 7023 M = 1 , NPLWKP
            COVERL = COVERL + CPTWFP(M+NINDW) * CONJG (CPTWFP(M+NINDD))
 7023     CONTINUE
          COVERL =  - COVERL
          DO 7024 M = 1 , NPLWKP
            CPTWFP(M+NINDW) = CPTWFP(M+NINDW) + COVERL * CPTWFP(M+NINDD)
 7024     CONTINUE
 7022   CONTINUE
C=======================================================================
C  NORMALIZE PRESENT WAVEFUNCTION
C=======================================================================
        ANORM = 0.0
        DO 7027 M = 1 , NPLWKP
          ANORM = ANORM +  CPTWFP(M+NINDW) * CONJG (CPTWFP(M+NINDW))
 7027   CONTINUE
        ANORM = 1.0 / SQRT(ANORM)
        DO 7028 M = 1 , NPLWKP
          CPTWFP(M+NINDW) = ANORM * CPTWFP(M+NINDW)
 7028   CONTINUE
C
 7026   CONTINUE
C=======================================================================
C  START ITERATIONS FOR THE PRESENT BAND.
C
C  WE ARE LOOKING FOR THE WAVEFUNCTION WHICH WILL MINIMIZE
C  CORRESPONDING EIGENVALUE UNDER CONSTRAINT OF BEING ORTHOGONAL
C  TO ALL LOWER BANDS.
C  START WITH |Fo> AND ADD TRIAL FUNCTION |W>,
C  THEN MINIMIZE:
C
C      <Fo+W|H|Fo+W>
C  E = -------------
C       <Fo+W|Fo+W>
C=======================================================================
C
C       WRITE(*,*)' BANDS CYCLE EIGEN EIGENNEW XL1 XL2 R1 R2 XLAMDA'
        DO 5000 NIT = 1 , NITMAX
C ......................................................................
C    CALCULATE THE PRODUCT OF HAMILTONIAN AND WAVEFUNCTION
C    H|Fo> WHICH WILL GIVE ALSO AFTER PRECONDITIONING AND 
C    CONJUGATION THE SEARCH DIRECTION, |W>.
C ......................................................................
C
C=======================================================================
C    1 ) CALCULATE THE PRODUCT OF THE POTENTIAL AND THE WAVEFUNCTION IN
C         real SPACE AND FOURIER TRANSFORM TO reciprocal SPACE
C=======================================================================
C INITIALISE THE ARRAYS USED IN THE FOURIER TRANSFORM TO ZERO. THIS STEP
C MUST BE PERFORMED TO ENSURE THAT THE COEFFICIENTS OF THE PLANE WAVES
C BEYOND THE CUT-OFF ENERGY ARE ZERO.
C=======================================================================
          DO 5561 M = 1 , MPLWV
            CPTOWR(M) = (0.0,0.0)
            CPTNWR(M)=(0.0,0.0)
            CWORK(M) = (0.0,0.0)
 5561     CONTINUE
          DO 5510 M = 1 , NPLWKP
            CPTOWR(NINDPW(M)) = CPTWFP(M+NINDW)
 5510     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                                                         
          DO 2002 NNN = 1 , NPLWV
            CPTNWR(NNN)=(CPTNWR(NNN)+CV(NNN)*CPTOWR(NNN))*RINPLW
 2002     CONTINUE
C=======================================================================
C TRANSFORM (WAVEFUNCTION*POTENTIAL) INTO RECIPROCAL SPACE
C=======================================================================
          CALL FFT3D(CPTNWR,CWORK,NGPTAR,-1)
C=======================================================================
C    2) ADDED THE KINETIC ENERGY TERM, WE OBTAIN THE GRADIENT
C=======================================================================
          DO 2010 M=1,NPLWKP
            CGRA(M) = CPTNWR(NINDPW(M)) + 
     &                CPTWFP(M+NINDW) * DATAKE(1+(7*(M-1)))
 2010     CONTINUE
C=======================================================================
C    IF THE PSEUDOPOTENTIAL IS NON-LOCAL, WE ADD THE NON-LOCAL
C    CONTRIBUTION. (MUST BE IN THE KLEIMAN-BYLANDER FORM)
C    THIS IS DONE IN SUBROUTINE VNLWAV (Vnl*WAV)     16-MAR-90 X.WENG
C=======================================================================
          IF (IVPTYP.NE.0 .AND.NLPOT.EQ.0) THEN
            IF (ICLOCK.EQ.1) CALL PCLOCK(203)
            CALL VNLWAV(NRPLWV,NIONS,NSPEC,NIONSP,NPLWKP,
     &                  DNLKG(1,0,NPKPT),VOLC,PSCALE,CPTWFP,
     &                  CWRK20,CWRK21,CWRK22,CWRK23,CPHSGR,VGNL,CELFRC,
     &                  CDUM,CORGR,NN,IVPTYN)
            DO 2015 M = 1 , NPLWKP
              CGRA(M) = CGRA(M) +  CELFRC(M)
 2015       CONTINUE
            IF (ICLOCK.EQ.1) CALL PCLOCK(204)
          END IF
C
          DO 2017 M = 1 , NPLWKP
            CGRA(M)  = - CGRA(M)
            CORGR(M) = CGRA(M)
 2017     CONTINUE
C=======================================================================
C    3) ORTHOGONALISE GRADIENT TO PRESENT BAND VECTOR
C       changed ->>         ** LOWER AND PRESENT BANDS ONLY **
C=======================================================================
          DO 2022 NB = 1 , NN
            NINDD = NRPLWV * ( NB - 1 )
            COVERL = (0.0,0.0)
            DO 2023 M = 1 , NPLWKP
              COVERL = COVERL + CORGR(M) * CONJG (CPTWFP(M+NINDD))
 2023       CONTINUE
            COVERL = - COVERL
            DO 2024 M = 1 , NPLWKP
              CORGR(M) = CORGR(M) + COVERL * CPTWFP(M+NINDD)
 2024       CONTINUE
 2022     CONTINUE
C=======================================================================
C    4) PRE-CONDITION THE GRADIENT
C=======================================================================
          NCELIN = NN + NBANDS * ( NPKPT - 1 )
          ENBAKE = REAL (CELEN(NCELIN))
          IF (ENBAKE.LT.1.0) THEN
            ENKEIN = 1.0
          ELSE
            ENKEIN = 1.0 / ENBAKE
          ENDIF
          DO 2031 M = 1 , NPLWKP
            X = DATAKE(1+7*(M-1)) * ENKEIN
            PCNUM = 27.0 + ( 18.0 + ( 12.0 + 8.0 * X ) * X ) * X
            PCDEN = PCNUM + 16.0 * ( X * X )**2
            PRECON(M) = PCNUM / PCDEN
C            PRECON(M) = (1.0,0.0)
            CDUM(M) = CORGR(M) * PRECON(M)
 2031     CONTINUE
C=======================================================================
C    5) ORTHOGONALISE PRECONDITIONED RESIDUAL TO LOWER AND PRESENT BANDS
C=======================================================================
          DO 2032 NB = 1 , NN
            NINDD = NRPLWV * ( NB - 1 )
C
C     COVERL=CDOTC(NPLWKP,CPTWFP(NINDD+1),1,CORGR,1)
C     COVERL=-COVERL
C     CALL CAXPY(NPLWKP,COVERL,CPTWFP(NINDD+1),1,CORGR,1)
C
            COVERL = (0.0,0.0)
            DO 2033 M = 1 , NPLWKP
              COVERL = COVERL + CDUM(M) * CONJG (CPTWFP(NINDD+M))
 2033       CONTINUE
            DO 2034 M = 1 , NPLWKP
              CDUM(M) = CDUM(M) - COVERL * CPTWFP(NINDD+M)
 2034       CONTINUE
 2032     CONTINUE
C=======================================================================
C AFTER FIRST ITERATION USE CONJUGATE DIRECTIONS
C=======================================================================
          IF (NIT.EQ.1) THEN
            DO 2035 M = 1 , NPLWKP
C              CDIR(M)=CORGR(M)
              CDIR(M) = CDUM(M)
 2035       CONTINUE
          ELSE
            CDOT1 = (0.0,0.0)
            CDOT2 = (0.0,0.0)
            DO 2040 M = 1 , NPLWKP
C
C     From Numerical Recipes p304, eq(10.6.5):
C     1) Fletcher-Reeves method  GAMMA = G(i+1)*G(i+1)/G(i)*G(i)
C     2) Poloak-Ribiere method   GAMMA = [G(i+1)-G(i)]*G(i+1)/G(i)*G(i)
C
              CDOT1 = CDOT1 + CORGR(M) * CONJG (CDUM(M))
              CDOT2 = CDOT2 + COGRPI(M) * CONJG (COGRPI(M)) * PRECON(M)
C
C      CDOT1=CDOT1+(CORGR(M)-COGRPI(M)) *CONJG(CDUM(M))
C      CDOT2=CDOT2+COGRPI(M)*CONJG(COGRPI(M))
 2040       CONTINUE
            CGAMMA = CDOT1 / CDOT2
            DO 2050 M = 1 , NPLWKP
              CDIR(M) = CDUM(M) + CGAMMA * CDIRPI(M)
 2050       CONTINUE
          ENDIF
          DO 2060 M = 1 , NPLWKP
            CDIRPI(M) = CDIR(M)
            COGRPI(M) = CORGR(M)
 2060     CONTINUE
C=======================================================================
C FINALLY ORTHOGONALISE SEARCH DIRECTION TO PRESENT BAND 
C=======================================================================
          DO 2071 NB = 1 , NN
            NINDD = NRPLWV * ( NB - 1 )
          COVERL = (0.0,0.0)
          DO 2070 M = 1 , NPLWKP
            COVERL = COVERL + CDIR(M) * CONJG (CPTWFP(M+NINDD))
 2070     CONTINUE
          DO 2080 M = 1 , NPLWKP
            CDIR(M) = CDIR(M) - COVERL * CPTWFP(M+NINDD)
 2080     CONTINUE
 2071     CONTINUE
C=======================================================================
C CALCULATE THE PRODUCT OF THE HAMILTONIAN MATRIX WITH THE NEW DIRECTION
C IN REAL SPACE, H|W>
C=======================================================================
          DO 3000 M = 1 , NPLWV
            CPTOWR(M) = (0.0,0.0) 
            CVD(M) = (0.0,0.0) 
 3000     CONTINUE
          DO 3010 M = 1 , NPLWKP
            CPTOWR(NINDPW(M)) = CDIR(M)
 3010     CONTINUE
C=======================================================================
C TRANSFORM THE DIRECTION 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)                                                      
            DO 3020 NNN = 1 , NPLWV
              CPTOWR(NNN) = ( CPTNWR(NNN) + CV(NNN) * CPTOWR(NNN) )
     &                      * RINPLW
 3020       CONTINUE
          ELSE
            DO 3021 NNN = 1 , NPLWV                                
              CPTOWR(NNN) = CV(NNN) * CPTOWR(NNN) * RINPLW
 3021       CONTINUE                                               
          ENDIF                                                         
C=======================================================================
C TRANSFORM (WAVEFUNCTION*POTENTIAL) INTO RECIPROCAL SPACE
C=======================================================================
          CALL FFT3D(CPTOWR,CWORK,NGPTAR,-1)
C=======================================================================
C    2) ADDED THE KINETIC ENERGY TERM, WE OBTAIN THE GRADIENT
C=======================================================================
          DO 3030 M = 1 , NPLWKP
            CVD(M) = CPTOWR(NINDPW(M)) + CDIR(M) * DATAKE(1+(7*(M-1)))
 3030     CONTINUE
C=======================================================================
C    IF THE PSEUDOPOTENTIAL IS NON-LOCAL, WE ADD THE NON-LOCAL
C    CONTRIBUTION. (MUST BE IN THE KLEIMAN-BYLANDER FORM)
C    THIS IS DONE IN SUBROUTINE VNLWAV (Vnl*WAV)     16-MAR-90 X.WENG
C=======================================================================
          IF (IVPTYP.NE.0 .AND. NLPOT.EQ.0) THEN
            IF (ICLOCK.EQ.1) CALL PCLOCK(203)
            CALL VNLWAV(NRPLWV,NIONS,NSPEC,NIONSP,NPLWKP,
     &                  DNLKG(1,0,NPKPT),VOLC,PSCALE,CDIR,
     &                  CWRK20,CWRK21,CWRK22,CWRK23,CPHSGR,VGNL,CELFRC,
     &                  CDUM,CORGR,1,IVPTYP)
            DO 3040 M = 1 , NPLWKP
              CVD(M) = CVD(M) + CELFRC(M)
 3040       CONTINUE
          END IF
C=======================================================================
C  CALCULATE A NEW STEP ALONG THE SEARCH DIRECTION:
C  |F> = |Fo> + Lambda |W>, WHERE LAMBDA COMES FROM THE
C  MINIMIZATION OF EIGENVALUE - SECOND ORDER EQUATION.
C
C   XA = <Fo|H|Fo>
C   XB = <W|H|Fo> + c.c.
C   XC = <W|H|W>
C   XD = <W|W>
C
C (ASSUME THAT <W|Fo>=0)
C
C                    XA + Lambda * XB + Lambda**2 * XC
C   eigenvalue = min ---------------------------------
C                          1  + Lambda**2 * XD
C=======================================================================
          XA = 0.0
          XB = 0.0
          XC = 0.0
          XD = 0.0
          DO 3060 M = 1 , NPLWKP
            XA = XA - CPTWFP(M+NINDW) * CONJG (CGRA(M))
            XB = XB + CPTWFP(M+NINDW) * CONJG (CVD(M)) 
     &              - CDIR(M) * CONJG (CGRA(M))
            XC = XC + CDIR(M) * CONJG (CVD(M))
            XD = XD +  CDIR(M) * CONJG (CDIR(M))
 3060     CONTINUE
          A1 = XB * XD
          A2 = XA * XD - XC 
          A3 = - XB
          IF (ABS(A1).GT.1.E-8) THEN
            R1 = - (  A2 + SQRT ( A2**2 - A3 * A1 ) ) / A1
            R2 = - (  A2 - SQRT ( A2**2 - A3 * A1 ) ) / A1
          ELSE
            R1 = - A3 / ( 2.0 * A2 )
            R2 = R1
          END IF
C
          XL1 = ( XA + XB * R1 + XC * R1**2 ) / 
     &         ( 1.0 + XD * R1**2 )
          XL2 = ( XA + XB * R2 + XC * R2**2 ) / 
     &         ( 1.0 + XD * R2**2 )
C
          AANEW = MIN(XL1,XL2)
          IF (XL1.LT.XL2) THEN
            XLAMDA = R1
          ELSE
            XLAMDA = R2
          END IF
C         WRITE(6,101) NN,NIT,XA,AANEW,A1,A2,A3,XLAMDA
101    FORMAT(1X,2I3,7F10.5)
C=======================================================================
C UPDATE THE WAVEFUNCTION AND NORMALIZE IT
C=======================================================================
          ANORM = 0.0
          DO 3070 M = 1 , NPLWKP
            CPTWFP(M+NINDW) = CPTWFP(M+NINDW) + XLAMDA * CDIR(M)
            ANORM = ANORM +  CPTWFP(M+NINDW) * CONJG (CPTWFP(M+NINDW))
 3070     CONTINUE
          ANORM = 1.0 / SQRT(ANORM)
          DO 3080 M = 1 , NPLWKP
            CPTWFP(M+NINDW) = ANORM * CPTWFP(M+NINDW)
 3080     CONTINUE
C======================================================================
C     MOVE ONTO THE NEXT ITERATION
C======================================================================
 5000   CONTINUE
        EIGVAL(NN,NPKPT) = AANEW      
C======================================================================
C     MOVE ONTO THE NEXT BAND
C======================================================================
 2001 CONTINUE
C=======================================================================
C SUB-SPACE ROTATION (FINDS OUT EIGEN STATES AND EIGEN VALUES)
C THE CELEN AND VNL ARE UPDATED
C=======================================================================
      DO 2800 NN = 1,NBANDS 
        NINDX = NRPLWV * ( NN - 1 )
        NNN = NN + NBANDS * ( NPKPT - 1 )
        CELEN(NNN) = (0.0,0.0)
        DO 2810 M = 1,NPLWKP
          CELEN(NNN) = CELEN(NNN) +
     &    DATAKE(1+7*(M-1)) * CPTWFP(M+NINDX) * CONJG(CPTWFP(M+NINDX))
 2810   CONTINUE
C        WRITE(*,*)'  KINETIC ENERGY IS: ', CELEN(NN)
 2800 CONTINUE
C
      IF (ISBROT.NE.0 .AND. N.GE.NITER) THEN
        CALL SUBROT(NBANDS,NKPTS,NPLWV,MPLWV, NRPLWV,NINDPW,NPLWKP,CV,
     &       CPTWFP,CPTWFL,VOLC, CELEN,VNL,NGPTAR,DATAKE,
     &       CPTOWR,CPTNWR,CWORK,CGRA,NPKPT,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)            
        DO 6006 NB = 1 , NBANDS
          EIGVAL(NB,NPKPT) = AUX(NB)
 6006   CONTINUE
      END IF
C
      IF (IPRINT.GE.1) WRITE(*,*)' BYE FROM CONGRA-BS'
      RETURN
      END
