C**************************************************************************
C     CONJUGATE GRADIENT METHOD FOR TOTAL ENERGY MINIMISATION 
C     WITH REGARD TO ELECTRONIC BAND COEFFICIENTS
C
C     M.P. TETER, M.C. PAYNE, AND D.C. ALLAN, PHYS. REV. B40, 12255 (1989)
C
C     VERSION OF 27-FEB-1990
C    
C     MODIFICATIONS MADE BY X WENG:    
C
C     28-Feb-90  Tested on Si 8 atoms unit cell. Seems work.
C     05-Mar-90  Correction made as suggested by MC Payne (/PRECON)
C                ORTHOGONALISE GRADIENT TO *ALL* BAND VECTORS 
C     16-Mar-90  Non-Local pseudoptenial used
C     31-MAR-90  Tested on Si and MgO
C     14-APR-90  VNL replaces AUX ( for NKPTS.GT.1)
C     11-JUN-90  PRECON PART IS CHANGED,
C                was: ENKE.LE.10   etc.
C                now: ENBAKE.LE.1  etc.
C     11-JUN-90  was: GAMMA = CORGR*CONJG(CORGR)/ COGRPI*CONJG(COGRPI)  
C                now: GAMMA = CORGR*CONJG(CDUM(M)) / COGRPI*CONJG(COGRPI)     
C
C     14-JUN-90  IVPTYP for local/nonlocal Vps added  
C                IPRINT for print-out added  
C
C     16-JUN-90  Subspace rotation attached (ISBROT=0 skip it)
C
C     17-JUN-90  CV0, CV1, AND CV2 ARE REDUCED TO ONE K-POINT A TIME
C                TO REDUCE THE DEMAND ON MEMORY SPACE. THIS WILL NOT 
C                RESULT AN INCREASE IN COMPTUTIONAL TIME FOR DYNAMIC 
C                SIMULATION, BUT WILL DO SO FOR THE STATIC CALCULATION,
C                WHERE CV0 ETC REMAIN THE SAME.
C
C     26-JUN-90  RESTRUCTURED THE NON-LOCAL PART. THE USAGE OF CV0,CV1 AND
C                CV2 IS ABANDONED. WE NOW USE
C 	         1) CPHSGR(NRPLWV,NIONS,NSPEC) = EXP(iq.Rn)
C	         2) VGNL(NRPLWV,0:2,NSPEC) = Vnl(l,q)
C	         IN ORDER TO REDUCE THE REQUIRED SPACE. TESTED OK
C
C
C     18-JUL-90  TIME STATEMENTS ADDED. TO PRINT OUT TIME, SET ICLOCK=1
C
C
C***********************************************************************
C IN VALUES
C
C CV(r)    ... THE REAL SPACE TOTAL KOHN-SHAM POTENTIAL
C CHDENR(r) .. THE REAL SPACE CHARGE DENSITY
C CVION(r) ... THE REAL SPACE IONIC POTENTIAL
C DIRDAT(G) .. THE RECIPROCAL SPACE GRID OF 1/G**2 USED FOR HARTREE POT 
C CHDENG() ... EMPTY
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 DENC   .... THE HARTREE ENERGY CORRECTION
C XCENC  ..... THE EXCHANGE CORRELATION ENERGY CORRECTION
C PSCENC  ..... THE CORE ENERGY
C TEWEN   .... THE EWALD ENERGY
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 CONGRA(NGX,NGY,NGZ,NBANDS,NKPTS,NPLWV,MPLWV,
     &    NRPLWV,ENMAX,NINDPW,NPLWKP,WTKPT,CV, CPTWFP,
     &    CPTWFL,DIRC,RECC,VOLC,CELEN,NGPTAR,SUMWEI,DATAKE,
     &    CPTOWR,CPTNWR,CWORK,
     &    NBANOC,CGRA,CORGR,CDIR,COGRPI,CDIRPI,
     &    PRECON,CVION,CVD,CHDENR,DENC,
     &    XCENC,TEWEN,PSCENC,NPKPT,CHDENG,XCFDAT,XCPDAT,DIRDAT,CDUM,
     &    HR, HI, AUX, FV1, FV2, FV3, CH0, NITMAX, TOTEN1,
     &    NSPEC, NIONS, NIONSP,PSCALE, 
     &    DNLKG, CPHSGR,VGNL, CELFRC, CWRK20,
     &    CWRK21,CWRK22,CWRK23,VNL,IVPTYP,IPRINT,ISBROT,
     &    IOCCUP,OCC,EIGVAL, ISYMM,NLPOT,
     &    SIGKE, SIGXC,EXCDAT, ICLOCK,IVPTYN,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 CVION(*)
      DIMENSION CHDENR(*)
      DIMENSION CVD(*)
      DIMENSION CHDENG(*)
      DIMENSION WTKPT(*)
      DIMENSION DIRDAT(*)
      DIMENSION NIONSP(*)
      DIMENSION OCC(NBANDS,NKPTS), EIGVAL(NBANDS,NKPTS)
      DIMENSION EXCDAT(*), XCFDAT(*), XCPDAT(*)
      DIMENSION DIRC(3,3),RECC(3,3) 
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
C SIGKE(6,NKPTS) = THE FORCE ON THE UNIT CELL DUE TO THE CHANGE IN THE
C          KINETIC ENERGY ON CHANGING THE SIZE OF THE CELL FOR EACH
C          K POINT
C
C=======================================================================
      DIMENSION CPTOWR(*)
      DIMENSION CPTNWR(*)
      DIMENSION CWORK(*)
      DIMENSION PRECON(*)
      DIMENSION CGRA(*),CORGR(*),CDIR(*),COGRPI(*),CDIRPI(*)
      DIMENSION CDUM (*)
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 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 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)
      DIMENSION IVPTYN(NSPEC)
C=======================================================================
C     HAMMER'S
C=======================================================================
      DIMENSION TAU(6), SIGKE(6), SIGXC(6)
C
      DATA ALPHA0/0.05/
C
      WTKPTT = WTKPT(NPKPT)
      IF (ISYMM.EQ.1) THEN
        WTKPT(NPKPT) = SUMWEI
        IF (IPRINT.GE.1) WRITE (*,*) 'SYMM:',WTKPTT, 'NATURAL:',SUMWEI
      END IF
C
      RINPLW = 1.0 / NPLWV
C======================================================================
C CALCULATE THE TOTAL ENERGY
C=========================================================================
      ENPOT=0.0
      DO 1000 M=1,NPLWV
        ENPOT=ENPOT+REAL(CHDENR(M)*CV(M))
 1000 CONTINUE
      ENPOT=ENPOT*RINPLW
C=========================================================================
C     THE CONTRIBUTION TO TOTAL ENERGY DUE TO THE NON-LOCAL PARTS
C     OF PSEUDOPOTENTAIL, STORED IN VNL(NBANDS,NKPTS)
C     CPTWFP IS IN THE RECIPROCAL SPACE NOW.
C
C     ENVNL: IS TOTAL ENERGY OF DUE TO NON-LOCAL POTENTIAL, SUM OVER ALL  
C     BANDS.
C=========================================================================
      IF (ICLOCK.EQ.1) CALL PCLOCK(201)
      IF (IVPTYP.EQ.0 .OR. NLPOT.EQ.1) GOTO 1003
      DO 1001 NB=1,NBANDS
        CALL VNLWAV(NRPLWV,NIONS,NSPEC,NIONSP,NPLWKP,DNLKG(1,0,NPKPT),
     &              VOLC,PSCALE,CPTWFP,CWRK20,CWRK21,CWRK22,
     &              CWRK23,CPHSGR,VGNL,CELFRC,CDUM,CORGR,NB,
     &              IVPTYN)
        VNL(NB,NPKPT)=0.0
        NINDW=(NB-1)*NRPLWV
        DO 1002 M=1,NPLWKP
          VNL(NB,NPKPT)=VNL(NB,NPKPT)+
     &                  REAL(CONJG(CPTWFP(M+NINDW))*CELFRC(M))
1002    CONTINUE
1001  CONTINUE 
      IF (IPRINT.GE.3) WRITE(*,4000) (VNL(NB,NPKPT),NB=1,NBANDS)
      IF(ICLOCK.EQ.1) CALL PCLOCK(202)
C
1003  CONTINUE
      ENKE=0.0
      ENVNL=0.0
      DO 1010 NK=1,NKPTS
        IF (NK.EQ.NPKPT) GO TO 1010
        DO 1020 NB=1,NBANDS
          ENKE =ENKE +(2.0*WTKPT(NK)*REAL(CELEN((NK-1)*NBANDS+NB)))*
     &          OCC(NB,NK)
          IF(IVPTYP.EQ.1)ENVNL=ENVNL+2.0*WTKPT(NK)*VNL(NB,NK)*OCC(NB,NK)
 1020   CONTINUE
 1010 CONTINUE
C
      NK = NPKPT
      ENKE1 = ENKE
      ENVNL1 = ENVNL
      DO 1021 NB=1,NBANDS
        ENKE1 =ENKE1 +(2.0*WTKPTT*REAL(CELEN((NK-1)*NBANDS+NB)))*
     &        OCC(NB,NK)
        ENKE =ENKE +(2.0*WTKPT(NK)*REAL(CELEN((NK-1)*NBANDS+NB)))*
     &        OCC(NB,NK)
        IF(IVPTYP.EQ.1) THEN
          ENVNL1=ENVNL1+2.0*WTKPTT*VNL(NB,NK)*OCC(NB,NK)
          ENVNL =ENVNL +2.0*WTKPT(NK)*VNL(NB,NK)*OCC(NB,NK)
        END IF
 1021 CONTINUE
C
      TOTEN=ENPOT+ENKE+DENC+XCENC+TEWEN+PSCENC + ENVNL
      TOTEN1 = TOTEN
      TOTEN2=ENPOT+ENKE1+DENC+XCENC+TEWEN+PSCENC + ENVNL1
      IF (IPRINT.GE.1) THEN
        WRITE(*,*)' Welcome to CONGRA: NITMAX =', NITMAX
        IF (IPRINT.GE.2) THEN
          WRITE(*,4001) ENKE1,ENPOT
          IF (IVPTYP.NE.0) WRITE(*,4002) ENVNL1
          WRITE(*,4003) DENC,XCENC,TEWEN,PSCENC
        END IF
        IF (ISYMM.NE.0) WRITE(*,1) TOTEN2
        WRITE(*,1) TOTEN
      END IF
 4000 FORMAT(' Vnl',6F12.6)
 4001 FORMAT('  ENKE  = ', F16.7, /, '  ENPOT = ', F16.7)
 4002 FORMAT('  ENVNL = ', F16.7)
 4003 FORMAT('  DENC  = ', F16.7,
     &   /,  '  XCENC = ', F16.7, /, '  TEWEN = ', F16.7,
     &   /,  '  PSCENC= ', F16.7)
C=======================================================================
C CALCULATE THE STRESS ON THE UNIT CELL
C
C                                                      2
C                  [   h**2]  --  !           i(k+G)r !
C  Sigma      = -2 [2*-----]  >   ! c(n,k+G) e        !  (k+G)  (k+G)
C       Ai,Aj    ! [   2 m ]  --  !                   !       Ai     Aj
C                !           k,G,n
C               spin
C=======================================================================
      DO 2009 NN=1,NBANOC
        NINDW=NRPLWV*(NN-1)
        NINDKE=0
        DO 2005 M=1,6
          TAU(M)=0.0
 2005   CONTINUE
        DO 2008 N=1,NPLWKP
          NINDKE=7*(N-1)
          WFMAG=REAL(CPTWFP(NINDW+N)*CONJG(CPTWFP(NINDW+N)))
          TAU(1)=TAU(1)+DATAKE(2+NINDKE)*WFMAG
          TAU(2)=TAU(2)+DATAKE(3+NINDKE)*WFMAG
          TAU(3)=TAU(3)+DATAKE(4+NINDKE)*WFMAG
          TAU(4)=TAU(4)+DATAKE(5+NINDKE)*WFMAG
          TAU(5)=TAU(5)+DATAKE(6+NINDKE)*WFMAG
          TAU(6)=TAU(6)+DATAKE(7+NINDKE)*WFMAG
 2008   CONTINUE
C=======================================================================
C SUM THE CONTRIBUTIONS TO THE FORCE ON THE UNIT CELL FROM EACH K POINT
C SCALING THE STRESSES BY 2/VOLC, MULTIPLYING BY THE K POINT WEIGHT
C AND BY OCC FOR ELECTRON SPINS AND OCCUPANCIES
C=======================================================================
        WGHT=WTKPTT*4.0/VOLC*OCC(NN,NPKPT) 
        SIGKE(1)=SIGKE(1)-TAU(1)*WGHT
        SIGKE(2)=SIGKE(2)-TAU(2)*WGHT
        SIGKE(3)=SIGKE(3)-TAU(3)*WGHT
        SIGKE(5)=SIGKE(5)-TAU(5)*WGHT
        SIGKE(4)=SIGKE(4)-TAU(4)*WGHT
        SIGKE(6)=SIGKE(6)-TAU(6)*WGHT
 2009 CONTINUE
C=======================================================================
C
C START CONJUGATE GRADIENTS LOOP
C
C=======================================================================
       IF(ICLOCK.EQ.1) CALL PCLOCK(203)
       DO 2001 NN=1,NBANDS
        NINDEN=NN
        NINDKE=0
        NINDW=NRPLWV*(NN-1)
        OCCNN=OCC(NN,NPKPT)
        FACTOR=2.0*WTKPT(NPKPT)*OCCNN
C 
C ......................................................................
C    CALCULATE THE PRODUCT OF HAMILTONIAN AND WAVEFUNCTION
C    
C    H SI
C
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)
          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)
        DO 5000 NIT=1,NITMAX
          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 2002 NNN=1,NPLWV
              CPTNWR(NNN)=(CPTNWR(NNN)+CV(NNN)*CPTOWR(NNN))*RINPLW
 2002       CONTINUE
          ELSE 
            DO 2003 NNN=1,NPLWV                                   
              CPTNWR(NNN)=CV(NNN)*CPTOWR(NNN)*RINPLW
 2003       CONTINUE                                              
          END IF
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.EQ.0 .OR. NLPOT.EQ.1) GOTO 2016
C          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
C          IF(ICLOCK.EQ.1) CALL PCLOCK(204)
C
 2016     CONTINUE
          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 ->>         ** TO ALL BAND VECTORS **    X.WENG 05-MAR-90
C=======================================================================
          DO 2022 NB=1,NBANDS
            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 ALL BANDS
C=======================================================================
          DO 2032 NB=1,NBANDS
            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
              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))
C             CDOT1=CDOT1+(CORGR(M)-COGRPI(M)) *CONJG(CDUM(M))
C
C             CDOT2=CDOT2+COGRPI(M)*CONJG(COGRPI(M))
              CDOT2=CDOT2+COGRPI(M)*CONJG(COGRPI(M))*PRECON(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 AND NORMALISE
C=======================================================================
          COVERL=(0.0,0.0)
          DO 2070 M=1,NPLWKP
            COVERL=COVERL+CDIR(M)*CONJG(CPTWFP(M+NINDW))
 2070     CONTINUE
          DO 2080 M=1,NPLWKP
            CDIR(M)=CDIR(M)-COVERL*CPTWFP(M+NINDW)
 2080     CONTINUE
          ANORM=0.0
          DO 2090 M=1,NPLWKP
            ANORM=ANORM+REAL(CDIR(M)*CONJG(CDIR(M)))
 2090     CONTINUE
          FNORM=1.0/(SQRT(ANORM))
          DO 2095 M=1,NPLWKP
            CDIR(M)=CDIR(M)*FNORM
 2095     CONTINUE
C=======================================================================
C NOW DETERMINE OPTIMUM STEP LENGTH ALONG SEARCH DIRECTION
C STEP 1 DETERMINE THE COEFFICIENT OF THE 'SIN' TERM IN THE ENERGY
C EXPRESSION E = E(AV) + ES SIN(2*THETA) + EC COS(2*THETA)
C=======================================================================
          EN1=0.0
          DO 2100 M=1,NPLWKP
            EN1=EN1+REAL(CDIR(M)*CONJG(CGRA(M)))
 2100     CONTINUE
          ESIN=-EN1*2.0*WTKPT(NPKPT)*OCCNN
C=======================================================================
C CALCULATE ENERY AT THE BEGINNING OF THE STEP
C=======================================================================
          EN2=TOTEN
C=======================================================================
C CALCULATE A NEW VECTOR AT A TRIAL DISTANCE ALONG THE SEARCH DIRECTION
C=======================================================================
          ALPHA=ALPHA0
          ALCOS=COS(ALPHA)
          ALSIN=SIN(ALPHA)
          TALCOS=COS(2.0*ALPHA)
          TALSIN=SIN(2.0*ALPHA)
          DO 2110 M=1,NPLWKP
            CPTWFL(M)=ALCOS*CPTWFP(M+NINDW)+ALSIN*CDIR(M)
 2110     CONTINUE
C=======================================================================
C REPEAT PROCESS OF CALCULATING THE ENERGY BUT WITH THE NEW VECTOR
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 2112 M=1,MPLWV
            CPTNWR(M)=(0.0,0.0)
            CWORK(M)=(0.0,0.0)
 2112     CONTINUE
          DO 2113 M=1,NPLWKP
            CPTNWR(NINDPW(M))=CPTWFL(M)
 2113     CONTINUE
          NCELIN=NN+(NBANDS*(NPKPT-1))
          CELEN(NCELIN)=(0.0,0.0)
          DO 2114 M=1,NPLWKP
            CELEN(NCELIN)=CELEN(NCELIN)+DATAKE(1+7*(M-1))*CPTWFL(M)*
     &                    CONJG(CPTWFL(M))
 2114     CONTINUE
C=======================================================================
C TRANSFORM THE WAVEFUNCTION INTO REAL SPACE, AND FIND THE NEW CHDENR(r)
C=======================================================================
          CALL FFT3D(CPTNWR,CWORK,NGPTAR,1)
          DO 3000 M=1,NPLWV
            CHDENR(M)=CHDENR(M)+FACTOR*
     &         (CPTNWR(M)*CONJG(CPTNWR(M))-CPTOWR(M)*CONJG(CPTOWR(M)))
 3000     CONTINUE
          DO 3010 M=1,NPLWV
            CPTOWR(M)=CPTNWR(M)
 3010     CONTINUE
C=======================================================================
C CALL ROUTINE TO CALCULATE NON-LOCAL ENERGY IN REAL SPACE
C=======================================================================
          IF(IVPTYP.EQ.1 .AND. NLPOT.EQ.1) THEN
            CALL ENRLNL(VOLC,NGX,NGY,NGZ,VNL,NN,NPKPT,MPLWV,NRPLWV,
     &       NRGRPT,NIONSP,NRLPPI,NBANDS,NKPTS,NSPEC,CWORK,CPTNWR,          
     &       CPHGRD,NRLNL,NIONST,MXRLNL,IRLNL,PRLSCA,VRLGRD,NADGRD,
     &       MXRLSH,CESAVE)                      
          ENDIF                                                             
          DO 3020 M=1,NPLWV
            CPTNWR(M)=CHDENR(M)*RINPLW
 3020     CONTINUE
C=======================================================================
C TRANSFORM THE NEW CHARGE DENSITY INTO RECIPROCAL SPACE
C=======================================================================
          CALL FFT3D(CPTNWR,CWORK,NGPTAR,-1)
C======================================================================
C CALC. HARTREE POTENTIAL AND ENERGY CORRECTION IN RECIPROCAL SPACE
C======================================================================
          DENC=0.0
          DO 9122 M=1,NPLWV
            CVD(M)=CPTNWR(M)*DIRDAT(M)
            DENC=DENC+REAL(CVD(M)*CONJG(CPTNWR(M)))
 9122     CONTINUE
          DENC=-DENC/2
          DO 3030 M=1,NPLWV
            CPTNWR(M)=CVD(M)
 3030     CONTINUE
C======================================================================
C TRANSFORM HARTREE POTENTIAL INTO REAL SPACE
C======================================================================
          CALL FFT3D(CPTNWR,CWORK,NGPTAR,1)
          DO 3035 M=1,NPLWV
            CVD(M)=CPTNWR(M)
 3035     CONTINUE
          DO 3040 M=1,NPLWV
            CV(M)=CHDENR(M)
 3040     CONTINUE
C======================================================================
C CALC. EXCHANGE CORRELATION ENERGY CORRECTION USING CHDENR(r)
C======================================================================
          CALL FEXCTR(NPLWV,CV,VOLC,SIGXC,XCENC,XCENER,EXCDAT,XCFDAT,
     &                XCPDAT)
          DO 3050 M=1,NPLWV
            CV(M)=CV(M)+CVION(M)
            CV(M)=CV(M)+CVD(M)
 3050     CONTINUE
C======================================================================
C CALCULATE THE TOTAL ENERGY (IN REAL SPACE)
C======================================================================
          ENPOT=0.0
          DO 1100 M=1,NPLWV
            ENPOT=ENPOT+REAL(CHDENR(M)*CV(M))
 1100     CONTINUE
          ENPOT=ENPOT*RINPLW
C======================================================================
C     THE CONTRIBUTION TO TOTAL ENERGY FROM THE NON-LOCAL PSEUDOPOTENTIAL.
C     THIS TIME, ONLY THAT OF THE NNth BAND NEED TO BE CALCULATED
C     THE TRIAL WAVE FUNCTION IS IN CPTWPL, IT IS IN THE RECIPROCAL SPACE
C======================================================================
          IF(IVPTYP.EQ.0 .OR. NLPOT.EQ.1) GOTO 1103 
C          IF(ICLOCK.EQ.1) CALL PCLOCK(205)
          CALL VNLWAV(NRPLWV,NIONS,NSPEC,NIONSP,NPLWKP,DNLKG(1,0,NPKPT),
     &                VOLC, PSCALE,CPTWFL,CWRK20,CWRK21,CWRK22,
     &                CWRK23,CPHSGR,VGNL,CELFRC,CDUM,CORGR,1,
     &                IVPTYN)
          VNL(NN,NPKPT)=0.0
          DO 1102 M=1,NPLWKP
            VNL(NN,NPKPT)=VNL(NN,NPKPT)+REAL(CONJG(CPTWFL(M))*CELFRC(M))
 1102     CONTINUE
C          IF(ICLOCK.EQ.1) CALL PCLOCK(206)
C
1103      CONTINUE
          ENKE=0.0
          ENVNL=0.0
          DO 1110 NK=1,NKPTS
            DO 1120 NB=1,NBANDS
              ENKE=ENKE+(2.0*WTKPT(NK)*REAL(CELEN((NK-1)*NBANDS+NB)))*
     &             OCC(NB,NK)
              IF(IVPTYP.EQ.1) ENVNL=ENVNL+ 2.0*WTKPT(NK)*VNL(NB,NK)*
     &                              OCC(NB,NK)
 1120       CONTINUE
 1110     CONTINUE
C
          IF(IPRINT.GE.3) THEN 
            WRITE(*,4000) (VNL(NB,NPKPT),NB=1,NBANDS)
            WRITE(*,4001) ENKE,ENPOT
            IF (IVPTYP.EQ.1) WRITE(*,4002) ENVNL
            WRITE(*,4003) DENC,XCENC,TEWEN,PSCENC
          END IF
          TOTEN=ENPOT+ENKE+DENC+XCENC+TEWEN+PSCENC +ENVNL
          EN3=TOTEN
          ECOS=(EN2-EN3+ESIN*TALSIN)/(1.0-TALCOS)
          EAVE=EN2-ECOS
          IF(IPRINT.GE.3) THEN 
            WRITE(*,*)' EN2, EN3:', EN2, EN3
            WRITE(*,*)' ESIN, ECOS, EAVE:', ESIN, ECOS, EAVE
          END IF
C=======================================================================
C CALCULATE STEPLENGTH
C=======================================================================
          IF(ECOS.EQ.0.0) THEN
            TTHETA=0.0
          ELSE
            TTHETA=ATAN(ESIN/ECOS)
          END IF
          THETA=TTHETA/2.0
          IF(THETA.LT.-0.1) THETA=THETA+1.570796327
          IF(IPRINT.GE.3) WRITE(*,*)' THETA:', THETA
C=======================================================================
C UPDATE WAVEFUNCTION
C=======================================================================
          THCOS=COS(THETA)
          THSIN=SIN(THETA)
          DO 2200 M=1,NPLWKP
            CPTWFP(M+NINDW)=THCOS*CPTWFP(M+NINDW)+THSIN*CDIR(M)
 2200     CONTINUE
          DO 3012 M=1,MPLWV
            CPTNWR(M)=(0.0,0.0)
            CWORK(M)=(0.0,0.0)
 3012     CONTINUE
          DO 3013 M=1,NPLWKP
            CPTNWR(NINDPW(M))=CPTWFP(M+NINDW)
 3013     CONTINUE
          NCELIN=NN+(NBANDS*(NPKPT-1))
          CELEN(NCELIN)=(0.0,0.0)
          DO 3014 M=1,NPLWKP
            CELEN(NCELIN)=CELEN(NCELIN)+DATAKE(1+7*(M-1))*
     &                    CPTWFP(M+NINDW)*CONJG(CPTWFP(M+NINDW))
 3014     CONTINUE
C======================================================================
C     UPDATE THE NON-LOCAL ENERGY OF NNth BAND
C======================================================================
          IF(IVPTYP.EQ.0 .OR. NLPOT.EQ.1) GOTO 3016
C          IF(ICLOCK.EQ.1) CALL PCLOCK(207)
          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)
          VNL(NN,NPKPT)=0.0
          DO 3015 M=1,NPLWKP
            VNL(NN,NPKPT)=VNL(NN,NPKPT)+REAL(CONJG(CPTWFP(M+NINDW))*
     &                    CELFRC(M))
 3015    CONTINUE
C          IF(ICLOCK.EQ.1) CALL PCLOCK(208)
C======================================================================
C     UPDATE THE TOTAL CHARGE DENSITY
C======================================================================
3016      CONTINUE
C======================================================================
C     TRANSFORM WAVEFUNCTION INTO REAL SPACE
C======================================================================
          CALL FFT3D(CPTNWR,CWORK,NGPTAR,1)
C======================================================================
C     CALCULATE NEW CHARGE DENSITY
C======================================================================
          DO 3100 M=1,NPLWV
            CHDENR(M)=CHDENR(M)+FACTOR*
     +        ( CPTNWR(M)*CONJG(CPTNWR(M))-CPTOWR(M)*CONJG(CPTOWR(M)) )
 3100     CONTINUE
          DO 3110 M=1,NPLWV
            CPTOWR(M)=CPTNWR(M)
 3110     CONTINUE
C=======================================================================
C CALL ROUTINE TO CALCULATE NON-LOCAL ENERGY IN REAL SPACE
C=======================================================================
         IF(IVPTYP.EQ.1 .AND. NLPOT.EQ.1) THEN                         
            CALL ENRLNL(VOLC,NGX,NGY,NGZ,VNL,NN,NPKPT,MPLWV,NRPLWV,
     &       NRGRPT,NIONSP,NRLPPI,NBANDS,NKPTS,NSPEC,CWORK,CPTNWR,      
     &       CPHGRD,NRLNL,NIONST,MXRLNL,IRLNL,PRLSCA,VRLGRD,NADGRD,
     &       MXRLSH,CESAVE)                                             
          ENDIF
C======================================================================
C     TRANSFORM CHARGE DENSITY INTO RECIPROCAL SPACE
C======================================================================
          DO 3120 M=1,NPLWV
            CPTNWR(M)=CHDENR(M)*RINPLW
 3120     CONTINUE
          CALL FFT3D(CPTNWR,CWORK,NGPTAR,-1)
          DO 3121 M=1,NPLWV
            CHDENG(M)=CPTNWR(M)
 3121     CONTINUE
C======================================================================
C     CALCULATE COULOMB POTENTIAL AND CORRECTION TO TOTAL ENERGY 
C     IN RECIPROCAL SPACE
C======================================================================
          DENC=0.0
          DO 9123 M=1,NPLWV
            CVD(M)=CPTNWR(M)*DIRDAT(M)
            DENC=DENC+REAL(CVD(M)*CONJG(CPTNWR(M)))
 9123     CONTINUE
          DENC=-DENC/2
C======================================================================
C     FFT THE HARTREE POTENTIAL INTO REAL SPACE
C======================================================================
          DO 3130 M=1,NPLWV
            CPTNWR(M)=CVD(M)
 3130     CONTINUE
          CALL FFT3D(CPTNWR,CWORK,NGPTAR,1)
          DO 3135 M=1,NPLWV
            CVD(M)=CPTNWR(M)
 3135     CONTINUE
C======================================================================
C     CALCULATE EXCHANGE-CORRELATION ENERGY IN REAL SPACE
C======================================================================
          DO 3140 M=1,NPLWV
            CV(M)=CHDENR(M)
 3140     CONTINUE
          CALL FEXCTR(NPLWV,CV,VOLC,SIGXC,XCENC,XCENER,EXCDAT,
     &                                             XCFDAT,XCPDAT)
          DO 3150 M=1,NPLWV
            CV(M)=CV(M)+CVION(M)
            CV(M)=CV(M)+CVD(M)
 3150     CONTINUE
C======================================================================
C
C     CV(M) IS THE CHARGE DENSITY IN REAL SPACE
C
C======================================================================
C CALCULATE THE TOTAL ENERGY
C======================================================================
          ENPOT=0.0
          ENKE=0.0
          ENVNL=0.0
          DO 1200 M=1,NPLWV
            ENPOT=ENPOT+REAL(CHDENR(M)*CV(M))
 1200     CONTINUE
          ENPOT=ENPOT*RINPLW
          DO 1210 NK=1,NKPTS
            DO 1220 NB=1,NBANDS
            ENKE=ENKE+(2.0*WTKPT(NK)*REAL(CELEN((NK-1)*NBANDS+NB)))*
     &           OCC(NB,NK)
            IF(IVPTYP.EQ.1) ENVNL=ENVNL+ 2.0*WTKPT(NK)*VNL(NB,NK)*
     &                      OCC(NB,NK)
 1220       CONTINUE
 1210     CONTINUE
          TOTEN=ENPOT+ENKE+DENC+XCENC+TEWEN+PSCENC +ENVNL
C======================================================================
C     MOVE ONTO THE NEXT ITERATION
C======================================================================
 5000   CONTINUE
C======================================================================
C     MOVE ONTO THE NEXT BAND
C======================================================================
 2001 CONTINUE
      IF(ICLOCK.EQ.1) CALL PCLOCK(203)
      IF(IPRINT.GE.3) WRITE (*,4000) (VNL(NB,NPKPT),NB=1,NBANDS)
      IF (IPRINT.GE.2) THEN
        WRITE (*,4001) ENKE,ENPOT
        IF (IVPTYP.EQ.1) WRITE (*,4002) ENVNL
        WRITE (*,4003) DENC,XCENC,TEWEN,PSCENC
      END IF
      TOTEN=ENPOT+ENKE+DENC+XCENC+TEWEN+PSCENC +ENVNL
      IF(IPRINT.GE.1) WRITE(*,1) TOTEN
      IF (ISYMM.NE.0) THEN
        ENVNL1=0.0
        ENKE1 =0.0
        NK = NPKPT
        DO 6011 NB=1,NBANDS
          ENKE1=ENKE1+(2.0*REAL(CELEN((NK-1)*NBANDS+NB)))*
     &         OCC(NB,NK)*(WTKPTT-WTKPT(NK))
          IF(IVPTYP.EQ.1) ENVNL1=ENVNL1+ 2.0*VNL(NB,NK)*
     &                    OCC(NB,NK)*(WTKPTT-WTKPT(NK))
 6011   CONTINUE
        TOTEN2 = TOTEN + ENKE1 + ENVNL1
        IF(IPRINT.GE.1) WRITE(*,1) TOTEN2
      END IF
 1    FORMAT(1X,'THE TOTAL ENERGY IS   ',F16.8)
C=======================================================================
C SUB-SPACE ROTATION (FINDS OUT EIGEN STATES AND EIGEN VALUES) 
C THE CELEN AND VNL ARE UPDATED
C=======================================================================
      IF(ISBROT.EQ.0) GOTO 9999
         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
      ENVNL=0.0
      ENKE =0.0
      DO 6010 NK=1,NKPTS
        DO 6010 NB=1,NBANDS
          ENKE=ENKE+(2.0*WTKPT(NK)*REAL(CELEN((NK-1)*NBANDS+NB)))*
     &         OCC(NB,NK)
          IF(IVPTYP.EQ.1) ENVNL=ENVNL+ 2.0*WTKPT(NK)*VNL(NB,NK)*
     &                    OCC(NB,NK)
 6010 CONTINUE
      IF (IPRINT.GE.1) THEN
        WRITE(*,*)' AFTER SUBROT: ENKE =', ENKE
        IF (IVPTYP.EQ.1) WRITE(*,*)'               ENVNL=', ENVNL
        TOTEN=ENPOT+ENKE+DENC+XCENC+TEWEN+PSCENC +ENVNL
        WRITE (*,1) TOTEN 
      END IF
C======================================================================
C
C======================================================================
 9999 CONTINUE
      IF(IPRINT.GE.2)WRITE(*,*)' BYE FROM CONGRA, VERSION OF 16-JUN-92'
      IF (ISYMM.NE.0) WTKPT(NPKPT)=WTKPTT
      RETURN
      END 
