C***********************************************************************
      SUBROUTINE WFTRI(NGX,NGY,NGZ,NBANDS,NKPTS,NPLWV,
     &                 MPLWV,NRPLWV,NINDPW,NPLWKP,
     &                 WTKPT,CV,CPTWFP,CPTWFL,VOLC,CELEN,NGPTAR,DATAKE,
     &                 CPTOWR,CPTNWR,CWORK,NUNIT,NBANOC, CGRA,
     &                 CORGR,CDIR,COGRPI,CDIRPI,PRECON,
     &                 CHDENR,NPKPT,CDUM,NSPEC,
     &                 NIONS,NIONSP,PSCALE,
     &                 DNLKG,CPHSGR,VGNL,CELFRC,
     &                 CWRK20,CWRK21,CWRK22,CWRK23,VNL,IVPTYP,OCC,
     &                 ENGRSI,ALPHA,CVTRUE,IVPTYN,NLPOT,NRGRPT,NRLPPI,      
     &                 CPHGRD,NRLNL,NIONST,MXRLNL,IRLNL,PRLSCA,VRLGRD,
     &                 NADGRD,MXRLSH,CESAVE)       
C=======================================================================
      IMPLICIT COMPLEX (C)
      DIMENSION NINDPW(*)
      DIMENSION CV(*)
      DIMENSION CPTWFP(*)
      DIMENSION CPTWFL(*)
      DIMENSION CELEN(*)
      DIMENSION NGPTAR(*)
      DIMENSION DATAKE(*)
      DIMENSION CHDENR(*)
      DIMENSION WTKPT(*)
      DIMENSION NIONSP(*)
      DIMENSION OCC(NBANDS,NKPTS)
C=======================================================================
C
C                      DIMENSION STATEMENTS
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 (*)
      DIMENSION CPHSGR(NRPLWV,NIONS,NSPEC),VGNL(NRPLWV,0:2,NSPEC)
      DIMENSION DNLKG(NRPLWV,0:3,NKPTS)
      DIMENSION PSCALE(0:2,NSPEC)
      DIMENSION VNL(NBANDS,NKPTS),CELFRC(NRPLWV)
      DIMENSION CWRK20(NIONS),CWRK21(3,NIONS)
      DIMENSION CWRK22(NIONS),CWRK23(3,3,NIONS)
      DIMENSION CVTRUE(*)
      DIMENSION IVPTYN(*)
      DIMENSION  IRLNL(MXRLNL,NSPEC),PRLSCA(MXRLNL,NSPEC)
      DIMENSION  NRLNL(NSPEC),VRLGRD(NRGRPT,MXRLSH,NIONST)
      DIMENSION  NADGRD(NRGRPT,NIONST),CPHGRD(NRGRPT,NIONST)
      DIMENSION CESAVE(NIONST,20),NRLPPI(NIONST)
C=======================================================================
C     DIMENSION FOR SUBROUTINE ROTATION (WENG, 27-FEB-90)
C=======================================================================
C
      RINPLW = 1.0 / NPLWV
C=======================================================================
C START CONJUGATE GRADIENTS LOOP
C=======================================================================
      DO 2001 NN = 1 , NBANDS
        NINDW = NRPLWV * ( NN - 1 )
C=======================================================================
C CALCULATE THE PRODUCT OF THE POTENTIAL AND THE WAVEFUNCTION IN REAL
C 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 5562 M = 1 , NRPLWV
          CELFRC(M) = (0.0,0.0)
 5562   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 THIS PRODUCT IS SAVED IN CWORK AND TRANSFORMED INTO RECIPROCAL SPACE
C TO BE FURTHER SENT INTO CELFRC ARRAY - AS FOR RECIPROCAL SPACE
C VERSION OF NONLOCAL POTENTIALS.
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) * RINPLW                            
 2002     CONTINUE                                                  
          CALL FFT3D(CPTNWR,CWORK,NGPTAR,-1)
          DO 5511 M = 1 , NPLWKP               
            CELFRC(M) = CPTNWR(NINDPW(M))
 5511     CONTINUE                             
        END IF                                     
        DO 2005 NNN=1,NPLWV                      
          CPTNWR(NNN)=CVTRUE(NNN)*CPTOWR(NNN)*RINPLW 
 2005   CONTINUE                                 
C=======================================================================
C TRANSFORM THE WAVEFUNCTION*POTENTIAL INTO RECIPROCAL SPACE
C=======================================================================
        CALL FFT3D(CPTNWR,CWORK,NGPTAR,-1)
C=======================================================================
C CALCULATE HAMILTONIAN TIMES PSI - STORED IN CPTWFL
C=======================================================================
        DO 2003 M = 1 , NPLWKP
          CPTWFL(M) = CPTWFP(M+NINDW) * DATAKE(1+(7*(M-1))) +
     &               CPTNWR(NINDPW(M))
 2003   CONTINUE
C=======================================================================
C  FOR NON-LOCAL PSEUDOPOTENTIAL WE ADD THE NON-LOCAL CONTRIBUTION
C  ALSO VNL ARE CALCULATED TO BE USED IN THE MAIN PROGRAM WHILE
C  EVALUATING ENVNL FOR NEW IONIC POSITIONS - OLD ELECTRONIC FUNCTIONS
C=======================================================================
        IF (IVPTYP.NE.0) THEN
          IF (NLPOT.EQ.0)
     &    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 1102 M = 1 , NPLWKP
            CPTWFL(M) = CPTWFL(M) + CELFRC(M)
 1102     CONTINUE
        END IF
C=======================================================================
C  H * WF FOR THE SCREENED POTENTIAL
C=======================================================================
        DO 2004 NNN = 1 , NPLWV
          CPTNWR(NNN) = CV(NNN) * CPTOWR(NNN) * RINPLW
 2004   CONTINUE
        CALL FFT3D(CPTNWR,CWORK,NGPTAR,-1)
C=======================================================================
C CALCULATE THE GRADIENT
C=======================================================================
        DO 2010 M = 1 , NPLWKP
          CGRA(M) = CPTWFP(M+NINDW) * DATAKE(1+(7*(M-1))) +
     &              CPTNWR(NINDPW(M)) 
 2010   CONTINUE
        IF (IVPTYP.NE.0) THEN
          DO 2011 M = 1 , NPLWKP
            CGRA(M) = CGRA(M) + CELFRC(M)
 2011     CONTINUE
        END IF
C
        DO 2017 M = 1 , NPLWKP
          CORGR(M) = - CGRA(M)
 2017   CONTINUE
C=======================================================================
C ORTHOGONALISE GRADIENT TO ALL BAND VECTOR
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
          DO 2024 M = 1 , NPLWKP
            CORGR(M) = CORGR(M) - COVERL * CPTWFP(M+NINDD)
 2024     CONTINUE
 2022   CONTINUE
C=======================================================================
C PRE-CONDITION THE GRADIENT
C=======================================================================
        NCELIN = NN + NBANDS * ( NPKPT - 1 ) 
        ENKE = REAL ( CELEN(NCELIN) )
        IF (ENKE.LT.1.0) THEN
          ENKEIN = 1.0
        ELSE
          ENKEIN = 1.0 / ENKE
        ENDIF
        DO 2031 M = 1 , NPLWKP
          X = DATAKE(1+7*(M-1)) * ENKEIN
          PCNUM = 27 + ( 18 * X ) + ( 12 * X * X ) + ( 8 * X * X * X )
          PCDEN = PCNUM + ( 16 * X * X * X * X )
          PRECON(M) = PCNUM / PCDEN
          CDUM(M) = CORGR(M) * PRECON(M)
 2031   CONTINUE
C=======================================================================
C ORTHOGONALISE PRECONDITIONED RESIDUAL TO ALL BANDS
C=======================================================================
        DO 2032 NB = 1 , NBANDS
          NINDD = NRPLWV * ( NB - 1 )
          COVERL = 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=======================================================================
        DO 2035 M = 1 , NPLWKP
          CDIR(M) = CDUM(M)
 2035    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 STORE THE GRADIENT VECTOR
C=======================================================================
        WRITE (NUNIT) ( CDIR(I) , I = 1 , NRPLWV )
C=======================================================================
C CALCULATE THE CONTRIBUTION TO THE ELECTRONIC GRADIENT TERM
C=======================================================================
        EN1 = 0.0
        DO 2100 M = 1 , NPLWKP
          EN1 = EN1 + REAL ( CDIR(M) * CONJG ( CPTWFL(M) ) )
 2100   CONTINUE
        OCCNN = 2.0 * WTKPT(NPKPT) * OCC(NN,NPKPT)
        ENGRSI = ENGRSI + EN1 * OCCNN
C=======================================================================
C CALCULATE A NEW VECTOR AT A TRIAL DISTANCE ALONG THE SEARCH DIRECTION
C=======================================================================
        ALCOS = COS ( ALPHA )
        ALSIN = SIN ( 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
        IF (IVPTYP.NE.0 .AND. NLPOT.EQ.0) THEN
          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 1103 M = 1 , NPLWKP
            VNL(NN,NPKPT) = VNL(NN,NPKPT) + 
     &                      REAL ( CONJG (CPTWFL(M)) * CELFRC(M) ) 
 1103     CONTINUE
        END IF
C=======================================================================
C TRANSFORM THE WAVEFUNCTION INTO REAL SPACE
C=======================================================================
        CALL FFT3D(CPTNWR,CWORK,NGPTAR,1)
        DO 3000 M = 1 , NPLWV
          CHDENR(M) = CHDENR(M) + OCCNN *
     &              ( CPTNWR(M) * CONJG ( CPTNWR(M) ) -
     &                CPTOWR(M) * CONJG ( CPTOWR(M) ) ) 
 3000   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 END LOOPS
C======================================================================
 2001 CONTINUE
      RETURN
      END
